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      isortByPt(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}