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
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/FastJets.hh"

namespace Rivet {

  /// @brief Isolated triphotons at 8 TeV
  class ATLAS_2017_I1644367 : public Analysis {
  public:

    // Constructor
    DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1644367);
      
    // Book histograms and initialise projections before the run
    void init() {
      
      const FinalState fs;
	    declare(fs, "FS");

      FastJets fj(fs, FastJets::KT, 0.5);
	    fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
	    declare(fj, "KtJetsD05");

	    IdentifiedFinalState photonfs(Cuts::abspid == PID::PHOTON && Cuts::abseta < 2.37 && Cuts::pT > 15*GeV);
      declare(photonfs, "Photon");
	
	    // Histograms
	    _h["etg1"]     = bookHisto1D( 1, 1, 1);
	    _h["etg2"]     = bookHisto1D( 2, 1, 1);
	    _h["etg3"]     = bookHisto1D( 3, 1, 1);
	    _h["dphig1g2"] = bookHisto1D( 4, 1, 1);
	    _h["dphig1g3"] = bookHisto1D( 5, 1, 1);
	    _h["dphig2g3"] = bookHisto1D( 6, 1, 1);
	    _h["detag1g2"] = bookHisto1D( 7, 1, 1);
	    _h["detag1g3"] = bookHisto1D( 8, 1, 1);
	    _h["detag2g3"] = bookHisto1D( 9, 1, 1);
	    _h["mg1g2"]    = bookHisto1D(10, 1, 1);
	    _h["mg1g3"]    = bookHisto1D(11, 1, 1);
	    _h["mg2g3"]    = bookHisto1D(12, 1, 1);
	    _h["mg1g2g3"]  = bookHisto1D(13, 1, 1);

    }
     
