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}