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 D+ -> K+K-pi+
10  class CLEO_2008_I791716 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_2008_I791716);
15
16
17    /// @name Analysis methods
18    /// @{
19
20    /// Book histograms and initialise projections before the run
21    void init() {
22      // Initialise and register projections
23      UnstableParticles ufs = UnstableParticles(Cuts::abspid==411);
24      declare(ufs, "UFS");
25      DecayedParticles DP(ufs);
26      DP.addStable(PID::PI0);
27      DP.addStable(PID::K0S);
28      DP.addStable(PID::ETA);
29      DP.addStable(PID::ETAPRIME);
30      declare(DP, "DP");
31      // histos
32      book(_h_Kppi,1,1,1);
33      book(_h_Kmpi,1,1,2);
34      book(_h_KK  ,1,1,3);
35      book(_dalitz, "dalitz",50,0.,2.,50,0.9,3.1);
36    }
37
38    /// Perform the per-event analysis
39    void analyze(const Event& event) {
40      static const map<PdgId,unsigned int> & mode   = { { 321,1},{-321,1}, { 211,1}};
41      static const map<PdgId,unsigned int> & modeCC = { { 321,1},{-321,1}, {-211,1}};
42      DecayedParticles DP = apply<DecayedParticles>(event, "DP");
43      // loop over particles
44      for(unsigned int ix=0;ix<DP.decaying().size();++ix) {
45	int sign = 1;
46	if (DP.decaying()[ix].pid()>0 && DP.modeMatches(ix,3,mode)) {
47	  sign=1;
48	}
49	else if  (DP.decaying()[ix].pid()<0 && DP.modeMatches(ix,3,modeCC)) {
50	  sign=-1;
51	}
52	else
53	  continue;
54	const Particles & Kp = DP.decayProducts()[ix].at( sign*321);
55	const Particles & Km = DP.decayProducts()[ix].at(-sign*321);
56	const Particles & pip= DP.decayProducts()[ix].at( sign*211);
57	double mminus = (Km[0].momentum()+pip[0].momentum() ).mass2();
58	double mplus  = (Kp[0].momentum()+pip[0].momentum() ).mass2();
59	double mKK    = (Kp[0].momentum()+Km[0].momentum()).mass2();
60	_h_Kppi->fill(mplus);
61	_h_Kmpi->fill(mminus);
62	_h_KK  ->fill(mKK);
63	_dalitz->fill(mminus,mKK); 
64      }
65    }
66
67
68    /// Normalise histograms etc., after the run
69    void finalize() {
70      normalize(_h_Kppi);
71      normalize(_h_Kmpi);
72      normalize(_h_KK  );
73      normalize(_dalitz);
74    }
75
76    /// @}
77
78
79    /// @name Histograms
80    /// @{
81    Histo1DPtr _h_Kppi,_h_Kmpi,_h_KK;
82    Histo2DPtr _dalitz;
83    /// @}
84
85
86  };
87
88
89  RIVET_DECLARE_PLUGIN(CLEO_2008_I791716);
90
91}