Project
Loading...
Searching...
No Matches
GeneratorPythia8.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
13
17#include <fairlogger/Logger.h>
18#include "TParticle.h"
19#include "TF1.h"
20#include "TRandom.h"
24#if PYTHIA_VERSION_INTEGER >= 8310
25#include "Pythia8/HIInfo.h"
26#else
27#include "Pythia8/HIUserHooks.h"
28#endif
29#include "Pythia8Plugins/PowhegHooks.h"
30#include "TSystem.h"
33#include <filesystem>
35
36#include <iostream>
37#include <unordered_map>
38#include <numeric>
39
40namespace o2
41{
42namespace eventgen
43{
44
46
47/*****************************************************************/
48/*****************************************************************/
49
50// the default construct uses the GeneratorPythia8Param singleton to extract a config and delegates
51// to the proper constructor
53{
54 LOG(info) << "GeneratorPythia8 constructed from GeneratorPythia8Param ConfigurableParam";
55}
56
57/*****************************************************************/
58
59GeneratorPythia8::GeneratorPythia8(Pythia8GenConfig const& config) : Generator("ALICEo2", "ALICEo2 Pythia8 Generator")
60{
64
65 mInterface = reinterpret_cast<void*>(&mPythia);
66 mInterfaceName = "pythia8";
67
68 LOG(info) << "Instance \'Pythia8\' generator with following parameters";
69 LOG(info) << "config: " << config.config;
70 LOG(info) << "hooksFileName: " << config.hooksFileName;
71 LOG(info) << "hooksFuncName: " << config.hooksFuncName;
72
73 mGenConfig = config;
74
78}
79
80/*****************************************************************/
81
82GeneratorPythia8::GeneratorPythia8(const Char_t* name, const Char_t* title) : Generator(name, title)
83{
86 mInterface = reinterpret_cast<void*>(&mPythia);
87 mInterfaceName = "pythia8";
88}
89
91{
92 // check first of all if Init not yet called and seed not <0
93 if (mIsInitialized) {
94 return false;
95 }
96 if (seed < 0) {
97 return false;
98 }
99 // sets the initial seed and applies the correct Pythia
100 // range
101 mInitialRNGSeed = seed % (MAX_SEED + 1);
102 LOG(info) << "GeneratorPythia8: Setting initial seed to " << mInitialRNGSeed;
103 return true;
104}
105
106/*****************************************************************/
108{
113
116
117 auto seed = mInitialRNGSeed;
118 if (seed == -1) {
119 // Will use the mInitialRNGSeed if it was set.
120 // Otherwise will seed the generator with the state of
121 // TRandom::GetSeed. This is the seed that is influenced from
122 // SimConfig --seed command line options options.
123 seed = gRandom->TRandom::GetSeed(); // this uses the "original" seed
124 // we advance the seed by one so that the next Pythia8 generator gets a different value
125 if (mThisPythia8InstanceID > 0) {
126 gRandom->Rndm();
127 LOG(info) << "Multiple Pythia8 generator instances detected .. automatically adjusting seed further to avoid overlap ";
128 seed = seed ^ gRandom->GetSeed(); // this uses the "current" seed
129 }
130 // apply max seed cuttof
131 seed = seed % (MAX_SEED + 1);
132 LOG(info) << "GeneratorPythia8: Using random seed from gRandom % 900000001: " << seed;
133 }
134 mPythia.readString("Random:setSeed on");
135 mPythia.readString("Random:seed " + std::to_string(seed));
136}
137
138/*****************************************************************/
139
141{
146
152
154 if (!mConfig.empty()) {
155 std::stringstream ss(mConfig);
156 std::string config;
157 while (getline(ss, config, ' ')) {
158 config = gSystem->ExpandPathName(config.c_str());
159 LOG(info) << "Reading configuration from file: " << config;
160 if (!mPythia.readFile(config, true)) {
161 LOG(fatal) << "Failed to init \'Pythia8\': problems with configuration file "
162 << config;
163 return false;
164 }
165 }
166 }
167
169 if (!mHooksFileName.empty()) {
170 LOG(info) << "Applying \'Pythia8\' user hooks: " << mHooksFileName << " -> " << mHooksFuncName;
171 auto hooks = o2::conf::GetFromMacro<Pythia8::UserHooks*>(mHooksFileName, mHooksFuncName, "Pythia8::UserHooks*", "pythia8_user_hooks");
172 if (!hooks) {
173 LOG(fatal) << "Failed to init \'Pythia8\': problem with user hooks configuration ";
174 return false;
175 }
176 setUserHooks(hooks);
177 }
178
179#if PYTHIA_VERSION_INTEGER < 8300
188 mPythia.readString("HadronLevel:Decay off");
189#endif
190 if (mPythia.settings.mode("Beams:frameType") == 4) {
191 // Hook for POWHEG
192 // Read in key POWHEG merging settings
193 int vetoMode = mPythia.settings.mode("POWHEG:veto");
194 int MPIvetoMode = mPythia.settings.mode("POWHEG:MPIveto");
195 bool loadHooks = (vetoMode > 0 || MPIvetoMode > 0);
196 // Add in user hooks for shower vetoing
197 std::shared_ptr<Pythia8::PowhegHooks> powhegHooks;
198 if (loadHooks) {
199 // Set ISR and FSR to start at the kinematical limit
200 if (vetoMode > 0) {
201 mPythia.readString("SpaceShower:pTmaxMatch = 2");
202 mPythia.readString("TimeShower:pTmaxMatch = 2");
203 }
204 // Set MPI to start at the kinematical limit
205 if (MPIvetoMode > 0) {
206 mPythia.readString("MultipartonInteractions:pTmaxMatch = 2");
207 }
208 powhegHooks = std::make_shared<Pythia8::PowhegHooks>();
209 mPythia.setUserHooksPtr((Pythia8::UserHooksPtr)powhegHooks);
210 }
211 }
213 mPythia.particleData.addParticle(1000100200, "20Ne", 6, 30, 0, 19.992440);
215 if (!mPythia.init()) {
216 LOG(fatal) << "Failed to init \'Pythia8\': init returned with error";
217 return false;
218 }
219
221
222 mIsInitialized = true;
223
225 return true;
226}
227
228/*****************************************************************/
229void GeneratorPythia8::setUserHooks(Pythia8::UserHooks* hooks)
230{
231#if PYTHIA_VERSION_INTEGER < 8300
232 mPythia.setUserHooksPtr(hooks);
233#else
234 mPythia.setUserHooksPtr(std::shared_ptr<Pythia8::UserHooks>(hooks));
235#endif
236}
237
238/*****************************************************************/
239
240Bool_t
242{
244 if (!mPythia.next()) {
245 return false;
246 }
247
248#if PYTHIA_VERSION_INTEGER < 8300
263 auto nParticles = mPythia.event.size();
264 for (int iparticle = 0; iparticle < nParticles; iparticle++) {
265 auto& aParticle = mPythia.event[iparticle];
266 aParticle.xProd(0.);
267 aParticle.yProd(0.);
268 aParticle.zProd(0.);
269 aParticle.tProd(0.);
270 }
271
273 if (!mPythia.moreDecays())
274 return false;
275#endif
276
278 return true;
279}
280
281/*****************************************************************/
283 const std::vector<int>& old2New,
284 size_t index,
285 std::vector<bool>& done,
286 GetRelatives getter,
287 SetRelatives setter,
288 FirstLastRelative firstLast,
289 const std::string& what,
290 const std::string& ind)
291{
292 // Utility to find new index, or -1 if not found
293 auto findNew = [old2New](size_t old) -> int {
294 return old2New[old];
295 };
296 int newIdx = findNew(index);
297 int hepmc = event[index].statusHepMC();
298
299 LOG(debug) << ind
300 << index << " -> "
301 << newIdx << " (" << hepmc << ") ";
302 if (done[index]) {
303 LOG(debug) << ind << " already done";
304 return;
305 }
306
307 // Our list of new relatives
308 using IdList = std::pair<int, int>;
309 constexpr int invalid = 0xFFFFFFF;
310 IdList newRelatives = std::make_pair(invalid, -invalid);
311
312 // Utility to add id
313 auto addId = [](IdList& l, size_t id) {
314 l.first = std::min(int(id), l.first);
315 l.second = std::max(int(id), l.second);
316 };
317
318 // Get particle and relatives
319 auto& particle = event[index];
320 auto relatives = getter(particle);
321
322 LOG(debug) << ind << " Check " << what << "s ["
323 << std::setw(3) << firstLast(particle).first << ","
324 << std::setw(3) << firstLast(particle).second << "] "
325 << relatives.size();
326
327 for (auto relativeIdx : relatives) {
328 int newRelative = findNew(relativeIdx);
329 if (newRelative >= 0) {
330 // If this relative is to be kept, then append to list of new
331 // relatives.
332 LOG(debug) << ind << " "
333 << what << " "
334 << relativeIdx << " -> "
335 << newRelative << " to be kept" << std::endl;
336 addId(newRelatives, newRelative);
337 continue;
338 }
339 LOG(debug) << ind << " "
340 << what << " "
341 << relativeIdx << " not to be kept "
342 << (done[relativeIdx] ? "already done" : "to be done")
343 << std::endl;
344
345 // Below is code for when the relative is not to be kept
346 auto& relative = event[relativeIdx];
347 if (not done[relativeIdx]) {
348 // IF the relative hasn't been processed yet, do so now
350 old2New, // Map from old to new
351 relativeIdx, // current particle index
352 done, // cache flag
353 getter, // get mother relatives
354 setter, // set mother relatives
355 firstLast, // get first and last
356 what, // what we're looking at
357 ind + " "); // Logging indent
358 }
359
360 // If this relative was already done, then get its relatives and
361 // add them to the list of new relatives.
362 auto grandRelatives = firstLast(relative);
363 int grandRelative1 = grandRelatives.first;
364 int grandRelative2 = grandRelatives.second;
365 assert(grandRelative1 != invalid);
366 assert(grandRelative2 != -invalid);
367 if (grandRelative1 > 0) {
368 addId(newRelatives, grandRelative1);
369 }
370 if (grandRelative2 > 0) {
371 addId(newRelatives, grandRelative2);
372 }
373 LOG(debug) << ind << " "
374 << what << " "
375 << relativeIdx << " gave new relatives "
376 << grandRelative1 << " -> " << grandRelative2;
377 }
378 LOG(debug) << ind << " Got "
379 << (newRelatives.second - newRelatives.first + 1) << " new "
380 << what << "s ";
381
382 if (newRelatives.first != invalid) {
383 // If the first relative is not invalid, then the second isn't
384 // either (possibly the same though).
385 int newRelative1 = newRelatives.first;
386 int newRelative2 = newRelatives.second;
387 setter(particle, newRelative1, newRelative2);
388 LOG(debug) << ind << " " << what << "s: "
389 << firstLast(particle).first << " ("
390 << newRelative1 << "),"
391 << firstLast(particle).second << " ("
392 << newRelative2 << ")";
393
394 } else {
395 setter(particle, 0, 0);
396 }
397 done[index] = true;
398}
399
400/*****************************************************************/
402{
403 // Mapping from old to new index.
404 std::vector<int> old2new(event.size(), -1);
405
406 // Particle 0 is a system particle, and we will skip that in the
407 // following.
408 size_t newId = 0;
409
410 // Loop over particles and store those we need
411 for (size_t i = 1; i < event.size(); i++) {
412 auto& particle = event[i];
413 if (select(particle)) {
414 ++newId;
415 old2new[i] = newId;
416 }
417 }
418 // Utility to find new index, or -1 if not found
419 auto findNew = [old2new](size_t old) -> int {
420 return old2new[old];
421 };
422
423 // First loop, investigate mothers - from the bottom
424 auto getMothers = [](const Pythia8::Particle& particle) { return particle.motherList(); };
425 auto setMothers = [](Pythia8::Particle& particle, int m1, int m2) { particle.mothers(m1, m2); };
426 auto firstLastMothers = [](const Pythia8::Particle& particle) { return std::make_pair(particle.mother1(), particle.mother2()); };
427
428 std::vector<bool> motherDone(event.size(), false);
429 for (size_t i = 1; i < event.size(); ++i) {
431 old2new, // Map from old to new
432 i, // current particle index
433 motherDone, // cache flag
434 getMothers, // get mother relatives
435 setMothers, // set mother relatives
436 firstLastMothers, // get first and last
437 "mother"); // what we're looking at
438 }
439
440 // Second loop, investigate daughters - from the top
441 auto getDaughters = [](const Pythia8::Particle& particle) {
442 // In case of |status|==13 (diffractive), we cannot use
443 // Pythia8::Particle::daughterList as it will give more than
444 // just the immediate daughters. In that cae, we do it
445 // ourselves.
446 if (std::abs(particle.status()) == 13) {
447 int d1 = particle.daughter1();
448 int d2 = particle.daughter2();
449 if (d1 == 0 and d2 == 0) {
450 return std::vector<int>();
451 }
452 if (d2 == 0) {
453 return std::vector<int>{d1};
454 }
455 if (d2 > d1) {
456 std::vector<int> ret(d2-d1+1);
457 std::iota(ret.begin(), ret.end(), d1);
458 return ret;
459 }
460 return std::vector<int>{d2,d1};
461 }
462 return particle.daughterList(); };
463 auto setDaughters = [](Pythia8::Particle& particle, int d1, int d2) { particle.daughters(d1, d2); };
464 auto firstLastDaughters = [](const Pythia8::Particle& particle) { return std::make_pair(particle.daughter1(), particle.daughter2()); };
465
466 std::vector<bool> daughterDone(event.size(), false);
467 for (size_t i = event.size() - 1; i > 0; --i) {
469 old2new, // Map from old to new
470 i, // current particle index
471 daughterDone, // cache flag
472 getDaughters, // get mother relatives
473 setDaughters, // set mother relatives
474 firstLastDaughters, // get first and last
475 "daughter"); // what we're looking at
476 }
477
478 // Make a pruned event
479 Pythia8::Event pruned;
480 pruned.init("Pruned event", &mPythia.particleData);
481 pruned.reset();
482
483 for (size_t i = 1; i < event.size(); i++) {
484 int newIdx = findNew(i);
485 if (newIdx < 0) {
486 continue;
487 }
488
489 auto particle = event[i];
490 int realIdx = pruned.append(particle);
491 assert(realIdx == newIdx);
492 }
493
494 // We may have that two or more mothers share some daughters, but
495 // that one or more mothers have more daughters than the other
496 // mothers, and hence not all daughters point back to all mothers.
497 // This can happen, for example, if a beam particle radiates
498 // on-shell particles before an interaction with any daughters
499 // from the other mothers. Thus, we need to take care of that or
500 // the event record will be corrupted.
501 //
502 // What we do is that for all particles, we look up the daughters.
503 // Then for each daughter, we check the mothers of those
504 // daughters. If this list of mothers include other mothers than
505 // the currently investigated mother, we must change the mothers
506 // of the currently investigated daughters.
507 using IdList = std::pair<int, int>;
508 // Utility to add id
509 auto addId = [](IdList& l, size_t id) {
510 l.first = std::min(int(id), l.first);
511 l.second = std::max(int(id), l.second);
512 };
513 constexpr int invalid = 0xFFFFFFF;
514
515 std::vector<bool> shareDone(pruned.size(), false);
516 for (size_t i = 1; i < pruned.size(); i++) {
517 if (shareDone[i]) {
518 continue;
519 }
520
521 auto& particle = pruned[i];
522 auto daughters = particle.daughterList();
523 IdList allDaughters = std::make_pair(invalid, -invalid);
524 IdList allMothers = std::make_pair(invalid, -invalid);
525 addId(allMothers, i);
526 for (auto daughterIdx : daughters) {
527 // Add this daughter to set of all daughters
528 addId(allDaughters, daughterIdx);
529 auto& daughter = pruned[daughterIdx];
530 auto otherMothers = daughter.motherList();
531 for (auto otherMotherIdx : otherMothers) {
532 // Add this mother to set of all mothers. That is, take all
533 // mothers of the current daughter of the current particle
534 // and store that. In this way, we register mothers that
535 // share a daughter with the current particle.
536 addId(allMothers, otherMotherIdx);
537 // We also need to take all the daughters of this shared
538 // mother and reister those.
539 auto& otherMother = pruned[otherMotherIdx];
540 int otherDaughter1 = otherMother.daughter1();
541 int otherDaughter2 = otherMother.daughter2();
542 if (otherDaughter1 > 0) {
543 addId(allDaughters, otherDaughter1);
544 }
545 if (otherDaughter2 > 0) {
546 addId(allDaughters, otherDaughter2);
547 }
548 }
549 // At this point, we have added all mothers of current
550 // daughter, and all daughters of those mothers.
551 }
552 // At this point, we have all mothers that share daughters with
553 // the current particle, and we have all of the daughters
554 // too.
555 //
556 // We can now update the daughter information on all mothers
557 int minDaughter = allDaughters.first;
558 int maxDaughter = allDaughters.second;
559 int minMother = allMothers.first;
560 int maxMother = allMothers.second;
561 if (minMother != invalid) {
562 // If first mother isn't invalid, then second isn't either
563 for (size_t motherIdx = minMother; motherIdx <= maxMother; //
564 motherIdx++) {
565 shareDone[motherIdx] = true;
566 if (minDaughter == invalid) {
567 pruned[motherIdx].daughters(0, 0);
568 } else {
569 pruned[motherIdx].daughters(minDaughter, maxDaughter);
570 }
571 }
572 }
573 if (minDaughter != invalid) {
574 // If least mother isn't invalid, then largest mother will not
575 // be invalid either.
576 for (size_t daughterIdx = minDaughter; daughterIdx <= maxDaughter; //
577 daughterIdx++) {
578 if (minMother == invalid) {
579 pruned[daughterIdx].mothers(0, 0);
580 } else {
581 pruned[daughterIdx].mothers(minMother, maxMother);
582 }
583 }
584 }
585 }
586 if (mGenConfig.verbose) {
587 LOG(info) << "Pythia event was pruned from " << event.size()
588 << " to " << pruned.size() << " particles";
589 }
590 // Assign our pruned event to the event passed in
591 event = pruned;
592}
593
594/*****************************************************************/
596{
597 mUserFilterFcn = [](Pythia8::Particle const&) -> bool { return true; };
598
599 std::string filter = mGenConfig.particleFilter;
600 if (filter.size() > 0) {
601 LOG(info) << "Initializing the callback for user-based particle pruning " << filter;
602 auto expandedFileName = o2::utils::expandShellVarsInFileName(filter);
603 if (std::filesystem::exists(expandedFileName)) {
604 // if the filter is in a file we will compile the hook on the fly
605 mUserFilterFcn = o2::conf::GetFromMacro<UserFilterFcn>(expandedFileName, "filterPythia()", "o2::eventgen::GeneratorPythia8::UserFilterFcn", "o2mc_pythia8_userfilter_hook");
606 LOG(info) << "Hook initialized from file " << expandedFileName;
607 } else {
608 // if it's not a file we interpret it as a C++ lambda string and JIT it directly;
609 LOG(error) << "Did not find a file " << expandedFileName << " ; Will not execute hook";
610 }
611 mApplyPruning = true;
612 }
613}
614
615/*****************************************************************/
616
617Bool_t
619{
622 // The right moment to filter out unwanted stuff (like parton-level
623 // event information) Here, we aim to filter out everything before
624 // hadronization with the motivation to reduce the size of the MC
625 // event record in the AOD.
626
627 std::function<bool(const Pythia8::Particle&)> partonSelect = [](const Pythia8::Particle&) { return true; };
628 bool includeParton = mGenConfig.includePartonEvent;
629 if (not includeParton) {
630
631 // Select pythia particles
632 partonSelect = [](const Pythia8::Particle& particle) {
633 switch (particle.statusHepMC()) {
634 case 1: // Final st
635 case 2: // Decayed
636 case 4: // Beam
637 return true;
638 }
639 // For example to keep diffractive particles
640 // if (particle.id() == 9902210) return true;
641 return false;
642 };
643 mApplyPruning = true;
644 }
645
646 if (mApplyPruning) {
647 auto finalSelect = [partonSelect, this](const Pythia8::Particle& p) { return partonSelect(p) && mUserFilterFcn(p); };
648 pruneEvent(event, finalSelect);
649 }
650
651 /* loop over particles */
652 auto nParticles = event.size();
653 for (Int_t iparticle = 1; iparticle < nParticles; iparticle++) {
654 // first particle is system
655 auto particle = event[iparticle];
656 auto pdg = particle.id();
657 auto st = o2::mcgenstatus::MCGenStatusEncoding(particle.statusHepMC(), //
658 particle.status()) //
660 mParticles.push_back(TParticle(particle.id(), // Particle type
661 st, // status
662 particle.mother1() - 1, // first mother
663 particle.mother2() - 1, // second mother
664 particle.daughter1() - 1, // first daughter
665 particle.daughter2() - 1, // second daughter
666 particle.px(), // X-momentum
667 particle.py(), // Y-momentum
668 particle.pz(), // Z-momentum
669 particle.e(), // Energy
670 particle.xProd(), // Production X
671 particle.yProd(), // Production Y
672 particle.zProd(), // Production Z
673 particle.tProd())); // Production t
674 mParticles.back().SetBit(ParticleStatus::kToBeDone, //
675 particle.statusHepMC() == 1);
676 }
677
679 return kTRUE;
680}
681
682/*****************************************************************/
683
685{
688
689 eventHeader->putInfo<std::string>(Key::generator, "pythia8");
690 eventHeader->putInfo<int>(Key::generatorVersion, PYTHIA_VERSION_INTEGER);
691 eventHeader->putInfo<std::string>(Key::processName, mPythia.info.name());
692 eventHeader->putInfo<int>(Key::processCode, mPythia.info.code());
693 eventHeader->putInfo<float>(Key::weight, mPythia.info.weight());
694
695 auto& info = mPythia.info;
696
697 eventHeader->putInfo<int>(Key::acceptedEvents, info.nAccepted());
698 eventHeader->putInfo<int>(Key::attemptedEvents, info.nTried());
699
700 // Set PDF information
701 eventHeader->putInfo<int>(Key::pdfParton1Id, info.id1pdf());
702 eventHeader->putInfo<int>(Key::pdfParton2Id, info.id2pdf());
703 eventHeader->putInfo<float>(Key::pdfX1, info.x1pdf());
704 eventHeader->putInfo<float>(Key::pdfX2, info.x2pdf());
705 eventHeader->putInfo<float>(Key::pdfScale, info.QFac());
706 eventHeader->putInfo<float>(Key::pdfXF1, info.pdf1());
707 eventHeader->putInfo<float>(Key::pdfXF2, info.pdf2());
708
709 // Set cross section
710 eventHeader->putInfo<float>(Key::xSection, info.sigmaGen() * 1e9);
711 eventHeader->putInfo<float>(Key::xSectionError, info.sigmaErr() * 1e9);
712
713 // Set event scale and nMPI
714 eventHeader->putInfo<float>(Key::eventScale, info.QRen());
715 eventHeader->putInfo<int>(Key::mpi, info.nMPI());
716
717 // Set weights (overrides cross-section for each weight)
718 size_t iw = 0;
719 auto xsecErr = info.weightContainerPtr->getTotalXsecErr();
720 for (auto w : info.weightContainerPtr->getTotalXsec()) {
721 std::string post = (iw == 0 ? "" : "_" + std::to_string(iw));
722 eventHeader->putInfo<float>(Key::weight + post, info.weightValueByIndex(iw));
723 eventHeader->putInfo<float>(Key::xSection + post, w * 1e9);
724 eventHeader->putInfo<float>(Key::xSectionError + post, xsecErr[iw] * 1e9);
725 iw++;
726 }
727
728#if PYTHIA_VERSION_INTEGER < 8300
729 auto hiinfo = mPythia.info.hiinfo;
730#else
731 auto hiinfo = mPythia.info.hiInfo;
732#endif
733
734 if (hiinfo) {
736 eventHeader->SetB(hiinfo->b());
737 eventHeader->putInfo<double>(Key::impactParameter, hiinfo->b());
739#if PYTHIA_VERSION_INTEGER >= 8310
740 eventHeader->putInfo<double>(Key::planeAngle, hiinfo->phi());
741#endif
742 auto bImp = hiinfo->b();
744 int nColl, nPart;
745 int nPartProtonProj, nPartNeutronProj, nPartProtonTarg, nPartNeutronTarg;
746 int nRemnProtonProj, nRemnNeutronProj, nRemnProtonTarg, nRemnNeutronTarg;
747 int nFreeNeutronProj, nFreeProtonProj, nFreeNeutronTarg, nFreeProtonTarg;
748 getNcoll(nColl);
749 getNpart(nPart);
750 getNpart(nPartProtonProj, nPartNeutronProj, nPartProtonTarg, nPartNeutronTarg);
751 getNremn(nRemnProtonProj, nRemnNeutronProj, nRemnProtonTarg, nRemnNeutronTarg);
752 getNfreeSpec(nFreeNeutronProj, nFreeProtonProj, nFreeNeutronTarg, nFreeProtonTarg);
753 eventHeader->putInfo<int>(Key::nColl, nColl);
754 // These are all non-HepMC3 fields - of limited use
755 eventHeader->putInfo<int>("Npart", nPart);
756 eventHeader->putInfo<int>("Npart_proj_p", nPartProtonProj);
757 eventHeader->putInfo<int>("Npart_proj_n", nPartNeutronProj);
758 eventHeader->putInfo<int>("Npart_targ_p", nPartProtonTarg);
759 eventHeader->putInfo<int>("Npart_targ_n", nPartNeutronTarg);
760 eventHeader->putInfo<int>("Nremn_proj_p", nRemnProtonProj);
761 eventHeader->putInfo<int>("Nremn_proj_n", nRemnNeutronProj);
762 eventHeader->putInfo<int>("Nremn_targ_p", nRemnProtonTarg);
763 eventHeader->putInfo<int>("Nremn_targ_n", nRemnNeutronTarg);
764 eventHeader->putInfo<int>("Nfree_proj_n", nFreeNeutronProj);
765 eventHeader->putInfo<int>("Nfree_proj_p", nFreeProtonProj);
766 eventHeader->putInfo<int>("Nfree_targ_n", nFreeNeutronTarg);
767 eventHeader->putInfo<int>("Nfree_targ_p", nFreeProtonTarg);
768
769 // --- HepMC3 conforming information ---
770 // This is how the Pythia authors define Ncoll
771 // eventHeader->putInfo<int>(Key::nColl,
772 // hiinfo->nAbsProj() + hiinfo->nDiffProj() +
773 // hiinfo->nAbsTarg() + hiinfo->nDiffTarg() -
774 // hiiinfo->nCollND() - hiinfo->nCollDD());
775 eventHeader->putInfo<int>(Key::nPartProjectile,
776 hiinfo->nAbsProj() + hiinfo->nDiffProj());
777 eventHeader->putInfo<int>(Key::nPartTarget,
778 hiinfo->nAbsTarg() + hiinfo->nDiffTarg());
779 eventHeader->putInfo<int>(Key::nCollHard, hiinfo->nCollNDTot());
780 }
781}
782
783/*****************************************************************/
784
785void GeneratorPythia8::selectFromAncestor(int ancestor, Pythia8::Event& inputEvent, Pythia8::Event& outputEvent)
786{
787
792 // recursive selection via lambda function
793 std::set<int> selected;
794 std::function<void(int)> select;
795 select = [&](int i) {
796 selected.insert(i);
797 auto dl = inputEvent[i].daughterList();
798 for (auto j : dl) {
799 select(j);
800 }
801 };
802 select(ancestor);
803
804 // map selected particle index to output index
805 std::map<int, int> indexMap;
806 int index = outputEvent.size();
807 for (auto i : selected) {
808 indexMap[i] = index++;
809 }
810
811 // adjust mother/daughter indices and append to output event
812 for (auto i : selected) {
813 auto p = mPythia.event[i];
814 auto m1 = indexMap[p.mother1()];
815 auto m2 = indexMap[p.mother2()];
816 auto d1 = indexMap[p.daughter1()];
817 auto d2 = indexMap[p.daughter2()];
818 p.mothers(m1, m2);
819 p.daughters(d1, d2);
820
821 outputEvent.append(p);
822 }
823}
824
825/*****************************************************************/
826
827void GeneratorPythia8::getNcoll(const Pythia8::Info& info, int& nColl)
828{
829
832#if PYTHIA_VERSION_INTEGER < 8300
833 auto hiinfo = info.hiinfo;
834#else
835 auto hiinfo = info.hiInfo;
836#endif
837 nColl = 0;
838 if (!hiinfo) {
839 LOG(warn) << "No heavy-ion information from Pythia";
840 return;
841 }
842
843 // This is how the Pythia authors define Ncoll
844 nColl = (hiinfo->nAbsProj() + hiinfo->nDiffProj() +
845 hiinfo->nAbsTarg() + hiinfo->nDiffTarg() -
846 hiinfo->nCollND() - hiinfo->nCollDD());
847
848 if (not hiinfo->subCollisionsPtr()) {
849#if PYTHIA_VERSION_INTEGER < 8310
850 LOG(fatal) << "No sub-collision pointer from Pythia";
851#endif
852 return;
853 }
854
855 // loop over sub-collisions
856 auto scptr = hiinfo->subCollisionsPtr();
857 for (auto sc : *scptr) {
858
859 // wounded nucleon flag in projectile/target
860 auto pW = sc.proj->status() == Pythia8::Nucleon::ABS; // according to C.Bierlich this should be == Nucleon::ABS
861 auto tW = sc.targ->status() == Pythia8::Nucleon::ABS;
862
863 // increase number of collisions if both are wounded
864 if (pW && tW) {
865 nColl++;
866 }
867 }
868}
869
870/*****************************************************************/
871
872void GeneratorPythia8::getNpart(const Pythia8::Info& info, int& nPart)
873{
874
877#if PYTHIA_VERSION_INTEGER < 8300
878 auto hiinfo = info.hiinfo;
879#else
880 auto hiinfo = info.hiInfo;
881#endif
882 nPart = 0;
883 if (not hiinfo) {
884 return;
885 }
886
887 // This is how the Pythia authors calculate Npart
888 nPart = (hiinfo->nAbsProj() + hiinfo->nDiffProj() +
889 hiinfo->nAbsTarg() + hiinfo->nDiffTarg());
890 if (not hiinfo->subCollisionsPtr()) {
891 return;
892 }
893
894 int nProtonProj, nNeutronProj, nProtonTarg, nNeutronTarg;
895 getNpart(info, nProtonProj, nNeutronProj, nProtonTarg, nNeutronTarg);
896 nPart = nProtonProj + nNeutronProj + nProtonTarg + nNeutronTarg;
897}
898
899/*****************************************************************/
900
901void GeneratorPythia8::getNpart(const Pythia8::Info& info, int& nProtonProj, int& nNeutronProj, int& nProtonTarg, int& nNeutronTarg)
902{
903
906#if PYTHIA_VERSION_INTEGER < 8300
907 auto hiinfo = info.hiinfo;
908#else
909 auto hiinfo = info.hiInfo;
910#endif
911
912 nProtonProj = nNeutronProj = nProtonTarg = nNeutronTarg = 0;
913 if (!hiinfo) {
914 return;
915 }
916
917 nProtonProj = hiinfo->nAbsProj() + hiinfo->nDiffProj();
918 nProtonTarg = hiinfo->nAbsTarg() + hiinfo->nDiffTarg();
919 if (not hiinfo->subCollisionsPtr()) {
920#if PYTHIA_VERSION_INTEGER < 8310
921 LOG(fatal) << "No sub-collision pointer from Pythia";
922#endif
923 return;
924 }
925
926 // keep track of wounded nucleons
927 std::vector<Pythia8::Nucleon*> projW;
928 std::vector<Pythia8::Nucleon*> targW;
929
930 // loop over sub-collisions
931 auto scptr = hiinfo->subCollisionsPtr();
932 for (auto sc : *scptr) {
933
934 // wounded nucleon flag in projectile/target
935 auto pW = sc.proj->status() == Pythia8::Nucleon::ABS || sc.proj->status() == Pythia8::Nucleon::DIFF; // according to C.Bierlich this should be == Nucleon::ABS || Nucleon::DIFF
936 auto tW = sc.targ->status() == Pythia8::Nucleon::ABS || sc.targ->status() == Pythia8::Nucleon::DIFF;
937
938 // increase number of wounded projectile nucleons if not yet in the wounded vector
939 if (pW && std::find(projW.begin(), projW.end(), sc.proj) == projW.end()) {
940 projW.push_back(sc.proj);
941 if (sc.proj->id() == 2212) {
942 nProtonProj++;
943 } else if (sc.proj->id() == 2112) {
944 nNeutronProj++;
945 }
946 }
947
948 // increase number of wounded target nucleons if not yet in the wounded vector
949 if (tW && std::find(targW.begin(), targW.end(), sc.targ) == targW.end()) {
950 targW.push_back(sc.targ);
951 if (sc.targ->id() == 2212) {
952 nProtonTarg++;
953 } else if (sc.targ->id() == 2112) {
954 nNeutronTarg++;
955 }
956 }
957 }
958}
959
960/*****************************************************************/
961
962void GeneratorPythia8::getNremn(const Pythia8::Event& event, int& nProtonProj, int& nNeutronProj, int& nProtonTarg, int& nNeutronTarg)
963{
964
967 // reset
968 nProtonProj = nNeutronProj = nProtonTarg = nNeutronTarg = 0;
969 auto nNucRem = 0;
970
971 // particle loop
972 auto nparticles = event.size();
973 for (int ipa = 0; ipa < nparticles; ++ipa) {
974 const auto particle = event[ipa];
975 auto pdg = particle.id();
976
977 // nuclear remnants have pdg code = ±10LZZZAAA9
978 if (pdg < 1000000000) {
979 continue; // must be nucleus
980 }
981 if (pdg % 10 != 9) {
982 continue; // first digit must be 9
983 }
984 nNucRem++;
985
986 // extract A, Z and L from pdg code
987 pdg /= 10;
988 auto A = pdg % 1000;
989 pdg /= 1000;
990 auto Z = pdg % 1000;
991 pdg /= 1000;
992 auto L = pdg % 10;
993
994 if (particle.pz() > 0.) {
995 nProtonProj = Z;
996 nNeutronProj = A - Z;
997 }
998 if (particle.pz() < 0.) {
999 nProtonTarg = Z;
1000 nNeutronTarg = A - Z;
1001 }
1002
1003 } // end of particle loop
1004
1005 if (nNucRem > 2) {
1006 LOG(warning) << " GeneratorPythia8: found more than two nuclear remnants (weird)";
1007 }
1008}
1009/*****************************************************************/
1010
1011/*****************************************************************/
1012
1013void GeneratorPythia8::getNfreeSpec(const Pythia8::Info& info, int& nFreenProj, int& nFreepProj, int& nFreenTarg, int& nFreepTarg)
1014{
1017#if PYTHIA_VERSION_INTEGER < 8300
1018 auto hiinfo = info.hiinfo;
1019#else
1020 auto hiinfo = info.hiInfo;
1021#endif
1022
1023 if (!hiinfo) {
1024 return;
1025 }
1026
1027 double b = hiinfo->b();
1028
1029 static o2::zdc::FragmentParam frag; // data-driven model to get free spectators given impact parameter
1030
1031 TF1 const& fneutrons = frag.getfNeutrons();
1032 TF1 const& fsigman = frag.getsigmaNeutrons();
1033 TF1 const& fprotons = frag.getfProtons();
1034 TF1 const& fsigmap = frag.getsigmaProtons();
1035
1036 // Calculating no. of free spectators from parametrization
1037 int nneu[2] = {0, 0};
1038 for (int i = 0; i < 2; i++) {
1039 float nave = fneutrons.Eval(b);
1040 float sigman = fsigman.Eval(b);
1041 float nfree = gRandom->Gaus(nave, 0.68 * sigman * nave);
1042 nneu[i] = (int)nfree;
1043 if (nave < 0 || nneu[i] < 0) {
1044 nneu[i] = 0;
1045 }
1046 if (nneu[i] > 126) {
1047 nneu[i] = 126;
1048 }
1049 }
1050 //
1051 int npro[2] = {0, 0};
1052 for (int i = 0; i < 2; i++) {
1053 float pave = fprotons.Eval(b);
1054 float sigmap = fsigman.Eval(b);
1055 float pfree = gRandom->Gaus(pave, 0.68 * sigmap * pave) / 0.7;
1056 npro[i] = (int)pfree;
1057 if (pave < 0 || npro[i] < 0) {
1058 npro[i] = 0;
1059 }
1060 if (npro[i] > 82) {
1061 npro[i] = 82;
1062 }
1063 }
1064
1065 nFreenProj = nneu[0];
1066 nFreenTarg = nneu[1];
1067 nFreepProj = npro[0];
1068 nFreepTarg = npro[1];
1069 /*****************************************************************/
1070}
1071
1072} /* namespace eventgen */
1073} /* namespace o2 */
o2::monitoring::tags::Key Key
int32_t i
bool done
@ kToBeDone
uint32_t j
Definition RawData.h:0
std::ostringstream debug
benchmark::State & st
Definition A.h:16
void putInfo(std::string const &key, T const &value)
void setHooksFuncName(std::string val)
void selectFromAncestor(int ancestor, Pythia8::Event &inputEvent, Pythia8::Event &outputEvent)
void pruneEvent(Pythia8::Event &event, Select select)
static std::atomic< int > Pythia8InstanceCounter
void setHooksFileName(std::string val)
void investigateRelatives(Pythia8::Event &event, const std::vector< int > &old2New, size_t index, std::vector< bool > &done, GetRelatives getter, SetRelatives setter, FirstLastRelative firstLast, const std::string &what, const std::string &ind="")
void getNremn(int &nProtonProj, int &nNeutronProj, int &nProtonTarg, int &nNeutronTarg)
void seedGenerator()
performs seeding of the random state of Pythia (called from Init)
void setUserHooks(Pythia8::UserHooks *hooks)
void updateHeader(o2::dataformats::MCEventHeader *eventHeader) override
void getNfreeSpec(int &nFreenProj, int &nFreepProj, int &nFreenTarg, int &nFreepTarg)
void setConfig(std::string val)
std::string mInterfaceName
Definition Generator.h:124
std::vector< TParticle > mParticles
Definition Generator.h:144
Bool_t Init() override
Definition Generator.cxx:57
TF1 const & getfProtons() const
TF1 const & getfNeutrons() const
TF1 const & getsigmaNeutrons() const
TF1 const & getsigmaProtons() const
struct _cl_event * event
Definition glcorearb.h:2982
GLuint index
Definition glcorearb.h:781
GLuint const GLchar * name
Definition glcorearb.h:781
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLint GLint GLint GLint GLint GLint GLint GLbitfield GLenum filter
Definition glcorearb.h:1308
GLubyte GLubyte GLubyte GLubyte w
Definition glcorearb.h:852
GLuint id
Definition glcorearb.h:650
std::vector< InputSpec > select(char const *matcher="")
std::string expandShellVarsInFileName(std::string const &input)
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
std::string to_string(gsl::span< T, Size > span)
Definition common.h:52
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"