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