Rivet analyses referenceATLAS_2017_I1591327Inclusive diphoton cross-sections at 8 TeVExperiment: ATLAS (LHC) Inspire ID: 1591327 Status: VALIDATED Authors:
Beam energies: (4000.0, 4000.0) GeV Run details:
A measurement of the production cross section for two isolated photons in proton--proton collisions at a center-of-mass energy of $\sqrt{s}=8$ TeV is presented. The results are based on an integrated luminosity of 20.2 fb$^{-1}$ recorded by the ATLAS detector at the Large Hadron Collider. The measurement considers photons with pseudorapidities satisfying $|\eta^\gamma| < 1.37$ or $1.56 < |\eta^\gamma| < 2.37$ and transverse energies of respectively $E^\gamma_\text{T,1}>40$ GeV and $E^\gamma_\text{T,2}> 30$ GeV for the two leading photons ordered in transverse energy produced in the interaction. The background due to hadronic jets and electrons is subtracted using data-driven techniques. The fiducial cross sections are corrected for detector effects and measured differentially as a function of six kinematic observables. The measured cross section integrated within the fiducial volume is $16.8\pm0.8$ pb. The data are compared to fixed-order QCD calculations at next-to-leading-order and next-to-next-to-leading-order accuracy as well as next-to-leading-order computations including resummation of initial-state gluon radiation at next-to-next-to-leading logarithm or matched to a parton shower, with relative uncertainties varying from 5% to 20%. Source code: ATLAS_2017_I1591327.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/IdentifiedFinalState.hh"
5#include "Rivet/Projections/FastJets.hh"
6
7namespace Rivet {
8
9
10 /// Isolated diphoton + X differential cross-sections
11 class ATLAS_2017_I1591327 : public Analysis {
12 public:
13
14 // Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1591327);
16
17
18 // Book histograms and initialise projections before the run
19 void init() {
20
21 const FinalState fs;
22 declare(fs, "FS");
23
24 FastJets fj(fs, JetAlg::KT, 0.5);
25 _area_def = new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec());
26 fj.useJetArea(_area_def);
27 declare(fj, "KtJetsD05");
28
29 IdentifiedFinalState photonfs(Cuts::abseta < 2.37 && Cuts::pT > 30*GeV);
30 photonfs.acceptId(PID::PHOTON);
31 declare(photonfs, "Photon");
32
33 // Histograms
34 book(_h_M , 2, 1, 1);
35 book(_h_pT , 3, 1, 1);
36 book(_h_at , 4, 1, 1);
37 book(_h_phistar, 5, 1, 1);
38 book(_h_costh , 6, 1, 1);
39 book(_h_dPhi , 7, 1, 1);
40 }
41
42
43 // Perform the per-event analysis
44 void analyze(const Event& event) {
45
46 // Require at least 2 photons in final state
47 const Particles photons = apply<IdentifiedFinalState>(event, "Photon").particlesByPt();
48 if (photons.size() < 2) vetoEvent;
49
50 // Compute the median energy density
51 _ptDensity.clear();
52 _sigma.clear();
53 _Njets.clear();
54 vector<vector<double> > ptDensities;
55 vector<double> emptyVec;
56 ptDensities.assign(ETA_BINS.size()-1, emptyVec);
57
58 // Get jets, and corresponding jet areas
59 const shared_ptr<fastjet::ClusterSequenceArea> clust_seq_area = apply<FastJets>(event, "KtJetsD05").clusterSeqArea();
60 for (const fastjet::PseudoJet& jet : apply<FastJets>(event, "KtJetsD05").pseudojets(0.0*GeV)) {
61 const double aeta = fabs(jet.eta());
62 const double pt = jet.perp();
63 const double area = clust_seq_area->area(jet);
64 if (area < 1e-3) continue;
65 const int ieta = binIndex(aeta, ETA_BINS);
66 if (ieta != -1) ptDensities[ieta].push_back(pt/area);
67 }
68
69 // Compute median jet properties over the jets in the event
70 for (size_t b = 0; b < ETA_BINS.size()-1; ++b) {
71 double median = 0.0, sigma = 0.0;
72 int Njets = 0;
73 if (ptDensities[b].size() > 0) {
74 std::sort(ptDensities[b].begin(), ptDensities[b].end());
75 int nDens = ptDensities[b].size();
76 median = (nDens % 2 == 0) ? (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2 : ptDensities[b][(nDens-1)/2];
77 sigma = ptDensities[b][(int)(.15865*nDens)];
78 Njets = nDens;
79 }
80 _ptDensity.push_back(median);
81 _sigma.push_back(sigma);
82 _Njets.push_back(Njets);
83 }
84
85 // Loop over photons and fill vector of isolated ones
86 Particles isolated_photons;
87 for (const Particle& photon : photons) {
88 // Check if it's a prompt photon (needed for SHERPA 2->5 sample, otherwise I also get photons from hadron decays in jets)
89 if (!photon.isPrompt()) continue;
90
91 // Remove photons in ECAL crack region
92 if (inRange(photon.abseta(), 1.37, 1.56)) continue;
93 const double eta_P = photon.eta();
94 const double phi_P = photon.phi();
95
96 // Compute isolation via particles within an R=0.4 cone of the photon
97 const Particles fs = apply<FinalState>(event, "FS").particles();
98 FourMomentum mom_in_EtCone;
99 for (const Particle& p : fs) {
100 // Reject if not in cone
101 if (deltaR(photon.momentum(), p.momentum()) > 0.4) continue;
102 // Reject if in the 5x7 cell central core
103 if (fabs(eta_P - p.eta()) < 0.025 * 5 * 0.5 &&
104 fabs(phi_P - p.phi()) < PI/128. * 7 * 0.5) continue;
105 // Sum momentum
106 mom_in_EtCone += p.momentum();
107 }
108 // Now figure out the correction (area*density)
109 const double EtCone_area = M_PI*sqr(0.4) - (7*.025)*(5*M_PI/128.); // cone area - central core rectangle
110 const double correction = _ptDensity[binIndex(fabs(eta_P), ETA_BINS)] * EtCone_area;
111
112 // Discard the photon if there is more than 11 GeV of cone activity
113 // NOTE: Shouldn't need to subtract photon itself (it's in the central core)
114 if (mom_in_EtCone.Et() - correction > 11*GeV) continue;
115 // Add isolated photon to list
116 isolated_photons.push_back(photon);
117 }
118
119 // Require at least two isolated photons
120 if (isolated_photons.size() < 2) vetoEvent;
121
122 // Select leading pT pair
123 isortByPt(isolated_photons);
124 const FourMomentum y1 = isolated_photons[0];
125 const FourMomentum y2 = isolated_photons[1];
126
127 // Leading photon should have pT > 40 GeV, subleading > 30 GeV
128 if (y1.pT() < 40.*GeV) vetoEvent;
129 if (y2.pT() < 30.*GeV) vetoEvent;
130
131 // Require the two photons to be separated (dR>0.4)
132 if (deltaR(y1,y2) < 0.4) vetoEvent;
133
134 const FourMomentum yy = y1 + y2;
135 const double Myy = yy.mass();
136 const double pTyy = yy.pT();
137 const double dPhiyy = mapAngle0ToPi(y1.phi() - y2.phi());
138
139 // phi*
140 const double costhetastar_ = fabs(tanh(( y1.eta() - y2.eta() ) / 2.));
141 const double sinthetastar_ = sqrt(1. - pow(costhetastar_, 2));
142 const double phistar = tan(0.5 * (PI - dPhiyy)) * sinthetastar_;
143
144 // a_t
145 const Vector3 t_hat(y1.x()-y2.x(), y1.y()-y2.y(), 0.);
146 const double factor = t_hat.mod();
147 const Vector3 t_hatx(t_hat.x()/factor, t_hat.y()/factor, t_hat.z()/factor);
148 const Vector3 At(y1.x()+y2.x(), y1.y()+y2.y(), 0.);
149 // Compute a_t transverse component with respect to t_hat
150 const double at = At.cross(t_hatx).mod();
151
152 // Fill histograms
153 _h_M->fill(Myy);
154 _h_pT->fill(pTyy);
155 _h_dPhi->fill(dPhiyy);
156 _h_costh->fill(costhetastar_);
157 _h_phistar->fill(phistar);
158 _h_at->fill(at);
159 }
160
161
162 // Normalise histograms etc., after the run
163 void finalize() {
164 const double sf = crossSection()/femtobarn / sumOfWeights();
165 scale(_h_M, sf); scale(_h_pT, sf); scale(_h_dPhi, sf);
166 scale(_h_costh, sf); scale(_h_phistar, sf); scale(_h_at, sf);
167 }
168
169
170 private:
171
172 Histo1DPtr _h_M, _h_pT, _h_dPhi, _h_costh, _h_phistar, _h_at;
173
174 fastjet::AreaDefinition* _area_def;
175
176 const vector<double> ETA_BINS = {0.0, 1.5, 3.0};
177 vector<double> _ptDensity, _sigma, _Njets;
178
179 };
180
181
182 RIVET_DECLARE_PLUGIN(ATLAS_2017_I1591327);
183
184}
|
