PWGHF: add semimuonic process with acceptance -4.3 to -2.2

MC/config/PWGHF/ini/GeneratorHF_mu_bbbar_gap5_Mode2_accSmall.ini

@@ -0,0 +1,7 @@
+### The external generator derives from GeneratorPythia8.
+[GeneratorExternal]
+fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
+funcName=GeneratorPythia8GapTriggeredBeauty(5, -4.3, -2.2)
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_beauty_with_mudecays_Mode2.cfg
+includePartonEvent=true

MC/config/PWGHF/ini/GeneratorHF_mu_ccbar_gap5_Mode2_accSmall.ini

@@ -0,0 +1,7 @@
+### The external generator derives from GeneratorPythia8.
+[GeneratorExternal]
+fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
+funcName=GeneratorPythia8GapTriggeredCharm(5, -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_mu_bbbar_gap5_Mode2_accSmall.C

@@ -0,0 +1,103 @@
+int External() {
+
+  int checkPdgDecayMuon = 13;
+  int checkPdgQuark = 5;
+
+  float ratioTrigger = 1. / 5; // 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{};
+  int nEventsInj{};
+  int nQuarks{};
+  int nMuons{};
+
+  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 == checkPdgQuark) {
+        nEventsInj++;
+      }
+    } // if event header
+
+    int nmuonsev = 0;
+    int nmuonsevinacc = 0;
+
+    for (auto &track : *tracks) {
+      auto pdg = track.GetPdgCode();
+      if (std::abs(pdg) == checkPdgQuark) {
+        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.) == 5 ||
+            int(std::abs(gmpdg) / 1000.) == 5) {
+          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: ", checkPdgQuark)
+            << nEventsInj << "\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;
+  }
+  std::cout << "#muons: " << nMuons << "\n";
+  std::cout << "#muons in acceptance: " << nMuonsInAcceptance << "\n";
+
+  return 0;
+} // external

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

@@ -0,0 +1,103 @@
+int External() {
+
+  int checkPdgDecayMuon = 13;
+  int checkPdgQuark = 4;
+
+  float ratioTrigger = 1. / 5; // 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{};
+  int nEventsInj{};
+  int nQuarks{};
+  int nMuons{};
+
+  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 == checkPdgQuark) {
+        nEventsInj++;
+      }
+    } // if event header
+
+    int nmuonsev = 0;
+    int nmuonsevinacc = 0;
+
+    for (auto &track : *tracks) {
+      auto pdg = track.GetPdgCode();
+      if (std::abs(pdg) == checkPdgQuark) {
+        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.) == 5 ||
+            int(std::abs(gmpdg) / 1000.) == 5) {
+          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: ", checkPdgQuark)
+            << nEventsInj << "\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;
+  }
+  std::cout << "#muons: " << nMuons << "\n";
+  std::cout << "#muons in acceptance: " << nMuonsInAcceptance << "\n";
+
+  return 0;
+} // external

MC/config/PWGHF/pythia8/generator/pythia8_beauty_with_mudecays_Mode2.cfg

