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 BELLE pi0 spectrum at Upsilon(4S)
 10  ///
 11  /// @author Peter Richardson
 12  class BELLE_2001_I554520 : public Analysis {
 13  public:
 14
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2001_I554520);
 16
 17
 18    void init() {
 19      declare(UnstableParticles(), "UFS");
 20      book(_histdSigDp ,1, 1, 1); // spectrum
 21      book(_histMult   ,2, 1, 1); // multiplicity
 22      book(_weightSum, "TMP/weightSum");
 23    }
 24
 25
 26    void analyze(const Event& e) {
 27      // Find the upsilons
 28      Particles upsilons;
 29      // First in unstable final state
 30      const UnstableParticles& ufs = apply<UnstableParticles>(e, "UFS");
 31      for (const Particle& p : ufs.particles())
 32        if (p.pid()==300553) upsilons.push_back(p);
 33      // Then in whole event if fails
 34      if (upsilons.empty()) {
 35        for(ConstGenParticlePtr p: HepMCUtils::particles(e.genEvent())) {
 36          if (p->pdg_id() != 300553) continue;
 37          ConstGenVertexPtr pv = p->production_vertex();
 38          bool passed = true;
 39          if (pv) {
 40            for (ConstGenParticlePtr pp: HepMCUtils::particles(pv, Relatives::PARENTS)){
 41              if ( p->pdg_id() == pp->pdg_id() ) {
 42                passed = false;
 43                break;
 44              }
 45            }
 46          }
 47          if (passed) upsilons.push_back(Particle(p));
 48        }
 49      }
 50
 51      // Find upsilons
 52      for (const Particle& p : upsilons) {
 53        _weightSum->fill();
 54        // Find the neutral pions from the decay
 55        vector<ConstGenParticlePtr> pions;
 56        findDecayProducts(p.genParticle(), pions);
 57        const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 58        for (size_t ix=0; ix<pions.size(); ++ix) {
 59          const double pcm = cms_boost.transform(FourMomentum(pions[ix]->momentum())).p();
 60          _histdSigDp->fill(pcm);
 61        }
 62        _histMult->fill(0, pions.size());
 63      }
 64    }
 65
 66
 67    void finalize() {
 68      scale(_histdSigDp, 1./ *_weightSum);
 69      scale(_histMult  , 1./ *_weightSum);
 70    }
 71
 72
 73  private:
 74
 75    /// @{
 76    // count of weights
 77    CounterPtr _weightSum;
 78    /// Histograms
 79    Histo1DPtr _histdSigDp;
 80    BinnedProfilePtr<int> _histMult;
 81    /// @}
 82
 83
 84    void findDecayProducts(ConstGenParticlePtr p, vector<ConstGenParticlePtr>& pions) {
 85      ConstGenVertexPtr dv = p->end_vertex();
 86      for (ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
 87        const int id = pp->pdg_id();
 88        if (id == 111) {
 89          pions.push_back(pp);
 90        } else if (pp->end_vertex())
 91          findDecayProducts(pp, pions);
 92      }
 93    }
 94
 95  };
 96
 97
 98  RIVET_DECLARE_ALIASED_PLUGIN(BELLE_2001_I554520, BELLE_2001_S4598261);
 99
100}