|
| 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 "fairlogger/Logger.h" |
| 12 | +#include <cmath> |
| 13 | +#include <fstream> |
| 14 | +#include <string> |
| 15 | +#include <vector> |
| 16 | +using namespace Pythia8; |
| 17 | +#endif |
| 18 | + |
| 19 | +/// Pythia8 event generator for pp collisions |
| 20 | +/// Selection of events with leading particle (pt>ptThreshold) and containing at |
| 21 | +/// least one (anti)deuteron produced by simple coalescence |
| 22 | + |
| 23 | +class GeneratorPythia8AntidAndHighPt : public o2::eventgen::GeneratorPythia8 { |
| 24 | +public: |
| 25 | + /// Constructor |
| 26 | + GeneratorPythia8AntidAndHighPt(double p0 = 0.3, double pt_leading = 5.0) |
| 27 | + : o2::eventgen::GeneratorPythia8() { |
| 28 | + mP0 = p0; |
| 29 | + mPt_leading = pt_leading; |
| 30 | + } |
| 31 | + /// Destructor |
| 32 | + ~GeneratorPythia8AntidAndHighPt() = default; |
| 33 | + |
| 34 | + bool Init() override { |
| 35 | + addSubGenerator(0, "Pythia8 with (anti)deuterons and high pt particle"); |
| 36 | + return o2::eventgen::GeneratorPythia8::Init(); |
| 37 | + } |
| 38 | + |
| 39 | +protected: |
| 40 | + bool generateEvent() override { |
| 41 | + fmt::printf(">> Generating event %d\n", mGeneratedEvents); |
| 42 | + |
| 43 | + bool genOk = false; |
| 44 | + int localCounter{0}; |
| 45 | + while (!genOk) { |
| 46 | + if (GeneratorPythia8::generateEvent()) { |
| 47 | + genOk = selectEvent(mPythia.event); |
| 48 | + } |
| 49 | + localCounter++; |
| 50 | + } |
| 51 | + fmt::printf( |
| 52 | + ">> Generation of event of interest successful after %i iterations\n", |
| 53 | + localCounter); |
| 54 | + std::cout << std::endl << std::endl; |
| 55 | + notifySubGenerator(0); |
| 56 | + |
| 57 | + mGeneratedEvents++; |
| 58 | + |
| 59 | + return true; |
| 60 | + } |
| 61 | + |
| 62 | + bool selectEvent(Pythia8::Event &event) { |
| 63 | + |
| 64 | + bool has_particle_of_interest = false; |
| 65 | + |
| 66 | + // Deuteron Mass [GeV] |
| 67 | + double md = 1.87561294257; |
| 68 | + |
| 69 | + // Protons and Neutrons |
| 70 | + vector<int> proton_ID; |
| 71 | + vector<int> neutron_ID; |
| 72 | + vector<int> proton_status; |
| 73 | + vector<int> neutron_status; |
| 74 | + |
| 75 | + // ptMax |
| 76 | + double pt_max{0}; |
| 77 | + |
| 78 | + for (auto iPart{0}; iPart < event.size(); ++iPart) { |
| 79 | + |
| 80 | + // Only final-state particles |
| 81 | + if (event[iPart].status() <= 0) { |
| 82 | + continue; |
| 83 | + } |
| 84 | + |
| 85 | + //(Anti)Proton selection |
| 86 | + if (abs(event[iPart].id()) == 2212) { |
| 87 | + proton_ID.push_back(iPart); |
| 88 | + proton_status.push_back(0); |
| 89 | + } |
| 90 | + |
| 91 | + //(Anti)Neutron selection |
| 92 | + if (abs(event[iPart].id()) == 2112) { |
| 93 | + neutron_ID.push_back(iPart); |
| 94 | + neutron_status.push_back(0); |
| 95 | + } |
| 96 | + |
| 97 | + if (std::abs(event[iPart].eta()) < 0.8 && (!event[iPart].isNeutral()) && |
| 98 | + event[iPart].pT() > pt_max) |
| 99 | + pt_max = event[iPart].pT(); |
| 100 | + } |
| 101 | + |
| 102 | + // Skip Events with no leading particle |
| 103 | + if (pt_max < mPt_leading) |
| 104 | + return false; |
| 105 | + |
| 106 | + if (proton_ID.size() < 1 || neutron_ID.size() < 1) |
| 107 | + return false; |
| 108 | + |
| 109 | + for (uint32_t ip = 0; ip < proton_ID.size(); ++ip) { |
| 110 | + |
| 111 | + // Skip used protons |
| 112 | + if (proton_status[ip] < 0) { |
| 113 | + continue; |
| 114 | + } |
| 115 | + for (uint32_t in = 0; in < neutron_ID.size(); ++in) { |
| 116 | + |
| 117 | + // Skip used neutrons |
| 118 | + if (neutron_status[in] < 0) { |
| 119 | + continue; |
| 120 | + } |
| 121 | + |
| 122 | + auto sign_p = |
| 123 | + event[proton_ID[ip]].id() / abs(event[proton_ID[ip]].id()); |
| 124 | + auto sign_n = |
| 125 | + event[neutron_ID[in]].id() / abs(event[neutron_ID[in]].id()); |
| 126 | + |
| 127 | + auto p1 = event[proton_ID[ip]].p(); |
| 128 | + auto p2 = event[neutron_ID[in]].p(); |
| 129 | + auto p = p1 + p2; |
| 130 | + p1.bstback(p); |
| 131 | + p2.bstback(p); |
| 132 | + |
| 133 | + // Coalescence |
| 134 | + if (p1.pAbs() <= mP0 && p2.pAbs() <= mP0 && sign_p == sign_n) { |
| 135 | + p.e(std::hypot(p.pAbs(), md)); |
| 136 | + event.append(sign_p * 1000010020, 121, 0, 0, 0, 0, 0, 0, p.px(), |
| 137 | + p.py(), p.pz(), p.e(), md); |
| 138 | + event[proton_ID[ip]].statusNeg(); |
| 139 | + event[proton_ID[ip]].daughter1(event.size() - 1); |
| 140 | + event[neutron_ID[in]].statusNeg(); |
| 141 | + event[neutron_ID[in]].daughter1(event.size() - 1); |
| 142 | + proton_status[ip] = -1; |
| 143 | + neutron_status[in] = -1; |
| 144 | + has_particle_of_interest = true; |
| 145 | + } |
| 146 | + } |
| 147 | + } |
| 148 | + |
| 149 | + return has_particle_of_interest; |
| 150 | + } |
| 151 | + |
| 152 | +private: |
| 153 | + double mP0 = 0.3; |
| 154 | + double mPt_leading = 5.0; |
| 155 | + uint64_t mGeneratedEvents = 0; |
| 156 | +}; |
| 157 | + |
| 158 | +///___________________________________________________________ |
| 159 | +FairGenerator *generateAntidAndHighPt(double p0 = 0.3, |
| 160 | + double pt_leading = 5.0) { |
| 161 | + |
| 162 | + auto myGenerator = new GeneratorPythia8AntidAndHighPt(p0, pt_leading); |
| 163 | + auto seed = (gRandom->TRandom::GetSeed() % 900000000); |
| 164 | + myGenerator->readString("Random:setSeed on"); |
| 165 | + myGenerator->readString("Random:seed " + std::to_string(seed)); |
| 166 | + return myGenerator; |
| 167 | +} |
0 commit comments