1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4#include "Rivet/Projections/UnstableParticles.hh"
 5
 6namespace Rivet {
 7
 8
 9  /// @brief Add a short analysis description here
10  class BABAR_2006_I731865 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2006_I731865);
15
16
17    /// @name Analysis methods
18    /// @{
19
20    /// Book histograms and initialise projections before the run
21    void init() {
22
23      // Initialise and register projections
24      declare(FinalState(), "FS");
25      declare(UnstableParticles(), "UFS");
26
27      book(_nPhiEta, 1, 1, 2);
28    }
29
30    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
31      for (const Particle &child : p.children()) {
32        if(child.children().empty()) {
33          nRes[child.pid()]-=1;
34          --ncount;
35        }
36        else {
37          findChildren(child,nRes,ncount);
38        }
39      }
40    }
41
42    /// Perform the per-event analysis
43    void analyze(const Event& event) {
44      const FinalState& fs = apply<FinalState>(event, "FS");
45
46      map<long,int> nCount;
47      int ntotal(0);
48      for (const Particle& p : fs.particles()) {
49        nCount[p.pid()] += 1;
50        ++ntotal;
51      }
52      const FinalState& ufs = apply<FinalState>(event, "UFS");
53      for (const Particle& p : ufs.particles()) {
54        if(p.children().empty()) continue;
55        if(p.pid()!=333) continue;
56        map<long,int> nRes=nCount;
57        int ncount = ntotal;
58        findChildren(p,nRes,ncount);
59        for (const Particle& p2 : ufs.particles()) {
60          if (p2.pid()!=221) continue;
61          if (p2.parents()[0].isSame(p)) continue;
62          map<long,int> nResB = nRes;
63          int ncountB = ncount;
64          findChildren(p2,nResB,ncountB);
65          if (ncountB!=0) continue;
66          bool matched2 = true;
67          for (const auto& val : nResB) {
68            if (val.second!=0) {
69              matched2 = false;
70              break;
71            }
72          }
73          if (matched2) {
74            _nPhiEta->fill(Ecm);
75          }
76        }
77      }
78    }
79
80
81    /// Normalise histograms etc., after the run
82    void finalize() {
83      scale(_nPhiEta, crossSection()/ sumOfWeights() /femtobarn);
84    }
85
86    /// @}
87
88    /// @name Histograms
89    /// @{
90    BinnedHistoPtr<string> _nPhiEta;
91    const string Ecm = "10.58";
92    /// @}
93
94  };
95
96
97  RIVET_DECLARE_PLUGIN(BABAR_2006_I731865);
98
99}