Add Non-HFE–enhanced dataset

MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_nonhfe.C

@@ -0,0 +1,218 @@
+#include "FairGenerator.h"
+
+#include "TDatabasePDG.h"
+#include "TRandom.h"
+
+#include "Generators/GeneratorPythia8.h"
+#include "Generators/GeneratorPythia8Param.h"
+#include "Pythia8/Pythia.h"
+#include <fairlogger/Logger.h>
+
+#include <string>
+#include <vector>
+
+using namespace Pythia8;
+
+class GeneratorPythia8GapTriggeredPionEta : public o2::eventgen::GeneratorPythia8
+{
+ public:
+  /// default constructor
+  GeneratorPythia8GapTriggeredPionEta() = default;
+  GeneratorPythia8GapTriggeredPionEta(int inputTriggerRatio = 5, std::vector<int> quarkPdgList = {}, std::vector<int> hadronPdgList = {})
+  {
+    mGeneratedEvents = 0;
+    mInverseTriggerRatio = inputTriggerRatio;
+    mQuarkRapidityMin = -1.5;
+    mQuarkRapidityMax = 1.5;
+    mHadRapidityMin = -1.5;
+    mHadRapidityMax = 1.5;
+
+    mQuarkPdg = 0;
+    mHadronPdg = 0;
+    mQuarkPdgList = quarkPdgList;
+    mHadronPdgList = hadronPdgList;
+    Print();
+  }
+
+  ///  Destructor
+  GeneratorPythia8GapTriggeredPionEta() = default;
+  ///  Print the input
+  void Print()
+  {
+    LOG(info) << "********** GeneratorPythia8GapTriggeredHF configuration dump **********";
+    LOG(info) << Form("* Trigger ratio: %d", mInverseTriggerRatio);
+    LOG(info) << Form("* Quark pdg: %d", mQuarkPdg);
+    LOG(info) << Form("* Quark rapidity: %f - %f", mQuarkRapidityMin, mQuarkRapidityMax);
+    LOG(info) << Form("* Hadron pdg: %d", mHadronPdg);
+    LOG(info) << Form("* Hadron rapidity: %f - %f", mHadRapidityMin, mHadRapidityMax);
+    LOG(info) << Form("* Quark pdg list: ");
+    for (auto pdg : mQuarkPdgList) {
+      LOG(info) << Form("* %d ", pdg);
+    }
+    LOG(info) << Form("* Hadron pdg list: ");
+    for (auto pdg : mHadronPdgList) {
+      LOG(info) << Form("* %d ", pdg);
+    }
+    LOG(info) << "***********************************************************************";
+  }
+
+  bool Init() override
+  {
+    addSubGenerator(0, "Minimum bias");
+    addSubGenerator(1, "Down injected");
+    addSubGenerator(2, "Up injected");
+    addSubGenerator(3, "Strange injected");
+
+    return o2::eventgen::GeneratorPythia8::Init();
+  }
+  void setQuarkRapidity(float yMin, float yMax)
+  {
+    mQuarkRapidityMin = yMin;
+    mQuarkRapidityMax = yMax;
+  };
+  void setHadronRapidity(float yMin, float yMax)
+  {
+    mHadRapidityMin = yMin;
+    mHadRapidityMax = yMax;
+  };
+  void setUsedSeed(unsigned int seed)
+  {
+    mUsedSeed = seed;
+  };
+  unsigned int getUsedSeed() const
+  {
+    return mUsedSeed;
+  };
+
+ protected:
+  //__________________________________________________________________
+  bool generateEvent() override
+  {
+
+    // Simple straightforward check to alternate generators
+    if (mGeneratedEvents % mInverseTriggerRatio == 0) {
+      int nInjectedEvents = mGeneratedEvents / mInverseTriggerRatio;
+      // Alternate quarks if enabled (with the same ratio)
+      if (mQuarkPdgList.size() >= 1) {
+        int iQuark = nInjectedEvents % mQuarkPdgList.size();
+        mQuarkPdg = mQuarkPdgList[iQuark];
+        LOG(debug) << "SELECTED quark: " << mQuarkPdgList[iQuark];
+      }
+      // Alternate hadrons if enabled (with the same ratio)
+      if (mHadronPdgList.size() >= 1) {
+        int iHadron = (nInjectedEvents / std::max(mQuarkPdgList.size(), 1ul)) % mHadronPdgList.size();
+        mHadronPdg = mHadronPdgList[iHadron];
+        LOG(debug) << "SELECTED hadron: " << mHadronPdgList[iHadron];
+      }
+
+      // Generate event of interest
+      bool genOk = false;
+      while (!genOk) {
+        if (GeneratorPythia8::generateEvent()) {
+          genOk = selectEvent();
+        }
+      }
+      notifySubGenerator(mQuarkPdg);
+    } else {
+      // Generate minimum-bias event
+      bool genOk = false;
+      while (!genOk) {
+        genOk = GeneratorPythia8::generateEvent();
+      }
+      notifySubGenerator(0);
+    }
+
+    mGeneratedEvents++;
+
+    return true;
+  }
+  bool selectEvent()
+  {
+
+    bool isGoodAtPartonLevel{mQuarkPdgList.size() == 0};
+    bool isGoodAtHadronLevel{mHadronPdgList.size() == 0};
+
+    for (auto iPart{0}; iPart < mPythia.event.size(); ++iPart) {
+      // search for Q-Qbar mother with at least one Q in rapidity window
+      if (!isGoodAtPartonLevel) {
+        auto daughterList = mPythia.event[iPart].daughterList();
+        bool hasQ = false, hasQbar = false, atSelectedY = false;
+        for (auto iDau : daughterList) {
+          if (mPythia.event[iDau].id() == mQuarkPdg) {
+            hasQ = true;
+            if ((mPythia.event[iDau].y() > mQuarkRapidityMin) && (mPythia.event[iDau].y() < mQuarkRapidityMax)) {
+              atSelectedY = true;
+            }
+          }
+          if (mPythia.event[iDau].id() == -mQuarkPdg) {
+            hasQbar = true;
+            if ((mPythia.event[iDau].y() > mQuarkRapidityMin) && (mPythia.event[iDau].y() < mQuarkRapidityMax)) {
+              atSelectedY = true;
+            }
+          }
+        }
+        if (hasQ && hasQbar && atSelectedY) {
+          isGoodAtPartonLevel = true;
+        }
+      }
+      // search for hadron in rapidity window
+      if (!isGoodAtHadronLevel) {
+        int id = std::abs(mPythia.event[iPart].id());
+        float rap = mPythia.event[iPart].y();
+        if (id == mHadronPdg && rap > mHadRapidityMin && rap < mHadRapidityMax) {
+          isGoodAtHadronLevel = true;
+        }
+      }
+      // we send the trigger immediately (if there are no particles to replace, that can be different from the trigger ones)
+      if (isGoodAtPartonLevel && isGoodAtHadronLevel) {
+        LOG(debug) << "EVENT SELECTED: Found particle " << mPythia.event[iPart].id() << " at rapidity " << mPythia.event[iPart].y() << "\n";
+        return true;
+      }
+    }
+    // we send the trigger
+    if (isGoodAtPartonLevel && isGoodAtHadronLevel) {
+      return true;
+    }
+
+    return false;
+  };
+
+ private:
+  // Interface to override import particles
+  Pythia8::Event mOutputEvent;
+
+  // Properties of selection
+  int mQuarkPdg;
+  float mQuarkRapidityMin;
+  float mQuarkRapidityMax;
+  int mHadronPdg;
+  float mHadRapidityMin;
+  float mHadRapidityMax;
+  unsigned int mUsedSeed;
+
+  // Control gap-triggering
+  unsigned long long mGeneratedEvents;
+  int mInverseTriggerRatio;
+
+  // Control alternate trigger on charm and beauty quarks
+  std::vector<int> mQuarkPdgList = {};
+
+  // Control alternate trigger on different hadrons
+  std::vector<int> mHadronPdgList = {};
+};
+// Predefined generators:
+// Predefined generators:
+// Charm-enriched
+FairGenerator* GeneratorPythia8GapTriggeredPionAndEta(int inputTriggerRatio, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {}, std::vector<std::array<int, 2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {})
+{
+  auto myGen = new GeneratorPythia8GapTriggeredHF(inputTriggerRatio, std::vector<int>{1, 2, 3}, hadronPdgList, partPdgToReplaceList, freqReplaceList);
+  auto seed = (gRandom->TRandom::GetSeed() % 900000000);
+  myGen->setUsedSeed(seed);
+  myGen->readString("Random:setSeed on");
+  myGen->readString("Random:seed " + std::to_string(seed));
+  myGen->setQuarkRapidity(yQuarkMin, yQuarkMax);
+  if (hadronPdgList.size() != 0) {
+    myGen->setHadronRapidity(yHadronMin, yHadronMax);
+  }
+  return myGen;
+}

