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