1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/VisibleFinalState.hh"
  5#include "Rivet/Projections/VetoedFinalState.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Isolated photon + jets at 13 TeV
 13  class ATLAS_2017_I1645627 : public Analysis {
 14  public:
 15
 16    // Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1645627);
 18
 19    // Book histograms and initialise projections before the run
 20    void init() {
 21      const FinalState fs;
 22
 23      // calorimeter particles
 24      VisibleFinalState visFS(fs);
 25      VetoedFinalState calo_fs(visFS);
 26      calo_fs.addVetoPairId(PID::MUON);
 27      declare(calo_fs, "calo");
 28
 29      // Voronoi eta-phi tessellation with KT jets, for ambient energy density calculation
 30      FastJets fj(fs, JetAlg::KT, 0.5, JetMuons::NONE, JetInvisibles::NONE); // E-scheme used by default;
 31      fj.useJetArea(new fastjet::AreaDefinition(fastjet::voronoi_area, fastjet::VoronoiAreaSpec(1.0)));
 32      declare(fj, "KtJetsD05");
 33
 34      // photon
 35      PromptFinalState photonfs(Cuts::abspid == PID::PHOTON && Cuts::abseta < 2.37 && Cuts::pT > 125*GeV);
 36      declare(photonfs, "photons");
 37
 38      // Jets
 39      FastJets jetpro(fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::NONE);
 40      declare(jetpro, "Jets");
 41
 42      // Histograms
 43      book(_h_photon_pt     , 1, 1, 1);
 44      book(_h_jet_pt        , 2, 1, 1);
 45      book(_h_phjet_dphi    , 3, 1, 1);
 46      book(_h_phjet_mass    , 4, 1, 1);
 47      book(_h_phjet_costheta, 5, 1, 1);
 48
 49    }
 50
 51
 52    size_t getEtaBin(double eta) const {
 53      return binIndex(fabs(eta), _eta_bins_areaoffset);
 54    }
 55
 56
 57    // Perform the per-event analysis
 58    void analyze(const Event& event) {
 59
 60      // Get the photon
 61      const Particles& photons = apply<PromptFinalState>(event, "photons").particlesByPt(Cuts::abseta < 1.37 || Cuts::abseta > 1.56);
 62      if (photons.empty())  vetoEvent;
 63      const FourMomentum photon = photons[0].momentum();
 64
 65      // Get the jet
 66      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 100*GeV && Cuts::absrap < 2.37);
 67      idiscard(jets, deltaRLess(photon, 0.8));
 68      if (jets.empty())  vetoEvent;
 69      FourMomentum leadingJet = jets[0].momentum();
 70
 71      // Compute the jet pT densities
 72      vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
 73      FastJets fastjets = apply<FastJets>(event, "KtJetsD05");
 74      const auto clust_seq_area = fastjets.clusterSeqArea();
 75      for (const Jet& jet : fastjets.jets()) {
 76        const double area = clust_seq_area->area(jet); // Implicit call to pseudojet().
 77        //const double area2 = (clust_seq_area->area_4vector(jet)).perp(); // Area definition used in egammaTruthParticles.
 78        if (area > 1e-3 && jet.abseta() < _eta_bins_areaoffset.back()) {
 79          ptDensities.at(getEtaBin(jet.abseta())) += jet.pT()/area;
 80        }
 81      }
 82
 83      // Compute the median event energy density
 84      vector<double> ptDensity;
 85      for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
 86        ptDensity += ptDensities[b].empty() ? 0 : median(ptDensities[b]);
 87      }
 88
 89      // Compute photon isolation with a standard ET cone
 90      FourMomentum mom_in_EtCone;
 91      const Particles calo_fs = apply<VetoedFinalState>(event, "calo").particles();
 92      const double iso_dr = 0.4;
 93      for (const Particle& p : calo_fs) {
 94        // Check if it's in the cone of .4
 95        if (sqrt(2.0*(cosh(p.eta()-photon.eta()) - cos(p.phi()-photon.phi()))) >= iso_dr) continue;
 96        // Increment sum
 97        mom_in_EtCone += p.momentum();
 98      }
 99
100      // Remove the photon energy from the isolation
101      mom_in_EtCone -= photon;
102
103      // Figure out the correction (area*density)
104      const double etcone_area = PI*iso_dr*iso_dr;
105      const double correction = ptDensity[getEtaBin(photon.abseta())] * etcone_area;
106      // Require photon to be isolated
107      if ((mom_in_EtCone.Et()-correction) > (0.0042*photon.pT() + 10*GeV))  vetoEvent;
108
109      // Fill histos
110      const double photon_pt = photon.pT()/GeV;
111      const double jet_pt = leadingJet.pT()/GeV;
112      const double phjet_dphi = deltaPhi(photon, leadingJet);
113      const double photon_eta = photon.eta();
114      const double jet_y = leadingJet.rapidity();
115      _h_photon_pt->fill(photon_pt);
116      _h_jet_pt->fill(jet_pt);
117      _h_phjet_dphi->fill(phjet_dphi);
118
119      double dy = fabs(jet_y-photon_eta);
120      double phjet_costheta = tanh(dy/2.);
121      double phjet_mass= (photon+leadingJet).mass()/GeV;
122      if (phjet_mass <= 450.)  vetoEvent;
123      if (fabs(photon_eta + jet_y) >= 2.37)  vetoEvent;
124      if (phjet_costheta >= 0.83)  vetoEvent;
125      _h_phjet_costheta->fill(phjet_costheta);
126      _h_phjet_mass->fill(phjet_mass);
127    }
128
129
130    /// Normalise histograms etc., after the run
131    void finalize() {
132      const double sf = crossSection()/picobarn / sumOfWeights();
133      scale(_h_photon_pt, sf);
134      scale(_h_jet_pt, sf);
135      scale(_h_phjet_dphi, sf);
136      scale(_h_phjet_mass, sf);
137      scale(_h_phjet_costheta, sf);
138    }
139
140
141  private:
142
143    Histo1DPtr _h_photon_pt, _h_jet_pt;
144    Histo1DPtr _h_phjet_dphi, _h_phjet_mass, _h_phjet_costheta;
145
146    const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0, 4.0, 5.0};
147
148  };
149
150  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1645627);
151
152}