[PWGHF] Add possibility to use EVTGEN in HF generator + config for B2Jpsi

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

@@ -1,7 +1,12 @@
+R__LOAD_LIBRARY(EvtGen)
+R__ADD_INCLUDE_PATH($EVTGEN_ROOT/include)
+
 #include "FairGenerator.h"
 #include "Generators/GeneratorPythia8.h"
 #include "Generators/GeneratorPythia8Param.h"
+#include "EvtGen/EvtGen.hh"
 #include "Pythia8/Pythia.h"
+#include "Pythia8Plugins/EvtGen.h"
 #include "TRandom.h"
 #include "TDatabasePDG.h"
 #include <fairlogger/Logger.h>
@@ -33,6 +38,9 @@ public:
     mHadronPdgList = hadronPdgList;
     mPartPdgToReplaceList = partPdgToReplaceList;
     mFreqReplaceList = freqReplaceList;
+    mEvtGen = nullptr;
+    mUseEvtGen = false;
+    mEvtGenDecTable = "";
     // Ds1*(2700), Ds1*(2860), Ds3*(2860), Xic(3055)+, Xic(3080)+, Xic(3055)0, Xic(3080)0, LambdaC(2625), LambdaC(2595), LambdaC(2860), LambdaC(2880), LambdaC(2940), ThetaC(3100)
     mCustomPartPdgs = {30433, 40433, 437, 4315, 4316, 4325, 4326, 4124, 14122, 24124, 24126, 4125, 9422111};
     mCustomPartMasses[30433] = 2.714f;
@@ -114,6 +122,10 @@ public:
       pdgToReplace.push_back(mPartPdgToReplaceList[iRepl].at(0));
     }
 
+    if (mUseEvtGen) {
+      mEvtGen = new Pythia8::EvtGenDecays(&mPythia, mEvtGenDecTable.data(), gSystem->ExpandPathName("$EVTGEN_ROOT/share/EvtGen/evt.pdl"));
+    }
+
     return o2::eventgen::GeneratorPythia8::Init();
   }
 
@@ -135,6 +147,14 @@ public:
   {
     return mUsedSeed;
   };
+  void setUseEvtGenDecayer(std::string evtGenDecTable = "") {
+    mUseEvtGen = true;
+    if (evtGenDecTable.empty()) {
+      mEvtGenDecTable = gSystem->ExpandPathName("$EVTGEN_ROOT/share/EvtGen/DECAY.DEC");
+    } else {
+      mEvtGenDecTable = gSystem->ExpandPathName(evtGenDecTable.data());
+    }
+  }
 
 protected:
   //__________________________________________________________________
@@ -167,6 +187,10 @@ protected:
       {
         if (GeneratorPythia8::generateEvent())
         {
+          if (mUseEvtGen) {
+            mEvtGen->decay();
+            checkConsistency();
+          }
           genOk = selectEvent();
         }
       }
@@ -179,6 +203,12 @@ protected:
       while (!genOk)
       {
         genOk = GeneratorPythia8::generateEvent();
+        if (genOk) {
+          if (mUseEvtGen) {
+            mEvtGen->decay();
+            checkConsistency();
+          }
+        }
       }
       notifySubGenerator(0);
     }
@@ -197,6 +227,8 @@ protected:
 
     for (auto iPart{0}; iPart < mPythia.event.size(); ++iPart)
     {
+
+
       // search for Q-Qbar mother with at least one Q in rapidity window
       if (!isGoodAtPartonLevel)
       {
@@ -352,11 +384,22 @@ protected:
     return true;
   }
 
+  void checkConsistency() {
+    for (int iPart{1}; iPart<mPythia.event.size(); ++iPart) {
+      // verify if all particles of decay chain are in the TDatabasePDG
+      // taken from https://github.com/AliceO2Group/O2DPG/blob/master/MC/config/PWGDQ/EvtGen/GeneratorEvtGen.C
+      if (!TDatabasePDG::Instance()->GetParticle(abs(mPythia.event[iPart].id()))) {
+        // std::cout << "Particle code non known in TDatabasePDG " << mPythia.event[iPart].id() << " - set pdg = 89" << std::endl;
+        mPythia.event[iPart].id(89);
+      }
+    }
+  }
+
   private:
   // Interface to override import particles
   Pythia8::Event mOutputEvent;
 
-  // Properties of selection
+// Properties of selection
   int mQuarkPdg;
   float mQuarkRapidityMin;
   float mQuarkRapidityMax;
@@ -379,6 +422,12 @@ protected:
 
   // Control alternate trigger on different hadrons
   std::vector<int> mHadronPdgList = {};
+
+  // EVTGEN decayer from PYTHIA8 plugins
+  Pythia8::EvtGenDecays *mEvtGen;
+  bool mUseEvtGen;
+  std::string mEvtGenDecTable;
+
 };
 
 // Predefined generators:
@@ -446,3 +495,20 @@ FairGenerator *GeneratorPythia8GapHF(int inputTriggerRatio, float yQuarkMin = -1
 
   return myGen;
 }
