Add ini file for combined c->D->mu and b->B->D->mu

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

@@ -479,6 +479,22 @@ FairGenerator *GeneratorPythia8GapTriggeredCharmAndBeauty(int inputTriggerRatio,
   return myGen;
 }
 
+// Charm and beauty enriched (2/3 charm 1/3 beauty)
+FairGenerator *GeneratorPythia8GapTriggeredCharmBeautyCharm(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>{4, 5, 4}, 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;
+}
+
 FairGenerator *GeneratorPythia8GapHF(int inputTriggerRatio, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> quarkPdgList = {}, std::vector<int> hadronPdgList = {}, std::vector<std::array<int, 2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {})
 {
   if (hadronPdgList.size() == 0 && quarkPdgList.size() == 0)

MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkgSigmaC.ini

@@ -1,4 +1,4 @@
-#NEV_TEST> 100
+#NEV_TEST> 150
 ### 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/GeneratorHF_mu_ccbar_bbbar_charmdecays_gap3_Mode2_accSmall.ini

@@ -0,0 +1,8 @@
+#NEV_TEST> 150
+### The external generator derives from GeneratorPythia8
+[GeneratorExternal]
+fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
+funcName=GeneratorPythia8GapTriggeredCharmBeautyCharm(3, -4.3, -2.2)
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charm_with_mudecays_Mode2.cfg
+includePartonEvent=true

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

@@ -12,7 +12,7 @@ int External() {
     float ratioTrigger = 1./5.; // one event triggered out of 5
     std::array<std::array<int, 2>, 2> pdgReplParticles = {std::array{413, 14122}, std::array{413, 4124}};
     std::array<std::array<int, 2>, 2> pdgReplPartCounters = {std::array{0, 0}, std::array{0, 0}};
-    std::array<float, 2> freqRepl = {0.5, 0.5};
+    std::array<float, 2> freqRepl = {0.35, 0.5};
     std::map<int, int> sumOrigReplacedParticles = {{413, 0}};
 
     std::array<int, 3> checkPdgHadron{413, 14122, 4124};
@@ -120,7 +120,15 @@ int External() {
                         /// check that among the mothers there is a c/cbar quark
                         /// This means that the charm hadron comes from the c-quark hadronization, where the c/cbar quark
                         /// comes from a c-cbar pair present in the current event, tagged with a b-bbar (e.g. double-parton scattering)
-                        if(std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 4) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -4) == motherPdgCodes.end()) {
+                        /// we check also if it comes from a charmed diquark
+                        if(std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 4) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -4) == motherPdgCodes.end()
+                           && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 4101) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -4101) == motherPdgCodes.end()
+                           && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 4103) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -4103) == motherPdgCodes.end()
+                           && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 4201) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -4201) == motherPdgCodes.end()
+                           && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 4203) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -4203) == motherPdgCodes.end()
+                           && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 4301) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -4301) == motherPdgCodes.end()
+                           && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 4303) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -4303) == motherPdgCodes.end()
+                           && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 4403) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -4403) == motherPdgCodes.end()){
                             /// if we arrive here, it means that the hadron comes from the decay of a beauty hadron.
                             /// This can happen if and only if this is not a replaced one (i.e. native from Lambdab0 decay)
                             if (std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 5122) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -5122) == motherPdgCodes.end()) {

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

@@ -0,0 +1,114 @@
+int External() {
+
+  int checkPdgDecayMuon = 13;
+  int checkPdgQuarkOne{4};
+  int checkPdgQuarkTwo{5};
+  float ratioQuarkOneTwo{2.};
+
+  float ratioTrigger = 1./3; // one event triggered out of 5
+
+  std::string path{"o2sim_Kine.root"};
+
+  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{0};
+  int nEventsInjOne{0};
+  int nEventsInjTwo{0};
+  int nQuarks{0};
+  int nMuons{0};
+
+  int nMuonsInAcceptance{};
+
+  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++;
+      }
+    }
+
+    int nmuonsev = 0;
+    int nmuonsevinacc = 0;
+
+    for (auto &track : *tracks) {
+      auto pdg = track.GetPdgCode();
+      if (std::abs(pdg) == checkPdgQuarkOne || std::abs(pdg) == checkPdgQuarkTwo) {
+        nQuarks++;
+        continue;
+      } // pdgquark
+      auto y = track.GetRapidity();
+      if (std::abs(pdg) == checkPdgDecayMuon) {
+        int igmother = track.getMotherTrackId();
+        auto gmTrack = (*tracks)[igmother];
+        int gmpdg = gmTrack.GetPdgCode();
+        if (int(std::abs(gmpdg) / 100.) == 4 ||
+            int(std::abs(gmpdg) / 1000.) == 4) {
+          nMuons++;
+          nmuonsev++;
+          if (-4.3 < y && y < -2.2) {
+            nMuonsInAcceptance++;
+            nmuonsevinacc++;
+          }
+        } // gmpdg
+
+      } // pdgdecay
+
+    } // loop track
+    // std::cout << "#muons per event: " << nmuonsev << "\n";
+    // std::cout << "#muons in acceptance per event: " << nmuonsev << "\n";
+  } // events
+
+  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";
+  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 + nEventsInjTwo) < nEvents * ratioTrigger * 0.95 ||
+      (nEventsInjOne + nEventsInjTwo) > nEvents * ratioTrigger * 1.05) {
+    std::cerr << "Number of generated events injected with " << checkPdgQuarkOne << " and " << checkPdgQuarkTwo
+              << " different than expected\n";
+    return 1;
+  }
+  if (std::abs(float(nEventsInjOne) / nEventsInjTwo - ratioQuarkOneTwo) / ratioQuarkOneTwo > 0.1) {
+    std::cerr << "Ratio between quark " << checkPdgQuarkOne << " and quark " << checkPdgQuarkTwo
+              << " different than expected\n";
+    return 1;
+  }
+  std::cout << "#muons: " << nMuons << "\n";
+  std::cout << "#muons in acceptance: " << nMuonsInAcceptance << "\n";
+
+  return 0;
+} // external