1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/Beam.hh"
 4#include "Rivet/Projections/ChargedFinalState.hh"
 5#include "Rivet/Projections/UnstableParticles.hh"
 6
 7namespace Rivet {
 8
 9
10  /// @brief Charm meson spectra at 29 GeV
11  class HRS_1988_I23360 : public Analysis {
12  public:
13
14    /// Constructor
15    RIVET_DEFAULT_ANALYSIS_CTOR(HRS_1988_I23360);
16
17
18    /// @name Analysis methods
19    /// @{
20
21    /// Book histograms and initialise projections before the run
22    void init() {
23
24      // Initialise and register projections
25      declare(Beam(), "Beams");
26      declare(UnstableParticles(), "UFS");
27      const ChargedFinalState cfs;
28      declare(cfs, "CFS");
29
30      // Book histograms
31      book(_h_Dstar1,1,1,1);
32      book(_h_Dstar2,1,1,2);
33      book(_h_D0    ,2,1,1);
34      book(_h_Dp    ,2,1,2);
35    }
36
37
38    /// Perform the per-event analysis
39    void analyze(const Event& event) {
40      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
41      int nch = cfs.particles().size();
42      if(nch<5) vetoEvent;
43      // Get beams and average beam momentum
44      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
45      const double meanBeamMom = ( beams.first.p3().mod() +
46                                   beams.second.p3().mod() ) / 2.0;
47      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
48
49      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
50      for(const Particle & p : ufs.particles(Cuts::abspid==413 || Cuts::abspid==421 || Cuts::abspid==411 )) {
51         double xE = p.E()/meanBeamMom;
52	 int id = abs(p.pid());
53	 if(id==413) {
54	   _h_Dstar1->fill(xE);
55	   _h_Dstar2->fill(xE);
56	 }
57	 else if(id==421) {
58	   _h_D0->fill(xE);
59	 }
60	 else if(id==411) {
61	   _h_Dp->fill(xE);
62	 }
63      }
64    }
65
66
67    /// Normalise histograms etc., after the run
68    void finalize() {
69      double fact = sqr(sqrtS())*crossSection()/sumOfWeights()/microbarn;
70      scale(_h_Dstar1, fact);
71      normalize(_h_Dstar2);
72      scale(_h_D0,fact);
73      scale(_h_Dp,fact);
74    }
75
76    /// @}
77
78
79    /// @name Histograms
80    /// @{
81    Histo1DPtr _h_Dstar1,_h_Dstar2,_h_D0,_h_Dp;
82    /// @}
83
84
85  };
86
87
88  RIVET_DECLARE_PLUGIN(HRS_1988_I23360);
89
90
91}