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      // Book histos
22      for (double eVal : allowedEnergies()) {
23
24        const string en = toString(int(eVal));
25        if (isCompatibleWithSqrtS(eVal))  _sqs = en;
26
27        bool is900GeV(en == "900"s);
28        book(_h[en+"pi0"], 1 + is900GeV,1,1);
29      }
30      if (_sqs == "" && !merging()) {
31        throw BeamError("Invalid beam energy for " + name() + "\n");
32      }
33      book(_h["eta"], 3,1,1);
34      // Temporary plots with the binning of _h_etaToPion to construct the eta/pi0 ratio
35      book(_h["temp_eta"],  "TMP/h_eta" , refData(3,1,1));
36      book(_h["temp_pion"], "TMP/h_pion", refData(3,1,1));
37    }
38
39
40    /// Per-event analysis
41    void analyze(const Event& event) {
42
43      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
44      for (const Particle& p : ufs.particles()) {
45        const double normfactor = TWOPI*p.pT()/GeV*2*RAPMAX;
46        if (p.pid() == 111) {
47          // Neutral pion; ALICE corrects for pi0 feed-down from K_0_s and Lambda
48          if (p.hasAncestorWith(Cuts::pid == 310)  ||
49              p.hasAncestorWith(Cuts::pid == 3122) ||
50              p.hasAncestorWith(Cuts::pid == -3122)) continue; //< K_0_s, Lambda, Anti-Lambda
51          _h[_sqs+"pi0"]->fill(p.pT()/GeV, 1.0/normfactor);
52          if (_sqs == "7000"s)  _h["temp_pion"]->fill(p.pT()/GeV);
53        }
54        else if (p.pid() == 221 && _sqs == "7000"s) {
55          // eta meson (only for 7 TeV)
56          _h["eta"]->fill(p.pT()/GeV, 1.0/normfactor);
57          _h["temp_eta"]->fill(p.pT()/GeV);
58        }
59      }
60    }
61
62
63    /// Normalize histos and construct ratio
64    void finalize() {
65      const double sf = crossSection()/microbarn/sumOfWeights();
66      scale(_h, sf);
67
68      // first divide hists with binned axis
69      Estimate1DPtr tmp_ratio;
70      book(tmp_ratio, "TMP/ratio", _h["eta"]->xEdges());
71      divide(_h["temp_eta"], _h["temp_pion"], tmp_ratio);
72      // now convert to strings edges from ref data
73      BinnedEstimatePtr<string> ratio;
74      book(ratio, 4, 1, 1);
75      for (const auto& b : tmp_ratio->bins()) {
76        ratio->bin(b.index()).set(b.val(), b.err());
77      }
78    }
79
80
81  private:
82
83    const double RAPMAX = 0.8;
84    string _sqs = "";
85
86    map<string,Histo1DPtr> _h;
87
88  };
89
90
91  RIVET_DECLARE_PLUGIN(ALICE_2012_I1116147);
92
93}