1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief B -> Xs l+ l-
 9  class BELLE_2010_I899499 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2010_I899499);
14
15
16    /// @name Analysis methods
17    /// @{
18
19    /// Book histograms and initialise projections before the run
20    void init() {
21      // Initialise and register projections
22      declare(UnstableParticles(Cuts::abspid==521 || Cuts::abspid==511), "UFS");
23      // histos
24      for (unsigned int ix=0;ix<2;++ix) {
25        book(_h[ix],1,1,1+ix);
26      }
27      book(_c,"TMP/Nb");
28    }
29
30    void findDecayProducts(bool& charm, const Particle& mother,
31                           unsigned int& nK0, unsigned int& nKp,
32                           unsigned int& nKm, Particles & lp, Particles & lm) {
33      for (const Particle& p : mother.children()) {
34        int id = p.pid();
35        if (PID::isCharmHadron(p.pid())) charm = true;
36        else if ( id == PID::POSITRON || id == PID::ANTIMUON) lp.push_back(p);
37        else if ( id == PID::ELECTRON || id == PID::MUON    ) lm.push_back(p);
38        else if ( id == PID::KPLUS )    ++nKp;
39        else if (id == PID::KMINUS )    ++nKm;
40        else if (id == PID::K0S)        ++nK0;
41        else if (id == PID::PI0 || id == PID::PIPLUS || id == PID::PIMINUS) {
42          continue;
43        }
44        else if ( !p.children().empty() ) {
45          findDecayProducts(charm, p, nK0, nKp, nKm,lp,lm);
46        }
47      }
48    }
49
50    /// Perform the per-event analysis
51    void analyze(const Event& event) {
52      // Loop over bottoms
53      for (const Particle& bottom : apply<UnstableParticles>(event, "UFS").particles()) {
54       	// remove mixing entries etc
55        if (bottom.children()[0].abspid()==bottom.abspid()) continue;
56        _c->fill();
57        bool charm = false;
58        Particles lp,lm;
59        unsigned int nK0(0),nKp(0),nKm(0);
60        findDecayProducts(charm,bottom, nK0, nKp, nKm,lp,lm);
61        if (charm) continue;
62              unsigned int nk = nKp-nKm+nK0;
63        if ( nk % 2 == 0) continue;
64        if (lp.size()!=1 || lm.size()!=1 || lp[0].pid()!=-lm[0].pid()) continue;
65          FourMomentum qq = lp[0].mom()+lm[0].mom();
66          double q2 = qq.mass2();
67          if(q2<0.04) continue;
68          _h[0]->fill((bottom.mom()-qq).mass());
69          _h[1]->fill(q2);
70      }
71    }
72
73
74    /// Normalise histograms etc., after the run
75    void finalize() {
76      scale(_h, 0.5e6/ *_c);
77    }
78
79    /// @}
80
81
82    /// @name Histograms
83    /// @{
84    CounterPtr _c;
85    Histo1DPtr _h[2];
86    /// @}
87
88
89  };
90
91
92  RIVET_DECLARE_PLUGIN(BELLE_2010_I899499);
93
94}