1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief BaBar pion, kaon and proton production in the continuum
 10  /// @author Peter Richardson
 11  class BABAR_2013_I1238276 : public Analysis {
 12  public:
 13
 14    BABAR_2013_I1238276()
 15      : Analysis("BABAR_2013_I1238276")
 16    { }
 17
 18
 19    void init() {
 20      declare(Beam(), "Beams");
 21      declare(ChargedFinalState(), "FS");
 22
 23      book(_histPion_no_dec   ,1,1,1);
 24      book(_histKaon_no_dec   ,1,1,2);
 25      book(_histProton_no_dec ,1,1,3);
 26      book(_histPion_dec      ,2,1,1);
 27      book(_histKaon_dec      ,2,1,2);
 28      book(_histProton_dec    ,2,1,3);
 29    }
 30
 31
 32    void analyze(const Event& e) {
 33      // Loop through charged FS particles and look for charmed mesons/baryons
 34      const ChargedFinalState& fs = apply<ChargedFinalState>(e, "FS");
 35
 36      const Beam beamproj = apply<Beam>(e, "Beams");
 37      const ParticlePair& beams = beamproj.beams();
 38      const FourMomentum mom_tot = beams.first.momentum() + beams.second.momentum();
 39      const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(mom_tot.betaVec());
 40      MSG_DEBUG("CMS Energy sqrt s = " << beamproj.sqrtS());
 41
 42      for (const Particle& p : fs.particles()) {
 43        // check if prompt or not
 44        ConstGenParticlePtr pmother = p.genParticle();
 45        ConstGenVertexPtr ivertex = pmother->production_vertex();
 46        bool prompt = true;
 47        while (ivertex) {
 48          vector<ConstGenParticlePtr> inparts = HepMCUtils::particles(ivertex, Relatives::PARENTS);
 49          int n_inparts = inparts.size();
 50          if (n_inparts < 1) break;
 51          pmother = inparts[0]; // first mother particle
 52          int mother_pid = abs(pmother->pdg_id());
 53          if (mother_pid==PID::K0S || mother_pid==PID::LAMBDA) {
 54            prompt = false;
 55            break;
 56          }
 57          else if (mother_pid<6) {
 58            break;
 59          }
 60          ivertex = pmother->production_vertex();
 61        }
 62
 63        // momentum in CMS frame
 64        const double mom = cms_boost.transform(p.momentum()).vector3().mod();
 65        const int PdgId = p.abspid();
 66        MSG_DEBUG("pdgID = " << PdgId << " Momentum = " << mom);
 67        switch (PdgId) {
 68        case PID::PIPLUS:
 69          if(prompt) _histPion_no_dec->fill(mom);
 70          _histPion_dec   ->fill(mom);
 71          break;
 72        case PID::KPLUS:
 73          if(prompt) _histKaon_no_dec->fill(mom);
 74          _histKaon_dec   ->fill(mom);
 75          break;
 76        case PID::PROTON:
 77          if(prompt) _histProton_no_dec->fill(mom);
 78          _histProton_dec   ->fill(mom);
 79        default :
 80          break;
 81        }
 82      }
 83    }
 84
 85
 86    void finalize() {
 87      scale(_histPion_no_dec  ,1./sumOfWeights());
 88      scale(_histKaon_no_dec  ,1./sumOfWeights());
 89      scale(_histProton_no_dec,1./sumOfWeights());
 90      scale(_histPion_dec     ,1./sumOfWeights());
 91      scale(_histKaon_dec     ,1./sumOfWeights());
 92      scale(_histProton_dec   ,1./sumOfWeights());
 93    }
 94
 95
 96  private:
 97
 98    /// @{
 99    // Histograms for continuum data (sqrt(s) = 10.52 GeV)
100    // no K_S and Lambda decays
101    Histo1DPtr _histPion_no_dec;
102    Histo1DPtr _histKaon_no_dec;
103    Histo1DPtr _histProton_no_dec;
104    // including decays
105    Histo1DPtr _histPion_dec;
106    Histo1DPtr _histKaon_dec;
107    Histo1DPtr _histProton_dec;
108    /// @}
109
110  };
111
112
113  RIVET_DECLARE_PLUGIN(BABAR_2013_I1238276);
114
115}