1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Electron spectrum in D decays
  9  class CLEOC_2006_I715096 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEOC_2006_I715096);
 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(),"UFS");
 23      // Book histograms
 24      // specify custom binning
 25      book(_h_Dp, 1,1,1);
 26      book(_h_D0, 1,1,2);
 27    }
 28
 29    void findDecayProducts(Particle parent, Particles & em, Particles & ep,
 30			   Particles & nue, Particles & nueBar) {
 31      for(const Particle & p : parent.children()) {
 32	if(p.pid() == PID::EMINUS) {
 33	  em.push_back(p);
 34	}
 35	else if(p.pid() == PID::EPLUS) {
 36	  ep.push_back(p);
 37	}
 38	else if(p.pid() == PID::NU_E) {
 39	  nue.push_back(p);
 40	}
 41	else if(p.pid() == PID::NU_EBAR) {
 42	  nueBar.push_back(p);
 43	}
 44	else if(PID::isCharmHadron(p.pid())) {
 45	  findDecayProducts(p,em,ep,nue,nueBar);
 46	}
 47	else if(!PID::isHadron(p.pid())) {
 48	  findDecayProducts(p,em,ep,nue,nueBar);
 49	}
 50      }
 51    }
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55      // find and loop over psi(3770)
 56      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 57      for (const Particle& p : ufs.particles(Cuts::pid==30443)) {
 58	// boost to rest frame
 59	LorentzTransform boost;
 60	if (p.p3().mod() > 1*MeV)
 61	  boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 62	// loop over decay products
 63	for(const Particle & p2 : p.children()) {
 64	  if(p2.abspid()==411 || p2.abspid()==421) {
 65	    Particles em,ep,nue,nueBar;
 66	    findDecayProducts(p2,em,ep,nue,nueBar);
 67	    if(em.size()==1 && nueBar.size()==1) {
 68	      double pmod = boost.transform(em[0].momentum()).p3().mod();
 69	      if(p2.abspid()==411) _h_Dp->fill(pmod);
 70	      else                 _h_D0->fill(pmod);
 71	    }
 72	    else if(ep.size()==1 && nue.size()==1) {
 73	      double pmod = boost.transform(ep[0].momentum()).p3().mod();
 74	      if(p2.abspid()==411) _h_Dp->fill(pmod);
 75	      else                 _h_D0->fill(pmod);
 76	    }
 77	  }	
 78	}
 79      }
 80    }
 81
 82
 83    /// Normalise histograms etc., after the run
 84    void finalize() {
 85      normalize(_h_D0);
 86      normalize(_h_Dp);
 87    }
 88
 89    ///@}
 90
 91
 92    /// @name Histograms
 93    ///@{
 94    Histo1DPtr _h_Dp,_h_D0;
 95    ///@}
 96
 97
 98  };
 99
100
101  RIVET_DECLARE_PLUGIN(CLEOC_2006_I715096);
102
103}