1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4#include "Rivet/Projections/UnstableParticles.hh"
 5
 6namespace Rivet {
 7
 8
 9  /// Generic analysis looking at various distributions of final state particles
10  ///
11  /// @todo Rename as MC_HADRONS?
12  class MC_IDENTIFIED : public Analysis {
13  public:
14
15    /// Constructor
16    RIVET_DEFAULT_ANALYSIS_CTOR(MC_IDENTIFIED);
17
18
19    /// @name Analysis methods
20    /// @{
21
22    /// Book histograms and initialise projections before the run
23    void init() {
24
25      // Projections
26      const FinalState cnfs(Cuts::abseta < 5.0 && Cuts::pT > 500*MeV);
27      declare(cnfs, "FS");
28      declare(UnstableParticles(Cuts::abseta < 5.0 && Cuts::pT > 500*MeV), "UFS");
29
30      // Histograms
31      // @todo Choose E/pT ranged based on input energies... can't do anything about kin. cuts, though
32
33      book(_histStablePIDs  ,"MultsStablePIDs", 3335, -0.5, 3334.5);
34      book(_histDecayedPIDs ,"MultsDecayedPIDs", 3335, -0.5, 3334.5);
35      book(_histAllPIDs     ,"MultsAllPIDs", 3335, -0.5, 3334.5);
36
37      book(_histEtaPi       ,"EtaPi", 25, 0, 5);
38      book(_histEtaK        ,"EtaK", 25, 0, 5);
39      book(_histEtaLambda   ,"EtaLambda", 25, 0, 5);
40    }
41
42
43
44    /// Perform the per-event analysis
45    void analyze(const Event& event) {
46
47      // Unphysical (debug) plotting of all PIDs in the event, physical or otherwise
48      for(ConstGenParticlePtr gp: HepMCUtils::particles(event.genEvent())) {
49        _histAllPIDs->fill(abs(gp->pdg_id()));
50      }
51
52      // Charged + neutral final state PIDs
53      const FinalState& cnfs = apply<FinalState>(event, "FS");
54      for (const Particle& p : cnfs.particles()) {
55        _histStablePIDs->fill(p.abspid());
56      }
57
58      // Unstable PIDs and identified particle eta spectra
59      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
60      for (const Particle& p : ufs.particles()) {
61        _histDecayedPIDs->fill(p.pid());
62        const double eta_abs = p.abseta();
63        const PdgId pid = p.abspid(); //if (PID::isMeson(pid) && PID::hasStrange()) {
64        if (pid == 211 || pid == 111) _histEtaPi->fill(eta_abs);
65        else if (pid == 321 || pid == 130 || pid == 310) _histEtaK->fill(eta_abs);
66        else if (pid == 3122) _histEtaLambda->fill(eta_abs);
67      }
68
69    }
70
71
72
73    /// Finalize
74    void finalize() {
75      scale(_histStablePIDs, 1/sumOfWeights());
76      scale(_histDecayedPIDs, 1/sumOfWeights());
77      scale(_histAllPIDs, 1/sumOfWeights());
78      scale(_histEtaPi, 1/sumOfWeights());
79      scale(_histEtaK, 1/sumOfWeights());
80      scale(_histEtaLambda, 1/sumOfWeights());
81    }
82
83    /// @}
84
85
86  private:
87
88    /// @name Histograms
89    /// @{
90    Histo1DPtr _histStablePIDs, _histDecayedPIDs, _histAllPIDs;
91    Histo1DPtr _histEtaPi, _histEtaK, _histEtaLambda;
92    /// @}
93
94  };
95
96
97  RIVET_DECLARE_PLUGIN(MC_IDENTIFIED);
98
99}