1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Analyses/AtlasCommon.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Add a short analysis description here
 10  class ATLAS_2015_I1360290 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2015_I1360290);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      // Centrality projection.
 25      declareCentrality(ATLAS::SumET_PBPB_Centrality(), "ATLAS_PBPB_CENTRALITY", "sumETFwd","sumETFwd");
 26      // Trigger projection.
 27      declare(ATLAS::MinBiasTrigger(),"Trigger");
 28      // The measured final state.
 29      declare(ChargedFinalState (Cuts::abseta < 2. && Cuts::pT > 0.5*GeV && Cuts::pT < 150.0*GeV), "CFS");
 30
 31      taa = {26.3, 20.6, 14.4, 8.73, 5.05, 2.70, 1.34, 0.41};
 32      centData = {5., 10., 20., 30., 40., 50., 60., 80.};
 33
 34      for (int i = 0, N = centData.size(); i < N; ++i) {
 35        // eta hists starts from table 55 ( first 1.7 < pT < 2.0)
 36        book(histEta1[centData[i]], 55 + i, 1, 1);
 37        // From table 64, 6.7 < pT < 7.7
 38        book(histEta2[centData[i]], 64 + i, 1, 1 );
 39        // From table 73, 19.9 < pT < 22.8
 40        book(histEta3[centData[i]], 73 + i, 1, 1 );
 41        // From table 82, 59.8 < pT < 94.8
 42        book(histEta4[centData[i]], 82 + i, 1, 1 );
 43        // pt hists starts from table 2 on hepmc, |eta| < 2.0
 44        book(histpT[centData[i]], 2 + i, 1, 1);
 45        // keep track of sow in centrality bins.
 46        book(sow[centData[i]], "sow_" + toString(i));
 47      }
 48
 49    }
 50
 51
 52    /// Perform the per-event analysis
 53    void analyze(const Event& event) {
 54      if ( !apply<ATLAS::MinBiasTrigger>(event, "Trigger")() ) vetoEvent;
 55      const CentralityProjection& cent = apply<CentralityProjection>(event,"sumETFwd");
 56      double c = cent();
 57      // Find the correct centrality histograms
 58      auto hItr1 = histEta1.upper_bound(c);
 59      if (hItr1 == histEta1.end()) return;
 60      auto hItr2 = histEta2.upper_bound(c);
 61      if (hItr2 == histEta2.end()) return;
 62      auto hItr3 = histEta3.upper_bound(c);
 63      if (hItr3 == histEta3.end()) return;
 64      auto hItr4 = histEta4.upper_bound(c);
 65      if (hItr4 == histEta4.end()) return;
 66      auto hpTItr = histpT.upper_bound(c);
 67      if (hpTItr == histpT.end()) return;
 68      // Find the correct sow.
 69      auto sItr = sow.upper_bound(c);
 70      if (sItr == sow.end()) return;
 71      sItr->second->fill();
 72
 73      for (const auto& p : apply<ChargedFinalState>(event,"CFS").particles()) {
 74        const double pT = p.pT();
 75        const double eta = p.abseta();
 76        if (pT > 1.7 && pT < 2.0) hItr1->second->fill(eta, 0.5);
 77        else if (pT > 6.7 && pT < 7.7) hItr2->second->fill(eta, 0.5);
 78        else if (pT > 19.9 && pT < 22.8) hItr3->second->fill(eta, 0.5);
 79        else if (pT > 59.8 && pT < 94.8) hItr4->second->fill(eta, 0.5);
 80        if (eta < 2) hpTItr->second->fill(pT, 1.0/2./M_PI/pT/4.);
 81      }
 82
 83    }
 84
 85
 86    /// Normalise histograms etc., after the run
 87    void finalize() {
 88      for (int i = 0, N = centData.size(); i < N; ++i) {
 89        histEta1[centData[i]]->scaleW(1./sow[centData[i]]->sumW());
 90        histEta2[centData[i]]->scaleW(1./sow[centData[i]]->sumW());
 91        histEta3[centData[i]]->scaleW(1./sow[centData[i]]->sumW());
 92        histEta4[centData[i]]->scaleW(1./sow[centData[i]]->sumW());
 93        histpT[centData[i]]->scaleW(1./sow[centData[i]]->sumW()/taa[i]);
 94
 95      }
 96
 97    }
 98
 99    /// @}
100
101
102    /// @name Histograms
103    /// @{
104    // The centrality binned histograms
105    map<double, Histo1DPtr> histEta1;
106    map<double, Histo1DPtr> histEta2;
107    map<double, Histo1DPtr> histEta3;
108    map<double, Histo1DPtr> histEta4;
109    map<double, Histo1DPtr> histpT;
110    map<double, CounterPtr> sow;
111
112
113    vector<double> centData;
114    vector<double> taa;
115    /// @}
116
117
118  };
119
120
121  RIVET_DECLARE_PLUGIN(ATLAS_2015_I1360290);
122
123
124}