1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief Electron spectrum  in B decays
 9  class CRYSTAL_BALL_1989_I263581 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(CRYSTAL_BALL_1989_I263581);
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_all, 1, 1, 1);
26      book(_nB, "/TMP/nB");
27    }
28
29    void findDecayProducts(Particle parent, Particles& em) {
30      for (const Particle& p : parent.children()) {
31        if(p.abspid() == PID::EMINUS) {
32          em.push_back(p);
33        }
34        else  {
35          findDecayProducts(p,em);
36        }
37      }
38    }
39
40    /// Perform the per-event analysis
41    void analyze(const Event& event) {
42      if (_edges.empty())  _edges = _h_all->xEdges();
43      // find and loop over Upslion(4S)
44      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
45      for (const Particle& p : ufs.particles(Cuts::pid==300553)) {
46      	// boost to rest frame
47      	LorentzTransform cms_boost;
48      	if (p.p3().mod() > 1*MeV) {
49      	  cms_boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
50        }
51      	// loop over decay products
52      	for (const Particle& p2 : p.children()) {
53      	  if (p2.abspid()==511 || p2.abspid()==521) {
54            _nB->fill();
55      	    Particles em;
56            findDecayProducts(p2, em);
57            for (const Particle& electron : em) {
58              const double en = cms_boost.transform(electron.momentum()).E();
59              _h_all->fill(map2string(en));
60            }
61      	  }
62      	}
63      }
64    }
65
66    string map2string(const double val) const {
67      const size_t idx = _axis.index(val);
68      if (idx && idx <= _edges.size())  return _edges[idx-1];
69      return "OTHER";
70    }
71
72
73    /// Normalise histograms etc., after the run
74    void finalize() {
75      // normalize to number of B decays
76      scale(_h_all, 1./0.05/ *_nB);
77    }
78
79    /// @}
80
81
82    /// @name Histograms
83    /// @{
84    BinnedHistoPtr<string> _h_all;
85    CounterPtr _nB;
86    YODA::Axis<double> _axis = YODA::Axis<double>(48, 0.6, 3.0);
87    vector<string> _edges;
88    /// @}
89
90
91  };
92
93
94  RIVET_DECLARE_PLUGIN(CRYSTAL_BALL_1989_I263581);
95
96}