1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/WFinder.hh"
  6#include "Rivet/Projections/ZFinder.hh"
  7#include "fastjet/tools/Filter.hh"
  8#include "fastjet/tools/Pruner.hh"
  9#include "Rivet/Tools/AtlasCommon.hh"
 10
 11namespace Rivet {
 12
 13  class MC_JET_IN_HI : public Analysis {
 14  public:
 15
 16    /// @name Constructors etc.
 17    //@{
 18
 19    /// Constructor
 20    RIVET_DEFAULT_ANALYSIS_CTOR(MC_JET_IN_HI);
 21
 22    //@}
 23
 24
 25  public:
 26
 27    string ts(int in) {
 28      std::stringstream ss;
 29      ss << in;
 30      return ss.str();
 31    }
 32    /// @name Analysis methods
 33    //@{
 34
 35    /// Book histograms and initialise projections before the run
 36    void init() {
 37      // Declare centrality projection - we use the ATLAS PbPb definition
 38      // to be able to compare to data.
 39      declareCentrality(ATLAS::SumET_PBPB_Centrality(),"ATLAS_PBPB_CENTRALITY",
 40		      "sumETFwd","sumETFwd");
 41      // The final state where jets are found
 42      FinalState fs(Cuts::abseta < 2.5);
 43      declare(fs, "FS");
 44 
 45      ZFinder zfinder(fs, Cuts::abseta < 2.5 && Cuts::pT > 30*GeV, PID::MUON, 80*GeV, 100*GeV, 0.2,
 46                      ZFinder::ChargedLeptons::PROMPT, ZFinder::ClusterPhotons::NODECAY, ZFinder::AddPhotons::YES);
 47      declare(zfinder, "ZFinder");
 48
 49      // Z+jet jet collections
 50      declare(FastJets(zfinder.remainingFinalState(), FastJets::ANTIKT, 0.3), "JetsAK3");
 51      declare(FastJets(zfinder.remainingFinalState(), FastJets::ANTIKT, 0.5), "JetsAK5");
 52      declare(FastJets(zfinder.remainingFinalState(), FastJets::ANTIKT, 0.7), "JetsAK7");
 53      declare(FastJets(zfinder.remainingFinalState(), FastJets::ANTIKT, 0.9), "JetsAK9");
 54
 55      jetFinders = {"JetsAK3", "JetsAK5", "JetsAK7", "JetsAK9"};
 56      
 57      h_zpT.resize(jetFinders.size());
 58      h_jetpT.resize(jetFinders.size());
 59      for (size_t i = 0; i < jetFinders.size(); ++i) {
 60        string s = jetFinders[i];
 61        book(h_zpT[i], s + "zpT",logspace(50, 1.0,1000));
 62        book(h_jetpT[i], s + "jetpT",logspace(50, 1.0,1000));
 63      }
 64      book(incSow, "incSow");
 65
 66      centData = {0., 0.2, 0.4, 0.6, 0.8,};
 67      for (size_t i = 0; i < centData.size(); ++i) {
 68        book(c_jetpT[centData[i]], "cjetpT" + ts(i),logspace(100, 10.0,1000));
 69        book(c_zpT[centData[i]], "czpt" + ts(i),logspace(100, 10.0,1000));
 70	      book(sow[centData[i]], "sow_" + ts(i));
 71      }
 72    }
 73
 74    bool isBackToBack_zj(const ZFinder& zf, const fastjet::PseudoJet& psjet) {
 75      const FourMomentum& z = zf.bosons()[0].momentum();
 76      const FourMomentum jmom(psjet.e(), psjet.px(), psjet.py(), psjet.pz());
 77      return (deltaPhi(z, jmom) > 7.*M_PI/8. );
 78    }
 79
 80
 81    /// Perform the per-event analysis
 82    void analyze(const Event& event) {
 83      
 84      // Get the Z
 85      const ZFinder& zfinder = apply<ZFinder>(event, "ZFinder");
 86      if (zfinder.bosons().size() != 1) vetoEvent;
 87      Particle z = zfinder.bosons()[0];
 88      Particle l1 = zfinder.constituents()[0];
 89      Particle l2 = zfinder.constituents()[1];
 90
 91      // Require a high-pT Z (and constituents)
 92      if (l1.pT() < 10*GeV || l2.pT() < 10*GeV || z.pT() < 60*GeV) vetoEvent;
 93
 94      // Get the centrality
 95      const double c = apply<CentralityProjection>(event,"sumETFwd")();
 96      auto jetpTItr = c_jetpT.upper_bound(c);
 97      auto zpTItr = c_zpT.upper_bound(c);
 98      auto sowItr = sow.upper_bound(c);
 99      if (jetpTItr == c_jetpT.end() || zpTItr == c_zpT.end() || 
100        sowItr == sow.end()) vetoEvent;
101      sowItr->second->fill();
102      incSow->fill();
103      // Get the  jets
104      for (size_t i = 0; i < jetFinders.size(); ++i ) {
105        const PseudoJets& psjets = apply<FastJets>(event, 
106	  jetFinders[i]).pseudoJetsByPt(30.0*GeV);
107        if (!psjets.empty()) {
108        // Get the leading jet and make sure it's back-to-back with the Z
109        const fastjet::PseudoJet& j0 = psjets[0];
110	if (isBackToBack_zj(zfinder, j0)) {
111	    // Fill the centrality inclusive histograms
112	    h_zpT[i]->fill(z.pT());
113	    h_jetpT[i]->fill(j0.perp());
114	    // Fill centrality dept histograms only for R = 0.3
115	    if (i == 0) {
116              jetpTItr->second->fill(j0.perp());
117	      zpTItr->second->fill(z.pT());
118	    }
119	  }
120        }
121      }
122    }
123
124
125    /// Normalise histograms etc., after the run
126    void finalize() {
127       for(size_t i = 0; i < jetFinders.size(); ++i) {
128         h_zpT[i]->scaleW(1./incSow->sumW());
129         h_jetpT[i]->scaleW(1./incSow->sumW());
130       }
131       for (size_t i = 0; i < centData.size(); ++i) {
132         c_jetpT[centData[i]]->scaleW(1./sow[centData[i]]->sumW());
133         c_zpT[centData[i]]->scaleW(1./sow[centData[i]]->sumW());
134       }
135    }
136
137
138    //@}
139
140
141  private:
142    vector<string> jetFinders;
143    // Centrality inclusive histograms
144    vector<Histo1DPtr> h_zpT;
145    vector<Histo1DPtr> h_jetpT;
146    CounterPtr incSow;
147    // Centrality intervals
148    vector<double> centData;
149    // Centrality binned histograms
150    map<double, Histo1DPtr> c_jetpT;
151    map<double, Histo1DPtr> c_zpT;
152    map<double, CounterPtr> sow;
153
154  };
155
156
157  // The hook for the plugin system
158  RIVET_DECLARE_PLUGIN(MC_JET_IN_HI);
159
160}