    // Perform the per-event analysis
    void analyze(const Event& event) {

      const double weight = event.weight();
	    // Require at least 2 photons in final state
	    const Particles photons = apply<IdentifiedFinalState>(event, "Photon").particlesByPt(Cuts::abseta < 1.37 || Cuts::abseta > 1.5);
	    if (photons.size() < 3) vetoEvent;

	    // Get jets, and corresponding jet areas
   	  vector<vector<double> > ptDensities(ETA_BINS.size()-1);
      FastJets fastjets = apply<FastJets>(event, "KtJetsD05");
	    const auto clust_seq_area = fastjets.clusterSeqArea();
	    for (const Jet& jet : fastjets.jets()) {
	      const double area = clust_seq_area->area(jet);
	      if (area < 1e-3) continue;
	      const int ieta = binIndex(jet.abseta(), ETA_BINS);
	      if (ieta != -1) ptDensities[ieta].push_back(jet.pT()/area);
	    }

	    // Compute median jet properties over the jets in the event
	    vector<double> ptDensity;
	    for (size_t b = 0; b < ETA_BINS.size()-1; ++b) {
	      double median = 0.0;
	      if (ptDensities[b].size() > 0) {
	        std::sort(ptDensities[b].begin(), ptDensities[b].end());
	        int nDens = ptDensities[b].size();
	        median = (nDens % 2 == 0) ? (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2 : ptDensities[b][(nDens-1)/2];
        }

        ptDensity.push_back(median);
      }

      // Loop over photons and fill vector of isolated ones
      Particles isolated_photons;
      for (const Particle& photon : photons) {
	      if (photon.fromDecay()) continue;
	  
	      // Remove photons in ECAL crack region
	      const double eta_P = photon.eta();
	      const double phi_P = photon.phi();
	  
	      // Compute isolation via particles within an R=0.4 cone of the photon
	      const Particles fs = apply<FinalState>(event, "FS").particles();
	      FourMomentum mom_in_EtCone;
	      for (const Particle& p : fs) {
	        // Reject if not in cone
	        if (deltaR(photon.momentum(), p.momentum()) > 0.4)  continue;
	        // Reject if in the 5x7 cell central core
          if (fabs(eta_P - p.eta()) < 0.025 * 5 * 0.5 && fabs(phi_P - p.phi()) < PI/128. * 7 * 0.5)  continue;
          // Sum momentum
          mom_in_EtCone += p.momentum();
        }

        // Now figure out the correction (area*density)
        const double EtCone_area = M_PI*sqr(0.4) - (7*.025)*(5*M_PI/128.); // cone area - central core rectangle
        const double correction = ptDensity[binIndex(fabs(eta_P), ETA_BINS)] * EtCone_area;

        // Discard the photon if there is more than 11 GeV of cone activity
        // NOTE: Shouldn't need to subtract photon itself (it's in the central core)
        if (mom_in_EtCone.Et() - correction > 10*GeV)  continue;
        // Add isolated photon to list
        isolated_photons.push_back(photon);
	    }///loop over photons
	       
	    // Require at least two isolated photons
	    if (isolated_photons.size() < 3) vetoEvent;
	
	    // Select leading pT pair
	    sortByPt(isolated_photons);
	    const FourMomentum y1 = isolated_photons[0];
	    const FourMomentum y2 = isolated_photons[1];
	    const FourMomentum y3 = isolated_photons[2];
	
	    // Leading photon should have pT > 40 GeV, subleading > 30 GeV
	    if (y1.pT() < 27*GeV)  vetoEvent;
	    if (y2.pT() < 22*GeV)  vetoEvent;
	    if (y3.pT() < 15*GeV)  vetoEvent;
	      
      // Require the two photons to be separated (dR>0.4)
	    if (deltaR(y1,y2) < 0.45)  vetoEvent;
	    if (deltaR(y1,y3) < 0.45)  vetoEvent;
	    if (deltaR(y2,y3) < 0.45)  vetoEvent;
	      
       
	    const FourMomentum yyy = y1 + y2 + y3;
	    const FourMomentum y1y2 = y1 + y2;
	    const FourMomentum y1y3 = y1 + y3;
	    const FourMomentum y2y3 = y2 + y3;

	    const double Myyy = yyy.mass() / GeV;

	    const double dPhiy1y2 = mapAngle0ToPi(deltaPhi(y1, y2));
	    const double dPhiy1y3 = mapAngle0ToPi(deltaPhi(y1, y3));
	    const double dPhiy2y3 = mapAngle0ToPi(deltaPhi(y2, y3));

	    const double dEtay1y2 = fabs(y1.eta() - y2.eta());
	    const double dEtay1y3 = fabs(y1.eta() - y3.eta());
	    const double dEtay2y3 = fabs(y2.eta() - y3.eta());

      if(Myyy < 50.) vetoEvent;

      // Fill histograms

	    _h["etg1"]->fill(y1.pT() / GeV, weight);
	    _h["etg2"]->fill(y2.pT() / GeV, weight);
	    _h["etg3"]->fill(y3.pT() / GeV, weight);

	    _h["dphig1g2"]->fill(dPhiy1y2);
	    _h["dphig1g3"]->fill(dPhiy1y3);
	    _h["dphig2g3"]->fill(dPhiy2y3);

	    _h["detag1g2"]->fill(dEtay1y2);
	    _h["detag1g3"]->fill(dEtay1y3);
	    _h["detag2g3"]->fill(dEtay2y3);

	    _h["mg1g2"]->fill(y1y2.mass() / GeV, weight);
	    _h["mg1g3"]->fill(y1y3.mass() / GeV, weight);
	    _h["mg2g3"]->fill(y2y3.mass() / GeV, weight);

	    _h["mg1g2g3"]->fill(Myyy, weight);
	  }
	
	
	  // Normalise histograms etc., after the run
	  void finalize() {
	    const double sf = crossSection() / (femtobarn * sumOfWeights());
      for (auto &hist : _h) {  scale(hist.second, sf); }
	  }
	
	
	  private:
	
   	  map<string, Histo1DPtr> _h;
	    const vector<double> ETA_BINS = { 0.0, 1.5, 3.0 };
	
    };
  
    // The hook for the plugin system
    DECLARE_RIVET_PLUGIN(ATLAS_2017_I1644367);
}