Rivet analyses referenceL3_2004_I652683Jet rates and event shapes at LEP I+IIExperiment: L3 (LEP I+II) Inspire ID: 652683 Status: VALIDATED Authors:
Beam energies: (20.7, 20.7); (27.6, 27.6); (32.7, 32.7); (37.9, 37.9); (41.1, 41.1); (42.5, 42.5); (45.6, 45.6); (65.0, 65.0); (68.0, 68.0); (80.7, 80.7); (86.2, 86.2); (91.4, 91.4); (94.3, 94.3); (97.2, 97.2); (100.1, 100.1); (103.1, 103.1) GeV Run details:
Event shapes, charged particle multiplicity and scaled momentum distributions, flavour separated at the Z-boson peak. The analysis also includes non-flavour separated event shapes, charged particle multiplicity and scaled momentum distributions at events both above and below the Z pole. There are also jet distributions in the paper which are not currently implemented. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: L3_2004_I652683.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/FastJets.hh"
5#include "Rivet/Projections/FinalState.hh"
6#include "Rivet/Projections/ChargedFinalState.hh"
7#include "Rivet/Projections/Thrust.hh"
8#include "Rivet/Projections/ParisiTensor.hh"
9#include "Rivet/Projections/Hemispheres.hh"
10
11#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
12#include "Rivet/Projections/InitialQuarks.hh"
13#include "fastjet/EECambridgePlugin.hh"
14
15namespace Rivet {
16
17
18 /// Jet rates and event shapes at LEP I+II
19 class L3_2004_I652683 : public Analysis {
20 public:
21
22 /// Constructor
23 RIVET_DEFAULT_ANALYSIS_CTOR(L3_2004_I652683);
24
25 /// Book histograms and initialise projections before the run
26 void init() {
27 // Projections to use
28 const FinalState FS;
29 declare(FS, "FS");
30 declare(Beam(), "beams");
31 const ChargedFinalState CFS;
32 declare(CFS, "CFS");
33 const Thrust thrust(FS);
34 declare(thrust, "thrust");
35 declare(ParisiTensor(FS), "Parisi");
36 declare(Hemispheres(thrust), "Hemispheres");
37 declare(InitialQuarks(), "initialquarks");
38 FastJets jadeJets = FastJets(FS, JetAlg::JADE, 0.7, JetMuons::ALL, JetInvisibles::DECAY);
39 FastJets durhamJets = FastJets(FS, JetAlg::DURHAM, 0.7, JetMuons::ALL, JetInvisibles::DECAY);
40 FastJets cambridgeJets = FastJets(FS, JetAlg::CAM, 0.7, JetMuons::ALL, JetInvisibles::DECAY);
41 declare(jadeJets, "JadeJets");
42 declare(durhamJets, "DurhamJets");
43
44 // Book the histograms
45 if (isCompatibleWithSqrtS(91.2*GeV)) {
46 // z pole
47 book(_h_Thrust_udsc , 47, 1, 1);
48 book(_h_Thrust_bottom , 47, 1, 2);
49 book(_h_heavyJetmass_udsc , 48, 1, 1);
50 book(_h_heavyJetmass_bottom , 48, 1, 2);
51 book(_h_totalJetbroad_udsc , 49, 1, 1);
52 book(_h_totalJetbroad_bottom , 49, 1, 2);
53 book(_h_wideJetbroad_udsc , 50, 1, 1);
54 book(_h_wideJetbroad_bottom , 50, 1, 2);
55 book(_h_Cparameter_udsc , 51, 1, 1);
56 book(_h_Cparameter_bottom , 51, 1, 2);
57 book(_h_Dparameter_udsc , 52, 1, 1);
58 book(_h_Dparameter_bottom , 52, 1, 2);
59 book(_h_Ncharged , "/TMP/NCHARGED" , 28, 1, 57);
60 book(_h_Ncharged_udsc , "/TMP/NCHARGED_UDSC", 28, 1, 57);
61 book(_h_Ncharged_bottom , "/TMP/NCHARGED_B" , 27, 3, 57);
62 book(_h_scaledMomentum , 65, 1, 1);
63 book(_h_scaledMomentum_udsc , 65, 1, 2);
64 book(_h_scaledMomentum_bottom , 65, 1, 3);
65 }
66 else if(sqrtS()/GeV < 90) {
67 int i1(-1),i2(-1);
68 if(isCompatibleWithSqrtS(41.4*GeV)) {
69 i1=0; i2=1;
70 }
71 else if(isCompatibleWithSqrtS(55.3*GeV)) {
72 i1=0; i2=2;
73 }
74 else if(isCompatibleWithSqrtS(65.4*GeV)) {
75 i1=0; i2=3;
76 }
77 else if(isCompatibleWithSqrtS(75.7*GeV)) {
78 i1=1; i2=1;
79 }
80 else if(isCompatibleWithSqrtS(82.3*GeV)) {
81 i1=1; i2=2;
82 }
83 else if(isCompatibleWithSqrtS(85.1*GeV)) {
84 i1=1; i2=3;
85 }
86 else {
87 MSG_ERROR("Beam energy not supported!");
88 }
89
90 book(_h_thrust , 21+i1,1,i2);
91 book(_h_rho , 26+i1,1,i2);
92 book(_h_B_T , 31+i1,1,i2);
93 book(_h_B_W , 36+i1,1,i2);
94 }
95 else if(sqrtS()/GeV > 120) {
96 int i1(-1),i2(-1);
97 if(isCompatibleWithSqrtS(130.1*GeV)) {
98 i1=0; i2=1;
99 }
100 else if(isCompatibleWithSqrtS(136.1*GeV)) {
101 i1=0; i2=2;
102 }
103 else if(isCompatibleWithSqrtS(161.3*GeV)) {
104 i1=0; i2=3;
105 }
106 else if(isCompatibleWithSqrtS(172.3*GeV)) {
107 i1=1; i2=1;
108 }
109 else if(isCompatibleWithSqrtS(182.8*GeV)) {
110 i1=1; i2=2;
111 }
112 else if(isCompatibleWithSqrtS(188.6*GeV)) {
113 i1=1; i2=3;
114 }
115 else if(isCompatibleWithSqrtS(194.4*GeV)) {
116 i1=2; i2=1;
117 }
118 else if(isCompatibleWithSqrtS(200.2*GeV)) {
119 i1=2; i2=2;
120 }
121 else if(isCompatibleWithSqrtS(206.2*GeV)) {
122 i1=2; i2=3;
123 }
124 else {
125 MSG_ERROR("Beam energy not supported!");
126 }
127
128 book(_h_thrust , 23+i1,1,i2);
129 book(_h_rho , 28+i1,1,i2);
130 book(_h_B_T , 33+i1,1,i2);
131 book(_h_B_W , 38+i1,1,i2);
132 book(_h_C , 41+i1,1,i2);
133 book(_h_D , 44+i1,1,i2);
134 book(_h_N , "/TMP/NCHARGED", 22, 9, 53);
135 book(_h_xi , 66+i1,1,i2);
136 // and the jets
137 int i3 = 3*i1+i2;
138 book(_s["y_2_JADE"], i3,1,1);
139 book(_s["y_3_JADE"], i3,1,2);
140 book(_s["y_4_JADE"], i3,1,3);
141 book(_s["y_5_JADE"], i3,1,4);
142 book(_s["y_2_Durham"], 9+i3,1,1);
143 book(_s["y_3_Durham"], 9+i3,1,2);
144 book(_s["y_4_Durham"], 9+i3,1,3);
145 book(_s["y_5_Durham"], 9+i3,1,4);
146 if (i3==8 || i3==9) {
147 book(_s["y_2_Cambridge"], 11+i3,1,1);
148 book(_s["y_3_Cambridge"], 11+i3,1,2);
149 book(_s["y_4_Cambridge"], 11+i3,1,3);
150 book(_s["y_5_Cambridge"], 11+i3,1,4);
151 }
152 }
153
154 book(_sumW_udsc, "_sumW_udsc");
155 book(_sumW_b, "_sumW_b");
156 book(_sumW_ch, "_sumW_ch");
157 book(_sumW_ch_udsc, "_sumW_ch_udsc");
158 book(_sumW_ch_b, "_sumW_ch_b");
159
160 }
161
162
163 /// Perform the per-event analysis
164 void analyze(const Event& event) {
165 // Get beam average momentum
166 const ParticlePair& beams = apply<Beam>(event, "beams").beams();
167 const double beamMomentum = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
168
169 // InitialQuarks projection to have udsc events separated from b events
170 /// @todo Yuck!!! Eliminate when possible...
171 int iflav = 0;
172 // only need the flavour at Z pole
173 if (_h_Thrust_udsc) {
174 int flavour = 0;
175 const InitialQuarks& iqf = apply<InitialQuarks>(event, "initialquarks");
176 Particles quarks;
177 if ( iqf.particles().size() == 2 ) {
178 flavour = iqf.particles().front().abspid();
179 quarks = iqf.particles();
180 } else {
181 map<int, Particle> quarkmap;
182 for (const Particle& p : iqf.particles()) {
183 if (quarkmap.find(p.pid()) == quarkmap.end()) quarkmap[p.pid()] = p;
184 else if (quarkmap[p.pid()].E() < p.E()) quarkmap[p.pid()] = p;
185 }
186 double max_energy = 0.;
187 for (int i = 1; i <= 5; ++i) {
188 double energy = 0.;
189 if (quarkmap.find(i) != quarkmap.end())
190 energy += quarkmap[ i].E();
191 if (quarkmap.find(-i) != quarkmap.end())
192 energy += quarkmap[-i].E();
193 if (energy > max_energy)
194 flavour = i;
195 }
196 if (quarkmap.find(flavour) != quarkmap.end())
197 quarks.push_back(quarkmap[flavour]);
198 if (quarkmap.find(-flavour) != quarkmap.end())
199 quarks.push_back(quarkmap[-flavour]);
200 }
201 // Flavour label
202 /// @todo Change to a bool?
203 iflav = (flavour == PID::DQUARK || flavour == PID::UQUARK || flavour == PID::SQUARK || flavour == PID::CQUARK) ? 1 : (flavour == PID::BQUARK) ? 5 : 0;
204 }
205 // Update weight sums
206 if (iflav == 1) {
207 _sumW_udsc->fill();
208 } else if (iflav == 5) {
209 _sumW_b->fill();
210 }
211 _sumW_ch->fill();
212
213 // Charged multiplicity
214 const FinalState& cfs = apply<FinalState>(event, "CFS");
215 if(_h_Ncharged) _h_Ncharged->fill(cfs.size());
216 if (iflav == 1) {
217 _sumW_ch_udsc->fill();
218 _h_Ncharged_udsc->fill(cfs.size());
219 } else if (iflav == 5) {
220 _sumW_ch_b->fill();
221 _h_Ncharged_bottom->fill(cfs.size());
222 }
223 else if(_h_N) {
224 _h_N->fill(cfs.size());
225 }
226
227 // Scaled momentum
228 const Particles& chparticles = cfs.particlesByPt();
229 for (const Particle& p : chparticles) {
230 const Vector3 momentum3 = p.p3();
231 const double mom = momentum3.mod();
232 const double scaledMom = mom/beamMomentum;
233 const double logScaledMom = std::log(scaledMom);
234 if(_h_scaledMomentum) _h_scaledMomentum->fill(-logScaledMom);
235 if (iflav == 1) {
236 _h_scaledMomentum_udsc->fill(-logScaledMom);
237 } else if (iflav == 5) {
238 _h_scaledMomentum_bottom->fill(-logScaledMom);
239 }
240 else if(_h_xi) {
241 _h_xi->fill(-logScaledMom);
242 }
243 }
244
245 // Thrust
246 const Thrust& thrust = apply<Thrust>(event, "thrust");
247 if (iflav == 1) {
248 _h_Thrust_udsc->fill(thrust.thrust());
249 } else if (iflav == 5) {
250 _h_Thrust_bottom->fill(thrust.thrust());
251 }
252 else if(_h_thrust) {
253 _h_thrust->fill(1.-thrust.thrust());
254 }
255
256 // C and D Parisi parameters
257 const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
258 if (iflav == 1) {
259 _h_Cparameter_udsc->fill(parisi.C());
260 _h_Dparameter_udsc->fill(parisi.D());
261 } else if (iflav == 5) {
262 _h_Cparameter_bottom->fill(parisi.C());
263 _h_Dparameter_bottom->fill(parisi.D());
264 }
265 else if(_h_C) {
266 _h_C->fill(parisi.C());
267 _h_D->fill(parisi.D());
268 }
269
270 // The hemisphere variables
271 const Hemispheres& hemisphere = apply<Hemispheres>(event, "Hemispheres");
272 if (iflav == 1) {
273 _h_heavyJetmass_udsc->fill(hemisphere.scaledM2high());
274 _h_totalJetbroad_udsc->fill(hemisphere.Bsum());
275 _h_wideJetbroad_udsc->fill(hemisphere.Bmax());
276 } else if (iflav == 5) {
277 _h_heavyJetmass_bottom->fill(hemisphere.scaledM2high());
278 _h_totalJetbroad_bottom->fill(hemisphere.Bsum());
279 _h_wideJetbroad_bottom->fill(hemisphere.Bmax());
280 }
281 else if (_h_rho) {
282 _h_rho->fill(hemisphere.scaledM2high());
283 _h_B_T->fill(hemisphere.Bsum());
284 _h_B_W->fill(hemisphere.Bmax());
285 }
286 // jade jet rates
287 if (_s.count("y_2_JADE")) {
288 const FastJets& jadejet = apply<FastJets>(event, "JadeJets");
289 if (jadejet.clusterSeq()) {
290 const double y_23 = jadejet.clusterSeq()->exclusive_ymerge_max(2);
291 const double y_34 = jadejet.clusterSeq()->exclusive_ymerge_max(3);
292 const double y_45 = jadejet.clusterSeq()->exclusive_ymerge_max(4);
293 const double y_56 = jadejet.clusterSeq()->exclusive_ymerge_max(5);
294 for (auto& b : _s["y_2_JADE"]->bins()) {
295 const double ycut = std::stod(b.xEdge());
296 if (y_23 < ycut) _s["y_2_JADE"]->fill(b.xEdge());
297 }
298 for (auto& b : _s["y_3_JADE"]->bins()) {
299 const double ycut = std::stod(b.xEdge());
300 if (y_34 < ycut && y_23 > ycut) {
301 _s["y_3_JADE"]->fill(b.xEdge());
302 }
303 }
304 for (auto& b : _s["y_4_JADE"]->bins()) {
305 const double ycut = std::stod(b.xEdge());
306 if (y_45 < ycut && y_34 > ycut) {
307 _s["y_4_JADE"]->fill(b.xEdge());
308 }
309 }
310 for (auto& b : _s["y_5_JADE"]->bins()) {
311 const double ycut = std::stod(b.xEdge());
312 if (y_56 < ycut && y_45 > ycut) {
313 _s["y_5_JADE"]->fill(b.xEdge());
314 }
315 }
316 }
317 }
318 // Durham jet rates
319 if (_s.count("y_2_Durham")) {
320 const FastJets& durhamjet = apply<FastJets>(event, "DurhamJets");
321 if (durhamjet.clusterSeq()) {
322 const double y_23 = durhamjet.clusterSeq()->exclusive_ymerge_max(2);
323 const double y_34 = durhamjet.clusterSeq()->exclusive_ymerge_max(3);
324 const double y_45 = durhamjet.clusterSeq()->exclusive_ymerge_max(4);
325 const double y_56 = durhamjet.clusterSeq()->exclusive_ymerge_max(5);
326 for (auto& b : _s["y_2_Durham"]->bins()) {
327 const double ycut = std::stod(b.xEdge());
328 if (y_23 < ycut) _s["y_2_Durham"]->fill(b.xEdge());
329 }
330 for (auto& b : _s["y_3_Durham"]->bins()) {
331 const double ycut = std::stod(b.xEdge());
332 if (y_34 < ycut && y_23 > ycut) {
333 _s["y_3_Durham"]->fill(b.xEdge());
334 }
335 }
336 for (auto& b : _s["y_4_Durham"]->bins()) {
337 const double ycut = std::stod(b.xEdge());
338 if (y_45 < ycut && y_34 > ycut) {
339 _s["y_4_Durham"]->fill(b.xEdge());
340 }
341 }
342 for (auto& b : _s["y_5_Durham"]->bins()) {
343 const double ycut = std::stod(b.xEdge());
344 if (y_56 < ycut && y_45 > ycut) {
345 _s["y_5_Durham"]->fill(b.xEdge());
346 }
347 }
348 }
349 }
350 // Cambridge
351 if (_s.count("y_2_Cambridge")) {
352 PseudoJets pjs;
353 const FinalState& fs = apply<FinalState>(event, "FS");
354 for (size_t i = 0; i < fs.particles().size(); ++i) {
355 fastjet::PseudoJet pj = fs.particles()[i];
356 pjs.push_back(pj);
357 }
358 for (size_t i = 1; i < _s["y_2_Cambridge"]->numBins()+1; ++i) {
359 const auto& b = _s["y_2_Cambridge"]->bin(i);
360 const double ycut = std::stod(b.xEdge());
361 fastjet::EECambridgePlugin plugin(ycut);
362 fastjet::JetDefinition jdef(&plugin);
363 fastjet::ClusterSequence cseq(pjs, jdef);
364 size_t njet = cseq.inclusive_jets().size();
365 if (njet==2) {
366 _s["y_2_Cambridge"]->fill(b.xEdge());
367 }
368 else if (njet==3) {
369 if (i < _s["y_3_Cambridge"]->numBins()+1) {
370 _s["y_3_Cambridge"]->fill(b.xEdge());
371 }
372 }
373 else if(njet==4) {
374 if (i < _s["y_4_Cambridge"]->numBins()+1) {
375 _s["y_4_Cambridge"]->fill(b.xEdge());
376 }
377 }
378 else if(njet==5) {
379 if (i < _s["y_5_Cambridge"]->numBins()+1) {
380 _s["y_5_Cambridge"]->fill(b.xEdge());
381 }
382 }
383 }
384 }
385 }
386
387 Estimate1DPtr convertHisto(unsigned int ix,unsigned int iy, unsigned int iz, Histo1DPtr histo) {
388 Estimate1DPtr mult;
389 book(mult, ix, iy, iz);
390 barchart(histo, mult);
391 return mult;
392 }
393
394 /// Normalise histograms etc., after the run
395 void finalize() {
396 // Z pole plots
397 if(_h_Thrust_udsc) {
398 scale(_h_Thrust_udsc, 1/_sumW_udsc->sumW());
399 scale(_h_heavyJetmass_udsc, 1/_sumW_udsc->sumW());
400 scale(_h_totalJetbroad_udsc, 1/_sumW_udsc->sumW());
401 scale(_h_wideJetbroad_udsc, 1/_sumW_udsc->sumW());
402 scale(_h_Cparameter_udsc, 1/_sumW_udsc->sumW());
403 scale(_h_Dparameter_udsc, 1/_sumW_udsc->sumW());
404 scale(_h_Thrust_bottom, 1./_sumW_b->sumW());
405 scale(_h_heavyJetmass_bottom, 1./_sumW_b->sumW());
406 scale(_h_totalJetbroad_bottom, 1./_sumW_b->sumW());
407 scale(_h_wideJetbroad_bottom, 1./_sumW_b->sumW());
408 scale(_h_Cparameter_bottom, 1./_sumW_b->sumW());
409 scale(_h_Dparameter_bottom, 1./_sumW_b->sumW());
410 scale(_h_Ncharged, 1./_sumW_ch->sumW());
411 scale(_h_Ncharged_udsc, 1./_sumW_ch_udsc->sumW());
412 scale(_h_Ncharged_bottom, 1./_sumW_ch_b->sumW());
413 convertHisto(59,1,1,_h_Ncharged );
414 convertHisto(59,1,2,_h_Ncharged_udsc );
415 convertHisto(59,1,3,_h_Ncharged_bottom);
416
417 scale(_h_scaledMomentum, 1/_sumW_ch->sumW());
418 scale(_h_scaledMomentum_udsc, 1/_sumW_ch_udsc->sumW());
419 scale(_h_scaledMomentum_bottom, 1/_sumW_ch_b->sumW());
420 }
421 else {
422 if(_h_thrust) normalize(_h_thrust);
423 if(_h_rho) normalize(_h_rho);
424 if(_h_B_T) normalize(_h_B_T);
425 if(_h_B_W) normalize(_h_B_W);
426 if(_h_C) normalize(_h_C);
427 if(_h_D) normalize(_h_D);
428 if(_h_N) normalize(_h_N);
429 if(_h_xi) scale(_h_xi,1./sumOfWeights());
430
431
432 Estimate1DPtr mult;
433 if(_h_N) {
434 if(isCompatibleWithSqrtS(130.1*GeV)) {
435 convertHisto(60, 1, 1, _h_N);
436 }
437 else if(isCompatibleWithSqrtS(136.1*GeV)) {
438 convertHisto(60, 1, 2, _h_N);
439 }
440 else if(isCompatibleWithSqrtS(161.3*GeV)) {
441 convertHisto(60, 1, 3, _h_N);
442 }
443 else if(isCompatibleWithSqrtS(172.3*GeV)) {
444 convertHisto(61, 1, 1, _h_N);
445 }
446 else if(isCompatibleWithSqrtS(182.8*GeV)) {
447 convertHisto(61, 1, 2, _h_N);
448 }
449 else if(isCompatibleWithSqrtS(188.6*GeV)) {
450 convertHisto(61, 1, 3, _h_N);
451 }
452 else if(isCompatibleWithSqrtS(194.4*GeV)) {
453 convertHisto(62, 1, 1, _h_N);
454 }
455 else if(isCompatibleWithSqrtS(200.2*GeV)) {
456 convertHisto(62, 1, 2, _h_N);
457 }
458 else if(isCompatibleWithSqrtS(206.2*GeV)) {
459 convertHisto(62, 1, 3, _h_N);
460 }
461 }
462 // // the jets
463 scale(_s, 1./sumOfWeights());
464 }
465 }
466
467
468 /// Weight counters
469 CounterPtr _sumW_udsc, _sumW_b, _sumW_ch, _sumW_ch_udsc, _sumW_ch_b;
470
471 /// @name Histograms
472 /// @{
473 // at the Z pole
474 Histo1DPtr _h_Thrust_udsc, _h_Thrust_bottom;
475 Histo1DPtr _h_heavyJetmass_udsc, _h_heavyJetmass_bottom;
476 Histo1DPtr _h_totalJetbroad_udsc, _h_totalJetbroad_bottom;
477 Histo1DPtr _h_wideJetbroad_udsc, _h_wideJetbroad_bottom;
478 Histo1DPtr _h_Cparameter_udsc, _h_Cparameter_bottom;
479 Histo1DPtr _h_Dparameter_udsc, _h_Dparameter_bottom;
480 Histo1DPtr _h_Ncharged, _h_Ncharged_udsc, _h_Ncharged_bottom;
481 Histo1DPtr _h_scaledMomentum, _h_scaledMomentum_udsc, _h_scaledMomentum_bottom;
482 // at other enegies
483 Histo1DPtr _h_thrust, _h_rho, _h_B_T, _h_B_W, _h_C, _h_D, _h_N, _h_xi;
484 Histo1DPtr _h_y_2_JADE,_h_y_3_JADE,_h_y_4_JADE,_h_y_5_JADE;
485 map<string, BinnedHistoPtr<string>> _s;
486 /// @}
487
488 };
489
490
491 RIVET_DECLARE_PLUGIN(L3_2004_I652683);
492
493}
|