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