rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2017_I1644367

Isolated triphotons at 8 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1644367
Status: VALIDATED
Authors:
  • Josu Cantero Garcia
  • Christian Gutschow
References: Beams: p+ p+
Beam energies: (4000.0, 4000.0) GeV
Run details:
  • inclusive triphoton production

A measurement of the production of three isolated photons in proton-proton collisions at a centre-of-mass energy $\sqrt{s} = 8$ TeV is reported. The results are based on an integrated luminosity of 20.2 fb${}^{-1}$ collected with the ATLAS detector at the LHC. The differential cross sections are measured as functions of the transverse energy of each photon, the difference in azimuthal angle and in pseudorapidity between pairs of photons, the invariant mass of pairs of photons, and the invariant mass of the triphoton system. A measurement of the inclusive fiducial cross section is also reported. Next-to-leading-order perturbative QCD predictions are compared to the cross-section measurements. The predictions underestimate the measurement of the inclusive fiducial cross section and the differential measurements at low photon transverse energies and invariant masses. They provide adequate descriptions of the measurements at high values of the photon transverse energies, invariant mass of pairs of photons, and invariant mass of the triphoton system.

Source code: ATLAS_2017_I1644367.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief Isolated triphotons at 8 TeV
 11  class ATLAS_2017_I1644367 : public Analysis {
 12  public:
 13
 14    // Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1644367);
 16
 17    // Book histograms and initialise projections before the run
 18    void init() {
 19
 20      const FinalState fs;
 21      declare(fs, "FS");
 22
 23      FastJets fj(fs, JetAlg::KT, 0.5);
 24      fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
 25      declare(fj, "KtJetsD05");
 26
 27      PromptFinalState photonfs(Cuts::abspid == PID::PHOTON && Cuts::abseta < 2.37 && Cuts::pT > 15*GeV);
 28      declare(photonfs, "Photon");
 29
 30      // Histograms
 31      book(_h["etg1"]    ,  1, 1, 1);
 32      book(_h["etg2"]    ,  2, 1, 1);
 33      book(_h["etg3"]    ,  3, 1, 1);
 34      book(_h["dphig1g2"],  4, 1, 1);
 35      book(_h["dphig1g3"],  5, 1, 1);
 36      book(_h["dphig2g3"],  6, 1, 1);
 37      book(_h["detag1g2"],  7, 1, 1);
 38      book(_h["detag1g3"],  8, 1, 1);
 39      book(_h["detag2g3"],  9, 1, 1);
 40      book(_h["mg1g2"]   , 10, 1, 1);
 41      book(_h["mg1g3"]   , 11, 1, 1);
 42      book(_h["mg2g3"]   , 12, 1, 1);
 43      book(_h["mg1g2g3"] , 13, 1, 1);
 44
 45    }
 46
 47
 48    // Perform the per-event analysis
 49    void analyze(const Event& event) {
 50
 51      // Require at least 2 photons in final state
 52      const Particles photons = apply<PromptFinalState>(event, "Photon").particlesByPt(Cuts::abseta < 1.37 || Cuts::abseta > 1.5);
 53      if (photons.size() < 3) vetoEvent;
 54
 55      // Get jets, and corresponding jet areas
 56      vector<vector<double> > ptDensities(ETA_BINS.size()-1);
 57      FastJets fastjets = apply<FastJets>(event, "KtJetsD05");
 58      const auto clust_seq_area = fastjets.clusterSeqArea();
 59      for (const Jet& jet : fastjets.jets()) {
 60        const double area = clust_seq_area->area(jet);
 61        if (area < 1e-3) continue;
 62        const int ieta = binIndex(jet.abseta(), ETA_BINS);
 63        if (ieta != -1) ptDensities[ieta].push_back(jet.pT()/area);
 64      }
 65
 66      // Compute median jet properties over the jets in the event
 67      vector<double> ptDensity;
 68      for (size_t b = 0; b < ETA_BINS.size()-1; ++b) {
 69        double median = 0.0;
 70        if (ptDensities[b].size() > 0) {
 71          std::sort(ptDensities[b].begin(), ptDensities[b].end());
 72          int nDens = ptDensities[b].size();
 73          median = (nDens % 2 == 0) ? (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2 : ptDensities[b][(nDens-1)/2];
 74        }
 75
 76        ptDensity.push_back(median);
 77      }
 78
 79      // Loop over photons and fill vector of isolated ones
 80      Particles isolated_photons;
 81      for (const Particle& photon : photons) {
 82        if (!photon.isPrompt()) continue;
 83
 84        // Remove photons in ECAL crack region
 85        const double eta_P = photon.eta();
 86        const double phi_P = photon.phi();
 87
 88        // Compute isolation via particles within an R=0.4 cone of the photon
 89        const Particles fs = apply<FinalState>(event, "FS").particles();
 90        FourMomentum mom_in_EtCone;
 91        for (const Particle& p : fs) {
 92          // Reject if not in cone
 93          if (deltaR(photon.momentum(), p.momentum()) > 0.4)  continue;
 94          // Reject if in the 5x7 cell central core
 95          if (fabs(eta_P - p.eta()) < 0.025 * 5 * 0.5 && fabs(phi_P - p.phi()) < PI/128. * 7 * 0.5)  continue;
 96          // Sum momentum
 97          mom_in_EtCone += p.momentum();
 98        }
 99
100        // Now figure out the correction (area*density)
101        const double EtCone_area = M_PI*sqr(0.4) - (7*.025)*(5*M_PI/128.); // cone area - central core rectangle
102        const double correction = ptDensity[binIndex(fabs(eta_P), ETA_BINS)] * EtCone_area;
103
104        // Discard the photon if there is more than 11 GeV of cone activity
105        // NOTE: Shouldn't need to subtract photon itself (it's in the central core)
106        if (mom_in_EtCone.Et() - correction > 10*GeV)  continue;
107        // Add isolated photon to list
108        isolated_photons.push_back(photon);
109      }///loop over photons
110
111      // Require at least two isolated photons
112      if (isolated_photons.size() < 3) vetoEvent;
113
114      // Select leading pT pair
115      sortByPt(isolated_photons);
116      const FourMomentum y1 = isolated_photons[0];
117      const FourMomentum y2 = isolated_photons[1];
118      const FourMomentum y3 = isolated_photons[2];
119
120      // Leading photon should have pT > 40 GeV, subleading > 30 GeV
121      if (y1.pT() < 27*GeV)  vetoEvent;
122      if (y2.pT() < 22*GeV)  vetoEvent;
123      if (y3.pT() < 15*GeV)  vetoEvent;
124
125      // Require the two photons to be separated (dR>0.4)
126      if (deltaR(y1,y2) < 0.45)  vetoEvent;
127      if (deltaR(y1,y3) < 0.45)  vetoEvent;
128      if (deltaR(y2,y3) < 0.45)  vetoEvent;
129
130
131      const FourMomentum yyy = y1 + y2 + y3;
132      const FourMomentum y1y2 = y1 + y2;
133      const FourMomentum y1y3 = y1 + y3;
134      const FourMomentum y2y3 = y2 + y3;
135
136      const double Myyy = yyy.mass() / GeV;
137
138      const double dPhiy1y2 = mapAngle0ToPi(deltaPhi(y1, y2));
139      const double dPhiy1y3 = mapAngle0ToPi(deltaPhi(y1, y3));
140      const double dPhiy2y3 = mapAngle0ToPi(deltaPhi(y2, y3));
141
142      const double dEtay1y2 = fabs(y1.eta() - y2.eta());
143      const double dEtay1y3 = fabs(y1.eta() - y3.eta());
144      const double dEtay2y3 = fabs(y2.eta() - y3.eta());
145
146      if(Myyy < 50.) vetoEvent;
147
148      // Fill histograms
149
150      _h["etg1"]->fill(y1.pT() / GeV);
151      _h["etg2"]->fill(y2.pT() / GeV);
152      _h["etg3"]->fill(y3.pT() / GeV);
153
154      _h["dphig1g2"]->fill(dPhiy1y2);
155      _h["dphig1g3"]->fill(dPhiy1y3);
156      _h["dphig2g3"]->fill(dPhiy2y3);
157
158      _h["detag1g2"]->fill(dEtay1y2);
159      _h["detag1g3"]->fill(dEtay1y3);
160      _h["detag2g3"]->fill(dEtay2y3);
161
162      _h["mg1g2"]->fill(y1y2.mass() / GeV);
163      _h["mg1g3"]->fill(y1y3.mass() / GeV);
164      _h["mg2g3"]->fill(y2y3.mass() / GeV);
165
166      _h["mg1g2g3"]->fill(Myyy);
167    }
168
169
170    // Normalise histograms etc., after the run
171    void finalize() {
172      const double sf = crossSection() / (femtobarn * sumOfWeights());
173      for (auto &hist : _h) {  scale(hist.second, sf); }
174    }
175
176
177  private:
178
179    map<string, Histo1DPtr> _h;
180    const vector<double> ETA_BINS = { 0.0, 1.5, 3.0 };
181
182  };
183
184
185  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1644367);
186
187}