1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Babar charmonium spectra
  9  /// @author Peter Richardson
 10  class BABAR_2003_I593379 : public Analysis {
 11  public:
 12
 13    BABAR_2003_I593379()
 14      : Analysis("BABAR_2003_I593379")
 15    { }
 16
 17
 18    void analyze(const Event& e) {
 19      // Find the charmonia
 20      Particles upsilons;
 21      // First in unstable final state
 22      const UnstableParticles& ufs = apply<UnstableParticles>(e, "UFS");
 23      for (const Particle& p : ufs.particles())
 24        if (p.pid() == 300553) upsilons.push_back(p);
 25      // Then in whole event if fails
 26      if (upsilons.empty()) {
 27        for(ConstGenParticlePtr p: HepMCUtils::particles(e.genEvent())) {
 28          if (p->pdg_id() != 300553) continue;
 29          ConstGenVertexPtr pv = p->production_vertex();
 30          bool passed = true;
 31          if (pv) {
 32            for(ConstGenParticlePtr pp: HepMCUtils::particles(pv, Relatives::PARENTS)){
 33              if ( p->pdg_id() == pp->pdg_id() ) {
 34                passed = false;
 35                break;
 36              }
 37            }
 38          }
 39          if (passed) upsilons.push_back(Particle(p));
 40        }
 41      }
 42
 43      // Find upsilons
 44      for (const Particle& p : upsilons) {
 45        _weightSum->fill();
 46        // Find the charmonium resonances
 47        /// @todo Use Rivet::Particles
 48        vector<ConstGenParticlePtr> allJpsi, primaryJpsi, Psiprime, all_chi_c1, all_chi_c2, primary_chi_c1, primary_chi_c2;
 49        findDecayProducts(p.genParticle(), allJpsi, primaryJpsi, Psiprime,
 50                          all_chi_c1, all_chi_c2, primary_chi_c1, primary_chi_c2);
 51        const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(p.mom().betaVec());
 52        for (size_t i = 0; i < allJpsi.size(); i++) {
 53          const double pcm = cms_boost.transform(FourMomentum(allJpsi[i]->momentum())).p();
 54          _hist_all_Jpsi->fill(pcm);
 55        }
 56        _mult_JPsi->fill(double(allJpsi.size()));
 57        for (size_t i = 0; i < primaryJpsi.size(); i++) {
 58          const double pcm = cms_boost.transform(FourMomentum(primaryJpsi[i]->momentum())).p();
 59          _hist_primary_Jpsi->fill(pcm);
 60        }
 61        _mult_JPsi_direct->fill(double(primaryJpsi.size()));
 62        for (size_t i=0; i<Psiprime.size(); i++) {
 63          const double pcm = cms_boost.transform(FourMomentum(Psiprime[i]->momentum())).p();
 64          _hist_Psi_prime->fill(pcm);
 65        }
 66        _mult_Psi2S->fill(double(Psiprime.size()));
 67        for (size_t i = 0; i < all_chi_c1.size(); i++) {
 68          const double pcm = cms_boost.transform(FourMomentum(all_chi_c1[i]->momentum())).p();
 69          _hist_chi_c1->fill(pcm);
 70        }
 71        _mult_chi_c1->fill(double(all_chi_c1.size()));
 72        _mult_chi_c1_direct->fill(double(primary_chi_c1.size()));
 73        for (size_t i = 0; i < all_chi_c2.size(); i++) {
 74          const double pcm = cms_boost.transform(FourMomentum(all_chi_c2[i]->momentum())).p();
 75          _hist_chi_c2->fill(pcm);
 76        }
 77        _mult_chi_c2->fill(double(all_chi_c2.size()));
 78        _mult_chi_c2_direct->fill(double(primary_chi_c2.size()));
 79      }
 80    } // analyze
 81
 82
 83    void finalize() {
 84      scale(_hist_all_Jpsi     , 0.5*0.1 / *_weightSum);
 85      scale(_hist_chi_c1       , 0.5*0.1 / *_weightSum);
 86      scale(_hist_chi_c2       , 0.5*0.1 / *_weightSum);
 87      scale(_hist_Psi_prime    , 0.5*0.1 / *_weightSum);
 88      scale(_hist_primary_Jpsi , 0.5*0.1 / *_weightSum);
 89      scale(_mult_JPsi         , 0.5*100. / *_weightSum);
 90      scale(_mult_JPsi_direct  , 0.5*100. / *_weightSum);
 91      scale(_mult_chi_c1       , 0.5*100. / *_weightSum);
 92      scale(_mult_chi_c1_direct, 0.5*100. / *_weightSum);
 93      scale(_mult_chi_c2       , 0.5*100. / *_weightSum);
 94      scale(_mult_chi_c2_direct, 0.5*100. / *_weightSum);
 95      scale(_mult_Psi2S        , 0.5*100. / *_weightSum);
 96    } // finalize
 97
 98
 99    void init() {
100      declare(UnstableParticles(), "UFS");
101
102      book(_mult_JPsi          ,1, 1, 1);
103      book(_mult_JPsi_direct   ,1, 1, 2);
104      book(_mult_chi_c1        ,1, 1, 3);
105      book(_mult_chi_c1_direct ,1, 1, 4);
106      book(_mult_chi_c2        ,1, 1, 5);
107      book(_mult_chi_c2_direct ,1, 1, 6);
108      book(_mult_Psi2S         ,1, 1, 7);
109      book(_hist_all_Jpsi      ,2, 1, 1);
110      book(_hist_chi_c1        ,3, 1, 1);
111      book(_hist_chi_c2        ,3, 1, 2);
112      book(_hist_Psi_prime     ,4, 1, 1);
113      book(_hist_primary_Jpsi  ,6, 1, 1);
114
115      book(_weightSum, "TMP/weightSum");
116    } // init
117
118  private:
119
120    /// @{
121    // count of weights
122    CounterPtr _weightSum;
123    /// Histograms
124    Histo1DPtr _hist_all_Jpsi;
125    Histo1DPtr _hist_chi_c1;
126    Histo1DPtr _hist_chi_c2;
127    Histo1DPtr _hist_Psi_prime;
128    Histo1DPtr _hist_primary_Jpsi;
129
130    CounterPtr _mult_JPsi;
131    CounterPtr _mult_JPsi_direct;
132    CounterPtr _mult_chi_c1;
133    CounterPtr _mult_chi_c1_direct;
134    CounterPtr _mult_chi_c2;
135    CounterPtr _mult_chi_c2_direct;
136    CounterPtr _mult_Psi2S;
137    /// @}
138
139    void findDecayProducts(ConstGenParticlePtr p,
140                           vector<ConstGenParticlePtr>& allJpsi,
141                           vector<ConstGenParticlePtr>& primaryJpsi,
142                           vector<ConstGenParticlePtr>& Psiprime,
143                           vector<ConstGenParticlePtr>& all_chi_c1, vector<ConstGenParticlePtr>& all_chi_c2,
144                           vector<ConstGenParticlePtr>& primary_chi_c1, vector<ConstGenParticlePtr>& primary_chi_c2) {
145      ConstGenVertexPtr dv = p->end_vertex();
146      bool isOnium = false;
147      /// @todo Use better looping
148      for (ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::PARENTS)){
149        int id = pp->pdg_id();
150        id = id%1000;
151        id -= id%10;
152        id /= 10;
153        if (id==44) isOnium = true;
154      }
155      /// @todo Use better looping
156      for (ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
157        int id = pp->pdg_id();
158        if (id==100443) {
159          Psiprime.push_back(pp);
160        }
161        else if (id==20443) {
162          all_chi_c1.push_back(pp);
163          if (!isOnium) primary_chi_c1.push_back(pp);
164        }
165        else if (id==445) {
166          all_chi_c2.push_back(pp);
167          if (!isOnium) primary_chi_c2.push_back(pp);
168        }
169        else if (id==443) {
170          allJpsi.push_back(pp);
171          if (!isOnium) primaryJpsi.push_back(pp);
172        }
173        if (pp->end_vertex()) {
174          findDecayProducts(pp, allJpsi, primaryJpsi, Psiprime, all_chi_c1, all_chi_c2, primary_chi_c1, primary_chi_c2);
175        }
176      }
177    }
178
179  };
180
181
182  RIVET_DECLARE_PLUGIN(BABAR_2003_I593379);
183
184}