1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief q^2 in D0 and D+ semi-lepto
 9  class CLEOC_2009_I823313 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEOC_2009_I823313);
14
15
16    /// @name Analysis methods
17    /// @{
18
19    /// Book histograms and initialise projections before the run
20    void init() {
21
22      // Initialise and register projections
23      declare(UnstableParticles(), "UFS");
24      // histograms`
25      book(_h_q2_D0_pi,1,1,1);
26      book(_h_q2_D0_K ,1,1,2);
27      book(_h_q2_Dp_pi,1,1,3);
28      book(_h_q2_Dp_K ,1,1,4);
29    }
30
31    // Calculate the Q2 using mother and daugher meson
32    double q2(const Particle& B, int mesonID) {
33      FourMomentum q = B.mom() - select(B.children(), Cuts::abspid==abs(mesonID))[0];
34      return q*q;
35    }
36
37    // Check for explicit decay into pdgids
38    bool isSemileptonicDecay(const Particle& mother, vector<int> ids) {
39      // Trivial check to ignore any other decays but the one in question modulo photons
40      const Particles children = mother.children(Cuts::pid!=PID::PHOTON);
41      if (children.size()!=ids.size()) return false;
42      // Check for the explicit decay
43      return all(ids, [&](int i){return count(children, hasPID(i))==1;});
44    }
45
46    /// Perform the per-event analysis
47    void analyze(const Event& event) {
48      // Loop over D mesons
49      for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==PID::D0 or
50									       Cuts::abspid==PID::DPLUS )) {
51        if (p.abspid()==PID::D0) {
52	  if(isSemileptonicDecay(p, {PID::PIMINUS, PID::POSITRON, PID::NU_E}) ||
53	     isSemileptonicDecay(p, {PID::PIPLUS , PID::ELECTRON, PID::NU_EBAR}) )
54	    _h_q2_D0_pi->fill(q2(p, PID::PIMINUS));
55	  else if(isSemileptonicDecay(p, {PID::KMINUS, PID::POSITRON, PID::NU_E}) ||
56		  isSemileptonicDecay(p, {PID::KPLUS , PID::ELECTRON, PID::NU_EBAR}))
57	    _h_q2_D0_K ->fill(q2(p, PID::KMINUS));
58        }
59	else if(p.abspid()==PID::DPLUS) {
60	  if(isSemileptonicDecay(p, {PID::PI0, PID::POSITRON, PID::NU_E})  ||
61	     isSemileptonicDecay(p, {PID::PI0, PID::ELECTRON, PID::NU_EBAR}))
62	    _h_q2_Dp_pi->fill(q2(p, PID::PI0));
63	  else if(isSemileptonicDecay(p, {-311, PID::POSITRON, PID::NU_E}))
64	    _h_q2_Dp_K ->fill(q2(p, -311));
65	  else if(isSemileptonicDecay(p, { 311, PID::ELECTRON, PID::NU_EBAR}))
66	    _h_q2_Dp_K ->fill(q2(p, 311));
67	  else if(isSemileptonicDecay(p, {PID::K0S, PID::POSITRON, PID::NU_E}) ||
68		  isSemileptonicDecay(p, {PID::K0S, PID::ELECTRON, PID::NU_EBAR}))
69	    _h_q2_Dp_K ->fill(q2(p, PID::K0S));
70	  else if(isSemileptonicDecay(p, {PID::K0L, PID::POSITRON, PID::NU_E}) ||
71		  isSemileptonicDecay(p, {PID::K0L, PID::ELECTRON, PID::NU_EBAR}))
72	    _h_q2_Dp_K ->fill(q2(p, PID::K0L));
73	}
74      }
75    }
76
77    /// Normalise histograms etc., after the run
78    void finalize() {
79      normalize(_h_q2_D0_pi);
80      normalize(_h_q2_D0_K );
81      normalize(_h_q2_Dp_pi);
82      normalize(_h_q2_Dp_K );
83    }
84
85    /// @}
86
87
88    /// @name Histograms
89    /// @{
90    Histo1DPtr _h_q2_D0_pi, _h_q2_D0_K, _h_q2_Dp_pi, _h_q2_Dp_K;
91    /// @}
92
93
94  };
95
96
97  RIVET_DECLARE_PLUGIN(CLEOC_2009_I823313);
98
99}