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