|
| 1 | +#if !defined(__CLING__) || defined(__ROOTCLING__) |
| 2 | +#include "FairGenerator.h" |
| 3 | +#include "FairPrimaryGenerator.h" |
| 4 | +#include "Generators/GeneratorPythia8.h" |
| 5 | +#include "Pythia8/Pythia.h" |
| 6 | +#include "TDatabasePDG.h" |
| 7 | +#include "TMath.h" |
| 8 | +#include "TParticlePDG.h" |
| 9 | +#include "TRandom3.h" |
| 10 | +#include "TSystem.h" |
| 11 | +#include "TVector2.h" |
| 12 | +#include "fairlogger/Logger.h" |
| 13 | +#include <cmath> |
| 14 | +#include <fstream> |
| 15 | +#include <string> |
| 16 | +#include <vector> |
| 17 | +using namespace Pythia8; |
| 18 | +#endif |
| 19 | + |
| 20 | +/// Event generator using Pythia ropes |
| 21 | +/// One event of interest is triggered after n minimum-bias events. |
| 22 | +/// Event of interest : an event that contains at least two generated (anti-)Lambdas |
| 23 | + |
| 24 | +class GeneratorPythia8DoubleLambda : public o2::eventgen::GeneratorPythia8 |
| 25 | +{ |
| 26 | +public: |
| 27 | + /// Constructor |
| 28 | + GeneratorPythia8DoubleLambda(int gapSize = 4, double minPt = 0.2, double maxPt = 10, double maxEta = 0.8) |
| 29 | + : o2::eventgen::GeneratorPythia8(), |
| 30 | + mGapSize(gapSize), |
| 31 | + mMinPt(minPt), |
| 32 | + mMaxPt(maxPt), |
| 33 | + mMaxEta(maxEta) |
| 34 | + { |
| 35 | + fmt::printf( |
| 36 | + ">> Pythia8 generator: two (anti-)Lambdas, gap = %d, minPtLambda = %f, maxPtLambda = %f, |etaLambda| < %f\n", gapSize, minPt, maxPt, maxEta); |
| 37 | + } |
| 38 | + /// Destructor |
| 39 | + ~GeneratorPythia8DoubleLambda() = default; |
| 40 | + |
| 41 | + bool Init() override |
| 42 | + { |
| 43 | + addSubGenerator(0, "Pythia8 events with two (anti-)Lambdas"); |
| 44 | + return o2::eventgen::GeneratorPythia8::Init(); |
| 45 | + } |
| 46 | + |
| 47 | +protected: |
| 48 | + /// Check if particle is physical primary or from HF decay |
| 49 | + bool isLambdaPhysicalPrimaryOrFromHF(const Pythia8::Particle &p, const Pythia8::Event &event) |
| 50 | + { |
| 51 | + // Select only final-state particles |
| 52 | + if (!p.isFinal()) |
| 53 | + { |
| 54 | + return false; |
| 55 | + } |
| 56 | + |
| 57 | + // Walk up ancestry |
| 58 | + int motherId = p.mother1(); |
| 59 | + |
| 60 | + while (motherId > 0) |
| 61 | + { |
| 62 | + |
| 63 | + // Get mother |
| 64 | + const auto &mother = event[motherId]; |
| 65 | + const int absMotherPdg = std::abs(mother.id()); |
| 66 | + |
| 67 | + // Check if particle is from HF decay |
| 68 | + if ((absMotherPdg / 100 == 4) || (absMotherPdg / 100 == 5) || (absMotherPdg / 1000 == 4) || (absMotherPdg / 1000 == 5)) |
| 69 | + { |
| 70 | + return true; |
| 71 | + } |
| 72 | + |
| 73 | + // Reject non-physical primary hadrons |
| 74 | + if (mother.isHadron() && mother.tau0() > 1.0) |
| 75 | + { |
| 76 | + return false; |
| 77 | + } |
| 78 | + motherId = mother.mother1(); |
| 79 | + } |
| 80 | + return true; |
| 81 | + } |
| 82 | + |
| 83 | + bool generateEvent() override |
| 84 | + { |
| 85 | + fmt::printf(">> Generating event %d\n", mGeneratedEvents); |
| 86 | + |
| 87 | + bool genOk = false; |
| 88 | + int localCounter{0}; |
| 89 | + constexpr int kMaxTries{100000}; |
| 90 | + |
| 91 | + // Accept mGapSize events unconditionally, then one triggered event |
| 92 | + if (mGeneratedEvents % (mGapSize + 1) < mGapSize) |
| 93 | + { |
| 94 | + genOk = GeneratorPythia8::generateEvent(); |
| 95 | + fmt::printf(">> Gap-event (no strangeness check)\n"); |
| 96 | + } |
| 97 | + else |
| 98 | + { |
| 99 | + while (!genOk && localCounter < kMaxTries) |
| 100 | + { |
| 101 | + if (GeneratorPythia8::generateEvent()) |
| 102 | + { |
| 103 | + genOk = selectEvent(mPythia.event); |
| 104 | + } |
| 105 | + localCounter++; |
| 106 | + } |
| 107 | + if (!genOk) |
| 108 | + { |
| 109 | + fmt::printf("Failed to generate triggered event after %d tries\n", kMaxTries); |
| 110 | + return false; |
| 111 | + } |
| 112 | + fmt::printf(">> Triggered event: event accepted after %d iterations (double (anti-)Lambdas)\n", localCounter); |
| 113 | + } |
| 114 | + |
| 115 | + notifySubGenerator(0); |
| 116 | + mGeneratedEvents++; |
| 117 | + return true; |
| 118 | + } |
| 119 | + |
| 120 | + bool selectEvent(Pythia8::Event &event) |
| 121 | + { |
| 122 | + |
| 123 | + std::vector<int> particleID; |
| 124 | + int nLambda{0}; |
| 125 | + |
| 126 | + for (int i = 0; i < event.size(); i++) |
| 127 | + { |
| 128 | + const auto &p = event[i]; |
| 129 | + if (std::abs(p.eta()) > mMaxEta || p.pT() < mMinPt || p.pT() > mMaxPt) |
| 130 | + continue; |
| 131 | + if (!isLambdaPhysicalPrimaryOrFromHF(p, event)) |
| 132 | + continue; |
| 133 | + |
| 134 | + if (std::abs(p.id()) == 3122) |
| 135 | + nLambda++; |
| 136 | + |
| 137 | + particleID.emplace_back(i); |
| 138 | + } |
| 139 | + if (nLambda < 2) |
| 140 | + return false; |
| 141 | + |
| 142 | + return true; |
| 143 | + } |
| 144 | + |
| 145 | +private: |
| 146 | + int mGapSize{4}; |
| 147 | + double mMinPt{0.2}; |
| 148 | + double mMaxPt{10.0}; |
| 149 | + double mMaxEta{0.8}; |
| 150 | + uint64_t mGeneratedEvents{0}; |
| 151 | +}; |
| 152 | + |
| 153 | +///___________________________________________________________ |
| 154 | +FairGenerator *generateDoubleLambda(int gap = 4, double minPt = 0.2, double maxPt = 10, double maxEta = 0.8) |
| 155 | +{ |
| 156 | + |
| 157 | + auto myGenerator = new GeneratorPythia8DoubleLambda(gap, minPt, maxPt, maxEta); |
| 158 | + auto seed = (gRandom->TRandom::GetSeed() % 900000000); |
| 159 | + myGenerator->readString("Random:setSeed on"); |
| 160 | + myGenerator->readString("Random:seed " + std::to_string(seed)); |
| 161 | + return myGenerator; |
| 162 | +} |
0 commit comments