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  /// @brief Measurement of isolated diphoton + X differential cross-sections
 11  ///
 12  /// Inclusive isolated gamma gamma cross-sections, differential in M(gg), pT(gg), dphi(gg)
 13  ///
 14  /// @author Giovanni Marchiori
 15  class ATLAS_2011_I916832 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I916832);
 20
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24      FinalState fs;
 25      declare(fs, "FS");
 26
 27      FastJets fj(fs, JetAlg::KT, 0.5);
 28      fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
 29      declare(fj, "KtJetsD05");
 30
 31      IdentifiedFinalState photonfs(Cuts::abseta < 2.37 && Cuts::pT > 16*GeV);
 32      photonfs.acceptId(PID::PHOTON);
 33      declare(photonfs, "Photon");
 34
 35      book(_h_M    ,1, 1, 1);
 36      book(_h_pT   ,2, 1, 1);
 37      book(_h_dPhi ,3, 1, 1);
 38    }
 39
 40
 41    size_t getEtaBin(double eta) const {
 42      const double aeta = fabs(eta);
 43      return binIndex(aeta, _eta_bins_areaoffset);
 44    }
 45
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49      // Require at least 2 photons in final state
 50      Particles photons = apply<IdentifiedFinalState>(event, "Photon").particlesByPt();
 51      if (photons.size() < 2) vetoEvent;
 52
 53      // Compute jet pT densities
 54      vector<double> ptDensity;
 55      vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
 56      const shared_ptr<fastjet::ClusterSequenceArea> clust_seq_area = apply<FastJets>(event, "KtJetsD05").clusterSeqArea();
 57      for (const Jet& jet : apply<FastJets>(event, "KtJetsD05").jets()) {
 58        const double area = clust_seq_area->area(jet); // .pseudojet() called implicitly
 59        /// @todo Should be 1e-4?
 60        if (area > 10e-4 && jet.abseta() < _eta_bins_areaoffset.back()) {
 61          ptDensities.at(getEtaBin(jet.abseta())).push_back(jet.pT()/area);
 62        }
 63      }
 64
 65      // Compute the median energy density
 66      for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
 67        ptDensity += ptDensities[b].empty() ? 0 : median(ptDensities[b]);
 68      }
 69
 70      // Loop over photons and fill vector of isolated ones
 71      Particles isolated_photons;
 72      for (const Particle& photon : photons) {
 73
 74        // Remove photons in crack
 75        if (inRange(photon.abseta(), 1.37, 1.52)) continue;
 76
 77        // Standard ET cone isolation
 78        const Particles& fs = apply<FinalState>(event, "FS").particles();
 79        FourMomentum mom_in_EtCone;
 80        for (const Particle& p : fs) {
 81          // Check if it's in the cone of .4
 82          if (deltaR(photon, p) >= 0.4) continue;
 83          // Veto if it's in the 5x7 central core
 84          if (fabs(deltaEta(photon, p)) < 0.025*5.0*0.5 &&
 85              fabs(deltaPhi(photon, p)) < (M_PI/128.)*7.0*0.5) continue;
 86          // Increment isolation cone ET sum
 87          mom_in_EtCone += p.momentum();
 88        }
 89
 90        // Now figure out the correction (area*density)
 91        const double ETCONE_AREA = M_PI*.4*.4 - (7.0*.025)*(5.0*M_PI/128.);
 92        const double correction = ptDensity[getEtaBin(photon.abseta())] * ETCONE_AREA;
 93
 94        // Shouldn't need to subtract photon
 95        // NB. Using expected cut at hadron/particle level, not cut at reco level
 96        if (mom_in_EtCone.Et() - correction > 4.0*GeV) continue;
 97
 98        // Add to list of isolated photons
 99        isolated_photons.push_back(photon);
100      }
101
102      // Require at least two isolated photons
103      if (isolated_photons.size() < 2) vetoEvent;
104
105      // Select leading pT pair
106      std::sort(isolated_photons.begin(), isolated_photons.end(), cmpMomByPt);
107      const FourMomentum y1 = isolated_photons[0].momentum();
108      const FourMomentum y2 = isolated_photons[1].momentum();
109
110      // Require the two photons to be separated (dR>0.4)
111      if (deltaR(y1, y2) < 0.4) vetoEvent;
112
113      const FourMomentum yy = y1 + y2;
114      const double Myy = yy.mass()/GeV;
115      const double pTyy = yy.pT()/GeV;
116      const double dPhiyy = deltaPhi(y1.phi(), y2.phi());
117
118      _h_M->fill(Myy);
119      _h_pT->fill(pTyy);
120      _h_dPhi->fill(dPhiyy);
121    }
122
123
124    /// Normalise histograms etc., after the run
125    void finalize() {
126      scale(_h_M, crossSection()/picobarn/sumOfWeights());
127      scale(_h_pT, crossSection()/picobarn/sumOfWeights());
128      scale(_h_dPhi, crossSection()/picobarn/sumOfWeights());
129    }
130
131
132  private:
133
134    Histo1DPtr _h_M, _h_pT, _h_dPhi;
135    const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0};
136
137  };
138
139
140  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I916832, ATLAS_2011_S9120807);
141
142}