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  class ATLAS_2010_I882463 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2010_I882463);
 15
 16
 17    /// Book histograms and initialise projections before the run
 18    void init() {
 19      FinalState fs;
 20      declare(fs, "FS");
 21
 22      FastJets fj(fs, JetAlg::KT, 0.5);
 23      fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
 24      declare(fj, "KtJetsD05");
 25
 26      LeadingParticlesFinalState photonfs(FinalState(Cuts::abseta < 1.81 && Cuts::pT > 15*GeV));
 27      photonfs.addParticleId(PID::PHOTON);
 28      declare(photonfs, "LeadingPhoton");
 29
 30      size_t hist_bin = 0;
 31      for (size_t i = 0; i < _eta_bins.size()-1; ++i) {
 32        if (fabs(_eta_bins[i] - 1.37) < .0001) continue;
 33        book(_h_Et_photon[i] ,1+hist_bin, 1, 1);
 34        hist_bin += 1;
 35      }
 36    }
 37
 38
 39    size_t getEtaBin(double eta, bool area_eta) const {
 40      return (!area_eta) ? binIndex(fabs(eta), _eta_bins) : binIndex(fabs(eta), _eta_bins_areaoffset);
 41    }
 42
 43
 44    /// Perform the per-event analysis
 45    void analyze(const Event& event) {
 46      Particles photons = apply<LeadingParticlesFinalState>(event, "LeadingPhoton").particles();
 47      if (photons.size() != 1) vetoEvent;
 48
 49      const Particle& leadingPhoton = photons[0];
 50      if (inRange(leadingPhoton.abseta(), 1.37, 1.52)) vetoEvent;
 51      const int eta_bin = getEtaBin(leadingPhoton.abseta(), false);
 52
 53      const Particles& fs = apply<FinalState>(event, "FS").particles();
 54      FourMomentum mom_in_EtCone;
 55      for (const Particle& p : fs) {
 56        // Check if it's in the cone of .4
 57        if (deltaR(leadingPhoton, p) >= 0.4) continue;
 58        // Check if it's in the 5x7 central core
 59        if (fabs(deltaEta(leadingPhoton, p)) < .025*5.0*0.5 &&
 60            fabs(deltaPhi(leadingPhoton, p)) < (PI/128.)*7.0*0.5) continue;
 61        // Increment
 62        mom_in_EtCone += p.momentum();
 63      }
 64      MSG_DEBUG("Done with initial Et cone");
 65
 66      // Get the jet pT densities
 67      vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
 68      FastJets fastjets = apply<FastJets>(event, "KtJetsD05");
 69      const shared_ptr<fastjet::ClusterSequenceArea> clust_seq_area = fastjets.clusterSeqArea();
 70      for (const Jet& jet : fastjets.jets()) {
 71        const double area = clust_seq_area->area(jet); //< Implicit call to pseudojet()
 72        if (area > 1e-4 && jet.abseta() < _eta_bins_areaoffset.back()) {
 73          ptDensities.at(getEtaBin(jet.abseta(), true)) += jet.pT()/area;
 74        }
 75      }
 76
 77      // Now compute the median energy densities
 78      vector<double> ptDensity;
 79      for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
 80        ptDensity += ptDensities[b].empty() ? 0 : median(ptDensities[b]);
 81      }
 82
 83      // Now figure out the correction
 84      const double ETCONE_AREA = PI*.4*.4 - (7.0*.025)*(5.0*PI/128.);
 85      const double correction = ptDensity[getEtaBin(leadingPhoton.abseta(), true)]*ETCONE_AREA;
 86      MSG_DEBUG("Jet area correction done");
 87
 88      // Shouldn't need to subtract photon
 89      // NB. Using expected cut at hadron/particle level, not cut at reco level
 90      if (mom_in_EtCone.Et() - correction/*-leadingPhoton.Et()*/ > 4.0*GeV) vetoEvent;
 91      MSG_DEBUG("Passed isolation cut");
 92
 93      // Fill histogram
 94      _h_Et_photon[eta_bin]->fill(leadingPhoton.Et()/GeV);
 95    }
 96
 97
 98    /// Normalise histograms etc., after the run
 99    void finalize() {
100      for (size_t i = 0; i < _eta_bins.size()-1; ++i) {
101        if (fabs(_eta_bins[i] - 1.37) < .0001) continue;
102        scale(_h_Et_photon[i], crossSection()/picobarn/sumOfWeights());
103      }
104    }
105
106
107  private:
108
109    Histo1DPtr _h_Et_photon[6];
110
111    const vector<double> _eta_bins = {0.00, 0.60, 1.37, 1.52, 1.81};
112    const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0};
113
114  };
115
116
117
118  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2010_I882463, ATLAS_2010_S8914702);
119
120}