Add cc->mumu with natural decay

Author: motomioya

MC/config/PWGDQ/ini/GeneratorHF_ccbarToMuonsSemileptonic_fwdy_inel_TriggerGap_natural.ini

@@ -0,0 +1,14 @@
+### The setup uses an external event generator
+### This part sets the path of the file and the function call to retrieve it
+
+[GeneratorExternal]
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/GeneratorHFToMu_EvtGenFwdY_gaptriggered_dq.C 
+funcName = GeneratorHFToMu_EvtGenFwdY_gaptriggered_dq(-4.3, -2.3, false, false, false)
+
+### The external generator derives from GeneratorPythia8.
+### This part configures the bits of the interface: configuration and user hooks
+
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_cr2.cfg
+hooksFileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/hooks/pythia8_userhooks_qqbar.C
+hooksFuncName = pythia8_userhooks_ccbar(-4.3, -2.3)

MC/config/PWGDQ/ini/tests/GeneratorHF_ccbarToMuonsSemileptonic_fwdy_inel_TriggerGap_natural.C

@@ -0,0 +1,77 @@
+int External()
+{
+  int checkPdgDecay = 13;
+  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");
+  std::vector<o2::MCTrack>* tracks{};
+  tree->SetBranchAddress("MCTrack", &tracks);
+
+  int nLeptons{};
+  int nLeptonsInAcceptance{};
+  int nLeptonsToBeDone{};
+  int nSignalPairs{};
+  int nLeptonPairs{};
+  int nLeptonPairsInAcceptance{};
+  int nLeptonPairsToBeDone{};
+  auto nEvents = tree->GetEntries();
+
+  for (int i = 0; i < nEvents; i++) {
+    tree->GetEntry(i);
+    int nleptonseinacc = 0;
+    int nleptonse = 0;
+    int nleptonseToBeDone = 0;
+    int nopenHeavy = 0;
+    for (auto& track : *tracks) {
+      auto pdg = track.GetPdgCode();
+      auto y = track.GetRapidity();
+      if (std::abs(pdg) == checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if ( int(std::abs(gmpdg)/100.) == 4 || int(std::abs(gmpdg)/1000.) == 4 || int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5 ) {
+            nLeptons++;
+            nleptonse++;
+            if (-4.3 < y && y < -2.2) {
+              nleptonseinacc++;
+              nLeptonsInAcceptance++;
+            }
+            if (track.getToBeDone()) {
+              nLeptonsToBeDone++;
+              nleptonseToBeDone++;
+            }
+          }
+        }    
+      } else if (std::abs(pdg) == 411 || std::abs(pdg) == 421 || std::abs(pdg) == 431 || std::abs(pdg) == 4122 || std::abs(pdg) == 4132 || std::abs(pdg) == 4232 || std::abs(pdg) == 4332) {
+        nopenHeavy++;
+      }
+    }
+    if (nopenHeavy > 1) nSignalPairs++;
+    if (nleptonse > 1) nLeptonPairs++;
+    if (nleptonseToBeDone > 1) nLeptonPairsToBeDone++;
+    if (nleptonseinacc > 1) nLeptonPairsInAcceptance++;
+  }
+  std::cout << "#events: " << nEvents << "\n"
+  << "#muons in acceptance: " << nLeptonsInAcceptance << "\n"
+  << "#muon pairs in acceptance: " << nLeptonPairsInAcceptance << "\n"
+  << "#muons: " << nLeptons << "\n"
+  << "#muons to be done: " << nLeptonsToBeDone << "\n"
+  << "#signal pairs: " << nSignalPairs << "\n"
+  << "#muon pairs: " << nLeptonPairs << "\n"
+  << "#muon pairs to be done: " << nLeptonPairsToBeDone << "\n";
+  if (nLeptonPairs != nLeptonPairsToBeDone) {
+    std::cerr << "The number of muon pairs should be the same as the number of muon pairs which should be transported.\n";
+    return 1;
+  }
+  if (nLeptons != nLeptonsToBeDone) {
+    std::cerr << "The number of muons should be the same as the number of muons which should be transported.\n";
+    return 1;
+  }
+
+  return 0;
+}