+
+// Beauty-enriched with EvtGen decayer
+FairGenerator *GeneratorPythia8GapTriggeredBeautyWithEvtGen(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 = {}, std::string decayTable = "")
+{
+  auto myGen = new GeneratorPythia8GapTriggeredHF(inputTriggerRatio, std::vector<int>{5}, 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);
+  }
+  myGen->setUseEvtGenDecayer(decayTable);
+  return myGen;
+}

MC/config/PWGHF/ini/GeneratorHFTrigger_Bforced.ini

@@ -1,4 +1,4 @@
-#NEV_TEST> 20
+#NEV_TEST> 50
 [GeneratorExternal]
 fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
 funcName=GeneratorPythia8GapTriggeredBeauty(3, -5, 5)

MC/config/PWGHF/ini/GeneratorHF_D2H_bbbar_B2JPsi_gap5_Mode2.ini

@@ -0,0 +1,9 @@
+#NEV_TEST> 20
+### The external generator derives from GeneratorPythia8.
+[GeneratorExternal]
+fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
+funcName=GeneratorPythia8GapTriggeredBeautyWithEvtGen(5, -1.5, 1.5, -1.5, 1.5, {}, {}, {}, "${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/decayer/EVTGEN_DECAY_B2JPSI.DEC")
+
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_Mode2_noforcedecays.cfg
+includePartonEvent=false

MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkgSigmaC.ini

@@ -1,4 +1,4 @@
-#NEV_TEST> 40
+#NEV_TEST> 100
 ### The external generator derives from GeneratorPythia8.
 [GeneratorExternal]
 fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C

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

@@ -0,0 +1,113 @@
+int External() {
+    std::string path{"o2sim_Kine.root"};
+
+    int checkPdgQuark{5};
+    float ratioTrigger = 1./5; // one event triggered out of 5
+
+    std::vector<int> checkPdgHadron{443, 511, 521, 531, 5122};
+    std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted pdg of daughters
+        {443, {{-11, 11}, {-13, 13}}}, // J/psi
+        {511, {{313, 443}}}, // B0
+        {521, {{321, 443}}}, // B+
+        {531, {{333, 443}}}, // Bs0
+        {5122, {{321, 443, 2212}}} // Lb0
+    };
+
+    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");
+    std::vector<o2::MCTrack> *tracks{};
+    tree->SetBranchAddress("MCTrack", &tracks);
+    o2::dataformats::MCEventHeader *eventHeader = nullptr;
+    tree->SetBranchAddress("MCEventHeader.", &eventHeader);
+
+    int nEventsMB{}, nEventsInj{};
+    int nSignals{}, nSignalGoodDecay{};
+    int nJPsiToEE{}, nJPsiToMuMu{};
+    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 == checkPdgQuark) {
+                nEventsInj++;
+            }
+        }
+
+        for (auto &track : *tracks) {
+            auto pdg = track.GetPdgCode();
+            if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), std::abs(pdg)) != checkPdgHadron.end()) { // found signal
+                nSignals++; // count signal PDG
+
+                std::vector<int> pdgsDecay{};
+                std::vector<int> pdgsDecayAntiPart{};
+                for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
+                    auto pdgDau = tracks->at(j).GetPdgCode();
+                    if (pdg == 443 && pdgDau == 22) {
+                        continue;
+                    }
+                    pdgsDecay.push_back(pdgDau);
+                    if (pdgDau != 333) { // phi is antiparticle of itself
+                        pdgsDecayAntiPart.push_back(-pdgDau);
+                    } else {
+                        pdgsDecayAntiPart.push_back(pdgDau);
+                    }
+                }
+
+                std::sort(pdgsDecay.begin(), pdgsDecay.end());
+                std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end());
+
+                for (auto &decay : checkHadronDecays[std::abs(pdg)]) {
+                    if (pdgsDecay == decay || pdgsDecayAntiPart == decay) {
+                        nSignalGoodDecay++;
+                        if (pdg == 443) {
+                            if (decay == std::vector{-11, 11}) {
+                                nJPsiToEE++;
+                            } else if (decay == std::vector{-13, 13}) {
+                                nJPsiToMuMu++;
+                            }
+                        }
+                        break;
+                    }
+                }
+            }
+        }
+    }
+
+    std::cout << "--------------------------------\n";
+    std::cout << "# Events: " << nEvents << "\n";
+    std::cout << "# MB events: " << nEventsMB << "\n";
+    std::cout << Form("# events injected with %d quark pair: ", checkPdgQuark) << nEventsInj << "\n";
+    std::cout <<"# signal hadrons: " << nSignals << "\n";
+    std::cout <<"# signal hadrons decaying in the correct channel: " << nSignalGoodDecay << "\n";
+    std::cout <<"# J/Psi decaying to e+e-: " << nJPsiToEE << "\n";
+    std::cout <<"# J/Psi decaying to mu+mu-: " << nJPsiToMuMu << "\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 (nEventsInj < nEvents * ratioTrigger * 0.95 || nEventsInj > nEvents * ratioTrigger * 1.05) {
+        std::cerr << "Number of generated events injected with " << checkPdgQuark << " different than expected\n";
+        return 1;
+    }
+
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) : 1.0f;
+    if (1 - fracForcedDecays > 0.5 + uncFracForcedDecays) { // at least 50% in the main decay channels (we also have correlated backgrounds, mostly from other B decays)
+        std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
+        return 1;
+    }
+
+    return 0;
+}