@@ -0,0 +1,71 @@
+### authors: Fabrizio Grosa (fabrizio.grosa@cern.ch)
+###          Cristina Terrevoli (cristina.terrevoli@cern.ch)
+###          Fabio Catalano (fabio.catalano@cern.ch)
+###          Nicole Bastid (nicole.bastid@cern.ch)
+### Beauty decay muons 
+### last update:  November 2024
+
+### beams
+Beams:idA 2212			# proton
+Beams:idB 2212 			# proton
+Beams:eCM 13600. 		# GeV
+
+### processes
+SoftQCD:inelastic on		# all inelastic processes
+
+### decays
+ParticleDecays:limitTau0 on	
+ParticleDecays:tau0Max 10.
+
+### switching on Pythia Mode2
+ColourReconnection:mode 1
+ColourReconnection:allowDoubleJunRem off
+ColourReconnection:m0 0.3
+ColourReconnection:allowJunctions on
+ColourReconnection:junctionCorrection 1.20
+ColourReconnection:timeDilationMode 2
+ColourReconnection:timeDilationPar 0.18
+StringPT:sigma 0.335
+StringZ:aLund 0.36
+StringZ:bLund 0.56
+StringFlav:probQQtoQ 0.078
+StringFlav:ProbStoUD 0.2
+StringFlav:probQQ1toQQ0join 0.0275,0.0275,0.0275,0.0275
+MultiPartonInteractions:pT0Ref 2.15
+BeamRemnants:remnantMode 1
+BeamRemnants:saturation 5
+
+# Correct decay lengths (wrong in PYTHIA8 decay table)
+# Lb
+5122:tau0 = 0.4390
+# Xic0
+4132:tau0 = 0.0455
+# OmegaC
+4332:tau0 = 0.0803
+
+
+### switch off all decay channels for charm
+511:onMode = off
+521:onMode = off
+531:onMode = off
+5122:onMode = off
+5132:onMode = off
+5232:onMode = off
+5332:onMode = off
+
+###Semimuonic decays of charm
+
+### B0 -> mu X
+511:onIfAny = 13
+### B+/- -> mu X
+521:onIfAny = 13
+### B_s -> mu X
+531:onIfAny = 13
+### Lambda_b -> mu X
+5122:onIfAny = 13
+### Xsi_b -> mu X
+5132:onIfAny = 13
+### Xsi0_b -> mu X
+5232:onIfAny = 13
+### Omega_b -> mu X
+5332:onIfAny = 13

MC/config/PWGHF/pythia8/generator/pythia8_charm_with_mudecays_Mode2.cfg

@@ -0,0 +1,70 @@
+### authors: Fabrizio Grosa (fabrizio.grosa@cern.ch)
+###          Cristina Terrevoli (cristina.terrevoli@cern.ch)
+###          Fabio Catalano (fabio.catalano@cern.ch)
+###          Nicole Bastid (nicole.bastid@cern.ch)
+### charm decay muons 
+### last update:  November 2024
+
+### beams
+Beams:idA 2212			# proton
+Beams:idB 2212 			# proton
+Beams:eCM 13600. 		# GeV
+
+### processes
+SoftQCD:inelastic on		# all inelastic processes
+
+### decays
+ParticleDecays:limitTau0 on	
+ParticleDecays:tau0Max 10.
+
+### switching on Pythia Mode2
+ColourReconnection:mode 1
+ColourReconnection:allowDoubleJunRem off
+ColourReconnection:m0 0.3
+ColourReconnection:allowJunctions on
+ColourReconnection:junctionCorrection 1.20
+ColourReconnection:timeDilationMode 2
+ColourReconnection:timeDilationPar 0.18
+StringPT:sigma 0.335
+StringZ:aLund 0.36
+StringZ:bLund 0.56
+StringFlav:probQQtoQ 0.078
+StringFlav:ProbStoUD 0.2
+StringFlav:probQQ1toQQ0join 0.0275,0.0275,0.0275,0.0275
+MultiPartonInteractions:pT0Ref 2.15
+BeamRemnants:remnantMode 1
+BeamRemnants:saturation 5
+
+# Correct decay lengths (wrong in PYTHIA8 decay table)
+# Lb
+5122:tau0 = 0.4390
+# Xic0
+4132:tau0 = 0.0455
+# OmegaC
+4332:tau0 = 0.0803
+
+### switch off all decay channels for charm
+411:onMode = off
+421:onMode = off
+431:onMode = off
+4122:onMode = off
+4132:onMode = off
+4232:onMode = off
+4332:onMode = off
+
+###Semimuonic decays of charm
+
+### D+/- -> mu X
+411:onIfAny = 13
+### D0 -> mu X
+421:onIfAny = 13
+### D_s -> mu X
+431:onIfAny = 13
+### Lambda_c -> mu X
+4122:onIfAny = 13
+### Xsi0_c -> mu X
+4132:onIfAny = 13
+### Xsi+_c
+4232:onIfAny = 13
+### Omega_c -> mu X
+4332:onIfAny = 13