Rivet analyses referenceOPAL_1993_S2692198Measurement of photon production at LEP 1Experiment: OPAL (LEP Run 1) Inspire ID: 343181 Status: UNVALIDATED Authors:
Beam energies: (45.6, 45.6) GeV Run details:
Measurement of the production of photons in $e^+ e^-\to q \bar q$ events at LEP 1. Source code: OPAL_1993_S2692198.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/ChargedFinalState.hh"
6#include "Rivet/Projections/FastJets.hh"
7#include "fastjet/JadePlugin.hh"
8
9namespace Rivet {
10
11
12 /// @brief OPAL photon production
13 ///
14 /// @author Peter Richardson
15 class OPAL_1993_S2692198 : public Analysis {
16 public:
17
18 RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1993_S2692198);
19
20
21 /// @name Analysis methods
22 /// @{
23
24 void analyze(const Event& e) {
25 // Extract the photons
26 Particles photons;
27 Particles nonPhotons;
28 FourMomentum ptotal;
29 const FinalState& fs = apply<FinalState>(e, "FS");
30 for (const Particle& p : fs.particles()) {
31 ptotal+= p.momentum();
32 if (p.pid() == PID::PHOTON) {
33 photons.push_back(p);
34 } else {
35 nonPhotons.push_back(p);
36 }
37 }
38 // No photon return but still count for normalisation
39 if (photons.empty()) return;
40 // Definition of the Durham algorithm
41 fastjet::JetDefinition durham_def(fastjet::ee_kt_algorithm, fastjet::E_scheme, fastjet::Best);
42 // Definition of the JADE algorithm
43 fastjet::JadePlugin jade;
44 fastjet::JetDefinition jade_def = fastjet::JetDefinition(&jade);
45 // Now for the weird jet algorithm
46 double evis = ptotal.mass();
47 vector<fastjet::PseudoJet> input_particles;
48 // Pseudo-jets from the non photons
49 for (const Particle& p : nonPhotons) {
50 const FourMomentum p4 = p.momentum();
51 input_particles.push_back(fastjet::PseudoJet(p4.px(), p4.py(), p4.pz(), p4.E()));
52 }
53 // Pseudo-jets from all bar the first photon
54 for (size_t ix = 1; ix < photons.size(); ++ix) {
55 const FourMomentum p4 = photons[ix].momentum();
56 input_particles.push_back(fastjet::PseudoJet(p4.px(), p4.py(), p4.pz(), p4.E()));
57 }
58 // Now loop over the photons
59 for (size_t ix = 0; ix < photons.size(); ++ix) {
60 FourMomentum pgamma = photons[ix].momentum();
61 // Run the jet clustering DURHAM
62 fastjet::ClusterSequence clust_seq(input_particles, durham_def);
63 // Cluster the jets
64 for (size_t j = 0; j < _nPhotonDurham->numBins(); ++j) {
65 bool accept(true);
66 double ycut = _nPhotonDurham->bin(j).xMid(); ///< @todo Should this be xMin?
67 double dcut = sqr(evis)*ycut;
68 vector<fastjet::PseudoJet> exclusive_jets = sorted_by_E(clust_seq.exclusive_jets(dcut));
69 for (size_t iy = 0; iy < exclusive_jets.size(); ++iy) {
70 FourMomentum pjet(momentum(exclusive_jets[iy]));
71 double cost = pjet.p3().unit().dot(pgamma.p3().unit());
72 double ygamma = 2 * min(sqr(pjet.E()/evis), sqr(pgamma.E()/evis)) * (1 - cost);
73 if (ygamma < ycut) {
74 accept = false;
75 break;
76 }
77 }
78 if (!accept) continue;
79 _nPhotonDurham->fill(ycut, _nPhotonDurham->bin(j).xWidth());
80 size_t njet = min(size_t(4), exclusive_jets.size()) - 1;
81 if (j < _nPhotonJetDurham[njet]->numBins()) {
82 _nPhotonJetDurham[njet]->fillBin(j, _nPhotonJetDurham[njet]->bin(j).xWidth());
83 }
84 }
85 // Run the jet clustering JADE
86 fastjet::ClusterSequence clust_seq2(input_particles, jade_def);
87 for (size_t j = 0; j < _nPhotonJade->numBins(); ++j) {
88 bool accept(true);
89 double ycut = _nPhotonJade->bin(j).xMid(); ///< @todo Should this be xMin?
90 double dcut = sqr(evis)*ycut;
91 vector<fastjet::PseudoJet> exclusive_jets = sorted_by_E(clust_seq2.exclusive_jets(dcut));
92 for (size_t iy = 0; iy < exclusive_jets.size(); ++iy) {
93 FourMomentum pjet(momentum(exclusive_jets[iy]));
94 double cost = pjet.p3().unit().dot(pgamma.p3().unit());
95 double ygamma = 2.*pjet.E()*pgamma.E()/sqr(evis)*(1.-cost);
96 if (ygamma < ycut) {
97 accept = false;
98 break;
99 }
100 }
101 if (!accept) continue;
102 /// @todo Really want to use a "bar graph" here (i.e. ignore bin width)
103 _nPhotonJade->fill(ycut, _nPhotonJade->bin(j).xWidth());
104 size_t njet = min(size_t(4), exclusive_jets.size()) - 1;
105 if (j < _nPhotonJetJade[njet]->numBins()) {
106 _nPhotonJetJade[njet]->fillBin(j, _nPhotonJetJade[njet]->bin(j).xWidth());
107 }
108 }
109 // Add this photon back in for the next iteration of the loop
110 if (ix+1 != photons.size()) {
111 input_particles[nonPhotons.size()+ix] = fastjet::PseudoJet(pgamma.px(), pgamma.py(), pgamma.pz(), pgamma.E());
112 }
113 }
114 }
115
116
117 void init() {
118 // Projections
119 declare(FinalState(), "FS");
120
121 // Book datasets
122 book(_nPhotonJade ,1, 1, 1);
123 book(_nPhotonDurham ,2, 1, 1);
124 for (size_t ix = 0; ix < 4; ++ix) {
125 book(_nPhotonJetJade [ix] ,3 , 1, 1+ix);
126 book(_nPhotonJetDurham[ix] ,4 , 1, 1+ix);
127 }
128 }
129
130
131 /// Finalize
132 void finalize() {
133 const double fact = 1000/sumOfWeights();
134 scale(_nPhotonJade, fact);
135 scale(_nPhotonDurham, fact);
136 for (size_t ix = 0; ix < 4; ++ix) {
137 scale(_nPhotonJetJade [ix],fact);
138 scale(_nPhotonJetDurham[ix],fact);
139 }
140 }
141
142 /// @}
143
144
145 private:
146
147 Histo1DPtr _nPhotonJade;
148 Histo1DPtr _nPhotonDurham;
149 Histo1DPtr _nPhotonJetJade [4];
150 Histo1DPtr _nPhotonJetDurham[4];
151
152 };
153
154
155
156 RIVET_DECLARE_ALIASED_PLUGIN(OPAL_1993_S2692198, OPAL_1993_I343181);
157
158}
|