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 + 2 jets at 13 TeV
 13  class ATLAS_2019_I1772071 : public Analysis {
 14  public:
 15
 16    // Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2019_I1772071);
 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 > 150*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
 43      vector<string> observables = {"ETGamma", "pTjet", "RapJet","DeltaRapGammaJet",
 44                                    "DeltaPhiGammaJet", "DeltaRapJetJet", "DeltaPhiJetJet",
 45                                    "MassJetJet", "MassGammaJetJet"};
 46      vector<string> regions = {"Inclusive","Fragmentation", "Direct"};
 47
 48      int i=0;
 49      for (const string& region : regions){
 50        int j = 1;
 51	      for (const string& name : observables) {
 52	        book(_h[name+region], 9*i+j, 1, 1);
 53	        ++j;
 54	      }
 55	      ++i;
 56      }
 57    }
 58
 59
 60    size_t getEtaBin(double eta) const {
 61      return binIndex(fabs(eta), _eta_bins_areaoffset);
 62    }
 63
 64
 65    // Perform the per-event analysis
 66    void analyze(const Event& event) {
 67
 68      // Get the photon
 69      const Particles& photons = apply<PromptFinalState>(event, "photons").particlesByPt(Cuts::abseta < 1.37 || Cuts::abseta > 1.56);
 70      if (photons.empty())  vetoEvent;
 71      const FourMomentum photon = photons[0].momentum();
 72
 73      // Get the jet
 74      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 100*GeV && Cuts::absrap < 2.5);
 75      idiscard(jets, deltaRLess(photon, 0.8));
 76      if ( jets.size()<2 )  vetoEvent;
 77      FourMomentum leadingJet = jets[0].momentum();
 78      FourMomentum subleadingJet = jets[1].momentum();
 79
 80      // Compute the jet pT densities
 81      vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
 82      FastJets fastjets = apply<FastJets>(event, "KtJetsD05");
 83      const auto clust_seq_area = fastjets.clusterSeqArea();
 84      for (const Jet& jet : fastjets.jets()) {
 85        const double area = clust_seq_area->area(jet); // Implicit call to pseudojet().
 86        //const double area2 = (clust_seq_area->area_4vector(jet)).perp(); // Area definition used in egammaTruthParticles.
 87        if (area > 1e-3 && jet.abseta() < _eta_bins_areaoffset.back()) {
 88          ptDensities.at(getEtaBin(jet.abseta())) += jet.pT()/area;
 89        }
 90      }
 91
 92      // Compute the median event energy density
 93      vector<double> ptDensity;
 94      for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
 95        ptDensity += ptDensities[b].empty() ? 0 : median(ptDensities[b]);
 96      }
 97
 98      // Compute photon isolation with a standard ET cone
 99      FourMomentum mom_in_EtCone;
