Rivet analyses referenceBESIII_2022_I1864775Analysis of $J/\psi$ decays to $\Xi^-\bar{\Xi}^+$Experiment: BESIII (BEPC) Inspire ID: 1864775 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (1.6, 1.6) GeV Run details:
Analysis of the angular distribution of the baryons, and decay products, produced in $e^+e^-\to J/\psi \to \Xi^-\bar{\Xi}^+$. The decays $\Xi^-\to\Lambda^0\pi^-$, $\Lambda^0\to p\pi^-$ and their charged conjugates are used. Gives information about the decay and is useful for testing correlations in hadron decays. N.B. the data is not corrected and should only be used qualatively. Source code: BESIII_2022_I1864775.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 -> Xi- Xibar+
11 class BESIII_2022_I1864775 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2022_I1864775);
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 // Book histograms
29 book(_h_T1, "T1",20,-1.,1.);
30 book(_h_T2, "T2",20,-1.,1.);
31 book(_h_T3, "T3",20,-1.,1.);
32 book(_h_T4, "T4",20,-1.,1.);
33 book(_h_T5, "T5",20,-1.,1.);
34 book(_h_cTheta,"cTheta",20,-1.,1.);
35 book(_h_clam[0], "cthetaP",20,-1,1);
36 book(_h_clam[1], "cthetaM",20,-1,1);
37 book(_wsum,"TMP/wsum");
38 }
39
40 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
41 for( const Particle &child : p.children()) {
42 if(child.children().empty()) {
43 nRes[child.pid()]-=1;
44 --ncount;
45 }
46 else
47 findChildren(child,nRes,ncount);
48 }
49 }
50
51 /// Perform the per-event analysis
52 void analyze(const Event& event) {
53 // get the axis, direction of incoming electron
54 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
55 Vector3 axis;
56 if(beams.first.pid()>0)
57 axis = beams.first .momentum().p3().unit();
58 else
59 axis = beams.second.momentum().p3().unit();
60 // types of final state particles
61 const FinalState& fs = apply<FinalState>(event, "FS");
62 map<long,int> nCount;
63 int ntotal(0);
64 for (const Particle& p : fs.particles()) {
65 nCount[p.pid()] += 1;
66 ++ntotal;
67 }
68 // loop over lambda0 baryons
69 const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
70 Particle Xi,XiBar;
71 bool matched(false);
72 for (const Particle& p : ufs.particles(Cuts::abspid==3312)) {
73 if(p.children().empty()) continue;
74 map<long,int> nRes=nCount;
75 int ncount = ntotal;
76 findChildren(p,nRes,ncount);
77 matched=false;
78 // check for antiparticle
79 for (const Particle& p2 : ufs.particles(Cuts::pid==-p.pid())) {
80 if(p2.children().empty()) continue;
81 map<long,int> nRes2=nRes;
82 int ncount2 = ncount;
83 findChildren(p2,nRes2,ncount2);
84 if(ncount2==0) {
85 matched = true;
86 for(auto const & val : nRes2) {
87 if(val.second!=0) {
88 matched = false;
89 break;
90 }
91 }
92 // found baryon and antibaryon
93 if(matched) {
94 if(p.pid()>0) {
95 Xi = p;
96 XiBar = p2;
97 }
98 else {
99 Xi = p2;
100 XiBar = p;
101 }
102 break;
103 }
104 }
105 }
106 if(matched) break;
107 }
108 if(!matched) vetoEvent;
109 // find the lambda and antilambda
110 Particle Lambda,LamBar;
111 if ( Xi.children()[0].pid() ==3122 )
112 Lambda = Xi.children()[0];
113 else if ( Xi.children()[1].pid() ==3122 )
114 Lambda = Xi.children()[1];
115 else vetoEvent;
116 if ( XiBar.children()[0].pid() ==-3122 )
117 LamBar = XiBar.children()[0];
118 else if ( XiBar.children()[1].pid() ==-3122 )
119 LamBar = XiBar.children()[1];
120 else vetoEvent;
121 // boost to the Xi rest frame
122 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Xi.momentum().betaVec());
123 Vector3 e1z = Xi.momentum().p3().unit();
124 Vector3 e1y = e1z.cross(axis).unit();
125 Vector3 e1x = e1y.cross(e1z).unit();
126 FourMomentum pLambda = boost1.transform(Lambda.momentum());
127 Vector3 axis1 = pLambda.p3().unit();
128 double n1x(e1x.dot(axis1)),n1y(e1y.dot(axis1)),n1z(e1z.dot(axis1));
129 // boost to the Xi bar
130 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(XiBar.momentum().betaVec());
131 FourMomentum pLamBar = boost2.transform(LamBar.momentum());
132 Vector3 axis2 = pLamBar.p3().unit();
133 double n2x(e1x.dot(axis2)),n2y(e1y.dot(axis2)),n2z(e1z.dot(axis2));
134 double cosX = axis.dot(Xi.momentum().p3().unit());
135 double sinX = sqrt(1.-sqr(cosX));
136 double T1 = sqr(sinX)*n1x*n2x+sqr(cosX)*n1z*n2z;
137 double T2 = -sinX*cosX*(n1x*n2z+n1z*n2x);
138 double T3 = -sinX*cosX*n1y;
139 double T4 = -sinX*cosX*n2y;
140 double T5 = n1z*n2z-sqr(sinX)*n1y*n2y;
141 _h_T1->fill(cosX,T1);
142 _h_T2->fill(cosX,T2);
143 _h_T3->fill(cosX,T3);
144 _h_T4->fill(cosX,T4);
145 _h_T5->fill(cosX,T5);
146 _h_cTheta->fill(cosX);
147 _wsum->fill();
148 // finally for the lambda decay
149 Particle proton;
150 bool found(false);
151 if(Lambda.children()[0].pid()==2212 &&
152 Lambda.children()[1].pid()==-211) {
153 proton = Lambda.children()[0];
154 found = true;
155 }
156 else if(Lambda.children()[1].pid()==2212 &&
157 Lambda.children()[0].pid()==-211) {
158 proton = Lambda.children()[1];
159 found = true;
160 }
161 if(found) {
162 // first boost to Xi rest frame
163 FourMomentum pproton = boost1.transform(proton.momentum());
164 LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pLambda.betaVec());
165 Vector3 axis = pLambda.p3().unit();
166 FourMomentum pp = boost3.transform(pproton);
167 // calculate angle
168 double cTheta = pp.p3().unit().dot(axis);
169 _h_clam[0]->fill(cTheta);
170 }
171 // finally for the anti lambda decay
172 found = false;
173 if(LamBar.children()[0].pid()==-2212 &&
174 LamBar.children()[1].pid()==211) {
175 proton = LamBar.children()[0];
176 found = true;
177 }
178 else if(LamBar.children()[1].pid()==-2212 &&
179 LamBar.children()[0].pid()==211) {
180 proton = LamBar.children()[1];
181 found = true;
182 }
183 if(found) {
184 // first boost to Xi rest frame
185 FourMomentum pproton = boost2.transform(proton.momentum());
186 LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pLamBar.betaVec());
187 Vector3 axis = pLamBar.p3().unit();
188 FourMomentum pp = boost3.transform(pproton);
189 // calculate angle
190 double cTheta = pp.p3().unit().dot(axis);
191 _h_clam[1]->fill(cTheta);
192 }
193 }
194
195 pair<double,pair<double,double> > calcAlpha0(Histo1DPtr hist) {
196 if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
197 double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
198 double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
199 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
200 for (const auto& bin : hist->bins() ) {
201 double Oi = bin.sumW();
202 if(Oi==0.) continue;
203 double a = d*bin.xWidth();
204 double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
205 double Ei = bin.errW();
206 sum1 += a*Oi/sqr(Ei);
207 sum2 += b*Oi/sqr(Ei);
208 sum3 += sqr(a)/sqr(Ei);
209 sum4 += sqr(b)/sqr(Ei);
210 sum5 += a*b/sqr(Ei);
211 }
212 // calculate alpha
213 double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
214 // and error
215 double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
216 double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
217 double aa = sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
218 double dis = sqr(bb)-4.*aa*cc;
219 if(dis>0.) {
220 dis = sqrt(dis);
221 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
222 }
223 else {
224 return make_pair(alpha,make_pair(0.,0.));
225 }
226 }
227
228 pair<double,double> calcCoeff(unsigned int imode,Histo1DPtr hist) {
229 if(hist->numEntries()==0.) return make_pair(0.,0.);
230 double sum1(0.),sum2(0.);
231 for (const auto& bin : hist->bins() ) {
232 double Oi = bin.sumW();
233 if(Oi==0.) continue;
234 double ai(0.),bi(0.);
235 if(imode==0) {
236 bi = (pow(1.-sqr(bin.xMin()),1.5) - pow(1.-sqr(bin.xMax()),1.5))/3.;
237 }
238 else if(imode>=2 && imode<=4) {
239 bi = ( pow(bin.xMin(),3)*( -5. + 3.*sqr(bin.xMin())) +
240 pow(bin.xMax(),3)*( 5. - 3.*sqr(bin.xMax())))/15.;
241 }
242 else
243 assert(false);
244 double Ei = bin.errW();
245 sum1 += sqr(bi/Ei);
246 sum2 += bi/sqr(Ei)*(Oi-ai);
247 }
248 return make_pair(sum2/sum1,sqrt(1./sum1));
249 }
250
251 pair<double,double> calcAlpha(Histo1DPtr hist) {
252 if(hist->numEntries()==0.) return make_pair(0.,0.);
253 double sum1(0.),sum2(0.);
254 for (const auto& bin : hist->bins() ) {
255 double Oi = bin.sumW();
256 if(Oi==0.) continue;
257 double ai = 0.5*bin.xWidth();
258 double bi = 0.5*ai*(bin.xMax()+bin.xMin());
259 double Ei = bin.errW();
260 sum1 += sqr(bi/Ei);
261 sum2 += bi/sqr(Ei)*(Oi-ai);
262 }
263 return make_pair(sum2/sum1,sqrt(1./sum1));
264 }
265
266 /// Normalise histograms etc., after the run
267 void finalize() {
268 normalize(_h_cTheta);
269 scale(_h_T1,1./ *_wsum);
270 scale(_h_T2,1./ *_wsum);
271 scale(_h_T3,1./ *_wsum);
272 scale(_h_T4,1./ *_wsum);
273 scale(_h_T5,1./ *_wsum);
274 normalize(_h_clam[0]);
275 normalize(_h_clam[1]);
276 // calculate alpha0
277 pair<double,pair<double,double> > alpha0 = calcAlpha0(_h_cTheta);
278 Estimate0DPtr _h_alpha0;
279 book(_h_alpha0,2,1,1);
280 _h_alpha0->set(alpha0.first, alpha0.second);
281 double s2 = -1. + sqr(alpha0.first);
282 double s3 = 3 + alpha0.first;
283 double s1 = sqr(s3);
284 // alpha- and alpha+ from proton data
285 pair<double,double> c_T2 = calcCoeff(2,_h_T2);
286 pair<double,double> c_T3 = calcCoeff(3,_h_T3);
287 pair<double,double> c_T4 = calcCoeff(4,_h_T4);
288 double s4 = sqr(c_T2.first);
289 double s5 = sqr(c_T3.first);
290 double s6 = sqr(c_T4.first);
291 double disc = s1*s5*s6*(-9.*s2*s4 + 4.*s1*s5*s6);
292 if(disc<0.) return;
293 disc = sqrt(disc);
294 double aM = -sqrt(-1./s2/s6*(2.*s1*s5*s6+disc));
295 double aP = c_T4.first/c_T3.first*aM;
296 double aM_M = (2*(alpha0.first*c_T4.first*alpha0.second.first + c_T4.second*s2)*(disc + 2*s1*s5*s6)
297 - 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) +
298 (disc*(- 9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
299 + 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
300 - s2*s3*c_T3.first*c_T4.second)* s4
301 + 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))
302 /(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)))/
303 (2.*pow(c_T4.first,3)*pow(s2,2)*sqrt(-((disc + 2*s1*s5*s6)/(s2*s6))));
304 double aM_P = (2*(alpha0.first*c_T4.first*alpha0.second.second + c_T4.second*s2)*(disc + 2*s1*s5*s6)
305 - 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) +
306 (disc*(- 9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
307 + 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
308 - s2*s3*c_T3.first*c_T4.second)* s4
309 + 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))
310 /(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)))/
311 (2.*pow(c_T4.first,3)*pow(s2,2)*sqrt(-((disc + 2*s1*s5*s6)/(s2*s6))));
312 double aP_M = (c_T4.first*sqrt(-((disc + 2*s1*s5*s6)/ (s2*s6)))*
313 (-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)
314 + (disc*(-9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
315 + 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
316 - s2*s3*c_T3.first*c_T4.second)* s4 +
317 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))/
318 (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)))/
319 (disc + 2*s1*s5*s6)))/(2.*pow(c_T3.first,2));
320 double aP_P = (c_T4.first*sqrt(-((disc + 2*s1*s5*s6)/ (s2*s6)))*
321 (-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)
322 + (disc*(-9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
323 + 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
324 - s2*s3*c_T3.first*c_T4.second)* s4 +
325 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))/
326 (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)))/
327 (disc + 2*s1*s5*s6)))/(2.*pow(c_T3.first,2));
328 Estimate0DPtr _h_alphaM;
329 book(_h_alphaM,2,1,3);
330 _h_alphaM->set(aM, make_pair(-aM_M , -aM_P));
331
332 Estimate0DPtr _h_alphaP;
333 book(_h_alphaP,2,1,5);
334 _h_alphaP->set(aP, make_pair(-aP_M , -aP_P));
335 // now for Delta
336 double sDelta = (-2.*(3. + alpha0.first)*c_T3.first)/(aM*sqrt(1 - sqr(alpha0.first)));
337 double cDelta = (-3*(3 + alpha0.first)*c_T2.first)/(aM*aP*sqrt(1 - sqr(alpha0.first)));
338 double Delta = asin(sDelta);
339 if(cDelta<0.) Delta = M_PI-Delta;
340 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)
341 - (-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T2.first*c_T3.first*c_T4.second)*disc)/
342 (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));
343 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)
344 - (-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T2.first*c_T3.first*c_T4.second)*disc)/
345 (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));
346 ds_P /= sqrt(1.-sqr(sDelta));
347 ds_M /= sqrt(1.-sqr(sDelta));
348 Estimate0DPtr _h_sin;
349 book(_h_sin,2,1,2);
350 _h_sin->set(Delta, make_pair(-ds_P, -ds_M));
351 // final the lambdas
352 // xibar+
353 Estimate0DPtr _h_lamP;
354 book(_h_lamP, 2,1,7);
355 pair<double,double> alpha = calcAlpha(_h_clam[0]);
356 alpha.second = sqrt(sqr(alpha.second/alpha.first)+ 0.5*(sqr(aM_M)+sqr(aM_P))/sqr(aM));
357 alpha.first /= aM;
358 alpha.second *=alpha.first;
359 _h_lamP->set(alpha.first, alpha.second);
360 // xi-
361 Estimate0DPtr _h_lamM;
362 book(_h_lamM, 2,1,8);
363 alpha = calcAlpha(_h_clam[1]);
364 alpha.second = sqrt(sqr(alpha.second/alpha.first)+ 0.5*(sqr(aP_M)+sqr(aP_P))/sqr(aP));
365 alpha.first /= aP;
366 alpha.second *=alpha.first;
367 _h_lamM->set(alpha.first, alpha.second);
368 }
369
370 /// @}
371
372
373 /// @name Histograms
374 /// @{
375 Histo1DPtr _h_T1,_h_T2,_h_T3,_h_T4,_h_T5;
376 Histo1DPtr _h_cTheta,_h_clam[2];
377 CounterPtr _wsum;
378 /// @}
379
380
381 };
382
383
384 RIVET_DECLARE_PLUGIN(BESIII_2022_I1864775);
385
386}
|
