1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> e spectrum
  9  class ARGUS_1993_I357133 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1993_I357133);
 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      // histograms
 24      for (unsigned int ix=0; ix<2; ++ix) {
 25        book(_h[ix],1+ix,1,1);
 26      }
 27      book(_c,"TMP/nb");
 28    }
 29
 30    void findDecayProducts(const Particle& parent,
 31                           vector<Particles>&  em, vector<Particles>& ep,
 32                           vector<Particles>& nue, vector<Particles>& nueBar,
 33                           bool& charm, bool secondary) {
 34      for (const Particle & p : parent.children()) {
 35      	if (PID::isCharmHadron(p.pid())) {
 36      	  charm=true;
 37          findDecayProducts(p,em,ep,nue,nueBar,charm,true);
 38      	}
 39      	else if (p.pid() == PID::EMINUS) {
 40      	  em[secondary].push_back(p);
 41      	}
 42      	else if (p.pid() == PID::EPLUS) {
 43      	  ep[secondary].push_back(p);
 44      	}
 45      	else if (p.pid() == PID::NU_E ) {
 46      	  nue[secondary].push_back(p);
 47      	}
 48      	else if (p.pid() == PID::NU_EBAR ) {
 49      	  nueBar[secondary].push_back(p);
 50      	}
 51      	else if (PID::isBottomHadron(p.pid())) {
 52      	  findDecayProducts(p,em,ep,nue,nueBar,charm,false);
 53      	}
 54      	else if (!PID::isHadron(p.pid())) {
 55      	  findDecayProducts(p,em,ep,nue,nueBar,charm,secondary);
 56      	}
 57      }
 58    }
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 63      for (const Particle& p : ufs.particles(Cuts::pid==300553)) {
 64      	// boost to rest frame
 65      	LorentzTransform cms_boost;
 66      	if (p.p3().mod() > 1*MeV) {
 67      	  cms_boost = LorentzTransform::mkFrameTransformFromBeta(p.mom().betaVec());
 68        }
 69      	// loop over decay products
 70      	for (const Particle & p2 : p.children()) {
 71      	  if (p2.abspid()==511 || p2.abspid()==521) {
 72      	    _c->fill();
 73            // find decay products
 74            bool charm = false;
 75            vector<Particles> em(2),ep(2),nue(2),nueBar(2);
 76            findDecayProducts(p2,em,ep,nue,nueBar,charm,false);
 77            for (unsigned int isec=0; isec<2; ++isec) {
 78              FourMomentum pl,pnu;
 79              if (em[isec].size()==1 && nueBar[isec].size()==1 && em[isec][0].pid()+1==-nueBar[isec][0].pid()) {
 80                pl  = cms_boost.transform(em[isec][0].mom());
 81                pnu = cms_boost.transform(nueBar[isec][0].mom());
 82              }
 83              else if (ep[isec].size()==1 && nue[isec].size()==1 && nue[isec][0].pid()==-ep[isec][0].pid()+1) {
 84                pl  = cms_boost.transform(ep[isec][0].mom());
 85                pnu = cms_boost.transform(nue[isec][0].mom());
 86              }
 87              else {
 88                continue;
 89              }
 90              _h[isec]->fill(pl.p3().mod());
 91            }
 92          }
 93        }
 94      }
 95    }
 96
 97
 98    /// Normalise histograms etc., after the run
 99    void finalize() {
100      scale(_h, 1.0/ *_c);
101    }
102
103    /// @}
104
105
106    /// @name Histograms
107    /// @{
108    Histo1DPtr _h[2];
109    CounterPtr _c;
110    /// @}
111
112
113  };
114
115
116  RIVET_DECLARE_PLUGIN(ARGUS_1993_I357133);
117
118}