100      const Particles calo_fs = apply<VetoedFinalState>(event, "calo").particles();
101      const double iso_dr = 0.4;
102      for (const Particle& p : calo_fs) {
103        // Check if it's in the cone of .4
104        if (sqrt(2.0*(cosh(p.eta()-photon.eta()) - cos(p.phi()-photon.phi()))) >= iso_dr) continue;
105        // Increment sum
106        mom_in_EtCone += p.momentum();
107      }
108
109      // Remove the photon energy from the isolation
110      mom_in_EtCone -= photon;
111
112      // Figure out the correction (area*density)
113      const double etcone_area = PI*iso_dr*iso_dr;
114      const double correction = ptDensity[getEtaBin(photon.abseta())] * etcone_area;
115      // Require photon to be isolated
116      if ((mom_in_EtCone.Et()-correction) > (0.0042*photon.pT() + 10*GeV))  vetoEvent;
117
118      // Fill histos
119      const double photon_pt = photon.pT()/GeV;
120      const double jet_pt1 = leadingJet.pT()/GeV;
121      const double jet_pt2 = subleadingJet.pT()/GeV;
122      const double jet1_y = leadingJet.rapidity();
123      const double jet2_y = subleadingJet.rapidity();
124      const double phjet1_dphi = deltaPhi(photon, leadingJet);
125      const double phjet2_dphi = deltaPhi(photon, subleadingJet);
126      const double phjet1_drap = fabs(photon.eta()-leadingJet.rapidity());
127      const double phjet2_drap = fabs(photon.eta()-subleadingJet.rapidity());
128      const double jetjet_drap = fabs(leadingJet.rapidity()-subleadingJet.rapidity());
129      const FourMomentum jetjet = leadingJet+subleadingJet;
130      const double mjetjet = jetjet.mass()/GeV;
131      const FourMomentum phjet1 = photon+leadingJet;
132      const FourMomentum phjet2 = photon+subleadingJet;
133      const FourMomentum phjetjet = photon+leadingJet+subleadingJet;
134      const double mphjetjet = phjetjet.mass()/GeV;
135      const double jetjet_dphi = deltaPhi(subleadingJet, leadingJet);
136
137      _h["ETGammaInclusive"]->fill(photon_pt);
138      _h["pTjetInclusive"]->fill(jet_pt1);
139      _h["pTjetInclusive"]->fill(jet_pt2);
140      _h["RapJetInclusive"]->fill(fabs(jet1_y));
141      _h["RapJetInclusive"]->fill(fabs(jet2_y));
142      _h["DeltaRapGammaJetInclusive"]->fill(phjet1_drap);
143      _h["DeltaRapGammaJetInclusive"]->fill(phjet2_drap);
144      _h["DeltaPhiGammaJetInclusive"]->fill(phjet1_dphi);
145      _h["DeltaPhiGammaJetInclusive"]->fill(phjet2_dphi);
146      _h["MassJetJetInclusive"]->fill(mjetjet);
147      _h["DeltaPhiJetJetInclusive"]->fill(jetjet_dphi);
148      _h["DeltaRapJetJetInclusive"]->fill(jetjet_drap);
149      _h["MassGammaJetJetInclusive"]->fill(mphjetjet);
150
151      if (photon_pt>jet_pt1) {
152        _h["ETGammaDirect"]->fill(photon_pt);
153        _h["pTjetDirect"]->fill(jet_pt1);
154        _h["pTjetDirect"]->fill(jet_pt2);
155        _h["RapJetDirect"]->fill(fabs(jet1_y));
156        _h["RapJetDirect"]->fill(fabs(jet2_y));
157        _h["DeltaRapGammaJetDirect"]->fill(phjet1_drap);
158        _h["DeltaRapGammaJetDirect"]->fill(phjet2_drap);
159        _h["DeltaPhiGammaJetDirect"]->fill(phjet1_dphi);
160        _h["DeltaPhiGammaJetDirect"]->fill(phjet2_dphi);
161        _h["MassJetJetDirect"]->fill(mjetjet);
162        _h["DeltaPhiJetJetDirect"]->fill(jetjet_dphi);
163        _h["DeltaRapJetJetDirect"]->fill(jetjet_drap);
164        _h["MassGammaJetJetDirect"]->fill(mphjetjet);
165
166      }
167      else if (photon_pt < jet_pt2)	{
168        _h["ETGammaFragmentation"]->fill(photon_pt);
169        _h["pTjetFragmentation"]->fill(jet_pt1);
170        _h["pTjetFragmentation"]->fill(jet_pt2);
171        _h["RapJetFragmentation"]->fill(fabs(jet1_y));
172        _h["RapJetFragmentation"]->fill(fabs(jet2_y));
173        _h["DeltaRapGammaJetFragmentation"]->fill(phjet1_drap);
174        _h["DeltaRapGammaJetFragmentation"]->fill(phjet2_drap);
175        _h["DeltaPhiGammaJetFragmentation"]->fill(phjet1_dphi);
176        _h["DeltaPhiGammaJetFragmentation"]->fill(phjet2_dphi);
177        _h["MassJetJetFragmentation"]->fill(mjetjet);
178        _h["DeltaPhiJetJetFragmentation"]->fill(jetjet_dphi);
179        _h["DeltaRapJetJetFragmentation"]->fill(jetjet_drap);
180        _h["MassGammaJetJetFragmentation"]->fill(mphjetjet);
181      }
182    }
183
184
185    /// Normalise histograms etc., after the run
186    void finalize() {
187      const double sf = crossSection() / picobarn / sumOfWeights();
188	    scale(_h, sf);
189    }
190
191
192  private:
193
194    map<string,Histo1DPtr> _h;
195    const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0, 4.0, 5.0};
196
197  };
198
199  RIVET_DECLARE_PLUGIN(ATLAS_2019_I1772071);
200
201}