Rivet analyses referenceATLAS_2016_I1457605Inclusive prompt photons at 8 TeVExperiment: ATLAS (LHC) Inspire ID: 1457605 Status: VALIDATED Authors:
Beam energies: (4000.0, 4000.0) GeV Run details:
A measurement of the cross section for the inclusive production of isolated prompt photons in proton-proton collisions at a centre-of-mass energy of $\sqrt{s} = 8$ TeV is presented. The measurement covers the pseudorapidity ranges $|\eta^\gamma| < 1.37$ $1.56 < |\eta^\gamma| < 2.37$ in the transverse energy range $25 < E^\gamma_\text{T} < 1500$ GeV. The results are based on an integrated luminosity of 20.2 fb${}^{-1}$, recorded by the ATLAS detector at the LHC. Photon candidates are identified by combining information from the calorimeters and the inner tracker. The background is subtracted using a data-driven technique, based on the observed calorimeter shower-shape variables and the deposition of hadronic energy in a narrow cone around the photon candidate. Source code: ATLAS_2016_I1457605.cc 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/FastJets.hh"
7
8namespace Rivet {
9
10
11 /// Inclusive isolated prompt photon analysis with 2012 LHC data
12 class ATLAS_2016_I1457605 : public Analysis {
13 public:
14
15 /// Constructor
16 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2016_I1457605);
17
18 /// Book histograms and initialise projections before the run
19 void init() {
20
21 FinalState fs;
22 declare(fs, "FS");
23
24 // Consider the final state jets for the energy density calculation
25 FastJets fj(fs, JetAlg::KT, 0.5);
26 fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
27 declare(fj, "KtJetsD05");
28
29 // Consider the leading pt photon with |eta| < 2.37 and pT > 25 GeV
30 LeadingParticlesFinalState photonfs(PromptFinalState(FinalState(Cuts::abseta < 2.37 && Cuts::pT > 25*GeV)));
31 photonfs.addParticleId(PID::PHOTON);
32 declare(photonfs, "LeadingPhoton");
33
34 // Book the dsigma/dEt (in eta bins) histograms
35 for (size_t i = 0; i < _eta_bins.size() - 1; ++i) {
36 if (fuzzyEquals(_eta_bins[i], 1.37)) continue; // skip this bin
37 int offset = i > 2? 0 : 1;
38 book(_h_Et_photon[i] ,i + offset, 1, 1);
39 }
40
41 }
42
43
44 /// Return eta bin for either dsigma/dET histogram (area_eta=false) or energy density correction (area_eta=true)
45 size_t _getEtaBin(double eta_w, bool area_eta) const {
46 const double eta = fabs(eta_w);
47 if (!area_eta) {
48 return binIndex(eta, _eta_bins);
49 } else {
50 return binIndex(eta, _eta_bins_areaoffset);
51 }
52 }
53
54
55 /// Perform the per-event analysis
56 void analyze(const Event& event) {
57 // Retrieve leading photon
58 Particles photons = apply<LeadingParticlesFinalState>(event, "LeadingPhoton").particles();
59 if (photons.size() < 1) vetoEvent;
60 const Particle& leadingPhoton = photons[0];
61
62 // Veto events with photon in ECAL crack
63 if (inRange(leadingPhoton.abseta(), 1.37, 1.56)) vetoEvent;
64
65 // Compute isolation energy in cone of radius .4 around photon (all particles)
66 FourMomentum mom_in_EtCone;
67 Particles fs = apply<FinalState>(event, "FS").particles();
68 for (const Particle& p : fs) {
69 // Check if it's outside the cone of 0.4
70 if (deltaR(leadingPhoton, p) >= 0.4) continue;
71 // Except muons or neutrinos
72 if (PID::isNeutrino(p.abspid()) || p.abspid() == PID::MUON) continue;
73 // Increment isolation energy
74 mom_in_EtCone += p.momentum();
75 }
76 // Remove the photon energy from the isolation
77 mom_in_EtCone -= leadingPhoton.momentum();
78
79 // Get the area-filtered jet inputs for computing median energy density, etc.
80 vector<double> ptDensity;
81 vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
82 const FastJets& fast_jets = apply<FastJets>(event, "KtJetsD05");
83 const auto clust_seq_area = fast_jets.clusterSeqArea();
84 for (const Jet& jet : fast_jets.jets()) {
85 const double area = clust_seq_area->area(jet);
86 if (area > 1e-3 && jet.abseta() < _eta_bins_areaoffset.back())
87 ptDensities.at( _getEtaBin(jet.abseta(), true) ) += jet.pT()/area;
88 }
89 // Compute the median energy density, etc.
90 for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
91 const int njets = ptDensities[b].size();
92 ptDensity += (njets > 0) ? median(ptDensities[b]) : 0.0;
93 }
94 // Compute the isolation energy correction (cone area*energy density)
95 const double etCone_area = PI * sqr(0.4);
96 const double correction = ptDensity[_getEtaBin(leadingPhoton.abseta(), true)] * etCone_area;
97
98 // Apply isolation cut on area-corrected value
99 // cut is Etiso < 4.8GeV + 4.2E-03 * Et_gamma.
100 if (mom_in_EtCone.Et() - correction > 4.8*GeV + 0.0042*leadingPhoton.Et()) vetoEvent;
101
102 // Fill histograms
103 const size_t eta_bin = _getEtaBin(leadingPhoton.abseta(), false);
104 _h_Et_photon[eta_bin]->fill(leadingPhoton.Et());
105 }
106
107
108 /// Normalise histograms etc., after the run
109 void finalize() {
110 double sf = crossSection() / (picobarn * sumOfWeights());
111 for (size_t i = 0; i < _eta_bins.size()-1; ++i) {
112 if (fuzzyEquals(_eta_bins[i], 1.37)) continue;
113 scale(_h_Et_photon[i], sf);
114 }
115 }
116
117
118 private:
119
120 Histo1DPtr _h_Et_photon[5];
121
122 const vector<double> _eta_bins = {0.00, 0.60, 1.37, 1.56, 1.81, 2.37 };
123 const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0};
124
125 };
126
127
128 RIVET_DECLARE_PLUGIN(ATLAS_2016_I1457605);
129
130}
|