1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/DecayedParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Lambda_c+ -> Lambda0 e+ nu_e
 10  class BESIII_2022_I2127373 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2022_I2127373);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      UnstableParticles ufs = UnstableParticles(Cuts::pid==4122);
 25      declare(ufs, "UFS");
 26      DecayedParticles LAMBDAC(ufs);
 27      LAMBDAC.addStable(PID::PI0);
 28      LAMBDAC.addStable(PID::K0S);
 29      LAMBDAC.addStable(PID::ETA);
 30      LAMBDAC.addStable(PID::ETAPRIME);
 31      declare(LAMBDAC, "LAMBDAC");
 32      
 33      // Book histograms
 34      for(unsigned int ix=0;ix<4;++ix)
 35	book(_h[ix],1,1,1+ix);
 36    }
 37
 38
 39    /// Perform the per-event analysis
 40    void analyze(const Event& event) {
 41      static const map<PdgId,unsigned int> & mode = { { 2212,1}, {-211,1}, {-11,1}, { 12,1}};
 42      DecayedParticles LAMBDAC = apply<DecayedParticles>(event, "LAMBDAC");
 43      // loop over particles
 44      for(unsigned int ix=0;ix<LAMBDAC.decaying().size();++ix) {
 45      	if ( !LAMBDAC.modeMatches(ix,4,mode) ) continue;
 46	const Particle & pp = LAMBDAC.decayProducts()[ix].at(2212)[0];
 47	const Particle & pim= LAMBDAC.decayProducts()[ix].at(-211)[0];
 48	const Particle & ep = LAMBDAC.decayProducts()[ix].at( -11)[0];
 49	const Particle & nue= LAMBDAC.decayProducts()[ix].at(  12)[0];
 50	if(LAMBDAC.decaying()[ix].children(Cuts::pid==PID::LAMBDA).empty()) continue;
 51	FourMomentum pLambda = pp.momentum()+pim.momentum(); 
 52	FourMomentum qq = LAMBDAC.decaying()[ix].momentum()-pLambda;
 53	_h[0]->fill(qq.mass2());
 54	// boost momenta to LAMBDAC rest frame
 55	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(LAMBDAC.decaying()[ix].momentum().betaVec());
 56	FourMomentum pLam = boost.transform(pLambda);
 57	Matrix3 ptoz(-pLam.p3().unit(), Vector3(0,0,1));
 58	boost.preMult(ptoz);
 59	// the momenta in frane to W along z
 60	FourMomentum pD  = boost.transform(LAMBDAC.decaying()[ix].momentum());
 61	FourMomentum pP  = boost.transform(pp .momentum());
 62	FourMomentum ppi = boost.transform(pim.momentum());
 63	FourMomentum pe  = boost.transform(ep .momentum());
 64	FourMomentum pnu = boost.transform(nue.momentum());
 65	pLambda = pP+ppi;
 66	qq = pD-pLambda;
 67	LorentzTransform boostL = LorentzTransform::mkFrameTransformFromBeta(pLambda.betaVec());
 68	Vector3 axisP = boostL.transform(pP).p3().unit();
 69	_h[1]->fill(axisP.dot(pLambda.p3().unit()));
 70	LorentzTransform boostW = LorentzTransform::mkFrameTransformFromBeta(    qq.betaVec());
 71	Vector3 axisE = boostW.transform(pe).p3().unit();
 72	_h[2]->fill(-axisE.dot(qq.p3().unit()));
 73	axisP.setZ(0.);
 74	axisE.setZ(0.);
 75	double chi = atan2(axisE.cross(axisP).dot(qq.p3().unit()), axisE.dot(axisP));
 76	_h[3]->fill(chi);
 77      }
 78    }
 79
 80
 81    /// Normalise histograms etc., after the run
 82    void finalize() {
 83      for(unsigned int ix=0;ix<4;++ix)
 84	normalize(_h[ix]);
 85    }
 86
 87    /// @}
 88
 89
 90    /// @name Histograms
 91    /// @{
 92    Histo1DPtr _h[4];
 93    /// @}
 94
 95
 96  };
 97
 98
 99  RIVET_DECLARE_PLUGIN(BESIII_2022_I2127373);
100
101}