1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/LeadingParticlesFinalState.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// Inclusive isolated prompt photon analysis with 2012 LHC data
 12  class ATLAS_2016_I1457605 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2016_I1457605);
 17
 18    /// Book histograms and initialise projections before the run
 19    void init() {
 20
 21      FinalState fs;
 22      declare(fs, "FS");
 23
 24      // Consider the final state jets for the energy density calculation
 25      FastJets fj(fs, JetAlg::KT, 0.5);
 26      fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
 27      declare(fj, "KtJetsD05");
 28
 29      // Consider the leading pt photon with |eta| < 2.37 and pT > 25 GeV
 30      LeadingParticlesFinalState photonfs(PromptFinalState(FinalState(Cuts::abseta < 2.37 && Cuts::pT > 25*GeV)));
 31      photonfs.addParticleId(PID::PHOTON);
 32      declare(photonfs, "LeadingPhoton");
 33
 34      // Book the dsigma/dEt (in eta bins) histograms
 35      for (size_t i = 0; i < _eta_bins.size() - 1; ++i) {
 36        if (fuzzyEquals(_eta_bins[i], 1.37)) continue; // skip this bin
 37        int offset = i > 2? 0 : 1;
 38        book(_h_Et_photon[i] ,i + offset, 1, 1);
 39      }
 40
 41    }
 42
 43
 44    /// Return eta bin for either dsigma/dET histogram (area_eta=false) or energy density correction (area_eta=true)
 45    size_t _getEtaBin(double eta_w, bool area_eta) const {
 46      const double eta = fabs(eta_w);
 47      if (!area_eta) {
 48        return binIndex(eta, _eta_bins);
 49      } else {
 50        return binIndex(eta, _eta_bins_areaoffset);
 51      }
 52    }
 53
 54
 55    /// Perform the per-event analysis
 56    void analyze(const Event& event) {
 57      // Retrieve leading photon
 58      Particles photons = apply<LeadingParticlesFinalState>(event, "LeadingPhoton").particles();
 59      if (photons.size() < 1)  vetoEvent;
 60      const Particle& leadingPhoton = photons[0];
 61
 62      // Veto events with photon in ECAL crack
 63      if (inRange(leadingPhoton.abseta(), 1.37, 1.56)) vetoEvent;
 64
 65      // Compute isolation energy in cone of radius .4 around photon (all particles)
 66      FourMomentum mom_in_EtCone;
 67      Particles fs = apply<FinalState>(event, "FS").particles();
 68      for (const Particle& p : fs) {
 69        // Check if it's outside the cone of 0.4
 70        if (deltaR(leadingPhoton, p) >= 0.4) continue;
 71        // Except muons or neutrinos
 72        if (PID::isNeutrino(p.abspid()) || p.abspid() == PID::MUON) continue;
 73        // Increment isolation energy
 74        mom_in_EtCone += p.momentum();
 75      }
 76      // Remove the photon energy from the isolation
 77      mom_in_EtCone -= leadingPhoton.momentum();
 78
 79      // Get the area-filtered jet inputs for computing median energy density, etc.
 80      vector<double> ptDensity;
 81      vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
 82      const FastJets& fast_jets = apply<FastJets>(event, "KtJetsD05");
 83      const auto clust_seq_area = fast_jets.clusterSeqArea();
 84      for (const Jet& jet : fast_jets.jets()) {
 85        const double area = clust_seq_area->area(jet);
 86        if (area > 1e-3 && jet.abseta() < _eta_bins_areaoffset.back())
 87          ptDensities.at( _getEtaBin(jet.abseta(), true) ) += jet.pT()/area;
 88      }
 89      // Compute the median energy density, etc.
 90      for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
 91        const int njets = ptDensities[b].size();
 92        ptDensity += (njets > 0) ? median(ptDensities[b]) : 0.0;
 93      }
 94      // Compute the isolation energy correction (cone area*energy density)
 95      const double etCone_area = PI * sqr(0.4);
 96      const double correction = ptDensity[_getEtaBin(leadingPhoton.abseta(), true)] * etCone_area;
 97
 98      // Apply isolation cut on area-corrected value
 99      // cut is Etiso < 4.8GeV + 4.2E-03 * Et_gamma.
100      if (mom_in_EtCone.Et() - correction > 4.8*GeV + 0.0042*leadingPhoton.Et()) vetoEvent;
101
102      // Fill histograms
103      const size_t eta_bin = _getEtaBin(leadingPhoton.abseta(), false);
104      _h_Et_photon[eta_bin]->fill(leadingPhoton.Et());
105    }
106
107
108    /// Normalise histograms etc., after the run
109    void finalize() {
110      double sf = crossSection() / (picobarn * sumOfWeights());
111      for (size_t i = 0; i < _eta_bins.size()-1; ++i) {
112        if (fuzzyEquals(_eta_bins[i], 1.37)) continue;
113        scale(_h_Et_photon[i], sf);
114      }
115    }
116
117
118  private:
119
120    Histo1DPtr _h_Et_photon[5];
121
122    const vector<double> _eta_bins = {0.00, 0.60, 1.37, 1.56, 1.81, 2.37 };
123    const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0};
124
125  };
126
127
128  RIVET_DECLARE_PLUGIN(ATLAS_2016_I1457605);
129
130}