1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/LeadingParticlesFinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9
 10
 11  /// @brief MC validation analysis for single photon events
 12  class MC_PHOTONINC : public Analysis {
 13  public:
 14
 15    /// Default constructor
 16    MC_PHOTONINC()
 17      : Analysis("MC_PHOTONINC")
 18    {    }
 19
 20
 21    /// @name Analysis methods
 22    /// @{
 23
 24    /// Book histograms
 25    void init() {
 26      // General FS
 27      FinalState fs((Cuts::etaIn(-5.0, 5.0)));
 28      declare(fs, "FS");
 29
 30      // set photon cuts from input options
 31      const double etacut = getOption<double>("ABSETAGAMMAX", 2.5);
 32      const double ptcut = getOption<double>("PTGAMMIN", 30.);
 33      
 34      // Get leading photon
 35      LeadingParticlesFinalState photonfs(FinalState(Cuts::abseta < etacut && Cuts::pT >= ptcut*GeV));
 36      photonfs.addParticleId(PID::PHOTON);
 37      declare(photonfs, "LeadingPhoton");
 38
 39      // FS for isolation excludes the leading photon
 40      VetoedFinalState vfs(fs);
 41      vfs.addVetoOnThisFinalState(photonfs);
 42      declare(vfs, "JetFS");
 43
 44      book(_h_photon_pT ,"photon_pT", logspace(50, 30.0, 0.5*(sqrtS()>0.?sqrtS():14000.)));
 45      book(_h_photon_pT_lin ,"photon_pT_lin", 50, 0.0, 70.0);
 46      book(_h_photon_y ,"photon_y", 50, -5.0, 5.0);
 47    }
 48
 49
 50    /// Do the analysis
 51    void analyze(const Event& e) {
 52      // Get the photon
 53      const Particles photons = apply<FinalState>(e, "LeadingPhoton").particles();
 54      if (photons.size() != 1) {
 55        vetoEvent;
 56      }
 57      const FourMomentum photon = photons.front().momentum();
 58
 59      // Get all charged particles
 60      const FinalState& fs = apply<FinalState>(e, "JetFS");
 61      if (fs.empty()) {
 62        vetoEvent;
 63      }
 64
 65      // Passed cuts, so get the weight
 66
 67      // Isolate photon by ensuring that a 0.4 cone around it contains less than 7% of the photon's energy
 68      const double egamma = photon.E();
 69      double econe = 0.0;
 70      for (const Particle& p : fs.particles()) {
 71        if (deltaR(photon, p.momentum()) < 0.4) {
 72          econe += p.E();
 73          // Veto as soon as E_cone gets larger
 74          if (econe/egamma > 0.07) {
 75            vetoEvent;
 76          }
 77        }
 78      }
 79
 80      _h_photon_pT->fill(photon.pT());
 81      _h_photon_pT_lin->fill(photon.pT());
 82      _h_photon_y->fill(photon.rapidity());
 83    }
 84
 85
 86    // Finalize
 87    void finalize() {
 88      scale(_h_photon_pT, crossSectionPerEvent());
 89      scale(_h_photon_pT_lin, crossSectionPerEvent());
 90      scale(_h_photon_y, crossSectionPerEvent());
 91    }
 92
 93    /// @}
 94
 95
 96  private:
 97
 98    /// @name Histograms
 99    /// @{
100    Histo1DPtr _h_photon_pT;
101    Histo1DPtr _h_photon_pT_lin;
102    Histo1DPtr _h_photon_y;
103    /// @}
104
105  };
106
107
108
109  RIVET_DECLARE_PLUGIN(MC_PHOTONINC);
110
111}