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