1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/IdentifiedFinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Measurement of isolated diphoton + X differential cross-sections
 10  ///
 11  /// Inclusive isolated gamma gamma cross-sections, differential in M(gg), pT(gg),
 12  /// dphi(gg), cos(theta*)_CS
 13  ///
 14  /// @author Giovanni Marchiori
 15  ///
 16  class ATLAS_2012_I1199269 : public Analysis {
 17  public:
 18
 19    /// Constructor
 20    ATLAS_2012_I1199269()
 21      : Analysis("ATLAS_2012_I1199269")
 22    {    }
 23
 24
 25    /// Book histograms and initialise projections before the run
 26    void init() {
 27
 28      FinalState fs;
 29      declare(fs, "FS");
 30
 31      FastJets fj(fs, JetAlg::KT, 0.5);
 32      fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
 33      declare(fj, "KtJetsD05");
 34
 35      IdentifiedFinalState photonfs(Cuts::abseta < 2.37 && Cuts::pT > 22*GeV);
 36      photonfs.acceptId(PID::PHOTON);
 37      declare(photonfs, "Photon");
 38
 39      book(_h_M            ,1, 1, 1);
 40      book(_h_pT           ,2, 1, 1);
 41      book(_h_dPhi         ,3, 1, 1);
 42      book(_h_cosThetaStar ,4, 1, 1);
 43    }
 44
 45
 46    /// Perform the per-event analysis
 47    void analyze(const Event& event) {
 48
 49      // Require at least 2 photons in final state
 50      const Particles photons = apply<IdentifiedFinalState>(event, "Photon").particlesByPt();
 51      if (photons.size() < 2) vetoEvent;
 52
 53      // Get jets, and corresponding jet areas
 54      vector<vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
 55      const auto clust_seq_area = apply<FastJets>(event, "KtJetsD05").clusterSeqArea();
 56      for (const Jet& jet : apply<FastJets>(event, "KtJetsD05").jets()) {
 57        const double area = clust_seq_area->area(jet); // implicit .pseudojet()
 58        if (area < 1e-3) continue;
 59        const int ieta = binIndex(jet.abseta(), _eta_bins_areaoffset);
 60        if (ieta != -1) ptDensities[ieta].push_back(jet.pT()/area);
 61      }
 62
 63      // Compute median jet properties over the jets in the event
 64      vector<double> vptDensity; //, vsigma, vNjets;
 65      for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
 66        vptDensity += ptDensities[b].empty() ? 0 : median(ptDensities[b]);
 67      }
 68
 69
 70      // Loop over photons and fill vector of isolated ones
 71      Particles isolated_photons;
 72      for (const Particle& photon : photons) {
 73        /// Remove photons in ECAL crack region
 74        if (inRange(photon.abseta(), 1.37, 1.52)) continue;
 75        // Compute isolation via particles within an R=0.4 cone of the photon
 76        const Particles& fs = apply<FinalState>(event, "FS").particles();
 77        FourMomentum mom_in_EtCone;
 78        for (const Particle& p : fs) {
 79          // Reject if not in cone
 80          if (deltaR(photon, p) > 0.4) continue;
 81          // Reject if in the 5x7 cell central core
 82          if (fabs(deltaEta(photon, p)) < 0.025 * 5 * 0.5 &&
 83              fabs(deltaPhi(photon, p)) < PI/128. * 7 * 0.5) continue;
 84          // Sum momentum
 85          mom_in_EtCone += p.momentum();
 86        }
 87        // Now figure out the correction (area*density)
 88        const double ETCONE_AREA = PI*sqr(0.4) - (7*.025)*(5*PI/128.); // cone area - central core rectangle
 89        const double correction = vptDensity[binIndex(photon.abseta(), _eta_bins_areaoffset)] * ETCONE_AREA;
 90
 91        // Discard the photon if there is more than 4 GeV of cone activity
 92        // NOTE: Shouldn't need to subtract photon itself (it's in the central core)
 93        // NOTE: using expected cut at hadron/particle level, not at reco level
 94        if (mom_in_EtCone.Et() - correction > 4*GeV) continue;
 95        // Add isolated photon to list
 96        isolated_photons.push_back(photon);
 97      }
 98
 99      // Require at least two isolated photons and select leading pT pair
100      if (isolated_photons.size() < 2) vetoEvent;
101      sortByPt(isolated_photons);
102      const FourMomentum& y1 = isolated_photons[0].momentum();
103      const FourMomentum& y2 = isolated_photons[1].momentum();
104
105      // Leading photon should have pT > 25 GeV
106      if (y1.pT() < 25*GeV) vetoEvent;
107
108      // Require the two photons to be separated by dR > 0.4
109      if (deltaR(y1, y2) < 0.4) vetoEvent;
110
111      // Compute diphoton vector and fill histos
112      FourMomentum yy = y1 + y2;
113      const double costhetayy = 2 * y1.pT() * y2.pT() * sinh(y1.eta() - y2.eta()) / yy.mass() / add_quad(yy.mass(), yy.pT());
114      _h_M->fill(yy.mass()/GeV);
115      _h_pT->fill(yy.pT()/GeV);
116      _h_dPhi->fill(mapAngle0ToPi(y1.phi() - y2.phi()));
117      _h_cosThetaStar->fill(costhetayy);
118    }
119
120
121    /// Normalise histograms etc., after the run
122    void finalize() {
123      scale(_h_M, crossSection()/picobarn/sumOfWeights());
124      scale(_h_pT, crossSection()/picobarn/sumOfWeights());
125      scale(_h_dPhi, crossSection()/picobarn/sumOfWeights());
126      scale(_h_cosThetaStar, crossSection()/picobarn/sumOfWeights());
127    }
128
129
130  private:
131
132    Histo1DPtr _h_M, _h_pT, _h_dPhi, _h_cosThetaStar;
133
134    const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0};
135
136  };
137
138
139  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1199269);
140
141}