1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/LeadingParticlesFinalState.hh"
  5#include "Rivet/Projections/VetoedFinalState.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Measurement of isolated gamma + jet + X differential cross-sections
 12  ///
 13  /// Inclusive isolated gamma + jet cross-sections, differential in pT(gamma), for
 14  /// various photon and jet rapidity configurations.
 15  ///
 16  /// @author Giovanni Marchiori
 17  class ATLAS_2012_I1093738 : public Analysis {
 18  public:
 19
 20    // Constructor
 21    ATLAS_2012_I1093738()
 22      : Analysis("ATLAS_2012_I1093738")
 23    {    }
 24
 25
 26    // Book histograms and initialise projections before the run
 27    void init() {
 28      // Final state
 29      FinalState fs;
 30      declare(fs, "FS");
 31
 32      // Voronoi eta-phi tessellation with KT jets, for ambient energy density calculation
 33      FastJets fj(fs, JetAlg::KT, 0.5);
 34      fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
 35      declare(fj, "KtJetsD05");
 36
 37      // Leading photon
 38      LeadingParticlesFinalState photonfs(FinalState((Cuts::etaIn(-1.37, 1.37) && Cuts::pT >=  25.0*GeV)));
 39      photonfs.addParticleId(PID::PHOTON);
 40      declare(photonfs, "LeadingPhoton");
 41
 42      // FS excluding the leading photon
 43      VetoedFinalState vfs(fs);
 44      vfs.addVetoOnThisFinalState(photonfs);
 45      declare(vfs, "JetFS");
 46
 47      // Jets
 48      FastJets jetpro(vfs, JetAlg::ANTIKT, 0.4);
 49      jetpro.useInvisibles();
 50      declare(jetpro, "Jets");
 51
 52      book(_h_phbarrel_jetcentral_SS ,1, 1, 1);
 53      book(_h_phbarrel_jetmedium_SS  ,2, 1, 1);
 54      book(_h_phbarrel_jetforward_SS ,3, 1, 1);
 55
 56      book(_h_phbarrel_jetcentral_OS ,4, 1, 1);
 57      book(_h_phbarrel_jetmedium_OS  ,5, 1, 1);
 58      book(_h_phbarrel_jetforward_OS ,6, 1, 1);
 59    }
 60
 61
 62    int getEtaBin(double eta, int what) const {
 63      const double aeta = fabs(eta);
 64      if (what == 0) return binIndex(aeta, _eta_bins_ph);
 65      if (what == 1) return binIndex(aeta, _eta_bins_jet);
 66      return binIndex(aeta, _eta_bins_areaoffset);
 67    }
 68
 69
 70    // Perform the per-event analysis
 71    void analyze(const Event& event) {
 72
 73      // Get the photon
 74      const FinalState& photonfs = apply<FinalState>(event, "LeadingPhoton");
 75      if (photonfs.particles().size() < 1) vetoEvent;
 76      const FourMomentum photon = photonfs.particles().front().momentum();
 77
 78      // Get the jet
 79      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 20*GeV);
 80      if (jets.empty()) vetoEvent;
 81      FourMomentum leadingJet = jets[0].momentum();
 82
 83      // Require jet separated from photon
 84      if (deltaR(photon, leadingJet) < 1.0) vetoEvent;
 85
 86      // Veto if leading jet is outside plotted rapidity regions
 87      if (leadingJet.absrap() > 4.4) vetoEvent;
 88
 89      // Compute the jet pT densities
 90      vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
 91      FastJets fastjets = apply<FastJets>(event, "KtJetsD05");
 92      const shared_ptr<fastjet::ClusterSequenceArea> clust_seq_area = fastjets.clusterSeqArea();
 93      for (const Jet& jet : fastjets.jets()) {
 94        const double area = clust_seq_area->area(jet); //< Implicit call to pseudojet()
 95        if (area > 1e-4 && jet.abseta() < _eta_bins_areaoffset.back()) {
 96          ptDensities.at(getEtaBin(jet.abseta(), 2)) += jet.pT()/area;
 97        }
 98      }
 99
100      // Compute the median event energy density
101      /// @todo This looks equivalent to median(ptDensities[b]) -- isn't SKIPNHARDJETS meant to be used as an offset?
102      const unsigned int SKIPNHARDJETS = 0;
103      vector<double> ptDensity;
104      for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; b++) {
105        double median = 0.0;
106        if (ptDensities[b].size() > SKIPNHARDJETS) {
107          std::sort(ptDensities[b].begin(), ptDensities[b].end());
108          const int nDens = ptDensities[b].size() - SKIPNHARDJETS;
109          if (nDens % 2 == 0) {
110            median = (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2;
111          } else {
112            median = ptDensities[b][(nDens-1)/2];
113          }
114        }
115        ptDensity.push_back(median);
116      }
117
118      // Compute photon isolation with a standard ET cone
119      const Particles fs = apply<FinalState>(event, "JetFS").particles();
120      FourMomentum mom_in_EtCone;
121      const double ISO_DR = 0.4;
122      const double CLUSTER_ETA_WIDTH = 0.25*5.0;
123      const double CLUSTER_PHI_WIDTH = (PI/128.)*7.0;
124      for (const Particle& p : fs) {
125        // Check if it's in the cone of .4
126        if (deltaR(photon, p) >= ISO_DR) continue;
127        // Check if it's in the 5x7 central core
128        if (fabs(deltaEta(photon, p)) < CLUSTER_ETA_WIDTH*0.5 &&
129            fabs(deltaPhi(photon, p)) < CLUSTER_PHI_WIDTH*0.5) continue;
130        // Increment sum
131        mom_in_EtCone += p.momentum();
132      }
133
134      // Figure out the correction (area*density)
135      const double ETCONE_AREA = PI*ISO_DR*ISO_DR - CLUSTER_ETA_WIDTH*CLUSTER_PHI_WIDTH;
136      const double correction = ptDensity[getEtaBin(photon.abseta(),2)] * ETCONE_AREA;
137
138      // Require photon to be isolated
139      if (mom_in_EtCone.Et()-correction > 4.0*GeV) vetoEvent;
140
141      const int photon_jet_sign = sign( leadingJet.rapidity() * photon.rapidity() );
142
143      // Fill histos
144      const double abs_jet_rapidity = fabs(leadingJet.rapidity());
145      const double photon_pt = photon.pT()/GeV;
146      const double abs_photon_eta = fabs(photon.eta());
147      if (abs_photon_eta < 1.37) {
148        if (abs_jet_rapidity < 1.2) {
149          if (photon_jet_sign >= 1) {
150            _h_phbarrel_jetcentral_SS->fill(photon_pt);
151          } else {
152            _h_phbarrel_jetcentral_OS->fill(photon_pt);
153          }
154        } else if (abs_jet_rapidity < 2.8) {
155          if (photon_jet_sign >= 1) {
156            _h_phbarrel_jetmedium_SS->fill(photon_pt);
157          } else {
158            _h_phbarrel_jetmedium_OS->fill(photon_pt);
159          }
160        } else if (abs_jet_rapidity < 4.4) {
161          if (photon_jet_sign >= 1) {
162            _h_phbarrel_jetforward_SS->fill(photon_pt);
163          } else {
164            _h_phbarrel_jetforward_OS->fill(photon_pt);
165          }
166        }
167      }
168
169    }
170
171
172    /// Normalise histograms etc., after the run
173    void finalize() {
174      scale(_h_phbarrel_jetcentral_SS, crossSection()/picobarn/sumOfWeights());
175      scale(_h_phbarrel_jetcentral_OS, crossSection()/picobarn/sumOfWeights());
176      scale(_h_phbarrel_jetmedium_SS, crossSection()/picobarn/sumOfWeights());
177      scale(_h_phbarrel_jetmedium_OS, crossSection()/picobarn/sumOfWeights());
178      scale(_h_phbarrel_jetforward_SS, crossSection()/picobarn/sumOfWeights());
179      scale(_h_phbarrel_jetforward_OS, crossSection()/picobarn/sumOfWeights());
180    }
181
182
183  private:
184
185    Histo1DPtr _h_phbarrel_jetcentral_SS, _h_phbarrel_jetmedium_SS, _h_phbarrel_jetforward_SS;
186    Histo1DPtr _h_phbarrel_jetcentral_OS, _h_phbarrel_jetmedium_OS, _h_phbarrel_jetforward_OS;
187
188    const vector<double> _eta_bins_ph = {0.0, 1.37, 1.52, 2.37};
189    const vector<double> _eta_bins_jet = {0.0, 1.2, 2.8, 4.4};
190    const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0};
191
192  };
193
194
195  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1093738);
196
197
198}