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
150
151
152
153
154
155
156
157
158
| // -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "fastjet/JadePlugin.hh"
namespace Rivet {
/// @brief OPAL photon production
///
/// @author Peter Richardson
class OPAL_1993_S2692198 : public Analysis {
public:
RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1993_S2692198);
/// @name Analysis methods
/// @{
void analyze(const Event& e) {
// Extract the photons
Particles photons;
Particles nonPhotons;
FourMomentum ptotal;
const FinalState& fs = apply<FinalState>(e, "FS");
for (const Particle& p : fs.particles()) {
ptotal+= p.momentum();
if (p.pid() == PID::PHOTON) {
photons.push_back(p);
} else {
nonPhotons.push_back(p);
}
}
// No photon return but still count for normalisation
if (photons.empty()) return;
// Definition of the Durham algorithm
fastjet::JetDefinition durham_def(fastjet::ee_kt_algorithm, fastjet::E_scheme, fastjet::Best);
// Definition of the JADE algorithm
fastjet::JadePlugin jade;
fastjet::JetDefinition jade_def = fastjet::JetDefinition(&jade);
// Now for the weird jet algorithm
double evis = ptotal.mass();
vector<fastjet::PseudoJet> input_particles;
// Pseudo-jets from the non photons
for (const Particle& p : nonPhotons) {
const FourMomentum p4 = p.momentum();
input_particles.push_back(fastjet::PseudoJet(p4.px(), p4.py(), p4.pz(), p4.E()));
}
// Pseudo-jets from all bar the first photon
for (size_t ix = 1; ix < photons.size(); ++ix) {
const FourMomentum p4 = photons[ix].momentum();
input_particles.push_back(fastjet::PseudoJet(p4.px(), p4.py(), p4.pz(), p4.E()));
}
// Now loop over the photons
for (size_t ix = 0; ix < photons.size(); ++ix) {
FourMomentum pgamma = photons[ix].momentum();
// Run the jet clustering DURHAM
fastjet::ClusterSequence clust_seq(input_particles, durham_def);
// Cluster the jets
for (size_t j = 0; j < _nPhotonDurham->numBins(); ++j) {
bool accept(true);
double ycut = _nPhotonDurham->bin(j).xMid(); ///< @todo Should this be xMin?
double dcut = sqr(evis)*ycut;
vector<fastjet::PseudoJet> exclusive_jets = sorted_by_E(clust_seq.exclusive_jets(dcut));
for (size_t iy = 0; iy < exclusive_jets.size(); ++iy) {
FourMomentum pjet(momentum(exclusive_jets[iy]));
double cost = pjet.p3().unit().dot(pgamma.p3().unit());
double ygamma = 2 * min(sqr(pjet.E()/evis), sqr(pgamma.E()/evis)) * (1 - cost);
if (ygamma < ycut) {
accept = false;
break;
}
}
if (!accept) continue;
_nPhotonDurham->fill(ycut, _nPhotonDurham->bin(j).xWidth());
size_t njet = min(size_t(4), exclusive_jets.size()) - 1;
if (j < _nPhotonJetDurham[njet]->numBins()) {
_nPhotonJetDurham[njet]->fillBin(j, _nPhotonJetDurham[njet]->bin(j).xWidth());
}
}
// Run the jet clustering JADE
fastjet::ClusterSequence clust_seq2(input_particles, jade_def);
for (size_t j = 0; j < _nPhotonJade->numBins(); ++j) {
bool accept(true);
double ycut = _nPhotonJade->bin(j).xMid(); ///< @todo Should this be xMin?
double dcut = sqr(evis)*ycut;
vector<fastjet::PseudoJet> exclusive_jets = sorted_by_E(clust_seq2.exclusive_jets(dcut));
for (size_t iy = 0; iy < exclusive_jets.size(); ++iy) {
FourMomentum pjet(momentum(exclusive_jets[iy]));
double cost = pjet.p3().unit().dot(pgamma.p3().unit());
double ygamma = 2.*pjet.E()*pgamma.E()/sqr(evis)*(1.-cost);
if (ygamma < ycut) {
accept = false;
break;
}
}
if (!accept) continue;
/// @todo Really want to use a "bar graph" here (i.e. ignore bin width)
_nPhotonJade->fill(ycut, _nPhotonJade->bin(j).xWidth());
size_t njet = min(size_t(4), exclusive_jets.size()) - 1;
if (j < _nPhotonJetJade[njet]->numBins()) {
_nPhotonJetJade[njet]->fillBin(j, _nPhotonJetJade[njet]->bin(j).xWidth());
}
}
// Add this photon back in for the next iteration of the loop
if (ix+1 != photons.size()) {
input_particles[nonPhotons.size()+ix] = fastjet::PseudoJet(pgamma.px(), pgamma.py(), pgamma.pz(), pgamma.E());
}
}
}
void init() {
// Projections
declare(FinalState(), "FS");
// Book datasets
book(_nPhotonJade ,1, 1, 1);
book(_nPhotonDurham ,2, 1, 1);
for (size_t ix = 0; ix < 4; ++ix) {
book(_nPhotonJetJade [ix] ,3 , 1, 1+ix);
book(_nPhotonJetDurham[ix] ,4 , 1, 1+ix);
}
}
/// Finalize
void finalize() {
const double fact = 1000/sumOfWeights();
scale(_nPhotonJade, fact);
scale(_nPhotonDurham, fact);
for (size_t ix = 0; ix < 4; ++ix) {
scale(_nPhotonJetJade [ix],fact);
scale(_nPhotonJetDurham[ix],fact);
}
}
/// @}
private:
Histo1DPtr _nPhotonJade;
Histo1DPtr _nPhotonDurham;
Histo1DPtr _nPhotonJetJade [4];
Histo1DPtr _nPhotonJetDurham[4];
};
RIVET_DECLARE_ALIASED_PLUGIN(OPAL_1993_S2692198, OPAL_1993_I343181);
}
|