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/VetoedFinalState.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Measurement of isolated gamma + jet + X differential cross-sections
 13  class ATLAS_2013_I1244522 : public Analysis {
 14  public:
 15
 16    // Constructor
 17    ATLAS_2013_I1244522()
 18      : Analysis("ATLAS_2013_I1244522")
 19    {     }
 20
 21
 22    // Book histograms and initialise projections before the run
 23    void init() {
 24      FinalState fs;
 25
 26      // Voronoi eta-phi tassellation with KT jets, for ambient energy density calculation
 27      FastJets fj(fs, JetAlg::KT, 0.5);
 28      fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
 29      declare(fj, "KtJetsD05");
 30
 31      // Leading photon
 32      LeadingParticlesFinalState photonfs(PromptFinalState(FinalState((Cuts::etaIn(-2.37, 2.37) && Cuts::pT >=  45.0*GeV))));
 33      photonfs.addParticleId(PID::PHOTON);
 34      declare(photonfs, "LeadingPhoton");
 35
 36      // FS excluding the leading photon
 37      VetoedFinalState vfs(fs);
 38      vfs.addVetoOnThisFinalState(photonfs);
 39      declare(vfs, "JetFS");
 40
 41      // Jets
 42      FastJets jetpro(vfs, JetAlg::ANTIKT, 0.6);
 43      jetpro.useInvisibles();
 44      declare(jetpro, "Jets");
 45
 46      // Histograms
 47      book(_h_ph_pt      ,1, 1, 1);
 48      book(_h_jet_pt     ,2, 1, 1);
 49      book(_h_jet_rap    ,3, 1, 1);
 50      book(_h_dphi_phjet ,4, 1, 1);
 51      book(_h_costheta_biased_phjet ,5, 1, 1);
 52      book(_h_mass_phjet            ,6, 1, 1);
 53      book(_h_costheta_phjet        ,7, 1, 1);
 54
 55    }
 56
 57
 58    size_t getEtaBin(double eta) const {
 59      const double aeta = fabs(eta);
 60      return binIndex(aeta, _eta_bins_areaoffset);
 61    }
 62
 63
 64    // Perform the per-event analysis
 65    void analyze(const Event& event) {
 66
 67      // Get the photon
 68      Particles photons = apply<LeadingParticlesFinalState>(event, "LeadingPhoton").particles();
 69      if (photons.size() != 1 )  vetoEvent;
 70      const Particle& photon = photons[0];
 71
 72      if (inRange(photon.abseta(), 1.37, 1.52))  vetoEvent;
 73
 74      //Compute isolation energy in cone of radius .4 around photon (all particles)
 75      FourMomentum mom_in_EtCone;
 76      const Particles& fs = apply<VetoedFinalState>(event, "JetFS").particles();
 77      for (const Particle& p : fs) {
 78        // Check if it's outside the cone of 0.4
 79        if (deltaR(photon, p) >= 0.4) continue;
 80        // Increment isolation energy
 81        mom_in_EtCone += p.momentum();
 82      }
 83
 84      // Get the jets
 85      Jets alljets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 40*GeV);
 86      Jets jets;
 87      for (const Jet& jet : alljets)
 88        if (deltaR(photon, jet) > 1.0) jets += jet;
 89      if (jets.empty())  vetoEvent;
 90      Jet leadingJet = jets[0];
 91      if (leadingJet.absrap() > 2.37) vetoEvent;
 92
 93      // Get the area-filtered jet inputs for computing median energy density, etc.
 94      vector<double> ptDensity;
 95      vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
 96      FastJets fast_jets = apply<FastJets>(event, "KtJetsD05");
 97      const auto clust_seq_area = fast_jets.clusterSeqArea();
 98      for (const Jet& jet : fast_jets.jets()) {
 99        const double area = clust_seq_area->area(jet);
100        if (area > 1e-4 && jet.abseta() < _eta_bins_areaoffset.back())
101          ptDensities.at( getEtaBin(jet.abseta()) ).push_back(jet.pT()/area);
102      }
103
104      // Compute the median energy density, etc.
105      for (size_t b = 0; b < _eta_bins_areaoffset.size() - 1; ++b) {
106        const int njets = ptDensities[b].size();
107        ptDensity += (njets > 0) ? median(ptDensities[b]) : 0;
108      }
109
110      // Compute the isolation energy correction (cone area*energy density)
111      const double etCone_area = PI*sqr(0.4) - (5.0*.025)*(7.0*PI/128.);
112      const double correction = ptDensity[getEtaBin(photon.abseta())] * etCone_area;
113
114      // Apply isolation cut on area-corrected value
115      if (mom_in_EtCone.Et() - correction >= 4*GeV)  vetoEvent;
116
117      // Fill histos
118      const double dy = deltaRap(photon, leadingJet);
119      const double costheta_yj = tanh(dy/2);
120      _h_ph_pt->fill(photon.pT()/GeV);
121      _h_jet_pt->fill(leadingJet.pT()/GeV);
122      _h_jet_rap->fill(leadingJet.absrap());
123      _h_dphi_phjet->fill(deltaPhi(photon, leadingJet));
124      _h_costheta_biased_phjet->fill(costheta_yj);
125      if (costheta_yj < 0.829022) {
126        const FourMomentum yj = photon.momentum() + leadingJet.momentum();
127        if (yj.mass() > 160.939*GeV) {
128          if (fabs(photon.eta() + leadingJet.rap()) < 2.37) {
129            _h_mass_phjet->fill(yj.mass()/GeV);
130            _h_costheta_phjet->fill(costheta_yj);
131          }
132        }
133      }
134    }
135
136
137    /// Normalise histograms etc., after the run
138    void finalize() {
139      const double sf = crossSection() / picobarn / sumOfWeights();
140      scale(_h_ph_pt,                 sf);
141      scale(_h_jet_pt,                sf);
142      scale(_h_jet_rap,               sf);
143      scale(_h_dphi_phjet,            sf);
144      scale(_h_costheta_biased_phjet, sf);
145      scale(_h_mass_phjet,            sf);
146      scale(_h_costheta_phjet,        sf);
147    }
148
149
150  private:
151
152    Histo1DPtr _h_ph_pt, _h_jet_pt, _h_jet_rap, _h_dphi_phjet, _h_costheta_biased_phjet, _h_mass_phjet, _h_costheta_phjet;
153
154    const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0};
155
156  };
157
158
159  RIVET_DECLARE_PLUGIN(ATLAS_2013_I1244522);
160
161
162}