Rivet analyses referenceBABAR_2012_I1122034$X(3915)$ production in $\gamma\gamma\to J/\psi \omega$Experiment: BABAR (PEP-II) Inspire ID: 1122034 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (5.3, 5.3) GeV Run details:
Measurement of the mass and angle distributions in $\gamma\gamma\to J/\psi \omega$. Source code: BABAR_2012_I1122034.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5#include "Rivet/Projections/Beam.hh"
6
7namespace Rivet {
8
9
10 /// @brief gamma gamma -> X(3915) -> J/psi omega
11 class BABAR_2012_I1122034 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2012_I1122034);
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(FinalState(),"FS");
26 declare(UnstableParticles(Cuts::pid==223 or Cuts::pid==443), "UFS");
27 // histograms
28 book(_h_mass,1,1,1);
29 for (unsigned int ix=0;ix<4;++ix) {
30 book(_h_angle1[ix],2,1,1+ix);
31 if (ix<3) book(_h_angle2[ix],3,1,1+ix);
32 }
33 }
34
35 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
36 for (const Particle &child : p.children()) {
37 if (child.children().empty()) {
38 --nRes[child.pid()];
39 --ncount;
40 } else {
41 findChildren(child,nRes,ncount);
42 }
43 }
44 }
45
46 bool findScattered(Particle beam, double& q2) {
47 bool found = false;
48 Particle scat = beam;
49 while (!scat.children().empty()) {
50 found = false;
51 for (const Particle & p : scat.children()) {
52 if (p.pid()==scat.pid()) {
53 scat=p;
54 found=true;
55 break;
56 }
57 }
58 if (!found) break;
59 }
60 if (!found) return false;
61 q2 = -(beam.mom() - scat.mom()).mass2();
62 return true;
63 }
64
65 void findChildren(const Particle & p, Particles & pim, Particles & pip,
66 Particles & pi0, unsigned int &ncount) {
67 for (const Particle &child : p.children()) {
68 if(child.pid()==PID::PIPLUS) {
69 pip.push_back(child);
70 ncount+=1;
71 }
72 else if(child.pid()==PID::PIMINUS) {
73 pim.push_back(child);
74 ncount+=1;
75 }
76 else if(child.pid()==PID::PI0) {
77 pi0.push_back(child);
78 ncount+=1;
79 }
80 else if(child.children().empty()) {
81 ncount+=1;
82 }
83 else {
84 findChildren(child,pim,pip,pi0,ncount);
85 }
86 }
87 }
88
89 /// Perform the per-event analysis
90 void analyze(const Event& event) {
91 // find scattered leptons and calc Q2
92 const Beam& beams = apply<Beam>(event, "Beams");
93 double q12 = -1, q22 = -1;
94 if (!findScattered(beams.beams().first, q12)) vetoEvent;
95 if (!findScattered(beams.beams().second, q22)) vetoEvent;
96 // check the final state
97 const FinalState& fs = apply<FinalState>(event, "FS");
98 map<long,int> nCount;
99 int ntotal(0);
100 for (const Particle& p : fs.particles()) {
101 nCount[p.pid()] += 1;
102 ++ntotal;
103 }
104 // find the J/psi
105 const FinalState& ufs = apply<FinalState>(event, "UFS");
106 Particle omega,psi;
107 bool found=false;
108 for (const Particle &p1 : ufs.particles(Cuts::pid==443)) {
109 if (p1.children().empty()) continue;
110 map<long,int> nRes = nCount;
111 int ncount = ntotal;
112 findChildren(p1,nRes,ntotal);
113 for (const Particle& p2 : ufs.particles(Cuts::pid==223)) {
114 if (p2.children().empty()) continue;
115 map<long,int> nRes2 = nRes;
116 int ncount2 = ncount;
117 findChildren(p2,nRes2,ncount2);
118 found = true;
119 for (const auto& val : nRes2) {
120 if (abs(val.first)==11) {
121 if (val.second!=1) {
122 found = false;
123 break;
124 }
125 }
126 else if(val.second!=0) {
127 found = false;
128 break;
129 }
130 }
131 if (found) {
132 psi = p1;
133 omega = p2;
134 break;
135 }
136 }
137 }
138 if (!found) vetoEvent;
139 FourMomentum psum = omega.mom()+psi.mom();
140 if (psum.pT()>0.2) vetoEvent;
141 // mass distribution
142 _h_mass->fill(psum.mass());
143 // from now on we need specific decay modes of J/psi and omega
144 // first J/psi -> l+l-
145 if (psi.children().size()!=2) vetoEvent;
146 if (psi.children()[0].pid()!=-psi.children()[1].pid()) vetoEvent;
147 if (psi.children()[0].abspid()!=11 && psi.children()[0].abspid()!=13) vetoEvent;
148 Particle ep = psi.children()[0];
149 Particle em = psi.children()[1];
150 if (ep.pid()>0) swap(ep,em);
151 // omega decay
152 unsigned int ncount=0;
153 Particles pip,pim,pi0;
154 findChildren(omega,pim,pip,pi0,ncount);
155 if (ncount!=3 || !(pim.size()==1 && pip.size()==1 && pi0.size()==1)) vetoEvent;
156 // boost to gamma gamma frame
157 LorentzTransform boostCMS = LorentzTransform::mkFrameTransformFromBeta(psum.betaVec());
158 FourMomentum pPsi = boostCMS.transform(psi .mom());
159 FourMomentum pOmega = boostCMS.transform(omega .mom());
160 FourMomentum pLp = boostCMS.transform(ep .mom());
161 FourMomentum pPip = boostCMS.transform(pip[0].mom());
162 FourMomentum pPim = boostCMS.transform(pim[0].mom());
163 Vector3 axis(0.,0.,1.);
164 // lepton angle
165 double cosL = pLp.p3().unit().dot(axis);
166 _h_angle1[0]->fill(cosL);
167 LorentzTransform boostOmega = LorentzTransform::mkFrameTransformFromBeta(pOmega.betaVec());
168 pPip = boostOmega.transform(pPip);
169 pPim = boostOmega.transform(pPim);
170 // omega decay plane normal angle
171 Vector3 axisOmega = pPip.p3().cross(pPim.p3()).unit();
172 double cosN = axisOmega.dot(axis);
173 _h_angle1[1]->fill(cosN);
174 // angle lepton and omega
175 _h_angle1[2]->fill(pLp.p3().unit().dot(axisOmega));
176 // helicity angle
177 _h_angle1[3]->fill(pPsi.p3().unit().dot(psum.p3().unit()));
178 // now the new frame
179 Vector3 axisZ = pOmega.p3().unit();
180 Vector3 axisY = axisZ.cross(axisOmega);
181 Vector3 axisX = axisY.cross(axisZ);
182 // second set of angles
183 LorentzTransform boostPsi = LorentzTransform::mkFrameTransformFromBeta(pPsi.betaVec());
184 Vector3 axisL = boostPsi.transform(pLp).p3().unit();
185 _h_angle2[0]->fill(axisZ.dot(axisOmega));
186 _h_angle2[1]->fill(axisL.dot(pPsi.p3().unit()));
187 axisZ *=-1.;
188 axisX *=-1.;
189 Vector3 axisnp = axisL.cross(axisZ).unit();
190 double phiL = atan2(axisnp.dot(axisY),axisnp.dot(axisX));
191 double phiN = atan2(axisOmega.dot(axisY),axisOmega.dot(axisX));
192 _h_angle2[2]->fill(mapAngleMPiToPi(phiL-phiN)/M_PI*180.);
193 }
194
195
196 /// Normalise histograms etc., after the run
197 void finalize() {
198 normalize(_h_mass, 1.0, false);
199 normalize(_h_angle1, 1.0, false);
200 normalize(_h_angle2, 1.0, false);
201 }
202
203 /// @}
204
205
206 /// @name Histograms
207 /// @{
208 Histo1DPtr _h_mass,_h_angle1[4],_h_angle2[3];
209 /// @}
210
211
212 };
213
214
215 RIVET_DECLARE_PLUGIN(BABAR_2012_I1122034);
216
217}
|