1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief  J/Psi, Psi(2S), chi_c1 spectra at the Upsilon(4S)
  9  class CLEOII_1995_I382221 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEOII_1995_I382221);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // projections
 22      declare(UnstableParticles(Cuts::pid==300553), "UFS");
 23      // histos
 24      for (unsigned int ix=0; ix<4; ++ix) {
 25        book(_h[ix], 1+ix,1,1);
 26      }
 27      book(_c,"TMP/nUps");
 28    }
 29
 30    void findDecayProducts(Particle parent, Particles & primary,
 31			   Particles & secondary, bool second) {
 32      for (const Particle& p :parent.children()) {
 33        if (p.pid()==443) {
 34          if (second) secondary.push_back(p);
 35          else        primary.push_back(p);
 36       	  continue;
 37       	}
 38       	else if (p.pid()==100443 || p.pid()==20443) {
 39          if (second) secondary.push_back(p);
 40          else        primary.push_back(p);
 41          findDecayProducts(p,primary,secondary,true);
 42        }
 43        else if (p.pid()==10441||p.pid()==445) {
 44          findDecayProducts(p,primary,secondary,true);
 45        }
 46        else if(!p.children().empty()) {
 47          findDecayProducts(p,primary,secondary,second);
 48        }
 49      }
 50    }
 51
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 56      for (const Particle& p : ufs.particles()) {
 57        _c->fill();
 58        const LorentzTransform boost =
 59              LorentzTransform::mkFrameTransformFromBeta(p.mom().betaVec());
 60       	Particles primary,secondary;
 61       	findDecayProducts(p,primary,secondary,false);
 62       	for (const Particle& pp : primary) {
 63      	  const double pcm = boost.transform(pp.mom()).p3().mod();
 64          if (pp.pid()==443) {
 65            _h[0]->fill(pcm);
 66            _h[1]->fill(pcm);
 67          }
 68          else if(pp.pid()==100443) _h[2]->fill(pcm);
 69          else if(pp.pid()==20443)  _h[3]->fill(pcm);
 70        }
 71       	for (const Particle& pp : secondary) {
 72      	  const double pcm = boost.transform(pp.mom()).p3().mod();
 73          if (pp.pid()==443)        _h[0]->fill(pcm);
 74          else if (pp.pid()==100443) _h[2]->fill(pcm);
 75          else if (pp.pid()==20443)  _h[3]->fill(pcm);
 76        }
 77      }
 78    }
 79
 80
 81    /// Normalise histograms etc., after the run
 82    void finalize() {
 83      // 100 convert to % and 0.5 as 2 B decays per Upslion(4S)
 84	    scale(_h, 50./ *_c);
 85    }
 86
 87    /// @}
 88
 89
 90    /// @name Histograms
 91    /// @{
 92    Histo1DPtr _h[4];
 93    CounterPtr _c;
 94    /// @}
 95
 96
 97  };
 98
 99
100  RIVET_DECLARE_PLUGIN(CLEOII_1995_I382221);
101
102}