rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

BESIII_2023_I2689064

Kinematic distributions in $\Lambda_c^+\to \Lambda^0 \ell^+\nu_\ell$
Experiment: BESIII (BEPC)
Inspire ID: 2689064
Status: VALIDATED NOHEPDATA
Authors:
  • Peter Richardson
References:
  • Phys.Rev.D 108 (2023) 3, L031105
  • arXiv: 2306.02624
Beams: * *
Beam energies: ANY
Run details:
  • Any process producing Lambda_c+

Measurement of the kinematic distributions in $\Lambda_c^+\to \Lambda^0 \ell^+\nu_\ell$ by BES-III. N.B. The data were read from the paper, the decay for the asymmetries has been correct, however the kinematic distributions and may not have been corrected for acceptance, although the background given in the paper has been subtracted.

Source code: BESIII_2023_I2689064.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/DecayedParticles.hh"
  5#include "Rivet/Tools/HistoGroup.hh"
  6
  7
  8namespace Rivet {
  9
 10
 11  /// @brief  Lambda_c+ -> Lambda0 l+ nu_l
 12  class BESIII_2023_I2689064 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2023_I2689064);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      // Initialise and register projections
 26      UnstableParticles ufs = UnstableParticles(Cuts::pid==4122);
 27      declare(ufs, "UFS");
 28      DecayedParticles LAMBDAC(ufs);
 29      LAMBDAC.addStable(PID::PI0);
 30      LAMBDAC.addStable(PID::K0S);
 31      LAMBDAC.addStable(PID::ETA);
 32      LAMBDAC.addStable(PID::ETAPRIME);
 33      declare(LAMBDAC, "LAMBDAC");
 34      for (unsigned int ix=0; ix<2; ++ix) {
 35        book(_h_q[ix], 1, 1+ix, 1);
 36        for (unsigned int iy=0; iy<4; ++iy) {
 37          book(_h_kin[ix][iy], 3, ix+1, iy+1);
 38          if (iy<2) {
 39            book(_h_asym_l[ix][iy],"/TMP/h_asym_l_"+toString(ix+1)+"_"+toString(iy+1),refData(1,ix+1,2));
 40            book(_h_asym_p[ix][iy],"/TMP/h_asym_p_"+toString(ix+1)+"_"+toString(iy+1),refData(2,1,ix+1));
 41          }
 42        }
 43      }
 44      book(_b_ctheta, {0.,0.2,0.4,0.6,0.8,1.0,1.2,1.37});
 45      for (auto& b : _b_ctheta->bins()) {
 46        book(b, "/TMP/hLambda_"+toString(b.index()), 20, -1., 1.);
 47      }
 48      book(_c,"/TMP/n_lambda");
 49    }
 50
 51
 52    /// Perform the per-event analysis
 53    void analyze(const Event& event) {
 54      DecayedParticles LAMBDAC = apply<DecayedParticles>(event, "LAMBDAC");
 55      // loop over particles
 56      for (unsigned int ix=0; ix<LAMBDAC.decaying().size(); ++ix) {
 57        _c->fill();
 58        unsigned int il=0;
 59        if      (LAMBDAC.modeMatches(ix,4,mode1) ) il=0;
 60        else if (LAMBDAC.modeMatches(ix,4,mode2) ) il=1;
 61        else {
 62          continue;
 63        }
 64      	const Particle& pp = LAMBDAC.decayProducts()[ix].at(2212)[0];
 65      	const Particle& pim= LAMBDAC.decayProducts()[ix].at(-211)[0];
 66      	const Particle& ep = LAMBDAC.decayProducts()[ix].at( -11-2*il)[0];
 67      	const Particle& nue= LAMBDAC.decayProducts()[ix].at(  12+2*il)[0];
 68      	if (LAMBDAC.decaying()[ix].children(Cuts::pid==PID::LAMBDA).empty()) continue;
 69      	FourMomentum pLambda = pp.mom()+pim.mom();
 70       	FourMomentum qq = LAMBDAC.decaying()[ix].mom()-pLambda;
 71        double q2 = qq.mass2();
 72        _h_kin[il][0]->fill(q2);
 73        _h_q[il]->fill(q2);
 74        // boost momenta to LAMBDAC rest frame
 75        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(LAMBDAC.decaying()[ix].mom().betaVec());
 76        FourMomentum pLam = boost.transform(pLambda);
 77        Matrix3 ptoz(-pLam.p3().unit(), Vector3(0,0,1));
 78        boost.preMult(ptoz);
 79        // the momenta in frane to W along z
 80        FourMomentum pD  = boost.transform(LAMBDAC.decaying()[ix].mom());
 81        FourMomentum pP  = boost.transform(pp .mom());
 82        FourMomentum ppi = boost.transform(pim.mom());
 83        FourMomentum pe  = boost.transform(ep .mom());
 84        FourMomentum pnu = boost.transform(nue.mom());
 85        pLambda = pP+ppi;
 86        qq = pD-pLambda;
 87        LorentzTransform boostL = LorentzTransform::mkFrameTransformFromBeta(pLambda.betaVec());
 88        Vector3 axisP = boostL.transform(pP).p3().unit();
 89        const double cThetaP = axisP.dot(pLambda.p3().unit());
 90        _h_kin[il][1]->fill(cThetaP);
 91        _b_ctheta->fill(q2,cThetaP);
 92        if (cThetaP>0.) _h_asym_p[il][0]->fill(q2);
 93        else            _h_asym_p[il][1]->fill(q2);
 94        LorentzTransform boostW = LorentzTransform::mkFrameTransformFromBeta(    qq.betaVec());
 95        Vector3 axisE = boostW.transform(pe).p3().unit();
 96        double cThetaE = -axisE.dot(qq.p3().unit());
 97        _h_kin[il][2]->fill(cThetaE);
 98        if (cThetaE>0.) _h_asym_l[il][0]->fill(q2);
 99        else            _h_asym_l[il][1]->fill(q2);
100        axisP.setZ(0.);
101        axisE.setZ(0.);
102        double chi = atan2(axisE.cross(axisP).dot(qq.p3().unit()), axisE.dot(axisP));
103        _h_kin[il][3]->fill(chi);
104      }
105    }
106
107    pair<double,double> calcAlpha(Histo1DPtr hist) {
108      if (hist->numEntries()==0.)  return make_pair(0.,0.);
109      double sum1(0.),sum2(0.);
110      for (const auto& bin : hist->bins()) {
111        const double Oi = bin.sumW();
112        if (Oi==0.) continue;
113        const double ai = 0.5*(bin.xMax()-bin.xMin());
114        const double bi = 0.5*ai*(bin.xMax()+bin.xMin());
115        const double Ei = bin.errW();
116        sum1 += sqr(bi/Ei);
117        sum2 += bi/sqr(Ei)*(Oi-ai);
118      }
119      return make_pair(sum2/sum1,sqrt(1./sum1));
120    }
121
122    /// Normalise histograms etc., after the run
123    void finalize() {
124      const double tau = 201.5e-3; // lifetime in ps from PDG 2023
125      const double br = 0.641; // br for lambda -> p pi from PDG 2023
126      for (unsigned int ix=0; ix<2; ++ix) {
127        scale(_h_q[ix], 1.0/tau/br/ *_c);
128        normalize(_h_kin[ix], 1.0, false);
129      }
130      // asymmetries
131      Estimate1DPtr ap[2];
132      for (unsigned int ix=0; ix<2; ++ix) {
133        Estimate1DPtr as;
134        book(as, 1, 1+ix, 2);
135        asymm(_h_asym_l[ix][1], _h_asym_l[ix][0], as);
136        book(ap[ix], 2, 1, 1+ix);
137        asymm(_h_asym_p[ix][0], _h_asym_p[ix][1], ap[ix]);
138      }
139      Estimate1DPtr a2;
140      book(a2, 1, 3, 2);
141      divide(ap[1], ap[0], a2);
142      // ratios
143      Estimate1DPtr r1;
144      book(r1,1,3,1);
145      divide(_h_q[1], _h_q[0],r1);
146      // alpha parameters
147      Estimate1DPtr _h_alpha;
148      book(_h_alpha,2,2,1);
149      double lambda = 0.748; // PDG 2023
150      normalize(_b_ctheta);
151      for(unsigned int ix=0;ix<7;++ix) {
152	pair<double,double> alpha = calcAlpha(_b_ctheta->bin(ix+1));
153	alpha.first  /= lambda;
154	alpha.second /= lambda;
155	_h_alpha->bin(ix+1).set(alpha.first, alpha.second);
156      }
157    }
158
159    /// @}
160
161
162    /// @name Histograms
163    /// @{
164    Histo1DPtr _h_kin[2][4], _h_q[2] ,_h_asym_l[2][2], _h_asym_p[2][2];
165    Histo1DGroupPtr _b_ctheta;
166    CounterPtr _c;
167    const map<PdgId,unsigned int> mode1 = { { 2212,1}, {-211,1}, {-11,1}, { 12,1}};
168    const map<PdgId,unsigned int> mode2 = { { 2212,1}, {-211,1}, {-13,1}, { 14,1}};
169    /// @}
170
171
172  };
173
174
175  RIVET_DECLARE_PLUGIN(BESIII_2023_I2689064);
176
177}