1// -*- C++ -*-
  2#include "Rivet/Analyses/MC_JETS_BASE.hh"
  3#include "Rivet/Projections/LeadingParticlesFinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief MC validation analysis for photon + jets events
 11  class MC_PHOTONJETS : public MC_JETS_BASE {
 12  public:
 13
 14    /// Default constructor
 15    MC_PHOTONJETS()
 16      : MC_JETS_BASE("MC_PHOTONJETS", 4, "Jets")
 17    {    }
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms
 24    void init() {
 25      // General FS
 26      FinalState fs((Cuts::etaIn(-5.0, 5.0)));
 27      declare(fs, "FS");
 28
 29      // set ptcut from input option
 30      const double jetptcut = getOption<double>("PTJMIN", 20.0);
 31      _jetptcut = jetptcut * GeV;
 32
 33      // set clustering radius from input option
 34      const double R = getOption<double>("R", 0.4);
 35
 36      // set clustering algorithm from input option
 37      JetAlg clusterAlgo;
 38      const string algoopt = getOption("ALGO", "ANTIKT");
 39      if ( algoopt == "KT" ) {
 40	clusterAlgo = JetAlg::KT;
 41      } else if ( algoopt == "CA" ) {
 42	clusterAlgo = JetAlg::CA;
 43      } else if ( algoopt == "ANTIKT" ) {
 44	clusterAlgo = JetAlg::ANTIKT;
 45      } else {
 46	MSG_WARNING("Unknown jet clustering algorithm option " + algoopt + ". Defaulting to anti-kT");
 47	clusterAlgo = JetAlg::ANTIKT;
 48      }
 49      
 50      // set photon cuts from input options
 51      const double etacut = getOption<double>("ABSETAGAMMAX", 2.5);
 52      const double ptcut = getOption<double>("PTGAMMIN", 30.);
 53      
 54      // Get leading photon
 55      LeadingParticlesFinalState photonfs(FinalState(Cuts::abseta < etacut && Cuts::pT >= ptcut*GeV));
 56      photonfs.addParticleId(PID::PHOTON);
 57      declare(photonfs, "LeadingPhoton");
 58
 59      // FS for jets excludes the leading photon
 60      VetoedFinalState vfs(fs);
 61      vfs.addVetoOnThisFinalState(photonfs);
 62      declare(vfs, "JetFS");
 63      FastJets jetpro(vfs, clusterAlgo, R);
 64      declare(jetpro, "Jets");
 65
 66      book(_h_photon_jet1_deta ,"photon_jet1_deta", 50, -5.0, 5.0);
 67      book(_h_photon_jet1_dphi ,"photon_jet1_dphi", 20, 0.0, M_PI);
 68      book(_h_photon_jet1_dR ,"photon_jet1_dR", 25, 0.5, 7.0);
 69
 70      MC_JETS_BASE::init();
 71    }
 72
 73
 74    /// Do the analysis
 75    void analyze(const Event& e) {
 76      // Get the photon
 77      /// @todo share IsolatedPhoton projection between all MC_*PHOTON* analyses
 78      const Particles photons = apply<FinalState>(e, "LeadingPhoton").particles();
 79      if (photons.size() != 1) {
 80        vetoEvent;
 81      }
 82      const FourMomentum photon = photons.front().momentum();
 83
 84      // Get all charged particles
 85      const FinalState& fs = apply<FinalState>(e, "JetFS");
 86      if (fs.empty()) {
 87        vetoEvent;
 88      }
 89
 90      // Passed cuts, so get the weight
 91
 92      // Isolate photon by ensuring that a 0.4 cone around it contains less than 7% of the photon's energy
 93      const double egamma = photon.E();
 94      double econe = 0.0;
 95      for (const Particle& p : fs.particles()) {
 96        if (deltaR(photon, p.momentum()) < 0.4) {
 97          econe += p.E();
 98          // Veto as soon as E_cone gets larger
 99          if (econe/egamma > 0.07) {
100            vetoEvent;
101          }
102        }
103      }
104
105      const Jets& jets = apply<FastJets>(e, "Jets").jetsByPt(Cuts::pT > _jetptcut);
106      if (jets.size()>0) {
107        _h_photon_jet1_deta->fill(photon.eta()-jets[0].eta());
108        _h_photon_jet1_dphi->fill(mapAngle0ToPi(photon.phi()-jets[0].phi()));
109        _h_photon_jet1_dR->fill(deltaR(photon, jets[0].momentum()));
110      }
111
112      MC_JETS_BASE::analyze(e);
113    }
114
115
116    // Finalize
117    void finalize() {
118      scale(_h_photon_jet1_deta, crossSectionPerEvent());
119      scale(_h_photon_jet1_dphi, crossSectionPerEvent());
120      scale(_h_photon_jet1_dR, crossSectionPerEvent());
121
122      MC_JETS_BASE::finalize();
123    }
124
125    /// @}
126
127
128  private:
129
130    /// @name Histograms
131    /// @{
132    Histo1DPtr _h_photon_jet1_deta;
133    Histo1DPtr _h_photon_jet1_dphi;
134    Histo1DPtr _h_photon_jet1_dR;
135    /// @}
136
137  };
138
139
140
141  RIVET_DECLARE_PLUGIN(MC_PHOTONJETS);
142
143}