rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

BESIII_2022_I1864775

Analysis of $J/\psi$ decays to $\Xi^-\bar{\Xi}^+$
Experiment: BESIII (BEPC)
Inspire ID: 1864775
Status: VALIDATED NOHEPDATA
Authors:
  • Peter Richardson
References:
  • Nature 606 (2022) 7912, 64-69
Beams: e- e+
Beam energies: (1.6, 1.6) GeV
Run details:
  • e+e- > J/psi

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}