Rivet analyses referenceBESIII_2022_I2033855Analysis of $\psi(2S)\to\gamma\chi_{c(0,2)}$ decays with $\chi_{c(0,2)}\to \Xi^-\bar{\Xi}^+/\Xi^0\bar{\Xi}^0$Experiment: BESIII (BEPC) Inspire ID: 2033855 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (1.8, 1.8) GeV Run details:
Analysis of the angular distribution of the photons and baryons produced in $\psi(2S)\to\gamma\chi_{c(0,2)}$ decays with $\chi_{c(0,2)}\to \Xi^-\bar{\Xi}^+/\Xi^0\bar{\Xi}^0$ 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 $\alpha$ results are fully corrected. Source code: BESIII_2022_I2033855.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 -> Xi Xibar
11 class BESIII_2022_I2033855 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2022_I2033855);
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 // book hists
28 for(unsigned int ix=0;ix<3;++ix)
29 for(unsigned int iy=0;iy<2;++iy)
30 book(_h[ix][iy],4+ix,1,1+iy);
31 }
32
33 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
34 for( const Particle &child : p.children()) {
35 if(child.children().empty()) {
36 nRes[child.pid()]-=1;
37 --ncount;
38 }
39 else
40 findChildren(child,nRes,ncount);
41 }
42 }
43
44 /// Perform the per-event analysis
45 void analyze(const Event& event) {
46 // get the axis, direction of incoming electron
47 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
48 Vector3 axis;
49 if(beams.first.pid()>0)
50 axis = beams.first .momentum().p3().unit();
51 else
52 axis = beams.second.momentum().p3().unit();
53 // types of final state particles
54 const FinalState& fs = apply<FinalState>(event, "FS");
55 map<long,int> nCount;
56 int ntotal(0);
57 for (const Particle& p : fs.particles()) {
58 nCount[p.pid()] += 1;
59 ++ntotal;
60 }
61 // loop over chi_c states
62 Particle chi;
63 bool matched = false;
64 const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
65 for (const Particle& p : ufs.particles()) {
66 if(p.children().empty()) continue;
67 map<long,int> nRes=nCount;
68 int ncount = ntotal;
69 findChildren(p,nRes,ncount);
70 if(ncount==1) {
71 matched = true;
72 for(auto const & val : nRes) {
73 if(val.first==PID::PHOTON) {
74 if(val.second!=1) {
75 matched = false;
76 break;
77 }
78 }
79 else if(val.second!=0) {
80 matched = false;
81 break;
82 }
83 }
84 if(matched) {
85 chi=p;
86 break;
87 }
88 }
89 }
90 if(!matched) vetoEvent;
91 // have chi_c find psi2S
92 if(chi.parents().empty() || chi.children().size()!=2 ||
93 chi.children()[0].pid() != -chi.children()[1].pid()) vetoEvent;
94 Particle psi2S = chi.parents()[0];
95 if(psi2S.pid()!=100443 || psi2S.children().size()!=2) vetoEvent;
96 // then the first photon
97 Particle gamma1;
98 if(psi2S.children()[0].pid()==PID::PHOTON)
99 gamma1 = psi2S.children()[0];
100 else if(psi2S.children()[1].pid()==PID::PHOTON)
101 gamma1 = psi2S.children()[1];
102 else
103 vetoEvent;
104 // now the decay products of the chi_c
105 Particle bPlus,bMinus;
106 bool foundBaryon=false;
107 for(unsigned int ix=0;ix<2;++ix) {
108 if(chi.children()[ix].pid()==PID::XIMINUS ||
109 chi.children()[ix].pid()==PID::XI0 ) {
110 foundBaryon=true;
111 bPlus=chi.children()[ix];
112 }
113 else if(chi.children()[ix].pid()==-PID::XIMINUS ||
114 chi.children()[ix].pid()==-PID::XI0 ) {
115 bMinus=chi.children()[ix];
116 }
117 }
118 if(!foundBaryon) vetoEvent;
119 // type chi state
120 unsigned int ichi = 0;
121 if(chi.pid()==20443) ichi = 1;
122 else if(chi.pid()==445) ichi = 2;
123 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(chi.momentum().betaVec());
124 Vector3 e1z = gamma1.momentum().p3().unit();
125 FourMomentum pBaryon = boost1.transform(bPlus.momentum());
126 Vector3 axis1 = pBaryon.p3().unit();
127 double cBaryon = e1z.dot(axis1);
128 if(bPlus.pid()==PID::XIMINUS)
129 _h[ichi][0]->fill(cBaryon);
130 else
131 _h[ichi][1]->fill(cBaryon);
132 }
133
134 pair<double,pair<double,double> > calcAlpha0(Histo1DPtr hist) {
135 if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
136 double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
137 double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
138 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
139 for (const auto& bin : hist->bins() ) {
140 double Oi = bin.sumW();
141 if(Oi==0.) continue;
142 double a = d*(bin.xMax() - bin.xMin());
143 double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
144 double Ei = bin.errW();
145 sum1 += a*Oi/sqr(Ei);
146 sum2 += b*Oi/sqr(Ei);
147 sum3 += sqr(a)/sqr(Ei);
148 sum4 += sqr(b)/sqr(Ei);
149 sum5 += a*b/sqr(Ei);
150 }
151 // calculate alpha
152 double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
153 // and error
154 double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
155 double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
156 double aa = sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
157 double dis = sqr(bb)-4.*aa*cc;
158 if(dis>0.) {
159 dis = sqrt(dis);
160 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
161 }
162 else {
163 return make_pair(alpha,make_pair(0.,0.));
164 }
165 }
166
167 /// Normalise histograms etc., after the run
168 void finalize() {
169 for(unsigned int ix=0;ix<3;++ix) {
170 for(unsigned int iy=0;iy<2;++iy) {
171 normalize(_h[ix][iy],1.,false);
172 pair<double,pair<double,double> > alpha0 = calcAlpha0(_h[ix][iy]);
173 Estimate1DPtr _h_alpha0;
174 book(_h_alpha0,1+ix,1,1+iy);
175 _h_alpha0->bin(1).set(alpha0.first, alpha0.second);
176 }
177 }
178 }
179
180 /// @}
181
182
183 /// @name Histograms
184 /// @{
185 Histo1DPtr _h[3][2];
186 /// @}
187
188
189 };
190
191
192 RIVET_DECLARE_PLUGIN(BESIII_2022_I2033855);
193
194}
|