rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

BESIII_2019_I1726357

Measurement of $e^+e^-\to\Lambda^0\bar{\Lambda}^0$ at 2.396 GeV
Experiment: BESIII (BEPC)
Inspire ID: 1726357
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Phys.Rev.Lett. 123 (2019) no.12, 122003
Beams: e+ e-
Beam energies: (1.2, 1.2) GeV
Run details:
  • e+e- to hadrons

Measurement of the cross section, angular distribution and polarization for $e^+e^-\to\Lambda^0\bar{\Lambda}^0$ at 2.396 GeV by BESIII.

Source code: BESIII_2019_I1726357.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 Lambda Lambdabar cross section
 11  class BESIII_2019_I1726357 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2019_I1726357);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23      // Initialise and register projections
 24      declare(Beam(), "Beams");
 25      declare(FinalState(), "FS");
 26      declare(UnstableParticles(), "UFS");
 27      // histograms
 28      book(_h_sigma ,1,1,1);
 29      book(_h_cTheta,2,1,1);
 30      book(_h_pol, {-1, -0.8, -0.6, -0.4, -0.2, 0., 0.2, 0.4, 0.6, 0.8, 1.0});
 31      for (auto& b : _h_pol->bins()) {
 32        book(b, "/TMP/h_pol_"+to_string(b.index()+1), 20, -1.0, 1.0);
 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      const FinalState& fs = apply<FinalState>(event, "FS");
 57      // total hadronic and muonic cross sections
 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      // find the Lambdas
 65      bool matched = false;
 66      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
 67      Particle Lambda;
 68      for(unsigned int ix=0;ix<ufs.particles().size();++ix) {
 69        const Particle& p1 = ufs.particles()[ix];
 70        if(abs(p1.pid())!=3122) continue;
 71        // check fs
 72        bool fs = true;
 73        for (const Particle & child : p1.children()) {
 74          if(child.pid()==p1.pid()) {
 75            fs = false;
 76            break;
 77          }
 78        }
 79        if(!fs) continue;
 80        // find the children
 81        map<long,int> nRes = nCount;
 82        int ncount = ntotal;
 83        findChildren(p1,nRes,ncount);
 84        for(unsigned int iy=ix+1;iy<ufs.particles().size();++iy) {
 85          matched=false;
 86          const Particle& p2 = ufs.particles()[iy];
 87          if(abs(p2.pid())!=3122) continue;
 88          // check fs
 89          bool fs = true;
 90          for (const Particle & child : p2.children()) {
 91            if(child.pid()==p2.pid()) {
 92              fs = false;
 93              break;
 94            }
 95          }
 96          if (!fs) continue;
 97          map<long,int> nRes2 = nRes;
 98          int ncount2 = ncount;
 99          findChildren(p2, nRes2, ncount2);
100          if (ncount2!=0) continue;
101          matched=true;
102          for (const auto& val : nRes2) {
103            if(val.second!=0) {
104              matched = false;
105              break;
106            }
107          }
108          if (matched) {
109            _h_sigma->fill(2.396);
110            if(p1.pid()==PID::LAMBDA) {
111              Lambda=p1;
112            }
113            else {
114              Lambda=p2;
115            }
116            break;
117          }
118        }
119        if (matched) break;
120      }
121      // now for the polarization measurements
122      if (matched) {
123        double cTheta = Lambda.momentum().p3().unit().dot(axis);
124        _h_cTheta->fill(cTheta);
125        Particle proton;
126        if(Lambda.children().size()!=2) return;
127        if(Lambda.children()[0].pid()==PID::PROTON &&
128           Lambda.children()[1].pid()==PID::PIMINUS)
129          proton = Lambda.children()[0];
130        else if(Lambda.children()[1].pid()==PID::PROTON &&
131                Lambda.children()[0].pid()==PID::PIMINUS)
132          proton = Lambda.children()[1];
133        else return;
134        LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Lambda.momentum().betaVec());
135        Vector3 axis1 = boost1.transform(proton.momentum()).p3().unit();
136        double cPhi = axis1.dot(Lambda.momentum().p3().unit());
137        _h_pol->fill(cTheta, cPhi);
138      }
139    }
140
141    pair<double,double> calcAlpha(Histo1DPtr hist) {
142      if(hist->numEntries()==0.) return make_pair(0.,0.);
143      double sum1(0.),sum2(0.);
144      for (const auto& bin : hist->bins()) {
145        double Oi = bin.sumW();
146        if(Oi==0.) continue;
147        double ai = 0.5*bin.xWidth();
148        double bi = 0.5*ai*(bin.xMax()+bin.xMin());
149        double Ei = bin.errW();
150        sum1 += sqr(bi/Ei);
151        sum2 += bi/sqr(Ei)*(Oi-ai);
152      }
153      return make_pair(sum2/sum1,sqrt(1./sum1));
154    }
155
156    /// Normalise histograms etc., after the run
157    void finalize() {
158      scale(_h_sigma, crossSection()/sumOfWeights()/picobarn);
159      normalize(_h_cTheta);
160      Estimate1DPtr _h_alpha;
161      book(_h_alpha, 3, 1, 1);
162      for (auto& b : _h_pol->bins()) {
163        normalize(b);
164        pair<double,double> alpha = calcAlpha(b);
165        _h_alpha->bin(b.index()).set(alpha.first, alpha.second);
166      }
167    }
168    /// @}
169
170
171    /// @name Histograms
172    /// @{
173    Histo1DPtr _h_sigma,_h_cTheta;
174    Histo1DGroupPtr _h_pol;
175    /// @}
176
177
178  };
179
180
181  RIVET_DECLARE_PLUGIN(BESIII_2019_I1726357);
182
183}