1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief tau -> K phi
 9  class BELLE_2006_I725750 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2006_I725750);
14
15
16    /// @name Analysis methods
17    /// @{
18
19    /// Book histograms and initialise projections before the run
20    void init() {
21      declare(UnstableParticles(), "UFS");
22      book(_hist,  1, 1, 1);
23    }
24
25    void findDecayProducts(const Particle & mother,
26                           unsigned int & nstable,
27                           Particles& phi, Particles& K) {
28      for(const Particle & p : mother.children()) {
29	long id = p.abspid();
30        if (id == PID::PHI ) {
31          phi.push_back(p);
32          ++nstable;
33	}
34        else if (id == PID::KPLUS) {
35          K.push_back(p);
36          ++nstable;
37        }
38        else if ( !p.children().empty() ) {
39          findDecayProducts(p, nstable, phi,K);
40        }
41        else
42          ++nstable;
43      }
44    }
45
46    /// Perform the per-event analysis
47    void analyze(const Event& event) {
48      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
49      for (const Particle& p : ufs.particles(Cuts::abspid==PID::TAU)) {
50        Particles phi, K;
51        unsigned int nstable = 0;
52        // find the decay products we want
53        findDecayProducts(p, nstable, phi, K);
54        if (nstable != 3) continue;
55	// K phi
56        if (K.size() == 1 && phi.size() == 1)
57	  _hist->fill((phi[0].momentum()+K[0].momentum()).mass());
58      }
59    }
60
61
62    /// Normalise histograms etc., after the run
63    void finalize() {
64      normalize(_hist);
65    }
66
67    /// @}
68
69
70    /// @name Histograms
71    /// @{
72    Histo1DPtr _hist;
73    /// @}
74
75
76  };
77
78
79  RIVET_DECLARE_PLUGIN(BELLE_2006_I725750);
80
81}