1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4#include "Rivet/Projections/DecayedParticles.hh"
 5
 6namespace Rivet {
 7
 8
 9  /// @brief B0 -> D*- pi+ pi+ pi-
10  class BABAR_2016_I1487722 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2016_I1487722);
15
16
17    /// @name Analysis methods
18    /// @{
19
20    /// Book histograms and initialise projections before the run
21    void init() {
22      UnstableParticles ufs = UnstableParticles(Cuts::abspid==511);
23      declare(ufs, "UFS");
24      DecayedParticles B0(ufs);
25      B0.addStable(PID::PI0);
26      B0.addStable( 413);
27      B0.addStable(-413);
28      B0.addStable( 423);
29      B0.addStable(-423);
30      B0.addStable(PID::PI0);
31      declare(B0, "B0");
32      // histograms
33      book(_h,1,1,1);
34      // efficiency
35      const Estimate1D& ref = refData(2,1,1);
36      _edges = ref.xEdges();
37      _eff = ref.vals();
38    }
39
40    double eff(double mass) {
41      if(mass<_edges[0] || mass>_edges.back()) return 0.;
42      for(unsigned int ix=0;ix<_eff.size();++ix) {
43        if(mass<_edges[ix+1]) return _eff[ix];
44      }
45      return 0.;
46    }
47
48    /// Perform the per-event analysis
49    void analyze(const Event& event) {
50      static const map<PdgId,unsigned int> & mode1   = { { -413,1}, { 211,2}, {-211,1} };
51      static const map<PdgId,unsigned int> & mode1CC = { {  413,1}, {-211,2}, { 211,1} };
52      DecayedParticles B0 = apply<DecayedParticles>(event, "B0");
53      for(unsigned int ix=0;ix<B0.decaying().size();++ix) {
54        int sign = B0.decaying()[ix].pid()/B0.decaying()[ix].abspid();
55        if ( (sign== 1 && B0.modeMatches(ix,4,mode1) ) ||
56             (sign==-1 && B0.modeMatches(ix,4,mode1CC) ) ) {
57          FourMomentum ptotal;
58          for(const Particle & p : B0.decayProducts()[ix].at( sign*211) ) {
59            ptotal+=p.momentum();
60          }
61          for(const Particle & p : B0.decayProducts()[ix].at(-sign*211) ) {
62            ptotal+=p.momentum();
63          }
64          double m = ptotal.mass();
65          _h->fill(m, eff(m));
66        }
67      }
68    }
69
70
71    /// Normalise histograms etc., after the run
72    void finalize() {
73      normalize(_h);
74    }
75
76    /// @}
77
78
79    /// @name Histograms
80    /// @{
81    Histo1DPtr _h;
82    vector<double> _eff;
83    vector<double> _edges;
84    /// @}
85
86
87  };
88
89
90  RIVET_DECLARE_PLUGIN(BABAR_2016_I1487722);
91
92}