00001
00002 #include "Rivet/Analyses/MC_JetAnalysis.hh"
00003 #include "Rivet/Tools/Logging.hh"
00004 #include "Rivet/Projections/LeadingParticlesFinalState.hh"
00005 #include "Rivet/Projections/FastJets.hh"
00006 #include "Rivet/RivetAIDA.hh"
00007
00008 namespace Rivet {
00009
00010 class MC_PHOTONJETS : public MC_JetAnalysis {
00011 public:
00012
00013
00014 MC_PHOTONJETS()
00015 : MC_JetAnalysis("MC_PHOTONJETS", 4, "Jets")
00016 {
00017 setNeedsCrossSection(true);
00018 }
00019
00020
00021
00022
00023
00024
00025 void init() {
00026
00027 FinalState fs(-5.0, 5.0);
00028 addProjection(fs, "FS");
00029
00030
00031 LeadingParticlesFinalState photonfs(FinalState(-1.0, 1.0));
00032 photonfs.addParticleId(PHOTON);
00033 addProjection(photonfs, "LeadingPhoton");
00034
00035
00036 VetoedFinalState vfs(fs);
00037 vfs.addVetoOnThisFinalState(photonfs);
00038 addProjection(vfs, "JetFS");
00039 FastJets jetpro(vfs, FastJets::KT, 0.7);
00040 addProjection(jetpro, "Jets");
00041
00042 _h_photon_pT = bookHistogram1D("photon_pT", logBinEdges(50, 1.0, 0.5*sqrtS()));
00043 _h_photon_y = bookHistogram1D("photon_y", 20, -1.0, 1.0);
00044 _h_photon_jet1_deta = bookHistogram1D("photon_jet1_deta", 50, -5.0, 5.0);
00045 _h_photon_jet1_dphi = bookHistogram1D("photon_jet1_dphi", 20, 0.0, M_PI);
00046 _h_photon_jet1_dR = bookHistogram1D("photon_jet1_dR", 25, 0.5, 7.0);
00047
00048 MC_JetAnalysis::init();
00049 }
00050
00051
00052
00053 void analyze(const Event& e) {
00054
00055 const ParticleVector photons = applyProjection<FinalState>(e, "LeadingPhoton").particles();
00056 if (photons.size() != 1) {
00057 vetoEvent;
00058 }
00059 const FourMomentum photon = photons.front().momentum();
00060
00061
00062 const FinalState& fs = applyProjection<FinalState>(e, "JetFS");
00063 if (fs.empty()) {
00064 vetoEvent;
00065 }
00066
00067
00068 const double weight = e.weight();
00069
00070
00071 const double egamma = photon.E();
00072 double econe = 0.0;
00073 foreach (const Particle& p, fs.particles()) {
00074 if (deltaR(photon, p.momentum()) < 0.4) {
00075 econe += p.momentum().E();
00076
00077 if (econe/egamma > 0.07) {
00078 vetoEvent;
00079 }
00080 }
00081 }
00082
00083 _h_photon_pT->fill(photon.pT(),weight);
00084 _h_photon_y->fill(photon.rapidity(),weight);
00085
00086 const FastJets& jetpro = applyProjection<FastJets>(e, "Jets");
00087 const Jets& jets = jetpro.jetsByPt(20.0*GeV);
00088 if (jets.size()>0) {
00089 _h_photon_jet1_deta->fill(photon.eta()-jets[0].momentum().eta(), weight);
00090 _h_photon_jet1_dphi->fill(mapAngle0ToPi(photon.phi()-jets[0].momentum().phi()), weight);
00091 _h_photon_jet1_dR->fill(deltaR(photon, jets[0].momentum()), weight);
00092 }
00093
00094 MC_JetAnalysis::analyze(e);
00095 }
00096
00097
00098
00099 void finalize() {
00100 scale(_h_photon_pT, crossSectionPerEvent());
00101 scale(_h_photon_y, crossSectionPerEvent());
00102 scale(_h_photon_jet1_deta, crossSectionPerEvent());
00103 scale(_h_photon_jet1_dphi, crossSectionPerEvent());
00104 scale(_h_photon_jet1_dR, crossSectionPerEvent());
00105
00106 MC_JetAnalysis::finalize();
00107 }
00108
00109
00110
00111
00112 private:
00113
00114
00115
00116 AIDA::IHistogram1D * _h_photon_pT;
00117 AIDA::IHistogram1D * _h_photon_y;
00118 AIDA::IHistogram1D * _h_photon_jet1_deta;
00119 AIDA::IHistogram1D * _h_photon_jet1_dphi;
00120 AIDA::IHistogram1D * _h_photon_jet1_dR;
00121
00122
00123 };
00124
00125
00126
00127
00128 AnalysisBuilder<MC_PHOTONJETS> plugin_MC_PHOTONJETS;
00129
00130 }