|
| 1 | +R__ADD_INCLUDE_PATH($O2DPG_ROOT) |
| 2 | +///#include "FairGenerator.h" |
| 3 | +//#include "Generators/GeneratorPythia8.h" |
| 4 | +#include "Pythia8/Pythia.h" |
| 5 | +#include "MC/config/PWGGAJE/hooks/jets_hook.C" |
| 6 | +//#include "TRandom.h" |
| 7 | +//#include <fairlogger/Logger.h> |
| 8 | +// |
| 9 | +//#include <string> |
| 10 | +//#include <vector> |
| 11 | + |
| 12 | +// Jet-jet custom event generator |
| 13 | +// that alternates between 2 gun generators. |
| 14 | +// set up to inject MB events alongside jet-jet events |
| 15 | +// in 'MB-gap' mode. |
| 16 | +// The number of MB events injected, and the PYTHIA config |
| 17 | +// for each event type is defined by the user in the .ini |
| 18 | +// generator file (e.g. GeneratorJE_gapgen5_hook.ini) |
| 19 | +// |
| 20 | +// Author: Jaime Norman (jaime.norman@cern.ch) |
| 21 | + |
| 22 | +// o2-sim-dpl-eventgen --nEvents 10 --generator external --configKeyValues "GeneratorExternal.fileName=generator_pythia8_gaptrigger_jets_pythiabase.C;GeneratorExternal.funcName=getGeneratorPythia8GapGenJE()" |
| 23 | + |
| 24 | +namespace o2 |
| 25 | +{ |
| 26 | +namespace eventgen |
| 27 | +{ |
| 28 | + |
| 29 | +using namespace Pythia8; |
| 30 | + |
| 31 | + |
| 32 | +/// A very simple gap generator alternating between 2 different particle guns |
| 33 | +class GeneratorPythia8GapGenJE : public o2::eventgen::GeneratorPythia8 |
| 34 | +{ |
| 35 | +public: |
| 36 | + /// default constructor |
| 37 | + GeneratorPythia8GapGenJE(int inputTriggerRatio = 5,std::string pathMB = "",std::string pathSignal = "") { |
| 38 | + |
| 39 | + mGeneratedEvents = 0; |
| 40 | + mInverseTriggerRatio = inputTriggerRatio; |
| 41 | + |
| 42 | + auto seed = (gRandom->TRandom::GetSeed() % 900000000); |
| 43 | + |
| 44 | + cout << "Initalizing extra PYTHIA object used to generate min-bias events..." << endl; |
| 45 | + TString pathconfigMB = gSystem->ExpandPathName(TString(pathMB)); |
| 46 | + pythiaObjectMinimumBias.readFile(pathconfigMB.Data()); |
| 47 | + pythiaObjectMinimumBias.readString("Random:setSeed on"); |
| 48 | + pythiaObjectMinimumBias.readString("Random:seed " + std::to_string(seed)); |
| 49 | + pythiaObjectMinimumBias.init(); |
| 50 | + cout << "Initalization complete" << endl; |
| 51 | + cout << "Initalizing extra PYTHIA object used to generate signal events..." << endl; |
| 52 | + TString pathconfigSignal = gSystem->ExpandPathName(TString(pathSignal)); |
| 53 | + pythiaObjectSignal.readFile(pathconfigSignal.Data()); |
| 54 | + pythiaObjectSignal.readString("Random:setSeed on"); |
| 55 | + pythiaObjectSignal.readString("Random:seed " + std::to_string(seed)); |
| 56 | + // load jet hook to ensure at least one jet is within detector acceptance |
| 57 | + Pythia8::UserHooks *hook = pythia8_userhooks_jets(); |
| 58 | + pythiaObjectSignal.setUserHooksPtr(std::shared_ptr<Pythia8::UserHooks>(hook)); |
| 59 | + pythiaObjectSignal.init(); |
| 60 | + cout << "Initalization complete" << endl; |
| 61 | + // Add Sub generators |
| 62 | + addSubGenerator(0, "MB generator"); |
| 63 | + addSubGenerator(1, "jet-jet generator"); |
| 64 | + } |
| 65 | + |
| 66 | + |
| 67 | + /// Destructor |
| 68 | + ~GeneratorPythia8GapGenJE() = default; |
| 69 | + |
| 70 | + void setUsedSeed(unsigned int seed) |
| 71 | + { |
| 72 | + mUsedSeed = seed; |
| 73 | + }; |
| 74 | + unsigned int getUsedSeed() const |
| 75 | + { |
| 76 | + return mUsedSeed; |
| 77 | + }; |
| 78 | + |
| 79 | + bool generateEvent() override |
| 80 | + { |
| 81 | + |
| 82 | + // Simple straightforward check to alternate generators |
| 83 | + mPythia.event.reset(); |
| 84 | + |
| 85 | + if (mGeneratedEvents % mInverseTriggerRatio == 0) { |
| 86 | + LOG(info) << "Event " << mGeneratedEvents << ", generate signal event"; |
| 87 | + // Generate event of interest |
| 88 | + Bool_t mGenerationOK = kFALSE; |
| 89 | + while (!mGenerationOK) { |
| 90 | + mGenerationOK = pythiaObjectSignal.next(); |
| 91 | + } |
| 92 | + mPythia.event = pythiaObjectSignal.event; |
| 93 | + setEventHeaderProperties(pythiaObjectSignal); |
| 94 | + LOG(info) << "--- Print info properties custom..."; |
| 95 | + printEventHeaderProperties(pythiaObjectSignal); |
| 96 | + notifySubGenerator(1); |
| 97 | + } |
| 98 | + else { |
| 99 | + LOG(info) << "Event " << mGeneratedEvents << ", generate mb event"; |
| 100 | + // Generate minimum-bias event |
| 101 | + Bool_t mGenerationOK = kFALSE; |
| 102 | + while (!mGenerationOK) { |
| 103 | + mGenerationOK = pythiaObjectMinimumBias.next(); |
| 104 | + } |
| 105 | + mPythia.event = pythiaObjectMinimumBias.event; |
| 106 | + setEventHeaderProperties(pythiaObjectMinimumBias); |
| 107 | + LOG(info) << "--- Print info properties custom..."; |
| 108 | + printEventHeaderProperties(pythiaObjectMinimumBias); |
| 109 | + notifySubGenerator(0); |
| 110 | + } |
| 111 | + mGeneratedEvents++; |
| 112 | + return true; |
| 113 | + } |
| 114 | + |
| 115 | + // for testing |
| 116 | + void printEventHeaderProperties (Pythia8::Pythia &pythiaObject) { |
| 117 | + LOG(info) << "Info name = " << pythiaObject.info.name(); |
| 118 | + LOG(info) << "Info code = " << pythiaObject.info.code(); |
| 119 | + LOG(info) << "Info weight = " << pythiaObject.info.weight(); |
| 120 | + LOG(info) << "Info id1pdf = " << pythiaObject.info.id1pdf(); |
| 121 | + LOG(info) << "Info id2pdf = " << pythiaObject.info.id2pdf(); |
| 122 | + |
| 123 | + LOG(info) << "Info x1pdf = " << pythiaObject.info.x1pdf(); |
| 124 | + LOG(info) << "Info x2pdf = " << pythiaObject.info.x2pdf(); |
| 125 | + LOG(info) << "Info QFac = " << pythiaObject.info.QFac(); |
| 126 | + LOG(info) << "Info pdf1 = " << pythiaObject.info.pdf1(); |
| 127 | + LOG(info) << "Info pdf2 = " << pythiaObject.info.pdf2(); |
| 128 | + |
| 129 | + // Set cross section |
| 130 | + LOG(info) << "Info sigmaGen = " << pythiaObject.info.sigmaGen(); |
| 131 | + LOG(info) << "Info sigmaErr = " << pythiaObject.info.sigmaErr(); |
| 132 | + |
| 133 | + |
| 134 | + // Set event scale and nMPI |
| 135 | + LOG(info) << "Info QRen = " << pythiaObject.info.QRen(); |
| 136 | + LOG(info) << "Info nMPI = " << pythiaObject.info.nMPI(); |
| 137 | + |
| 138 | + // Set weights (overrides cross-section for each weight) |
| 139 | + size_t iw = 0; |
| 140 | + auto xsecErr = pythiaObject.info.weightContainerPtr->getTotalXsecErr(); |
| 141 | + for (auto w : pythiaObject.info.weightContainerPtr->getTotalXsec()) { |
| 142 | + std::string post = (iw == 0 ? "" : "_" + std::to_string(iw)); |
| 143 | + LOG(info) << "Info weight by index = " << pythiaObject.info.weightValueByIndex(iw); |
| 144 | + iw++; |
| 145 | + } |
| 146 | + |
| 147 | + } |
| 148 | + |
| 149 | + // in order to save the event weight we need to override the following function |
| 150 | + // from the inherited o2::eventgen::GeneratorPythia8 class. The event header properties |
| 151 | + // are created as members of this class, and are set using the active event generator |
| 152 | + // (MB or jet-jet), then propagated to the event header |
| 153 | + void updateHeader(o2::dataformats::MCEventHeader* eventHeader) override { |
| 154 | + /** update header **/ |
| 155 | + using Key = o2::dataformats::MCInfoKeys; |
| 156 | + |
| 157 | + eventHeader->putInfo<std::string>(Key::generator, "pythia8"); |
| 158 | + eventHeader->putInfo<int>(Key::generatorVersion, PYTHIA_VERSION_INTEGER); |
| 159 | + eventHeader->putInfo<std::string>(Key::processName, name); |
| 160 | + eventHeader->putInfo<int>(Key::processCode, code); |
| 161 | + eventHeader->putInfo<float>(Key::weight, weight); |
| 162 | + |
| 163 | + // Set PDF information |
| 164 | + eventHeader->putInfo<int>(Key::pdfParton1Id, id1pdf); |
| 165 | + eventHeader->putInfo<int>(Key::pdfParton2Id, id2pdf); |
| 166 | + eventHeader->putInfo<float>(Key::pdfX1, x1pdf); |
| 167 | + eventHeader->putInfo<float>(Key::pdfX2, x2pdf); |
| 168 | + eventHeader->putInfo<float>(Key::pdfScale, QFac); |
| 169 | + eventHeader->putInfo<float>(Key::pdfXF1, pdf1); |
| 170 | + eventHeader->putInfo<float>(Key::pdfXF2, pdf2); |
| 171 | + |
| 172 | + // Set cross section |
| 173 | + eventHeader->putInfo<float>(Key::xSection, sigmaGen * 1e9); |
| 174 | + eventHeader->putInfo<float>(Key::xSectionError, sigmaErr * 1e9); |
| 175 | + |
| 176 | + // Set event scale and nMPI |
| 177 | + eventHeader->putInfo<float>(Key::eventScale, QRen); |
| 178 | + eventHeader->putInfo<int>(Key::mpi, nMPI); |
| 179 | + LOG(info) << "--- updated header weight = " << weight; |
| 180 | + |
| 181 | + // The following is also set in the base class updateHeader function |
| 182 | + // but as far as I can tell, there is no Xsec weight in the default |
| 183 | + // header so this is not copied over for now |
| 184 | + |
| 185 | + //size_t iw = 0; |
| 186 | + //auto xsecErr = info.weightContainerPtr->getTotalXsecErr(); |
| 187 | + //for (auto w : info.weightContainerPtr->getTotalXsec()) { |
| 188 | + // std::string post = (iw == 0 ? "" : "_" + std::to_string(iw)); |
| 189 | + // eventHeader->putInfo<float>(Key::weight + post, info.weightValueByIndex(iw)); |
| 190 | + // eventHeader->putInfo<float>(Key::xSection + post, w * 1e9); |
| 191 | + // eventHeader->putInfo<float>(Key::xSectionError + post, xsecErr[iw] * 1e9); |
| 192 | + // iw++; |
| 193 | + //} |
| 194 | + } |
| 195 | + |
| 196 | + void setEventHeaderProperties (Pythia8::Pythia &pythiaObject) { |
| 197 | + |
| 198 | + auto& info = pythiaObject.info; |
| 199 | + |
| 200 | + name = info.name(); |
| 201 | + code = info.code(); |
| 202 | + weight = info.weight(); |
| 203 | + // Set PDF information |
| 204 | + id1pdf = info.id1pdf(); |
| 205 | + id2pdf = info.id2pdf(); |
| 206 | + x1pdf = info.x1pdf(); |
| 207 | + x2pdf = info.x2pdf(); |
| 208 | + QFac = info.QFac(); |
| 209 | + pdf1 = info.pdf1(); |
| 210 | + pdf2 = info.pdf2(); |
| 211 | + // Set cross section |
| 212 | + sigmaGen = info.sigmaGen(); |
| 213 | + sigmaErr = info.sigmaErr(); |
| 214 | + // Set event scale and nMPI |
| 215 | + QRen = info.QRen(); |
| 216 | + nMPI = info.nMPI(); |
| 217 | + } |
| 218 | + |
| 219 | +private: |
| 220 | + // Interface to override import particles |
| 221 | + Pythia8::Event mOutputEvent; |
| 222 | + |
| 223 | + // Properties of selection |
| 224 | + unsigned int mUsedSeed; |
| 225 | + |
| 226 | + // Control gap-triggering |
| 227 | + unsigned long long mGeneratedEvents; |
| 228 | + int mInverseTriggerRatio; |
| 229 | + |
| 230 | + // Handling generators |
| 231 | + Pythia8::Pythia pythiaObjectMinimumBias; |
| 232 | + Pythia8::Pythia pythiaObjectSignal; |
| 233 | + |
| 234 | + // header info - needed to save event properties |
| 235 | + std::string name; |
| 236 | + int code; |
| 237 | + float weight; |
| 238 | + // PDF information |
| 239 | + int id1pdf; |
| 240 | + int id2pdf; |
| 241 | + float x1pdf; |
| 242 | + float x2pdf; |
| 243 | + float QFac; |
| 244 | + float pdf1; |
| 245 | + float pdf2; |
| 246 | + // cross section |
| 247 | + float sigmaGen; |
| 248 | + float sigmaErr; |
| 249 | + // event scale and nMPI |
| 250 | + float QRen; |
| 251 | + int nMPI; |
| 252 | + |
| 253 | +}; |
| 254 | + |
| 255 | +} // namespace eventgen |
| 256 | +} // namespace o2 |
| 257 | + |
| 258 | +/** generator instance and settings **/ |
| 259 | + |
| 260 | +FairGenerator* getGeneratorPythia8GapGenJE(int inputTriggerRatio = 5, std::string pathMB = "",std::string pathSignal = "") { |
| 261 | + auto myGen = new o2::eventgen::GeneratorPythia8GapGenJE(inputTriggerRatio, pathMB, pathSignal); |
| 262 | + auto seed = (gRandom->TRandom::GetSeed() % 900000000); |
| 263 | + myGen->setUsedSeed(seed); |
| 264 | + myGen->readString("Random:setSeed on"); |
| 265 | + myGen->readString("Random:seed " + std::to_string(seed)); |
| 266 | + myGen->readString("HardQCD:all = on"); |
| 267 | + return myGen; |
| 268 | +} |
0 commit comments