Rivet analyses referenceBESIII_2022_I2071715Analysis of $J/\psi$ decays to $\Lambda^0\bar\Lambda^0$Experiment: BESIII () Inspire ID: 2071715 Status: VALIDATED NOHEPDATA No authors listed References:
Beam energies: (1.6, 1.6) GeV
Analysis of the angular distribution of the baryons, and decay products, produced in $e^+e^-\to J/\psi \to \Lambda^0\bar\Lambda^0$. This is an updated version of BESIII_2019_I1691850 with higher statistics. Gives information about the decay and is useful for testing correlations in hadron decays. N.B. the moment data is not corrected and should only be used qualatively. Source code: BESIII_2022_I2071715.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 JPsi > Lambda, Lambdabar
11 class BESIII_2022_I2071715 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2022_I2071715);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23
24 // Initialise and register projections
25 declare(Beam(), "Beams");
26 declare(UnstableParticles(), "UFS");
27 declare(FinalState(), "FS");
28
29 // Book histograms
30 book(_h_T1, "T1",20,-1.,1.);
31 book(_h_T2, "T2",20,-1.,1.);
32 book(_h_T3, "T3",20,-1.,1.);
33 book(_h_T4, "T4",20,-1.,1.);
34 book(_h_T5, "T5",20,-1.,1.);
35
36 book(_h_cThetaL,"cThetaL",20,-1.,1.);
37
38 book(_h_mu, 2,1,1);
39 book(_wsum,"TMP/wsum");
40 }
41
42 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
43 for( const Particle &child : p.children()) {
44 if(child.children().empty()) {
45 nRes[child.pid()]-=1;
46 --ncount;
47 }
48 else
49 findChildren(child,nRes,ncount);
50 }
51 }
52
53 /// Perform the per-event analysis
54 void analyze(const Event& event) {
55 // get the axis, direction of incoming electron
56 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
57 Vector3 axis;
58 if(beams.first.pid()>0)
59 axis = beams.first .momentum().p3().unit();
60 else
61 axis = beams.second.momentum().p3().unit();
62 // types of final state particles
63 const FinalState& fs = apply<FinalState>(event, "FS");
64 map<long,int> nCount;
65 int ntotal(0);
66 for (const Particle& p : fs.particles()) {
67 nCount[p.pid()] += 1;
68 ++ntotal;
69 }
70 // loop over lambda0 baryons
71 const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
72 Particle Lambda,LamBar;
73 bool matched(false);
74 for (const Particle& p : ufs.particles(Cuts::abspid==3122)) {
75 if(p.children().empty()) continue;
76 map<long,int> nRes=nCount;
77 int ncount = ntotal;
78 findChildren(p,nRes,ncount);
79 matched=false;
80 // check for antiparticle
81 for (const Particle& p2 : ufs.particles(Cuts::pid==-p.pid())) {
82 if(p2.children().empty()) continue;
83 map<long,int> nRes2=nRes;
84 int ncount2 = ncount;
85 findChildren(p2,nRes2,ncount2);
86 if(ncount2==0) {
87 matched = true;
88 for(auto const & val : nRes2) {
89 if(val.second!=0) {
90 matched = false;
91 break;
92 }
93 }
94 // found baryon and antibaryon
95 if(matched) {
96 if(p.pid()>0) {
97 Lambda = p;
98 LamBar = p2;
99 }
100 else {
101 Lambda = p2;
102 LamBar = p;
103 }
104 break;
105 }
106 }
107 }
108 if(matched) break;
109 }
110 if(!matched) vetoEvent;
111 Particle proton;
112 matched = false;
113 for (const Particle & p : Lambda.children()) {
114 if(p.pid()==2212) {
115 matched=true;
116 proton=p;
117 }
118 else if(p.pid()==PID::PHOTON)
119 vetoEvent;
120 }
121 if(!matched) vetoEvent;
122 Particle baryon;
123 matched = false;
124 for (const Particle & p : LamBar.children()) {
125 if(p.pid()==-2212) {
126 baryon=p;
127 matched=true;
128 }
129 else if(p.pid()==PID::PHOTON)
130 vetoEvent;
131 }
132 if(!matched) vetoEvent;
133 // boost to the Lambda rest frame
134 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Lambda.momentum().betaVec());
135 Vector3 e1z = Lambda.momentum().p3().unit();
136 Vector3 e1y = e1z.cross(axis).unit();
137 Vector3 e1x = e1y.cross(e1z).unit();
138 Vector3 axis1 = boost1.transform(proton.momentum()).p3().unit();
139 double n1x(e1x.dot(axis1)),n1y(e1y.dot(axis1)),n1z(e1z.dot(axis1));
140 // boost to the Lambda bar
141 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(LamBar.momentum().betaVec());
142 Vector3 axis2 = boost2.transform(baryon.momentum()).p3().unit();
143 double n2x(e1x.dot(axis2)),n2y(e1y.dot(axis2)),n2z(e1z.dot(axis2));
144 double cosL = axis.dot(Lambda.momentum().p3().unit());
145 double sinL = sqrt(1.-sqr(cosL));
146 double T1 = sqr(sinL)*n1x*n2x+sqr(cosL)*n1z*n2z;
147 double T2 = -sinL*cosL*(n1x*n2z+n1z*n2x);
148 double T3 = -sinL*cosL*n1y;
149 double T4 = -sinL*cosL*n2y;
150 double T5 = n1z*n2z-sqr(sinL)*n1y*n2y;
151 double mu = n1y-n2y;
152 _h_T1->fill(cosL,T1);
153 _h_T2->fill(cosL,T2);
154 _h_T3->fill(cosL,T3);
155 _h_T4->fill(cosL,T4);
156 _h_T5->fill(cosL,T5);
157 _h_mu->fill(cosL,mu);
158 _wsum->fill();
159 _h_cThetaL->fill(cosL);
160 }
161
162 pair<double,pair<double,double> > calcAlpha0(Histo1DPtr hist) {
163 if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
164 double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
165 double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
166 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
167 for (const auto& bin : hist->bins() ) {
168 double Oi = bin.sumW();
169 if(Oi==0.) continue;
170 double a = d*(bin.xMax() - bin.xMin());
171 double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
172 double Ei = bin.errW();
173 sum1 += a*Oi/sqr(Ei);
174 sum2 += b*Oi/sqr(Ei);
175 sum3 += sqr(a)/sqr(Ei);
176 sum4 += sqr(b)/sqr(Ei);
177 sum5 += a*b/sqr(Ei);
178 }
179 // calculate alpha
180 double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
181 // and error
182 double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
183 double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
184 double aa = sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
185 double dis = sqr(bb)-4.*aa*cc;
186 if(dis>0.) {
187 dis = sqrt(dis);
188 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
189 }
190 else {
191 return make_pair(alpha,make_pair(0.,0.));
192 }
193 }
194
195 pair<double,double> calcCoeff(unsigned int imode,Histo1DPtr hist) {
196 if(hist->numEntries()==0.) return make_pair(0.,0.);
197 double sum1(0.),sum2(0.);
198 for (const auto& bin : hist->bins() ) {
199 double Oi = bin.sumW();
200 if(Oi==0.) continue;
201 double ai(0.),bi(0.);
202 if(imode==0) {
203 bi = (pow(1.-sqr(bin.xMin()),1.5) - pow(1.-sqr(bin.xMax()),1.5))/3.;
204 }
205 else if(imode>=2 && imode<=4) {
206 bi = ( pow(bin.xMin(),3)*( -5. + 3.*sqr(bin.xMin())) +
207 pow(bin.xMax(),3)*( 5. - 3.*sqr(bin.xMax())))/15.;
208 }
209 else {
210 assert(false);
211 }
212 double Ei = bin.errW();
213 sum1 += sqr(bi/Ei);
214 sum2 += bi/sqr(Ei)*(Oi-ai);
215 }
216 return make_pair(sum2/sum1,sqrt(1./sum1));
217 }
218
219 /// Normalise histograms etc., after the run
220 void finalize() {
221 normalize(_h_cThetaL);
222 scale(_h_T1,1./ *_wsum);
223 scale(_h_T2,1./ *_wsum);
224 scale(_h_T3,1./ *_wsum);
225 scale(_h_T4,1./ *_wsum);
226 scale(_h_T5,1./ *_wsum);
227 scale(_h_mu,100.*0.02/ *_wsum);
228
229 // calculate alpha0
230 pair<double,pair<double,double> > alpha0 = calcAlpha0(_h_cThetaL);
231 Estimate1DPtr _h_alpha0;
232 book(_h_alpha0,1,1,1);
233 _h_alpha0->bin(1).set(alpha0.first, alpha0.second);
234 double s2 = -1. + sqr(alpha0.first);
235 double s3 = 3 + alpha0.first;
236 double s1 = sqr(s3);
237 // alpha- and alpha+ from proton data
238 pair<double,double> c_T2 = calcCoeff(2,_h_T2);
239 pair<double,double> c_T3 = calcCoeff(3,_h_T3);
240 pair<double,double> c_T4 = calcCoeff(4,_h_T4);
241 double s4 = sqr(c_T2.first);
242 double s5 = sqr(c_T3.first);
243 double s6 = sqr(c_T4.first);
244 double disc = s1*s5*s6*(-9.*s2*s4 + 4.*s1*s5*s6);
245 if(disc>=0.) {
246 disc = sqrt(disc);
247 double aM = sqrt(-1./s2/s6*(2.*s1*s5*s6+disc));
248 double aP = c_T4.first/c_T3.first*aM;
249 double aM_M = (2*(alpha0.first*c_T4.first*alpha0.second.first + c_T4.second*s2)*(disc + 2*s1*s5*s6)
250 - c_T4.first*s2*(4*s3*c_T3.first*c_T4.first*(c_T3.first*c_T4.first*alpha0.second.first +s3*c_T4.first*c_T3.second +s3*c_T3.first*c_T4.second) +
251 (disc*(- 9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
252 + 9*((1 - alpha0.first*(3 + 2*alpha0.first))* c_T3.first*c_T4.first*alpha0.second.first - s2*s3*c_T4.first*c_T3.second
253 - s2*s3*c_T3.first*c_T4.second)* s4
254 + 8*(c_T3.first*c_T4.first*alpha0.second.first + s3*c_T4.first*c_T3.second + s3*c_T3.first*c_T4.second)* s1*s5*s6))
255 /(4*pow(3 + alpha0.first,3)*pow(c_T3.first,3)*pow(c_T4.first,3) -9*s2*s3*c_T3.first*c_T4.first*s4)))/
256 (2.*pow(c_T4.first,3)*pow(s2,2)*sqrt(-((disc + 2*s1*s5*s6)/(s2*s6))));
257 double aM_P = (2*(alpha0.first*c_T4.first*alpha0.second.second + c_T4.second*s2)*(disc + 2*s1*s5*s6)
258 - c_T4.first*s2*(4*s3*c_T3.first*c_T4.first*(c_T3.first*c_T4.first*alpha0.second.second +s3*c_T4.first*c_T3.second +s3*c_T3.first*c_T4.second) +
259 (disc*(- 9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
260 + 9*((1 - alpha0.first*(3 + 2*alpha0.first))* c_T3.first*c_T4.first*alpha0.second.second - s2*s3*c_T4.first*c_T3.second
261 - s2*s3*c_T3.first*c_T4.second)* s4
262 + 8*(c_T3.first*c_T4.first*alpha0.second.second + s3*c_T4.first*c_T3.second + s3*c_T3.first*c_T4.second)* s1*s5*s6))
263 /(4*pow(3 + alpha0.first,3)*pow(c_T3.first,3)*pow(c_T4.first,3) -9*s2*s3*c_T3.first*c_T4.first*s4)))/
264 (2.*pow(c_T4.first,3)*pow(s2,2)*sqrt(-((disc + 2*s1*s5*s6)/(s2*s6))));
265 double aP_M = (c_T4.first*sqrt(-((disc + 2*s1*s5*s6)/ (s2*s6)))*
266 (-2*c_T3.second - (2*alpha0.first*c_T3.first*alpha0.second.first)/s2 + (c_T3.first*(4*s3*c_T3.first*c_T4.first*(c_T3.first*c_T4.first*alpha0.second.first + s3*c_T4.first*c_T3.second + s3*c_T3.first*c_T4.second)
267 + (disc*(-9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
268 + 9*((1 - alpha0.first*(3 + 2*alpha0.first))* c_T3.first*c_T4.first*alpha0.second.first - s2*s3*c_T4.first*c_T3.second
269 - s2*s3*c_T3.first*c_T4.second)* s4 +
270 8*(c_T3.first*c_T4.first*alpha0.second.first + s3*c_T4.first*c_T3.second + s3*c_T3.first*c_T4.second)* s1*s5*s6))/
271 (4* pow(3 + alpha0.first,3)* pow(c_T3.first,3)* pow(c_T4.first,3) - 9*s2*s3*c_T3.first*c_T4.first*s4)))/
272 (disc + 2*s1*s5*s6)))/(2.*pow(c_T3.first,2));
273 double aP_P = (c_T4.first*sqrt(-((disc + 2*s1*s5*s6)/ (s2*s6)))*
274 (-2*c_T3.second - (2*alpha0.first*c_T3.first*alpha0.second.second)/s2 + (c_T3.first*(4*s3*c_T3.first*c_T4.first*(c_T3.first*c_T4.first*alpha0.second.second + s3*c_T4.first*c_T3.second + s3*c_T3.first*c_T4.second)
275 + (disc*(-9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
276 + 9*((1 - alpha0.first*(3 + 2*alpha0.first))* c_T3.first*c_T4.first*alpha0.second.second - s2*s3*c_T4.first*c_T3.second
277 - s2*s3*c_T3.first*c_T4.second)* s4 +
278 8*(c_T3.first*c_T4.first*alpha0.second.second + s3*c_T4.first*c_T3.second + s3*c_T3.first*c_T4.second)* s1*s5*s6))/
279 (4* pow(3 + alpha0.first,3)* pow(c_T3.first,3)* pow(c_T4.first,3) - 9*s2*s3*c_T3.first*c_T4.first*s4)))/
280 (disc + 2*s1*s5*s6)))/(2.*pow(c_T3.first,2));
281 Estimate1DPtr _h_alphaM;
282 book(_h_alphaM,1,3,1);
283 _h_alphaM->bin(1).set(aM, make_pair(-aM_M , -aM_P));
284
285 Estimate1DPtr _h_alphaP;
286 book(_h_alphaP,1,4,1);
287 _h_alphaP->bin(1).set(aP, make_pair(-aP_M , -aP_P));
288 // now for Delta
289 double sDelta = (-2.*(3. + alpha0.first)*c_T3.first)/(aM*sqrt(1 - sqr(alpha0.first)));
290 double cDelta = (-3*(3 + alpha0.first)*c_T2.first)/(aM*aP*sqrt(1 - sqr(alpha0.first)));
291
292 double Delta = asin(sDelta);
293 if(cDelta<0.) Delta = M_PI-Delta;
294 double ds_P = (-9*c_T2.first*((-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T3.first*c_T4.first*c_T2.second + c_T2.first*c_T4.first*(c_T3.first*(alpha0.second.first + 3*alpha0.first*alpha0.second.first) -(-1 + alpha0.first)*(1 + alpha0.first)*(3 + alpha0.first)*c_T3.second)
295 - (-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T2.first*c_T3.first*c_T4.second)*disc)/
296 (pow(1 - pow(alpha0.first,2),1.5)*pow(c_T4.first,3)*pow(-((disc + 2*s1*s5*s6)/ (s2*s6)),1.5)*(-9*s2*s4 + 4*s1*s5*s6));
297 double ds_M = (-9*c_T2.first*((-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T3.first*c_T4.first*c_T2.second + c_T2.first*c_T4.first*(c_T3.first*(alpha0.second.second + 3*alpha0.first*alpha0.second.second) -(-1 + alpha0.first)*(1 + alpha0.first)*(3 + alpha0.first)*c_T3.second)
298 - (-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T2.first*c_T3.first*c_T4.second)*disc)/
299 (pow(1 - pow(alpha0.first,2),1.5)*pow(c_T4.first,3)*pow(-((disc + 2*s1*s5*s6)/ (s2*s6)),1.5)*(-9*s2*s4 + 4*s1*s5*s6));
300 ds_P /= sqrt(1.-sqr(sDelta));
301 ds_M /= sqrt(1.-sqr(sDelta));
302 Estimate1DPtr _h_sin;
303 book(_h_sin,1,2,1);
304 _h_sin->bin(1).set(Delta, make_pair(-ds_P, -ds_M));
305 }
306 }
307
308 /// @}
309
310
311 /// @name Histograms
312 /// @{
313 Histo1DPtr _h_T1,_h_T2,_h_T3,_h_T4,_h_T5;
314 Histo1DPtr _h_cThetaL;
315 Histo1DPtr _h_mu;
316 CounterPtr _wsum;
317 /// @}
318
319
320 };
321
322
323 RIVET_DECLARE_PLUGIN(BESIII_2022_I2071715);
324
325}
|