Rivet analyses referenceBABAR_2006_I719949$e^+e^-\to\rho^0\rho^0$ and $\rho^0\phi^0$ at $\sqrt{s}=10.58\,$GeVExperiment: BABAR (PEP-II) Inspire ID: 719949 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (5.3, 5.3) GeV Run details:
Measurement of the cross section, production angle and helicity angles in the decays for $e^+e^-\to\rho^0\rho^0$ and $\rho^0\phi^0$ at $\sqrt{s}=10.58\,$GeV. The cross section was taken from the tet of the paper and the corrected angular distributionss from figures 5 and 6. Source code: BABAR_2006_I719949.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/UnstableParticles.hh"
6
7namespace Rivet {
8
9
10 /// @brief e+e -> rho 0 rho and rho0 phi
11 class BABAR_2006_I719949 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2006_I719949);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23 // Initialise and register projections
24 declare(Beam(), "Beams");
25 declare(UnstableParticles(Cuts::pid==113 or
26 Cuts::pid==333), "UFS");
27 declare(FinalState(), "FS");
28 // histos
29 for(unsigned int ix=0;ix<3;++ix) {
30 book(_h_hel [ix],3,1,1+ix);
31 if(ix==2) continue;
32 book(_h_sigma[ix],1,1,1+ix);
33 book(_h_prod [ix],2,1,1+ix);
34 }
35 }
36
37 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
38 for( const Particle &child : p.children()) {
39 if(child.children().empty()) {
40 nRes[child.pid()]-=1;
41 --ncount;
42 }
43 else
44 findChildren(child,nRes,ncount);
45 }
46 }
47
48 /// Perform the per-event analysis
49 void analyze(const Event& event) {
50 // get the axis, direction of incoming electron
51 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
52 Vector3 axis;
53 if(beams.first.pid()>0)
54 axis = beams.first .momentum().p3().unit();
55 else
56 axis = beams.second.momentum().p3().unit();
57 // types of final state particles
58 const FinalState& fs = apply<FinalState>(event, "FS");
59 map<long,int> nCount;
60 int ntotal(0);
61 for (const Particle& p : fs.particles()) {
62 nCount[p.pid()] += 1;
63 ++ntotal;
64 }
65 // loop over rho mesons
66 const Particles vMesons = apply<UnstableParticles>(event, "UFS").particles();
67 Particle vectors[2];
68 bool matched(false);
69 for (unsigned int ix=0;ix<vMesons.size();++ix) {
70 if(vMesons[ix].children().empty()) continue;
71 map<long,int> nRes=nCount;
72 int ncount = ntotal;
73 findChildren(vMesons[ix],nRes,ncount);
74 matched=false;
75 for (unsigned int iy=ix+1;iy<vMesons.size();++iy) {
76 if(vMesons[iy].children().empty()) continue;
77 if(vMesons[ix].pid()==333 && vMesons[iy].pid()==333) continue;
78 map<long,int> nRes2=nRes;
79 int ncount2 = ncount;
80 findChildren(vMesons[iy],nRes2,ncount2);
81 if(ncount2==0) {
82 matched = true;
83 for(auto const & val : nRes2) {
84 if(val.second!=0) {
85 matched = false;
86 break;
87 }
88 }
89 if(matched) {
90 vectors[0] = vMesons[ix];
91 vectors[1] = vMesons[iy];
92 break;
93 }
94 }
95 }
96 if(matched) break;
97 }
98 if(!matched) vetoEvent;
99 if(vectors[0].pid()==333) swap(vectors[0],vectors[1]);
100 if (vectors[0].children().size()!=2) vetoEvent;
101 if (vectors[1].children().size()!=2) vetoEvent;
102 double cTheta = abs(axis.dot(vectors[0].momentum().p3().unit()));
103 if(cTheta>0.8) vetoEvent;
104 if(vectors[0].pid()==vectors[1].pid()) {
105 _h_sigma[0]->fill("10.58"s);
106 _h_prod [0]->fill(cTheta);
107 }
108 else {
109 _h_sigma[1]->fill("10.58"s);
110 _h_prod [1]->fill(cTheta);
111 }
112 // helicity angles
113 double cHel[2];
114 for(unsigned int ix=0;ix<2;++ix) {
115 int iMeson = vectors[ix].pid() == 113 ? 211 : 321;
116 Particle mP;
117 if(vectors[ix].children()[0].pid()== iMeson &&
118 vectors[ix].children()[1].pid()==-iMeson)
119 mP = vectors[ix].children()[0];
120 else if(vectors[ix].children()[1].pid()== iMeson &&
121 vectors[ix].children()[0].pid()==-iMeson)
122 mP = vectors[ix].children()[1];
123 else
124 vetoEvent;
125 // boost to the rho+ rest frame
126 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(vectors[ix].momentum().betaVec());
127 Vector3 e1z = vectors[ix].momentum().p3().unit();
128 Vector3 axis1 = boost1.transform(mP.momentum()).p3().unit();
129 cHel[ix] = e1z.dot(axis1);
130 }
131 if(vectors[0].pid()==vectors[1].pid()) {
132 _h_hel[0]->fill(cHel[0]);
133 _h_hel[0]->fill(cHel[1]);
134 }
135 else {
136 _h_hel[1]->fill(cHel[1]);
137 _h_hel[2]->fill(cHel[0]);
138 }
139 }
140
141
142 /// Normalise histograms etc., after the run
143 void finalize() {
144 double fact = crossSection()/sumOfWeights()/femtobarn;
145 for(unsigned int ix=0;ix<3;++ix) {
146 normalize(_h_hel [ix],1.,false);
147 if(ix==2) continue;
148 scale(_h_sigma[ix],fact);
149 normalize(_h_prod [ix],1.,false);
150 }
151 }
152
153 /// @}
154
155
156 /// @name Histograms
157 /// @{
158 BinnedHistoPtr<string> _h_sigma[2];
159 Histo1DPtr _h_prod[2],_h_hel[3];
160 /// @}
161
162
163 };
164
165
166 RIVET_DECLARE_PLUGIN(BABAR_2006_I719949);
167
168}
|