MC/config/PWGHF/ini/GeneratorHF_Non_Hfe.ini

@@ -0,0 +1,8 @@
+### The external generator derives from GeneratorPythia8.
+[GeneratorExternal]
+fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_nonhfe.C
+funcName=GeneratorPythia8GapTriggeredCharmAndBeauty(5, -1.5, 1.5, -1.0, 1.0, {11})
+
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_NonHfe.cfg
+includePartonEvent=true

MC/config/PWGHF/ini/tests/GeneratorHF_Non_Hfe.C

@@ -0,0 +1,157 @@
+int External() {
+  std::string path{"o2sim_Kine.root"};
+
+  int checkPdgDecayElectron = 11;
+  int checkPdgQuarkOne = 1;         // d quark
+  int checkPdgQuarkTwo = 2;         // u quark
+  int checkPdgQuarkThree = 3;       // s quark
+  float ratioTrigger = 1. / 5;      // one event triggered out of 5
+
+
+  TFile file(path.c_str(), "READ");
+  if (file.IsZombie()) {
+    std::cerr << "Cannot open ROOT file" << path << "\n";
+    return 1;
+  }
+
+  auto tree = (TTree *)file.Get("o2sim");
+  if (!tree) {
+    std::cerr << "Cannot find tree o2sim in file" << path << "\n";
+    return 1;
+  }
+
+  std::vector<o2::MCTrack> *tracks{};
+  tree->SetBranchAddress("MCTrack", &tracks);
+  o2::dataformats::MCEventHeader *eventHeader = nullptr;
+  tree->SetBranchAddress("MCEventHeader.", &eventHeader);
+
+  int nEventsMB{}, nEventsInjOne{}, nEventsInjTwo{}, nEventsInjThree{};
+  int nQuarksOne{}, nQuarksTwo{}, nQuarksThree{};
+  int nElectrons{};
+  auto nEvents = tree->GetEntries();
+
+  for (int i = 0; i < nEvents; i++) {
+    tree->GetEntry(i);
+
+    // check subgenerator information
+    if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) {
+      bool isValid = false;
+      int subGeneratorId = eventHeader->getInfo<int>(
+          o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid);
+      if (subGeneratorId == 0) {
+        nEventsMB++;
+      } else if (subGeneratorId == checkPdgQuarkOne) {
+        nEventsInjOne++;
+      } else if (subGeneratorId == checkPdgQuarkTwo) {
+        nEventsInjTwo++;
+      } else if (subGeneratorId == checkPdgQuarkThree) {
+        nEventsInjThree++;
+      }
+    } // if event header
+
+    int nelectronsev = 0;
+
+    for (auto &track : *tracks) {
+      auto pdg = track.GetPdgCode();
+      if (std::abs(pdg) == checkPdgQuarkOne) {
+        nQuarksOne++;
+        continue;
+      }
+      if (std::abs(pdg) == checkPdgQuarkTwo) {
+        nQuarksTwo++;
+        continue;
+      }
+      if (std::abs(pdg) == checkPdgQuarkThree) {
+             nQuarksThree++;
+              continue;
+       }
+
+
+      auto y = track.GetRapidity();
+      if (std::abs(pdg) == checkPdgDecayElectron) {
+        int igmother = track.getMotherTrackId();
+        auto gmTrack = (*tracks)[igmother];
+        int gmpdg = gmTrack.GetPdgCode();
+        if (int(std::abs(gmpdg) / 100.) == 1 ||
+            int(std::abs(gmpdg) / 1000.) == 1 ||
+            int(std::abs(gmpdg) / 100.) == 2 ||
+            int(std::abs(gmpdg) / 1000.) == 2 ||
+            int(std::abs(gmpdg) / 100.) == 3 ||
+            int(std::abs(gmpdg) / 1000.) == 3) {
+          nElectrons++;
+          nelectronsev++;
+        } // gmpdg
+      } // pdgdecay
+    } // loop track
+    // std::cout << "#electrons per event: " << nelectronsev << "\n";
+  }
+
+  std::cout << "--------------------------------\n";
+  std::cout << "# Events: " << nEvents << "\n";
+  std::cout << "# MB events: " << nEventsMB << "\n";
+  std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkOne)
+            << nEventsInjOne << "\n";
+  std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkTwo)
+            << nEventsInjTwo << "\n";
+  std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkThree)
+            << nEventsInjThree << "\n";
+  std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkOne) << nQuarksOne
+            << "\n";
+  std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkTwo) << nQuarksTwo
+            << "\n";
+  std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkThree) << nQuarksThree
+            << "\n";
+
+  if (nEventsMB < nEvents * (1 - ratioTrigger) * 0.95 ||
+      nEventsMB > nEvents * (1 - ratioTrigger) *
+                      1.05) { // we put some tolerance since the number of
+                              // generated events is small
+    std::cerr << "Number of generated MB events different than expected\n";
+    return 1;
+  }
+  if (nEventsInjOne < nEvents * ratioTrigger * 0.5 * 0.95 ||
+      nEventsInjOne > nEvents * ratioTrigger * 0.5 * 1.05) {
+    std::cerr << "Number of generated events injected with " << checkPdgQuarkOne
+              << " different than expected\n";
+    return 1;
+  }
+  if (nEventsInjTwo < nEvents * ratioTrigger * 0.5 * 0.95 ||
+      nEventsInjTwo > nEvents * ratioTrigger * 0.5 * 1.05) {
+    std::cerr << "Number of generated events injected with " << checkPdgQuarkTwo
+              << " different than expected\n";
+    return 1;
+  }
+  if (nEventsInjThree < nEvents * ratioTrigger * 0.5 * 0.95 ||
+      nEventsInjThree > nEvents * ratioTrigger * 0.5 * 1.05) {
+    std::cerr << "Number of generated events injected with " << checkPdgQuarkThree
+              << " different than expected\n";
+    return 1;
+  }
+  if (nQuarksOne <
+      nEvents *
+          ratioTrigger) { // we expect anyway more because the same quark is
+                          // repeated several time, after each gluon radiation
+    std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkOne
+              << " lower than expected\n";
+    return 1;
+  }
+  if (nQuarksTwo <
+      nEvents *
+          ratioTrigger) { // we expect anyway more because the same quark is
+                          // repeated several time, after each gluon radiation
+    std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkTwo
+              << " lower than expected\n";
+    return 1;
+  }
+    if (nQuarksThree <
+      nEvents *
+          ratioTrigger) { // we expect anyway more because the same quark is
+                          // repeated several time, after each gluon radiation
+    std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkThree
+              << " lower than expected\n";
+    return 1;
+  }
+  std::cout << "#electrons: " << nElectrons << "\n";
+
+  return 0;
+} // external