rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

BESIII_2022_I2071715

Analysis of $J/\psi$ decays to $\Lambda^0\bar\Lambda^0$
Experiment: BESIII ()
Inspire ID: 2071715
Status: VALIDATED NOHEPDATA
No authors listed References: Beams: e- e+
Beam energies: (1.6, 1.6) GeV
    No run details listed

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      Estimate0DPtr _h_alpha0;
232      book(_h_alpha0,1,1,1);
233      _h_alpha0->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	Estimate0DPtr _h_alphaM;
282	book(_h_alphaM,1,3,1);
283	_h_alphaM->set(aM, make_pair(-aM_M , -aM_P));
284
285	Estimate0DPtr _h_alphaP;
286	book(_h_alphaP,1,4,1);
287	_h_alphaP->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	Estimate0DPtr _h_sin;
303	book(_h_sin,1,2,1);
304	_h_sin->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}