Rivet analyses referenceBESIII_2011_I931195Analysis of $\psi(2S)\to\gamma\chi_{c(0,2)}$ decays using $\chi_{c(0,2)}\to \pi^+\pi^-/K^+K^-$Experiment: BESIII (BEPC) Inspire ID: 931195 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (1.8, 1.8) GeV Run details:
Analysis of the angular distribution of the photons and mesons produced in $e^+e^-\to \psi(2S) \to \gamma\chi_{c(0,2)}$ followed by $\chi_{c(0,2)}\to \pi^+\pi^-/K^+K^-$. Gives information about the decay and is useful for testing correlations in charmonium decays. N.B. the distributions were read from the figures in the paper and are not corrected and should only be used qualatively, however the $x$ and $y$ and multipole results are fully corrected. Source code: BESIII_2011_I931195.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_c0,2
11 class BESIII_2011_I931195 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2011_I931195);
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==10441 || Cuts::pid==445), "UFS");
26 declare(FinalState(), "FS");
27 for (unsigned int ichi=0;ichi<2;++ichi) {
28 for (unsigned int imeson=0;imeson<2;++imeson) {
29 book(_h_thy[ichi][imeson][0],"TMP/h_"+toString(ichi+1)+"_"+toString(imeson+1)+"_gamma",50,-1.,1.);
30 book(_h_thy[ichi][imeson][1],"TMP/h_"+toString(ichi+1)+"_"+toString(imeson+1)+"_meson",50,-1.,1.);
31 book(_h_thy[ichi][imeson][2],"TMP/h_"+toString(ichi+1)+"_"+toString(imeson+1)+"_phi" ,50,0.,2.*M_PI);
32 for(unsigned int iy=0;iy<3;++iy) {
33 book(_h_exp[ichi][imeson][iy],5+ichi,1+imeson,1+iy);
34 }
35 }
36 }
37 for (unsigned int ix=0;ix<3;++ix) {
38 for (unsigned int iy=0;iy<3;++iy) {
39 book(_c[ix][iy],"TMP/c_"+toString(ix+1)+"_"+toString(iy+1));
40 }
41 }
42 }
43
44 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
45 for (const Particle &child : p.children()) {
46 if (child.children().empty()) {
47 nRes[child.pid()]-=1;
48 --ncount;
49 }
50 else findChildren(child,nRes,ncount);
51 }
52 }
53
54 /// Perform the per-event analysis
55 void analyze(const Event& event) {
56 static const double cos20=0.9396926207859084;
57 // get the axis, direction of incoming electron
58 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
59 Vector3 axis;
60 if(beams.first.pid()>0)
61 axis = beams.first .momentum().p3().unit();
62 else
63 axis = beams.second.momentum().p3().unit();
64 // types of final state particles
65 const FinalState& fs = apply<FinalState>(event, "FS");
66 map<long,int> nCount;
67 int ntotal(0);
68 for (const Particle& p : fs.particles()) {
69 nCount[p.pid()] += 1;
70 ++ntotal;
71 }
72 // loop over chi_c states
73 Particle chi;
74 bool matched = false;
75 const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
76 for (const Particle& p : ufs.particles()) {
77 if(p.children().empty()) continue;
78 map<long,int> nRes=nCount;
79 int ncount = ntotal;
80 findChildren(p,nRes,ncount);
81 if(ncount==1) {
82 matched = true;
83 for(auto const & val : nRes) {
84 if(val.first==PID::PHOTON) {
85 if(val.second!=1) {
86 matched = false;
87 break;
88 }
89 }
90 else if(val.second!=0) {
91 matched = false;
92 break;
93 }
94 }
95 if(matched) {
96 chi=p;
97 break;
98 }
99 }
100 }
101 if(!matched) vetoEvent;
102 // have chi_c find psi2S
103 if(chi.parents().empty() || chi.children().size()!=2 ||
104 chi.children()[0].pid() != -chi.children()[1].pid()) vetoEvent;
105 Particle psi2S = chi.parents()[0];
106 if(psi2S.pid()!=100443 || psi2S.children().size()!=2) vetoEvent;
107 // then the first photon
108 Particle gamma1;
109 if(psi2S.children()[0].pid()==PID::PHOTON)
110 gamma1 = psi2S.children()[0];
111 else if(psi2S.children()[1].pid()==PID::PHOTON)
112 gamma1 = psi2S.children()[1];
113 else
114 vetoEvent;
115 // now the decay products of the chi_c
116 Particle mPlus,mMinus;
117 bool foundMeson=false;
118 for (unsigned int ix=0;ix<2;++ix) {
119 if(chi.children()[ix].pid()==PID::PIPLUS ||
120 chi.children()[ix].pid()==PID::KPLUS ) {
121 foundMeson=true;
122 mPlus=chi.children()[ix];
123 }
124 else if(chi.children()[ix].pid()==PID::PIMINUS ||
125 chi.children()[ix].pid()==PID::KMINUS ) {
126 mMinus=chi.children()[ix];
127 }
128 }
129 if(!foundMeson) vetoEvent;
130 // cut on photon angles
131 Vector3 aGamma = gamma1.p3().unit();
132 double cGammaCut = abs(axis.dot(aGamma));
133 // type chi state
134 unsigned int ichi= chi.pid()==445 ? 0 : 1;
135 // type of meson
136 unsigned int imeson = mPlus.pid()==PID::PIPLUS ? 0 : 1;
137 // first angle of gamma1 w.r.t beam
138 double cGamma = axis.dot(gamma1.momentum().p3().unit());
139 _h_thy[ichi][imeson][0]->fill(cGamma);
140 // axis in the chi frame
141 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(chi.momentum().betaVec());
142 Vector3 e1z = gamma1.momentum().p3().unit();
143 Vector3 e1y = e1z.cross(axis).unit();
144 Vector3 e1x = e1y.cross(e1z).unit();
145 FourMomentum pMeson = boost1.transform(mPlus.momentum());
146 Vector3 axis1 = pMeson.p3().unit();
147 double cMeson = e1z.dot(axis1);
148 _h_thy[ichi][imeson][1]->fill(cMeson);
149 double phi = atan2(e1y.dot(axis1),e1x.dot(axis1))+M_PI;
150 _h_thy[ichi][imeson][2]->fill(phi);
151 // moments to extract multipoles for chi_c2
152 if(ichi==0) {
153 double sGamma = sqrt(1.-sqr(cGamma));
154 double sMeson = sqrt(1.-sqr(cMeson));
155 double a3 = -3./sqrt(2.)*cos(phi)*sqr(sMeson)*2.*sMeson*cMeson*2.*sGamma*cGamma;
156 double a4 = sqrt(3.)*(3.*sqr(cMeson)-1.)*cos(phi)*2.*sMeson*cMeson*2.*sGamma*cGamma;
157 double a5 = sqrt(3./2.)*(3.*sqr(cMeson)-1.)*sqr(sGamma)*sqr(sMeson)*(2.*sqr(cos(phi))-1.);
158 _c[imeson][0]->fill(a3);
159 _c[imeson][1]->fill(a4);
160 _c[imeson][2]->fill(a5);
161 _c[2][0]->fill(a3);
162 _c[2][1]->fill(a4);
163 _c[2][2]->fill(a5);
164 }
165 // now fill experimental plots with cuts
166 if(cGammaCut>0.92 || (cGammaCut>0.8 && cGammaCut<0.86)) vetoEvent;
167 // cut on charged particles
168 if(abs(axis.dot(mPlus .p3().unit()))>0.93) vetoEvent;
169 if(abs(axis.dot(mMinus.p3().unit()))>0.93) vetoEvent;
170 // cut on angle of photon w.r.t. charged particles
171 if(abs(aGamma.dot(mPlus .p3().unit()))>cos20) vetoEvent;
172 if(abs(aGamma.dot(mMinus.p3().unit()))>cos20) vetoEvent;
173 // fill histos
174 _h_exp[ichi][imeson][0]->fill(cGamma);
175 _h_exp[ichi][imeson][1]->fill(cMeson);
176 _h_exp[ichi][imeson][2]->fill(phi);
177 }
178
179
180 /// Normalise histograms etc., after the run
181 void finalize() {
182 // first normalize the histograms
183 for (unsigned int ichi=0; ichi<2; ++ichi) {
184 for (unsigned int imeson=0; imeson<2; ++imeson) {
185 for (unsigned int iy=0; iy<3; ++iy) {
186 normalize(_h_thy[ichi][imeson][iy]);
187 normalize(_h_exp[ichi][imeson][iy]);
188 }
189 }
190 }
191 // extract the x and y values
192 double x, y, dx, dy;
193 for (unsigned int ix=0; ix<3; ++ix) {
194 Estimate0DPtr multX, multY;
195 book(multX, 1+ix, 1, 1);
196 book(multX, 1+ix, 1, 2);
197 *multX = (*_c[ix][0]/ *_c[ix][2]);
198 *multY = (*_c[ix][0]/ *_c[ix][2]);
199 x = multX->val();
200 dx = multX->errPos();
201 y = multY->val();
202 dy = multY->errPos();
203 }
204 // convert x and y to M1 and E2
205 double M1 = (3*sqrt(10) + sqrt(30)*x - 2*sqrt(15)*y)/(3.*(sqrt(2) + sqrt(6)*x + 2*sqrt(3)*y));
206 double E2 = (2*(2*sqrt(3) - 5*sqrt(2)*y - 2*sqrt(3)*sqr(y) + 4*x*(-1 + sqrt(6)*y)))/
207 ((sqrt(6) - 6*y)*(sqrt(2) + sqrt(6)*x + 2*sqrt(3)*y));
208 double e1 = (-4*sqrt(1.6666666666666667) + 4*sqrt(10)*y)/sqr(sqrt(2) + sqrt(6)*x + 2*sqrt(3)*y);
209 double e2 = (-4*sqrt(3.3333333333333335)*(2 + sqrt(3)*x))/sqr(sqrt(2) + sqrt(6)*x + 2*sqrt(3)*y);
210 double e3 = (20*(-sqrt(2) + y*(sqrt(3) + 3*sqrt(2)*y)))/((sqrt(6) - 6*y)*sqr(sqrt(2) + sqrt(6)*x + 2*sqrt(3)*y));
211 double e4 = (-10*(sqrt(6) - 3*sqrt(2)*x))/(3.*sqr(sqrt(2) + sqrt(6)*x + 2*sqrt(3)*y));
212 pair<double,double> dM1 = make_pair(sqrt(sqr(e1*dx)+sqr(e2*dy)),sqrt(sqr(e1*dx)+sqr(e2*dy)));
213 pair<double,double> dE2 = make_pair(sqrt(sqr(e3*dx)+sqr(e4*dy)),sqrt(sqr(e3*dx)+sqr(e4*dy)));
214 Estimate1DPtr multM1, multM2;
215 book(multM1, 4, 1, 1);
216 multM1->bin(1).set(M1, dM1);
217 book(multM1, 4, 1, 2);
218 multM2->bin(1).set(E2, dE2);
219 }
220
221 /// @}
222
223
224 /// @name Histograms
225 /// @{
226 Histo1DPtr _h_exp[2][2][3], _h_thy[2][2][3];
227 CounterPtr _c[3][3];
228 /// @}
229
230
231 };
232
233
234 RIVET_DECLARE_PLUGIN(BESIII_2011_I931195);
235
236}
|