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 (const auto& bin : hist->bins() ) {
108        double Oi = bin.sumW();
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.errW();
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      Estimate0DPtr h_alpha_lam;
143      book(h_alpha_lam, 5,1,1);
144      h_alpha_lam->set(alpha.first, alpha.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  RIVET_DECLARE_PLUGIN(BESIII_2012_I1121378);
161
162
163}