Rivet analyses referenceALEPH_2004_I636645Jet rates and event shapes at LEP I and IIExperiment: ALEPH (LEP Run 1 and 2) Inspire ID: 636645 Status: VALIDATED Authors:
Beam energies: (45.6, 45.6); (66.5, 66.5); (80.5, 80.5); (86.0, 86.0); (91.5, 91.5); (94.5, 94.5); (98.5, 98.5); (100.0, 100.0); (103.0, 103.0) GeV Run details:
Jet rates, event-shape variables and inclusive charged particle spectra are measured in $e^+ e^-$ collisions at CMS energies between 91 and 209 GeV. The previously published data at 91.2 GeV and 133 GeV have been re-processed and the higher energy data are presented here for the first time. Note that the data have been corrected to include neutrinos. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: ALEPH_2004_I636645.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/FastJets.hh"
5#include "Rivet/Projections/ChargedFinalState.hh"
6#include "Rivet/Projections/Thrust.hh"
7#include "Rivet/Projections/Sphericity.hh"
8#include "Rivet/Projections/ParisiTensor.hh"
9#include "Rivet/Projections/Hemispheres.hh"
10#include "Rivet/Projections/Beam.hh"
11
12namespace Rivet {
13
14
15 /// ALEPH jet rates and event shapes at LEP 1 and 2
16 class ALEPH_2004_I636645 : public Analysis {
17 public:
18
19 RIVET_DEFAULT_ANALYSIS_CTOR(ALEPH_2004_I636645);
20
21
22 void init() {
23 _initialisedJets = true;
24 _initialisedSpectra = true;
25 // TODO: According to the paper they seem to discard neutral particles
26 // between 1 and 2 GeV. That correction is included in the systematic
27 // uncertainties and overly complicated to program, so we ignore it.
28 const FinalState fs;
29 declare(fs, "FS");
30 FastJets durhamjets(fs, JetAlg::DURHAM, 0.7, JetMuons::ALL, JetInvisibles::ALL);
31 declare(durhamjets, "DurhamJets");
32
33 const Thrust thrust(fs);
34 declare(thrust, "Thrust");
35 declare(Sphericity(fs), "Sphericity");
36 declare(ParisiTensor(fs), "Parisi");
37 declare(Hemispheres(thrust), "Hemispheres");
38
39 const ChargedFinalState cfs;
40 declare(Beam(), "Beams");
41 declare(cfs, "CFS");
42
43 // Histos
44 // offset for the event shapes and jets
45 int offset = 0;
46 switch (int(sqrtS()/GeV + 0.5)) {
47 case 91: offset = 0; break;
48 case 133: offset = 1; break;
49 case 161: offset = 2; break;
50 case 172: offset = 3; break;
51 case 183: offset = 4; break;
52 case 189: offset = 5; break;
53 case 200: offset = 6; break;
54 case 206: offset = 7; break;
55 default:
56 _initialisedJets = false;
57 }
58 // event shapes
59 if(_initialisedJets) {
60 book(_h_thrust ,offset+54, 1, 1);
61 book(_h_heavyjetmass ,offset+62, 1, 1);
62 book(_h_totaljetbroadening ,offset+70, 1, 1);
63 book(_h_widejetbroadening ,offset+78, 1, 1);
64 book(_h_cparameter ,offset+86, 1, 1);
65 book(_h_thrustmajor ,offset+94, 1, 1);
66 book(_h_thrustminor ,offset+102, 1, 1);
67 book(_h_jetmassdifference ,offset+110, 1, 1);
68 book(_h_aplanarity ,offset+118, 1, 1);
69 if ( offset != 0 )
70 book(_h_planarity, offset+125, 1, 1);
71 book(_h_oblateness ,offset+133, 1, 1);
72 book(_h_sphericity ,offset+141, 1, 1);
73
74 // Durham n->m jet resolutions
75 book(_h_y_Durham[0] ,offset+149, 1, 1); // y12 d149 ... d156
76 book(_h_y_Durham[1] ,offset+157, 1, 1); // y23 d157 ... d164
77 if (offset < 6) { // there is no y34, y45 and y56 for 200 gev
78 book(_h_y_Durham[2] ,offset+165, 1, 1); // y34 d165 ... d172, but not 171
79 book(_h_y_Durham[3] ,offset+173, 1, 1); // y45 d173 ... d179
80 book(_h_y_Durham[4] ,offset+180, 1, 1); // y56 d180 ... d186
81 }
82 else if (offset == 6) {
83 _h_y_Durham[2] = Histo1DPtr();
84 _h_y_Durham[3] = Histo1DPtr();
85 _h_y_Durham[4] = Histo1DPtr();
86 }
87 else if (offset == 7) {
88 book(_h_y_Durham[2] ,172, 1, 1);
89 book(_h_y_Durham[3] ,179, 1, 1);
90 book(_h_y_Durham[4] ,186, 1, 1);
91 }
92
93 // Durham n-jet fractions
94 book(_h_R_Durham[0] ,offset+187, 1, 1); // R1 d187 ... d194
95 book(_h_R_Durham[1] ,offset+195, 1, 1); // R2 d195 ... d202
96 book(_h_R_Durham[2] ,offset+203, 1, 1); // R3 d203 ... d210
97 book(_h_R_Durham[3] ,offset+211, 1, 1); // R4 d211 ... d218
98 book(_h_R_Durham[4] ,offset+219, 1, 1); // R5 d219 ... d226
99 book(_h_R_Durham[5] ,offset+227, 1, 1); // R>=6 d227 ... d234
100 }
101 // offset for the charged particle distributions
102 offset = 0;
103 switch (int(sqrtS() + 0.5)) {
104 case 133: offset = 0; break;
105 case 161: offset = 1; break;
106 case 172: offset = 2; break;
107 case 183: offset = 3; break;
108 case 189: offset = 4; break;
109 case 196: offset = 5; break;
110 case 200: offset = 6; break;
111 case 206: offset = 7; break;
112 default:
113 _initialisedSpectra = false;
114 }
115 if (_initialisedSpectra) {
116 book(_h_xp , 2+offset, 1, 1);
117 book(_h_xi ,11+offset, 1, 1);
118 book(_h_xe ,19+offset, 1, 1);
119 book(_h_pTin ,27+offset, 1, 1);
120 if (offset == 7)
121 book(_h_pTout, 35, 1, 1);
122 book(_h_rapidityT ,36+offset, 1, 1);
123 book(_h_rapidityS ,44+offset, 1, 1);
124 }
125 book(_weightedTotalChargedPartNum, "_weightedTotalChargedPartNum");
126 if (!_initialisedSpectra && !_initialisedJets) {
127 MSG_WARNING("CoM energy of events sqrt(s) = " << sqrtS()/GeV
128 << " doesn't match any available analysis energy .");
129 }
130
131 book(_mult, 1, 1, 1);
132 book(_mean_totaljetbroadening, 53, 1, 1);
133 book(_mean_widejetbroadening, 53, 1, 2);
134 book(_mean_C, 53, 1, 3);
135 book(_mean_rho, 53, 1, 4);
136 book(_mean_thrust, 53, 1, 5);
137 for (const string& en : _mult.binning().edges<0>()) {
138 double end = std::stod(en)*GeV;
139 if(isCompatibleWithSqrtS(end)) {
140 _ecms = en;
141 break;
142 }
143 }
144 }
145
146
147 void analyze(const Event& e) {
148
149 const Thrust& thrust = apply<Thrust>(e, "Thrust");
150 _mean_thrust->fill(_ecms,1.- thrust.thrust());
151 const Sphericity& sphericity = apply<Sphericity>(e, "Sphericity");
152 const Hemispheres& hemi = apply<Hemispheres>(e, "Hemispheres");
153 _mean_totaljetbroadening->fill(_ecms,hemi.Bsum());
154 _mean_widejetbroadening->fill(_ecms,hemi.Bmax());
155 _mean_rho->fill(_ecms,hemi.scaledM2high());
156 const ParisiTensor& parisi = apply<ParisiTensor>(e, "Parisi");
157 _mean_C->fill(_ecms,parisi.C());
158
159 if(_initialisedJets) {
160 bool LEP1 = isCompatibleWithSqrtS(91.2*GeV,0.01);
161 // event shapes
162 double thr = LEP1 ? thrust.thrust() : 1.0 - thrust.thrust();
163 _h_thrust->fill(thr);
164 _h_thrustmajor->fill(thrust.thrustMajor());
165 if(LEP1)
166 _h_thrustminor->fill(log(thrust.thrustMinor()));
167 else
168 _h_thrustminor->fill(thrust.thrustMinor());
169 _h_oblateness->fill(thrust.oblateness());
170
171 _h_heavyjetmass->fill(hemi.scaledM2high());
172 _h_jetmassdifference->fill(hemi.scaledM2diff());
173 _h_totaljetbroadening->fill(hemi.Bsum());
174 _h_widejetbroadening->fill(hemi.Bmax());
175 _h_cparameter->fill(parisi.C());
176
177 _h_aplanarity->fill(sphericity.aplanarity());
178 if(_h_planarity)
179 _h_planarity->fill(sphericity.planarity());
180 _h_sphericity->fill(sphericity.sphericity());
181
182 // Jet rates
183 const FastJets& durjet = apply<FastJets>(e, "DurhamJets");
184 double log10e = log10(exp(1.));
185 if (durjet.clusterSeq()) {
186 double logynm1=0.;
187 double logyn;
188 for (size_t i=0; i<5; ++i) {
189 const double yn = durjet.clusterSeq()->exclusive_ymerge_max(i+1);
190 if (yn<=0.0) continue;
191 logyn = -log(yn);
192 if (_h_y_Durham[i]) {
193 _h_y_Durham[i]->fill(logyn);
194 }
195 if(!LEP1) logyn *= log10e;
196 for (const auto& b : _h_R_Durham[i]->bins()) {
197 const double xMin = b.xMin();
198 if (-xMin <= logynm1) break;
199 if (-xMin<logyn) _h_R_Durham[i]->fill(b.xMid(), b.xWidth());
200 }
201 logynm1 = logyn;
202 }
203 for (const auto& b : _h_R_Durham[5]->bins()) {
204 if (-b.xMin() <= logynm1) break;
205 _h_R_Durham[5]->fill(b.xMid(), b.xWidth());
206 }
207 }
208 if( !_initialisedSpectra) {
209 const ChargedFinalState& cfs = apply<ChargedFinalState>(e, "CFS");
210 const size_t numParticles = cfs.particles().size();
211 _weightedTotalChargedPartNum->fill(numParticles);
212 }
213 }
214
215 // charged particle distributions
216 if(_initialisedSpectra) {
217 const ChargedFinalState& cfs = apply<ChargedFinalState>(e, "CFS");
218 const size_t numParticles = cfs.particles().size();
219 _weightedTotalChargedPartNum->fill(numParticles);
220 const ParticlePair& beams = apply<Beam>(e, "Beams").beams();
221 const double meanBeamMom = ( beams.first.p3().mod() +
222 beams.second.p3().mod() ) / 2.0;
223 for (const Particle& p : cfs.particles()) {
224 const double xp = p.p3().mod()/meanBeamMom;
225 _h_xp->fill(xp);
226 const double logxp = -std::log(xp);
227 _h_xi->fill(logxp);
228 const double xe = p.E()/meanBeamMom;
229 _h_xe->fill(xe);
230 const double pTinT = dot(p.p3(), thrust.thrustMajorAxis());
231 const double pToutT = dot(p.p3(), thrust.thrustMinorAxis());
232 _h_pTin->fill(fabs(pTinT/GeV));
233 if(_h_pTout) _h_pTout->fill(fabs(pToutT/GeV));
234 const double momT = dot(thrust.thrustAxis() ,p.p3());
235 const double rapidityT = 0.5 * std::log((p.E() + momT) /
236 (p.E() - momT));
237 _h_rapidityT->fill(fabs(rapidityT));
238 const double momS = dot(sphericity.sphericityAxis(),p.p3());
239 const double rapidityS = 0.5 * std::log((p.E() + momS) /
240 (p.E() - momS));
241 _h_rapidityS->fill(fabs(rapidityS));
242 }
243 }
244 }
245
246 void finalize() {
247 if(!_initialisedJets && !_initialisedSpectra) return;
248
249 if (_initialisedJets) {
250 normalize(_h_thrust);
251 normalize(_h_heavyjetmass);
252 normalize(_h_totaljetbroadening);
253 normalize(_h_widejetbroadening);
254 normalize(_h_cparameter);
255 normalize(_h_thrustmajor);
256 normalize(_h_thrustminor);
257 normalize(_h_jetmassdifference);
258 normalize(_h_aplanarity);
259 if(_h_planarity) normalize(_h_planarity);
260 normalize(_h_oblateness);
261 normalize(_h_sphericity);
262
263 for (size_t n=0; n<6; ++n) {
264 scale(_h_R_Durham[n], 1./sumOfWeights());
265 }
266
267 for (size_t n = 0; n < 5; ++n) {
268 if (_h_y_Durham[n]) {
269 scale(_h_y_Durham[n], 1.0/sumOfWeights());
270 }
271 }
272 }
273
274 const double avgNumParts = dbl(*_weightedTotalChargedPartNum) / sumOfWeights();
275
276
277 _mult->fill(_ecms,avgNumParts);
278
279 if (_initialisedSpectra) {
280 normalize(_h_xp, avgNumParts);
281 normalize(_h_xi, avgNumParts);
282 normalize(_h_xe, avgNumParts);
283 normalize(_h_pTin , avgNumParts);
284 if (_h_pTout) normalize(_h_pTout, avgNumParts);
285 normalize(_h_rapidityT, avgNumParts);
286 normalize(_h_rapidityS, avgNumParts);
287 }
288 }
289
290
291 private:
292
293 bool _initialisedJets = false;
294 bool _initialisedSpectra = false;
295
296 BinnedProfilePtr<string> _mult;
297 BinnedProfilePtr<string> _mean_totaljetbroadening,_mean_widejetbroadening,
298 _mean_C,_mean_rho,_mean_thrust;
299 string _ecms;
300 Histo1DPtr _h_xp;
301 Histo1DPtr _h_xi;
302 Histo1DPtr _h_xe;
303 Histo1DPtr _h_pTin;
304 Histo1DPtr _h_pTout;
305 Histo1DPtr _h_rapidityT;
306 Histo1DPtr _h_rapidityS;
307 Histo1DPtr _h_thrust;
308 Histo1DPtr _h_heavyjetmass;
309 Histo1DPtr _h_totaljetbroadening;
310 Histo1DPtr _h_widejetbroadening;
311 Histo1DPtr _h_cparameter;
312 Histo1DPtr _h_thrustmajor;
313 Histo1DPtr _h_thrustminor;
314 Histo1DPtr _h_jetmassdifference;
315 Histo1DPtr _h_aplanarity;
316 Histo1DPtr _h_planarity;
317 Histo1DPtr _h_oblateness;
318 Histo1DPtr _h_sphericity;
319
320 Histo1DPtr _h_R_Durham[6];
321 Histo1DPtr _h_y_Durham[5];
322
323 CounterPtr _weightedTotalChargedPartNum;
324
325 };
326
327
328
329 RIVET_DECLARE_ALIASED_PLUGIN(ALEPH_2004_I636645, ALEPH_2004_S5765862);
330
331}
|