1//-*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  class ALICE_2012_I1116147 : public Analysis {
 9  public:
10
11    /// Constructor
12    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2012_I1116147);
13
14
15    /// Initialise projections and histograms
16    void init() {
17
18      const UnstableParticles ufs(Cuts::absrap < RAPMAX);
19      declare(ufs, "UFS");
20
21      // Check if cm energy is 7 TeV or 0.9 TeV
22      if (isCompatibleWithSqrtS(900*GeV))       _cm_energy_case = 1;
23      else if (isCompatibleWithSqrtS(7000*GeV)) _cm_energy_case = 2;
24      if (_cm_energy_case == 0)
25        throw UserError("Center of mass energy of the given input is neither 900 nor 7000 GeV.");
26
27      // Book histos
28      if (_cm_energy_case == 1) {
29        book(_h_pi0,       2,1,1);
30      } else {
31        book(_h_pi0,       1,1,1);
32        book(_h_eta,       3,1,1);
33        book(_h_etaToPion, 4,1,1);
34      }
35
36      // Temporary plots with the binning of _h_etaToPion to construct the eta/pi0 ratio
37      book(_temp_h_pion, "TMP/h_pion", refData(4,1,1));
38      book(_temp_h_eta , "TMP/h_eta",  refData(4,1,1));
39    }
40
41
42    /// Per-event analysis
43    void analyze(const Event& event) {
44
45      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
46      for (const Particle& p : ufs.particles()) {
47        const double normfactor = TWOPI*p.pT()/GeV*2*RAPMAX;
48        if (p.pid() == 111) {
49          // Neutral pion; ALICE corrects for pi0 feed-down from K_0_s and Lambda
50          if (p.hasAncestorWith(Cuts::pid == 310) || p.hasAncestorWith(Cuts::pid == 3122) || p.hasAncestorWith(Cuts::pid == -3122)) continue; //< K_0_s, Lambda, Anti-Lambda
51          _h_pi0->fill(p.pT()/GeV, 1.0/normfactor);
52          _temp_h_pion->fill(p.pT()/GeV);
53        } else if (p.pid() == 221 && _cm_energy_case == 2) {
54          // eta meson (only for 7 TeV)
55          _h_eta->fill(p.pT()/GeV, 1.0/normfactor);
56          _temp_h_eta->fill(p.pT()/GeV);
57        }
58      }
59    }
60
61
62    /// Normalize histos and construct ratio
63    void finalize() {
64      scale(_h_pi0, crossSection()/microbarn/sumOfWeights());
65      if (_cm_energy_case == 2) {
66        divide(_temp_h_eta, _temp_h_pion, _h_etaToPion);
67        scale(_h_eta, crossSection()/microbarn/sumOfWeights());
68      }
69    }
70
71
72  private:
73
74    const double RAPMAX = 0.8;
75    int _cm_energy_case = 0;
76
77    Histo1DPtr _h_pi0, _h_eta;
78    Histo1DPtr _temp_h_pion, _temp_h_eta;
79    Estimate1DPtr _h_etaToPion;
80
81  };
82
83
84  RIVET_DECLARE_PLUGIN(ALICE_2012_I1116147);
85
86}