Rivet analyses referenceOPAL_2004_I669402Event shape distributions and moments in $e^+ e^-$ -> hadrons at 91--209 GeVExperiment: OPAL (LEP 1 & 2) Inspire ID: 669402 Status: VALIDATED Authors:
Beam energies: (45.6, 45.6); (66.5, 66.5); (88.5, 88.5); (98.5, 98.5) GeV Run details:
Measurement of $e^+ e^-$ event shape variable distributions and their 1st to 5th moments in LEP running from the Z pole to the highest LEP 2 energy of 209 GeV.Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: OPAL_2004_I669402.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/Sphericity.hh"
7#include "Rivet/Projections/Thrust.hh"
8#include "Rivet/Projections/FastJets.hh"
9#include "Rivet/Projections/ParisiTensor.hh"
10#include "Rivet/Projections/Hemispheres.hh"
11#include <cmath>
12
13namespace Rivet {
14
15
16 /// @brief OPAL event shapes and moments at 91, 133, 177, and 197 GeV
17 ///
18 /// @author Andy Buckley
19 class OPAL_2004_I669402 : public Analysis {
20 public:
21
22 RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2004_I669402);
23
24
25 /// @name Analysis methods
26 /// @{
27
28 /// Energies: 91, 133, 177 (161-183), 197 (189-209) => index 0..4
29 int getHistIndex() {
30 int ih = -1;
31 if (inRange(sqrtS()/GeV, 89.9, 91.5)) {
32 ih = 0;
33 } else if (isCompatibleWithSqrtS(133*GeV)) {
34 ih = 1;
35 } else if (isCompatibleWithSqrtS(177*GeV)) { // (161-183)
36 ih = 2;
37 } else if (isCompatibleWithSqrtS(197*GeV)) { // (189-209)
38 ih = 3;
39 } else {
40 stringstream ss;
41 ss << "Invalid energy for OPAL_2004 analysis: "
42 << sqrtS() << " GeV != 91, 133, 177, or 197 GeV";
43 throw Error(ss.str());
44 }
45 assert(ih >= 0);
46 return ih;
47 }
48
49
50 void init() {
51 // Projections
52 declare(Beam(), "Beams");
53 const FinalState fs;
54 declare(fs, "FS");
55 const ChargedFinalState cfs;
56 declare(cfs, "CFS");
57 declare(FastJets(fs, JetAlg::DURHAM, 0.7), "DurhamJets");
58 declare(Sphericity(fs), "Sphericity");
59 declare(ParisiTensor(fs), "Parisi");
60 const Thrust thrust(fs);
61 declare(thrust, "Thrust");
62 declare(Hemispheres(thrust), "Hemispheres");
63 _isqrts = getHistIndex();
64
65 // Book histograms
66 book(_hist1MinusT[_isqrts] , 1, 1, _isqrts+1);
67 book(_histHemiMassH[_isqrts] , 2, 1, _isqrts+1);
68 book(_histCParam[_isqrts] , 3, 1, _isqrts+1);
69 book(_histHemiBroadT[_isqrts], 4, 1, _isqrts+1);
70 book(_histHemiBroadW[_isqrts], 5, 1, _isqrts+1);
71 book(_histY23Durham[_isqrts] , 6, 1, _isqrts+1);
72 book(_histTMajor[_isqrts] , 7, 1, _isqrts+1);
73 book(_histTMinor[_isqrts] , 8, 1, _isqrts+1);
74 book(_histAplanarity[_isqrts], 9, 1, _isqrts+1);
75 book(_histSphericity[_isqrts], 10, 1, _isqrts+1);
76 book(_histOblateness[_isqrts], 11, 1, _isqrts+1);
77 book(_histHemiMassL[_isqrts] , 12, 1, _isqrts+1);
78 book(_histHemiBroadN[_isqrts], 13, 1, _isqrts+1);
79 book(_histDParam[_isqrts] , 14, 1, _isqrts+1);
80 //
81 book(_hist1MinusTMom[_isqrts] , 15, 1, _isqrts+1);
82 book(_histHemiMassHMom[_isqrts] , 16, 1, _isqrts+1);
83 book(_histCParamMom[_isqrts] , 17, 1, _isqrts+1);
84 book(_histHemiBroadTMom[_isqrts], 18, 1, _isqrts+1);
85 book(_histHemiBroadWMom[_isqrts], 19, 1, _isqrts+1);
86 book(_histY23DurhamMom[_isqrts] , 20, 1, _isqrts+1);
87 book(_histTMajorMom[_isqrts] , 21, 1, _isqrts+1);
88 book(_histTMinorMom[_isqrts] , 22, 1, _isqrts+1);
89 book(_histSphericityMom[_isqrts], 23, 1, _isqrts+1);
90 book(_histOblatenessMom[_isqrts], 24, 1, _isqrts+1);
91 book(_histHemiMassLMom[_isqrts] , 25, 1, _isqrts+1);
92 book(_histHemiBroadNMom[_isqrts], 26, 1, _isqrts+1);
93
94 book(_sumWTrack2, "_sumWTrack2");
95 book(_sumWJet3, "_sumWJet3");
96 }
97
98
99 void analyze(const Event& event) {
100 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
101 const FinalState& cfs = apply<FinalState>(event, "CFS");
102 if (cfs.size() < 2) vetoEvent;
103
104 _sumWTrack2->fill();
105
106 // Thrusts
107 const Thrust& thrust = apply<Thrust>(event, "Thrust");
108 _hist1MinusT[_isqrts]->fill(1-thrust.thrust());
109 _histTMajor[_isqrts]->fill(thrust.thrustMajor());
110 _histTMinor[_isqrts]->fill(thrust.thrustMinor());
111 _histOblateness[_isqrts]->fill(thrust.oblateness());
112 for (int n = 1; n <= 5; ++n) {
113 _hist1MinusTMom[_isqrts]->fill(n, pow(1-thrust.thrust(), n));
114 _histTMajorMom[_isqrts]->fill(n, pow(thrust.thrustMajor(), n));
115 _histTMinorMom[_isqrts]->fill(n, pow(thrust.thrustMinor(), n));
116 _histOblatenessMom[_isqrts]->fill(n, pow(thrust.oblateness(), n));
117 }
118
119 // Jets
120 const FastJets& durjet = apply<FastJets>(event, "DurhamJets");
121 if (durjet.clusterSeq()) {
122 _sumWJet3->fill();
123 const double y23 = durjet.clusterSeq()->exclusive_ymerge_max(2);
124 if (y23>0.0) {
125 _histY23Durham[_isqrts]->fill(y23);
126 for (int n = 1; n <= 5; ++n) {
127 _histY23DurhamMom[_isqrts]->fill(n, pow(y23, n));
128 }
129 }
130 }
131
132 // Sphericities
133 const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
134 const double sph = sphericity.sphericity();
135 const double apl = sphericity.aplanarity();
136 _histSphericity[_isqrts]->fill(sph);
137 _histAplanarity[_isqrts]->fill(apl);
138 for (int n = 1; n <= 5; ++n) {
139 _histSphericityMom[_isqrts]->fill(n, pow(sph, n));
140 }
141
142 // C & D params
143 const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
144 const double cparam = parisi.C();
145 const double dparam = parisi.D();
146 _histCParam[_isqrts]->fill(cparam);
147 _histDParam[_isqrts]->fill(dparam);
148 for (int n = 1; n <= 5; ++n) {
149 _histCParamMom[_isqrts]->fill(n, pow(cparam, n));
150 }
151
152 // Hemispheres
153 const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
154 // The paper says that M_H/L are scaled by sqrt(s), but scaling by E_vis is the way that fits the data...
155 const double hemi_mh = hemi.scaledMhigh();
156 const double hemi_ml = hemi.scaledMlow();
157 /// @todo This shouldn't be necessary... what's going on? Memory corruption suspected :(
158 // if (std::isnan(hemi_ml)) {
159 // MSG_ERROR("NaN in HemiL! Event = " << numEvents());
160 // MSG_ERROR(hemi.M2low() << ", " << hemi.E2vis());
161 // }
162 if (!std::isnan(hemi_mh) && !std::isnan(hemi_ml)) {
163 const double hemi_bmax = hemi.Bmax();
164 const double hemi_bmin = hemi.Bmin();
165 const double hemi_bsum = hemi.Bsum();
166 _histHemiMassH[_isqrts]->fill(hemi_mh);
167 _histHemiMassL[_isqrts]->fill(hemi_ml);
168 _histHemiBroadW[_isqrts]->fill(hemi_bmax);
169 _histHemiBroadN[_isqrts]->fill(hemi_bmin);
170 _histHemiBroadT[_isqrts]->fill(hemi_bsum);
171 for (int n = 1; n <= 5; ++n) {
172 // if (std::isnan(pow(hemi_ml, n))) MSG_ERROR("NaN in HemiL moment! Event = " << numEvents());
173 _histHemiMassHMom[_isqrts]->fill(n, pow(hemi_mh, n));
174 _histHemiMassLMom[_isqrts]->fill(n, pow(hemi_ml, n));
175 _histHemiBroadWMom[_isqrts]->fill(n, pow(hemi_bmax, n));
176 _histHemiBroadNMom[_isqrts]->fill(n, pow(hemi_bmin, n));
177 _histHemiBroadTMom[_isqrts]->fill(n, pow(hemi_bsum, n));
178 }
179 }
180 }
181
182
183 void finalize() {
184 scale(_hist1MinusT[_isqrts], 1.0 / *_sumWTrack2);
185 scale(_histTMajor[_isqrts], 1.0 / *_sumWTrack2);
186 scale(_histTMinor[_isqrts], 1.0 / *_sumWTrack2);
187 scale(_histOblateness[_isqrts], 1.0 / *_sumWTrack2);
188 scale(_histSphericity[_isqrts], 1.0 / *_sumWTrack2);
189 scale(_histAplanarity[_isqrts], 1.0 / *_sumWTrack2);
190 scale(_histHemiMassH[_isqrts], 1.0 / *_sumWTrack2);
191 scale(_histHemiMassL[_isqrts], 1.0 / *_sumWTrack2);
192 scale(_histHemiBroadW[_isqrts], 1.0 / *_sumWTrack2);
193 scale(_histHemiBroadN[_isqrts], 1.0 / *_sumWTrack2);
194 scale(_histHemiBroadT[_isqrts], 1.0 / *_sumWTrack2);
195 scale(_histCParam[_isqrts], 1.0 / *_sumWTrack2);
196 scale(_histDParam[_isqrts], 1.0 / *_sumWTrack2);
197 scale(_histY23Durham[_isqrts], 1.0 / *_sumWJet3);
198 //
199 scale(_hist1MinusTMom[_isqrts], 1.0 / *_sumWTrack2);
200 scale(_histTMajorMom[_isqrts], 1.0 / *_sumWTrack2);
201 scale(_histTMinorMom[_isqrts], 1.0 / *_sumWTrack2);
202 scale(_histOblatenessMom[_isqrts], 1.0 / *_sumWTrack2);
203 scale(_histSphericityMom[_isqrts], 1.0 / *_sumWTrack2);
204 scale(_histHemiMassHMom[_isqrts], 1.0 / *_sumWTrack2);
205 scale(_histHemiMassLMom[_isqrts], 1.0 / *_sumWTrack2);
206 scale(_histHemiBroadWMom[_isqrts], 1.0 / *_sumWTrack2);
207 scale(_histHemiBroadNMom[_isqrts], 1.0 / *_sumWTrack2);
208 scale(_histHemiBroadTMom[_isqrts], 1.0 / *_sumWTrack2);
209 scale(_histCParamMom[_isqrts], 1.0 / *_sumWTrack2);
210 scale(_histY23DurhamMom[_isqrts], 1.0 / *_sumWJet3);
211 }
212
213 /// @}
214
215
216 private:
217
218 /// Beam energy index for histograms
219 int _isqrts = -1;
220
221 /// Counters of event weights passing the cuts
222 /// @{
223 CounterPtr _sumWTrack2, _sumWJet3;
224 /// @}
225
226 /// @name Event shape histos at 4 energies
227 /// @{
228 Histo1DPtr _hist1MinusT[4];
229 Histo1DPtr _histHemiMassH[4];
230 Histo1DPtr _histCParam[4];
231 Histo1DPtr _histHemiBroadT[4];
232 Histo1DPtr _histHemiBroadW[4];
233 Histo1DPtr _histY23Durham[4];
234 Histo1DPtr _histTMajor[4];
235 Histo1DPtr _histTMinor[4];
236 Histo1DPtr _histAplanarity[4];
237 Histo1DPtr _histSphericity[4];
238 Histo1DPtr _histOblateness[4];
239 Histo1DPtr _histHemiMassL[4];
240 Histo1DPtr _histHemiBroadN[4];
241 Histo1DPtr _histDParam[4];
242 /// @}
243
244 /// @name Event shape moment histos at 4 energies
245 /// @{
246 BinnedHistoPtr<int> _hist1MinusTMom[4];
247 BinnedHistoPtr<int> _histHemiMassHMom[4];
248 BinnedHistoPtr<int> _histCParamMom[4];
249 BinnedHistoPtr<int> _histHemiBroadTMom[4];
250 BinnedHistoPtr<int> _histHemiBroadWMom[4];
251 BinnedHistoPtr<int> _histY23DurhamMom[4];
252 BinnedHistoPtr<int> _histTMajorMom[4];
253 BinnedHistoPtr<int> _histTMinorMom[4];
254 BinnedHistoPtr<int> _histSphericityMom[4];
255 BinnedHistoPtr<int> _histOblatenessMom[4];
256 BinnedHistoPtr<int> _histHemiMassLMom[4];
257 BinnedHistoPtr<int> _histHemiBroadNMom[4];
258 /// @}
259
260 };
261
262
263
264 RIVET_DECLARE_ALIASED_PLUGIN(OPAL_2004_I669402, OPAL_2004_S6132243);
265
266}
|