Rivet analyses referenceBESIII_2013_I1203841Mass and angular distributions in $J/\psi\to\gamma\omega\phi$ decaysExperiment: BESIII (BEPC) Inspire ID: 1203841 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (1.6, 1.6) GeV Run details:
Measurement of mass and angular distributions in $J/\psi\to\gamma\omega\phi$ decays. Source code: BESIII_2013_I1203841.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5#include "Rivet/Projections/DecayedParticles.hh"
6
7namespace Rivet {
8
9
10 /// @brief J/psi -> gamma omega omega
11 class BESIII_2013_I1203841 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2013_I1203841);
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 UnstableParticles ufs = UnstableParticles(Cuts::abspid==443);
25 declare(ufs, "UFS");
26 DecayedParticles PSI(ufs);
27 PSI.addStable(PID::PHI);
28 PSI.addStable(PID::OMEGA);
29 declare(PSI, "PSI");
30 declare(Beam(), "Beams");
31 // histograms
32 for(unsigned int ix=0;ix<9;++ix)
33 book(_h[ix],1,1,1+ix);
34 }
35
36 // angle cuts due regions of BES calorimeter
37 bool vetoPhoton(const double & cTheta) {
38 return cTheta>0.92 || (cTheta>0.8 && cTheta<0.86);
39 }
40
41 void findChildren(const Particle & p, Particles & pim, Particles & pip,
42 Particles & pi0, unsigned int &ncount) {
43 for( const Particle &child : p.children()) {
44 if(child.pid()==PID::PIPLUS) {
45 pip.push_back(child);
46 ncount+=1;
47 }
48 else if(child.pid()==PID::PIMINUS) {
49 pim.push_back(child);
50 ncount+=1;
51 }
52 else if(child.pid()==PID::PI0) {
53 pi0.push_back(child);
54 ncount+=1;
55 }
56 else if(child.children().empty()) {
57 ncount+=1;
58 }
59 else
60 findChildren(child,pim,pip,pi0,ncount);
61 }
62 }
63
64 /// Perform the per-event analysis
65 void analyze(const Event& event) {
66 // get the axis, direction of incoming electron
67 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
68 Vector3 axis;
69 if(beams.first.pid()>0)
70 axis = beams.first .momentum().p3().unit();
71 else
72 axis = beams.second.momentum().p3().unit();
73 // find the J/psi decays
74 static const map<PdgId,unsigned int> & mode = { { 223,1}, { 333,1},{ 22,1}};
75 DecayedParticles PSI = apply<DecayedParticles>(event, "PSI");
76 if( PSI.decaying().size()!=1) vetoEvent;
77 if(!PSI.modeMatches(0,3,mode)) vetoEvent;
78 // particles
79 const Particle & phi = PSI.decayProducts()[0].at(333)[0];
80 const Particle & omega = PSI.decayProducts()[0].at(223)[0];
81 const Particle & gam = PSI.decayProducts()[0].at( 22)[0];
82 _h[0]->fill((omega.momentum()+phi.momentum()).mass());
83 _h[1]->fill((omega.momentum()+gam.momentum()).mass());
84 _h[2]->fill((phi .momentum()+gam.momentum()).mass());
85 double cTheta = axis.dot(gam.p3().unit());
86 if(vetoPhoton(abs(cTheta))) vetoEvent;
87 _h[3]->fill(cTheta);
88 // remaining angles
89 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(PSI.decaying()[0].momentum().betaVec());
90 FourMomentum pGamma = boost1.transform(gam.momentum());
91 FourMomentum pOmegaPhi = boost1.transform(omega.momentum()+phi.momentum());
92 Vector3 e1z = pGamma.p3().unit();
93 Vector3 e1y = e1z.cross(axis).unit();
94 Vector3 e1x = e1y.cross(e1z).unit();
95 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pOmegaPhi.betaVec());
96 Vector3 axis2 = boost2.transform(boost1.transform(phi.momentum())).p3().unit();
97 _h[5]->fill(e1z.dot(axis2));
98 double phiPhi = atan2(axis2.dot(e1y),axis2.dot(e1x));
99 if(phiPhi<0.) phiPhi+=2.*M_PI;
100 _h[7]->fill(phiPhi);
101 // now for the phi decays
102 if(phi.children().size()!=2|| phi.children()[0].pid()!=-phi.children()[1].pid() ||
103 phi.children()[0].abspid()!=321) vetoEvent;
104 Particle Km = phi.children()[0];
105 Particle Kp = phi.children()[1];
106 if(Kp.pid()<0) swap(Km,Kp);
107 FourMomentum pKp = boost2.transform(boost1.transform(Kp.momentum()));
108 FourMomentum pPhi = boost2.transform(boost1.transform(phi.momentum()));
109 LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pPhi.betaVec());
110 pKp = boost3.transform(pKp);
111 double cK = axis2.dot(pKp.p3().unit());
112 _h[6]->fill(cK);
113 // omega decay
114 unsigned int ncount=0;
115 Particles pip,pim,pi0;
116 findChildren(omega,pim,pip,pi0,ncount);
117 if( ncount!=3 || !(pim.size()==1 && pip.size()==1 && pi0.size()==1)) vetoEvent;
118 // boost to omega/phi frame
119 FourMomentum ppip = boost2.transform(boost1.transform(pip[0].momentum()));
120 FourMomentum ppim = boost2.transform(boost1.transform(pim[0].momentum()));
121 FourMomentum pOmega = boost2.transform(boost1.transform(omega.momentum()));
122 LorentzTransform boost4 = LorentzTransform::mkFrameTransformFromBeta(pOmega.betaVec());
123 Vector3 axisZ = pOmega.p3().unit();
124 ppip = boost4.transform(ppip);
125 ppim = boost4.transform(ppim);
126 Vector3 norm = ppip.p3().cross(ppim.p3()).unit();
127 double cOmega = norm.dot(axisZ);
128 _h[4]->fill(cOmega);
129 // angle between planes
130 Vector3 Trans1 = pKp.p3() - cK*pKp.p3().mod()*axis2;
131 Vector3 Trans2 = norm - cOmega*axisZ;
132 double chi = atan(Trans1.cross(Trans2).dot(axis2)/Trans1.dot(Trans2));
133 _h[8]->fill(abs(chi)/M_PI*180.);
134 }
135
136
137 /// Normalise histograms etc., after the run
138 void finalize() {
139 for(unsigned int ix=0;ix<9;++ix)
140 normalize(_h[ix],1.,false);
141 }
142
143 /// @}
144
145
146 /// @name Histograms
147 /// @{
148 Histo1DPtr _h[9];
149 /// @}
150
151
152 };
153
154
155 RIVET_DECLARE_PLUGIN(BESIII_2013_I1203841);
156
157}
|