Add EVTGEN in HF generator + config for B2Jpsi

Author: fgrosa

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,8 @@ public:
     mHadronPdgList = hadronPdgList;
     mPartPdgToReplaceList = partPdgToReplaceList;
     mFreqReplaceList = freqReplaceList;
+    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 +121,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 +146,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 +186,10 @@ protected:
       {
         if (GeneratorPythia8::generateEvent())
         {
+          if (mEvtGen) {
+            mEvtGen->decay();
+            checkConsistency();
+          }
           genOk = selectEvent();
         }
       }
@@ -179,6 +202,12 @@ protected:
       while (!genOk)
       {
         genOk = GeneratorPythia8::generateEvent();
+        if (genOk) {
+          if (mEvtGen) {
+            mEvtGen->decay();
+            checkConsistency();
+          }
+        }
       }
       notifySubGenerator(0);
     }
@@ -197,6 +226,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 +383,28 @@ 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
+      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);
+      }
+
+      // istat = mEvtstdhep->getIStat(i);
+
+      // if (istat != 1 && istat != 2)
+      //   std::cout << "ImportParticles: Attention unknown status code!" << std::endl;
+      // if (istat == 2)
+      //   istat = 11; // status decayed
+    }
+  }
+
   private:
   // Interface to override import particles
   Pythia8::Event mOutputEvent;
 
-  // Properties of selection
+// Properties of selection
   int mQuarkPdg;
   float mQuarkRapidityMin;
   float mQuarkRapidityMax;
@@ -379,6 +427,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 +500,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/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