rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

BESIII_2012_I1121378

Analysis of $J/psi\to \Lambda^0\bar\Sigma^0+\text{c.c.}$
Experiment: BESIII (BEPC)
Inspire ID: 1121378
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Phys.Rev. D86 (2012) 032008
Beams: e- e+
Beam energies: (1.6, 1.6) GeV
Run details:
  • e+e- > J/psi

Analysis of the angular distribution of the baryons produced in $e^+e^-\to J/\psi\to\Lambda^0\bar\Sigma^0+\text{c.c.}$. Gives information about the decay and is useful for testing correlations in hadron decays.

Source code: BESIII_2012_I1121378.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 -> sigma0 lambda0 analysis
 11  class BESIII_2012_I1121378 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2012_I1121378);
 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_lam, 1, 1, 2);
 31      book(_h_bar, 1, 1, 1);
 32      book(_h_all, "/TMP/h_all",20,-1.,1.);
 33
 34    }
 35
 36    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 37      for(const Particle &child : p.children()) {
 38	if(child.children().empty()) {
 39	  nRes[child.pid()]-=1;
 40	  --ncount;
 41	}
 42	else
 43	  findChildren(child,nRes,ncount);
 44      }
 45    }
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49      // get the axis, direction of incoming electron
 50      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 51      Vector3 axis;
 52      if(beams.first.pid()>0)
 53	axis = beams.first .momentum().p3().unit();
 54      else
 55	axis = beams.second.momentum().p3().unit();
 56      // types of final state particles
 57      const FinalState& fs = apply<FinalState>(event, "FS");
 58      map<long,int> nCount;
 59      int ntotal(0);
 60      for (const Particle& p : fs.particles()) {
 61	nCount[p.pid()] += 1;
 62	++ntotal;
 63      }
 64      // loop over lambda0 and sigma0 baryons
 65      const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
 66      for (const Particle& p : ufs.particles(Cuts::abspid==3122)) {
 67	int sign = p.pid()/3122;
 68       	if(p.children().empty()) continue;
 69       	map<long,int> nRes=nCount;
 70       	int ncount = ntotal;
 71       	findChildren(p,nRes,ncount);
 72	bool matched=false;
 73	// check for Sigma0 or Sigma0bar
 74	for (const Particle& p2 : ufs.particles(Cuts::pid==sign*3212)) {
 75	  if(p2.children().empty()) continue;
 76	  map<long,int> nRes2=nRes;
 77	  int ncount2 = ncount;
 78	  findChildren(p2,nRes2,ncount2);
 79	  if(ncount2==0) {
 80	    matched = true;
 81	    for(auto const & val : nRes2) {
 82	      if(val.second!=0) {
 83		matched = false;
 84		break;
 85	      }
 86	    }
 87	    // fond baryon and antibaryon
 88	    if(matched) {
 89	      // calc cosine
 90	      double ctheta = p .momentum().p3().unit().dot(axis);
 91	      if(p.pid()==3122)
 92		_h_lam->fill(ctheta);
 93	      else
 94		_h_bar->fill(ctheta);
 95	      _h_all->fill(ctheta);
 96	      break;
 97	    }
 98	  }
 99	}
100	if(matched) break;
101      }
102    }
103
104    pair<double,pair<double,double> > calcAlpha(Histo1DPtr hist) {
105      if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
106      double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
107      for (auto bin : hist->bins() ) {
108       	double Oi = bin.area();
109	if(Oi==0.) continue;
110	double a =  1.5*(bin.xMax() - bin.xMin());
111	double b = 0.5*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
112       	double Ei = bin.areaErr();
113	sum1 +=   a*Oi/sqr(Ei);
114	sum2 +=   b*Oi/sqr(Ei);
115	sum3 += sqr(a)/sqr(Ei);
116	sum4 += sqr(b)/sqr(Ei);
117	sum5 +=    a*b/sqr(Ei);
118      }
119      // calculate alpha
120      double alpha = (-3*sum1 + 9*sum2 + sum3 - 3*sum5)/(sum1 - 3*sum2 + 3*sum4 - sum5);
121      // and error
122      double cc = -pow((sum3 + 9*sum4 - 6*sum5),3);
123      double bb = -2*sqr(sum3 + 9*sum4 - 6*sum5)*(sum1 - 3*sum2 + 3*sum4 - sum5);
124      double aa =  sqr(sum1 - 3*sum2 + 3*sum4 - sum5)*(-sum3 - 9*sum4 + sqr(sum1 - 3*sum2 + 3*sum4 - sum5) + 6*sum5);      
125      double dis = sqr(bb)-4.*aa*cc;
126      if(dis>0.) {
127	dis = sqrt(dis);
128	return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
129      }
130      else {
131	return make_pair(alpha,make_pair(0.,0.));
132      }
133    }
134
135    /// Normalise histograms etc., after the run
136    void finalize() {
137
138      normalize(_h_lam);
139      normalize(_h_bar);
140      normalize(_h_all);
141      pair<double,pair<double,double> > alpha = calcAlpha(_h_all);
142      Scatter2DPtr h_alpha_lam;
143      book(h_alpha_lam, 5,1,1);
144      h_alpha_lam->addPoint(0.5, alpha.first, make_pair(0.5,0.5), make_pair(alpha.second.first,alpha.second.second) );
145
146    }
147
148    //@}
149
150
151    /// @name Histograms
152    //@{
153    Histo1DPtr _h_lam,_h_bar,_h_all;
154    //@}
155
156
157  };
158
159
160  // The hook for the plugin system
161  RIVET_DECLARE_PLUGIN(BESIII_2012_I1121378);
162
163
164}