1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief CLEO charmed mesons and baryons from fragmentation
 10  ///
 11  /// @author Peter Richardson
 12  class CLEO_2004_I645209 : public Analysis {
 13  public:
 14
 15    RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_2004_I645209);
 16
 17    void init() {
 18      declare(Beam(), "Beams");
 19      declare(UnstableParticles(), "UFS");
 20
 21      // continuum cross sections
 22      book(_sigmaDPlus      ,1,1,1);
 23      book(_sigmaD0A        ,1,1,2);
 24      book(_sigmaD0B        ,1,1,3);
 25      book(_sigmaDStarPlusA ,1,1,4);
 26      book(_sigmaDStarPlusB ,1,1,5);
 27      book(_sigmaDStar0A    ,1,1,6);
 28      book(_sigmaDStar0B    ,1,1,7);
 29
 30       // histograms for continuum data
 31      book(_histXpDplus      ,2, 1, 1);
 32      book(_histXpD0A        ,3, 1, 1);
 33      book(_histXpD0B        ,4, 1, 1);
 34      book(_histXpDStarPlusA ,5, 1, 1);
 35      book(_histXpDStarPlusB ,6, 1, 1);
 36      book(_histXpDStar0A    ,7, 1, 1);
 37      book(_histXpDStar0B    ,8, 1, 1);
 38      book(_histXpTotal      ,9, 1, 1);
 39
 40    }
 41
 42
 43    void analyze(const Event& e) {
 44      // Loop through unstable FS particles and look for charmed mesons/baryons
 45      const UnstableParticles& ufs = apply<UnstableParticles>(e, "UFS");
 46
 47      const Beam beamproj = apply<Beam>(e, "Beams");
 48      const ParticlePair& beams = beamproj.beams();
 49      const FourMomentum mom_tot = beams.first.momentum() + beams.second.momentum();
 50      LorentzTransform cms_boost;
 51      if (mom_tot.p3().mod() > 1*MeV)
 52        cms_boost = LorentzTransform::mkFrameTransformFromBeta(mom_tot.betaVec());
 53      const double s = sqr(beamproj.sqrtS());
 54
 55      // Particle masses from PDGlive (accessed online 16. Nov. 2009).
 56      for (const Particle& p : ufs.particles()) {
 57
 58        double xp = 0.0;
 59        double mH2 = 0.0;
 60        // 3-momentum in CMS frame
 61        const double mom = cms_boost.transform(p.momentum()).vector3().mod();
 62
 63        const int pdgid = p.abspid();
 64        MSG_DEBUG("pdgID = " << pdgid << "  mom = " << mom);
 65        switch (pdgid) {
 66
 67        case 421:
 68          MSG_DEBUG("D0 found");
 69          mH2 = 3.47763; // 1.86484^2
 70          xp = mom/sqrt(s/4.0 - mH2);
 71          _sigmaD0A->fill(Ecm);
 72          _sigmaD0B->fill(Ecm);
 73          _histXpD0A->fill(xp);
 74          _histXpD0B->fill(xp);
 75          _histXpTotal->fill(xp);
 76          break;
 77
 78        case 411:
 79          MSG_DEBUG("D+ found");
 80          mH2 = 3.49547; // 1.86962^2
 81          xp = mom/sqrt(s/4.0 - mH2);
 82          _sigmaDPlus->fill(Ecm);
 83          _histXpDplus->fill(xp);
 84          _histXpTotal->fill(xp);
 85          break;
 86
 87        case 413:
 88          MSG_DEBUG("D*+ found");
 89          mH2 = 4.04119; // 2.01027^2
 90          xp = mom/sqrt(s/4.0 - mH2);
 91          _sigmaDStarPlusA->fill(Ecm);
 92          _sigmaDStarPlusB->fill(Ecm);
 93          _histXpDStarPlusA->fill(xp);
 94          _histXpDStarPlusB->fill(xp);
 95          _histXpTotal->fill(xp);
 96          break;
 97
 98        case 423:
 99          MSG_DEBUG("D*0 found");
100          mH2 = 4.02793; // 2.00697**2
101          xp = mom/sqrt(s/4.0 - mH2);
102          _sigmaDStar0A->fill(Ecm);
103          _sigmaDStar0B->fill(Ecm);
104          _histXpDStar0A->fill(xp);
105          _histXpDStar0B->fill(xp);
106          _histXpTotal->fill(xp);
107          break;
108        }
109
110      }
111    }
112
113
114    void finalize() {
115
116      const double sf = crossSection()/picobarn/sumOfWeights();
117      scale(_sigmaDPlus     , sf);
118      scale(_sigmaD0A       , sf);
119      scale(_sigmaD0B       , sf);
120      scale(_sigmaDStarPlusA, sf);
121      scale(_sigmaDStarPlusB, sf);
122      scale(_sigmaDStar0A   , sf);
123      scale(_sigmaDStar0B   , sf);
124
125      scale(_histXpDplus     , sf);
126      scale(_histXpD0A       , sf);
127      scale(_histXpD0B       , sf);
128      scale(_histXpDStarPlusA, sf);
129      scale(_histXpDStarPlusB, sf);
130      scale(_histXpDStar0A   , sf);
131      scale(_histXpDStar0B   , sf);
132      scale(_histXpTotal     , sf/4.);
133    }
134
135
136  private:
137
138    // Histograms for the continuum cross sections
139    BinnedHistoPtr<string> _sigmaDPlus;
140    BinnedHistoPtr<string> _sigmaD0A;
141    BinnedHistoPtr<string> _sigmaD0B;
142    BinnedHistoPtr<string> _sigmaDStarPlusA;
143    BinnedHistoPtr<string> _sigmaDStarPlusB;
144    BinnedHistoPtr<string> _sigmaDStar0A;
145    BinnedHistoPtr<string> _sigmaDStar0B;
146
147    // Histograms for continuum data
148    Histo1DPtr _histXpDplus;
149    Histo1DPtr _histXpD0A;
150    Histo1DPtr _histXpD0B;
151    Histo1DPtr _histXpDStarPlusA;
152    Histo1DPtr _histXpDStarPlusB;
153    Histo1DPtr _histXpDStar0A;
154    Histo1DPtr _histXpDStar0B;
155    Histo1DPtr _histXpTotal;
156    const string Ecm = "10.56";
157
158  };
159
160
161  RIVET_DECLARE_ALIASED_PLUGIN(CLEO_2004_I645209, CLEO_2004_S5809304);
162
163}