1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/VisibleFinalState.hh"
  6#include "Rivet/Projections/VetoedFinalState.hh"
  7#include "Rivet/Projections/IdentifiedFinalState.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @author Peter Richardson
 14  class ATLAS_2012_I946427 : public Analysis {
 15  public:
 16
 17    /// @name Constructors etc.
 18    /// @{
 19
 20    /// Constructor
 21    ATLAS_2012_I946427()
 22      : Analysis("ATLAS_2012_I946427")
 23    {    }
 24
 25    /// @}
 26
 27
 28  public:
 29
 30    /// @name Analysis methods
 31    /// @{
 32
 33    /// Book histograms and initialise projections before the run
 34    void init() {
 35
 36      // photons
 37      IdentifiedFinalState photonfs(Cuts::abseta < 1.81 && Cuts::pT > 25*GeV);
 38      photonfs.acceptId(PID::PHOTON);
 39      declare(photonfs, "Photon");
 40
 41      //
 42      FinalState fs;
 43      declare(fs, "FS");
 44
 45      // Used for pTmiss
 46      declare(VisibleFinalState(Cuts::abseta < 4.9),"vfs");
 47
 48      // Book histograms
 49      book(_count_SR ,"count_SR", 1, 0., 1.);
 50
 51      book(_hist_ET_photon ,"hist_ET_photon", 48 , 20., 500.);
 52      book(_hist_met       ,"hist_met"      , 100,  0., 500.);
 53
 54    }
 55
 56
 57    /// Perform the per-event analysis
 58    void analyze(const Event& event) {
 59
 60
 61      // require at least 2 photons in final state
 62      Particles photons =
 63        apply<IdentifiedFinalState>(event, "Photon").particlesByPt();
 64      if (photons.size() < 2) {
 65        vetoEvent;
 66      }
 67
 68      // Loop over photons and fill vector of isolated ones
 69      Particles fs = apply<FinalState>(event, "FS").particles();
 70      Particles isolated_photons;
 71      for (const Particle& photon : photons) {
 72        // remove photons in crack
 73        double eta_P = photon.eta();
 74        if (fabs(eta_P)>=1.37 && fabs(eta_P)<1.52) continue;
 75
 76        double phi_P = photon.phi();
 77
 78        FourMomentum mom_in_EtCone = -photon.momentum();
 79        for (const Particle& p : fs) {
 80          // check if it's in the cone of .2
 81          if (deltaR(eta_P, phi_P, p.eta(),
 82                     p.phi()) >= 0.2) continue;
 83          mom_in_EtCone += p.momentum();
 84        }
 85        // apply isolation
 86        if(mom_in_EtCone.Et()>5.) continue;
 87
 88        // add photon to list of isolated ones
 89        isolated_photons.push_back(photon);
 90      }
 91
 92      // need two isolated photons
 93      if(isolated_photons.size() < 2 ) {
 94        vetoEvent;
 95      }
 96
 97      // pTmiss
 98      Particles vfs_particles =
 99        apply<VisibleFinalState>(event, "vfs").particles();
100      FourMomentum pTmiss;
101      for ( const Particle & p : vfs_particles ) {
102        pTmiss -= p.momentum();
103      }
104      double eTmiss = pTmiss.pT();
105
106      _hist_ET_photon->fill(isolated_photons[0].Et());
107      _hist_met      ->fill(eTmiss                             );
108
109      if(eTmiss>125.) _count_SR->fill(0.5);
110    }
111
112
113    void finalize() {
114
115      double norm = crossSection()/femtobarn*1.07/sumOfWeights();
116      // these are number of events at 1.07fb^-1 per 10 GeV
117      scale( _hist_ET_photon, 10. * norm );
118      // these are number of events at 1.07fb^-1 per  5 GeV
119      scale( _hist_met, 5. * norm );
120      // these are number of events at 1.07fb^-1
121      scale(_count_SR,norm);
122    }
123
124    /// @}
125
126
127  private:
128
129    Histo1DPtr _count_SR;
130    Histo1DPtr _hist_ET_photon;
131    Histo1DPtr _hist_met;
132
133  };
134
135
136  RIVET_DECLARE_PLUGIN(ATLAS_2012_I946427);
137
138}