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  class ATLAS_2014_I1307756 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    ATLAS_2014_I1307756()
 15      : Analysis("ATLAS_2014_I1307756")
 16    {    }
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      /// Initialise and register projections here
 26      FinalState fs;
 27      declare(fs, "FS");
 28
 29      FastJets fj(fs, JetAlg::KT, 0.5);
 30      fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
 31      declare(fj, "KtJetsD05");
 32
 33      IdentifiedFinalState photonfs(Cuts::abseta < 2.37 && Cuts::pT > 22*GeV);
 34      photonfs.acceptId(PID::PHOTON);
 35      declare(photonfs, "photons");
 36
 37      // Initialize event count here:
 38      book(_fidWeights, "_fidWeights");
 39    }
 40
 41
 42    int getEtaBin(double eta) const {
 43      double aeta = fabs(eta);
 44      return binIndex(aeta, _eta_bins_areaoffset);
 45    }
 46
 47
 48    /// Perform the per-event analysis
 49    void analyze(const Event& event) {
 50      /// Require at least 2 photons in final state
 51      Particles photons = apply<IdentifiedFinalState>(event, "photons").particlesByPt();
 52      if (photons.size() < 2) vetoEvent;
 53
 54      // Get jet pT densities
 55      vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
 56      const auto 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);
 59        if (area > 1e-4 && jet.abseta() < _eta_bins_areaoffset.back()) {
 60          ptDensities.at(getEtaBin(jet.abseta())) += jet.pT()/area;
 61        }
 62      }
 63
 64      /// Compute the median energy density per eta bin
 65      vector<double> ptDensity;
 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 find isolated ones
 71      Particles isolated_photons;
 72      for (const Particle& ph : photons) {
 73        Particles fs = apply<FinalState>(event, "FS").particles();
 74        FourMomentum mom_in_EtCone;
 75        for (const Particle& p : fs) {
 76
 77          // Reject if the particle is not in DR=0.4 cone
 78          if (deltaR(ph.momentum(), p.momentum()) > 0.4) continue;
 79
 80          // Reject if the particle falls in the photon core
 81          if (fabs(ph.eta() - p.eta()) < 0.025 * 7 * 0.5 &&
 82              fabs(ph.phi() - p.phi()) < PI/128. * 5 * 0.5) continue;
 83
 84          // Reject if the particle is a neutrino (muons are kept)
 85          if (p.isNeutrino()) continue;
 86
 87          // Sum momenta
 88          mom_in_EtCone += p.momentum();
 89        }
 90
 91        // Subtract the UE correction (area*density)
 92        const double ETCONE_AREA = M_PI*.4*.4 - (7.0*.025)*(5.0*M_PI/128.);
 93        const double correction = ptDensity[getEtaBin(ph.eta())] * ETCONE_AREA;
 94
 95        // Add isolated photon to list
 96        if (mom_in_EtCone.Et() - correction > 12*GeV) continue;
 97        isolated_photons.push_back(ph);
 98      }
 99
100      // Require at least two isolated photons
101      if (isolated_photons.size() < 2)  vetoEvent ;
102
103      // Select leading pT pair
104      std::sort(isolated_photons.begin(), isolated_photons.end(), cmpMomByPt);
105      const FourMomentum& y1 = isolated_photons[0].momentum();
106      const FourMomentum& y2 = isolated_photons[1].momentum();
107
108      // Compute invariant mass
109      const FourMomentum yy = y1 + y2;
110      const double Myy = yy.mass();
111
112      // If Myy >= 110 GeV, apply relative cuts
113      if (Myy >= 110*GeV && (y1.Et()/Myy < 0.4 || y2.Et()/Myy < 0.3) ) vetoEvent;
114
115      // Add to cross-section
116      _fidWeights->fill();
117    }
118
119    void finalize() {
120      scale(_fidWeights, crossSectionPerEvent()/femtobarn);
121    }
122
123    /// @}
124
125
126  private:
127
128    const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0};
129    CounterPtr _fidWeights;
130
131  };
132
133
134
135  RIVET_DECLARE_PLUGIN(ATLAS_2014_I1307756);
136
137}