Rivet analyses referenceBESIII_2017_I1507887Analysis of $\psi(2S)\to\gamma\chi_{c(1,2)}$ decays using $\chi_{c(1,2)}\to J/\psi\gamma$Experiment: BESIII (BEPC) Inspire ID: 1507887 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (1.8, 1.8) GeV Run details:
Analysis of the angular distribution of the photons and leptons produced in $e^+e^-\to \psi(2S) \to \gamma\chi_{c(1,2)}$ followed by $\chi_{c(1,2)}\to\gamma J/\psi$ and $J/\psi\to\ell^+\ell^-$ Gives information about the decay and is useful for testing correlations in charmonium decays. N.B. the data was read from the figures in the paper and is not corrected and should only be used qualatively. Source code: BESIII_2017_I1507887.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 psi(2S) -> gamma chi_c1,2
11 class BESIII_2017_I1507887 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2017_I1507887);
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==20443 || Cuts::pid==445), "UFS");
26 declare(FinalState(), "FS");
27 for(unsigned int ix=0;ix<10;++ix)
28 book(_h[ix],1,1,1+ix);
29 }
30
31 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
32 for( const Particle &child : p.children()) {
33 if(child.children().empty()) {
34 nRes[child.pid()]-=1;
35 --ncount;
36 }
37 else
38 findChildren(child,nRes,ncount);
39 }
40 }
41
42 // angle cuts due regions of BES calorimeter
43 bool vetoPhoton(const double & cTheta) {
44 return cTheta>0.92 || (cTheta>0.8 && cTheta<0.86);
45 }
46
47 /// Perform the per-event analysis
48 void analyze(const Event& event) {
49 // cos of 10 degress for cut
50 static const double cos10 = 0.984807753012208;
51 // get the axis, direction of incoming electron
52 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
53 Vector3 axis;
54 if(beams.first.pid()>0)
55 axis = beams.first .momentum().p3().unit();
56 else
57 axis = beams.second.momentum().p3().unit();
58 // types of final state particles
59 const FinalState& fs = apply<FinalState>(event, "FS");
60 map<long,int> nCount;
61 int ntotal(0);
62 for (const Particle& p : fs.particles()) {
63 nCount[p.pid()] += 1;
64 ++ntotal;
65 }
66 // loop over chi_c states
67 Particle chi;
68 bool matched = false;
69 const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
70 for (const Particle& p : ufs.particles()) {
71 if(p.children().empty()) continue;
72 map<long,int> nRes=nCount;
73 int ncount = ntotal;
74 findChildren(p,nRes,ncount);
75 if(ncount==1) {
76 matched = true;
77 for(auto const & val : nRes) {
78 if(val.first==PID::PHOTON) {
79 if(val.second!=1) {
80 matched = false;
81 break;
82 }
83 }
84 else if(val.second!=0) {
85 matched = false;
86 break;
87 }
88 }
89 if(matched) {
90 chi=p;
91 break;
92 }
93 }
94 }
95 if(!matched) vetoEvent;
96 // have chi_c find psi2S
97 if(chi.parents().empty() || chi.children().size()!=2) vetoEvent;
98 Particle psi2S = chi.parents()[0];
99 if(psi2S.pid()!=100443 || psi2S.children().size()!=2) vetoEvent;
100 // then the first photon
101 Particle gamma1;
102 if(psi2S.children()[0].pid()==PID::PHOTON)
103 gamma1 = psi2S.children()[0];
104 else if(psi2S.children()[1].pid()==PID::PHOTON)
105 gamma1 = psi2S.children()[1];
106 else
107 vetoEvent;
108 // cuts on the photon
109 if(vetoPhoton(abs(axis.dot(gamma1.p3().unit())))) vetoEvent;
110 // then the J/psi and second photon
111 Particle JPsi,gamma2;
112 if(chi.children()[0].pid()==PID::PHOTON &&
113 chi.children()[1].pid()==443) {
114 gamma2 = chi.children()[0];
115 JPsi = chi.children()[1];
116 }
117 else if(chi.children()[1].pid()==PID::PHOTON &&
118 chi.children()[0].pid()==443) {
119 gamma2 = chi.children()[1];
120 JPsi = chi.children()[0];
121 }
122 else
123 vetoEvent;
124 // cuts on the photon
125 if(vetoPhoton(abs(axis.dot(gamma2.p3().unit())))) vetoEvent;
126 // finally the leptons from J/psi decay
127 if(JPsi.children().size()!=2) vetoEvent;
128 if(JPsi.children()[0].pid()!=-JPsi.children()[1].pid()) vetoEvent;
129 if(JPsi.children()[0].abspid()!=PID::EMINUS &&
130 JPsi.children()[0].abspid()!=PID::MUON) vetoEvent;
131 Particle lm = JPsi.children()[0];
132 Particle lp = JPsi.children()[1];
133 if(lm.pid()<0) swap(lm,lp);
134 // cut between photons and charged tracks and on charged tracks
135 Vector3 dGamma[2] = {gamma1.momentum().p3().unit(),
136 gamma1.momentum().p3().unit()};
137 Vector3 dl [2] = {lm.momentum().p3().unit(),
138 lp.momentum().p3().unit()};
139 for(unsigned int ix=0;ix<2;++ix) {
140 // angle cut for charged tracks
141 if(abs(axis.dot(dl[ix]))>0.93) vetoEvent;
142 // angle between leptons and photons
143 for(unsigned int iy=0;iy<2;++iy)
144 if(abs(dGamma[ix].dot(dl[iy]))>cos10) vetoEvent;
145 }
146 // type chi state
147 unsigned int ichi= chi.pid()==445 ? 5 : 0;
148 // first angle of gamma1 w.r.t beam
149 _h[ichi]->fill(axis.dot(gamma1.momentum().p3().unit()));
150 // axis in the chi frame
151 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(chi.momentum().betaVec());
152 Vector3 e1z = gamma1.momentum().p3().unit();
153 Vector3 e1y = e1z.cross(axis).unit();
154 Vector3 e1x = e1y.cross(e1z).unit();
155 // cos theta_2 and phi 2 distributions
156 FourMomentum pGamma2 = boost1.transform(gamma2.momentum());
157 Vector3 axis1 = pGamma2.p3().unit();
158 _h[ichi+1]->fill(e1z.dot(axis1));
159 double phi2 = atan2(e1y.dot(axis1),e1x.dot(axis1));
160 if(phi2<-3.) phi2+=2.*M_PI;
161 _h[ichi+3]->fill(phi2);
162 // cos theta_3 and phi 3 distributions
163 FourMomentum pJpsi = boost1.transform(JPsi.momentum());
164 FourMomentum plp = boost1.transform( lp.momentum());
165 Vector3 axis3 = boost1.transform(gamma1.momentum()).p3().unit();
166 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pJpsi.betaVec());
167 Vector3 axis2 = boost2.transform(plp).p3().unit();
168 Vector3 e2z = gamma2.momentum().p3().unit();
169 Vector3 e2y = e2z.cross(axis3).unit();
170 Vector3 e2x = e2y.cross(e2z).unit();
171 _h[ichi+2]->fill(e2z.dot(axis2));
172 double phi3 = atan2(e2y.dot(axis2),e2x.dot(axis2));
173 if(phi3<-3.) phi3+=2.*M_PI;
174 _h[ichi+4]->fill(phi3);
175
176 }
177
178
179 /// Normalise histograms etc., after the run
180 void finalize() {
181 for(unsigned int ix=0;ix<10;++ix) {
182 normalize(_h[ix]);
183 }
184 }
185
186 /// @}
187
188
189 /// @name Histograms
190 /// @{
191 Histo1DPtr _h[10];
192 /// @}
193
194
195 };
196
197
198 RIVET_DECLARE_PLUGIN(BESIII_2017_I1507887);
199
200}
|