Non-prmopt psi2S to Jpsi+pipi

MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/BTOPSITOJPSIPIPI.DEC

@@ -0,0 +1,167 @@
+####
+Decay B0
+#       B -> cc= s
+0.000310000 psi(2S) K_S0                                    SVS; #[Reconstructed PDG2011]
+0.000310000 psi(2S) K_L0                                    SVS; #[Reconstructed PDG2011]
+#
+#
+0.000610000 psi(2S) K*0                                     SVV_HELAMP PKHplus PKphHplus PKHzero PKphHzero PKHminus PKphHminus; #[Reconstructed PDG2011]
+
+0.0004     psi(2S)  K+  pi-          PHSP;
+0.0002     psi(2S)  K0  pi0          PHSP;
+0.0002     psi(2S)  K0  pi-  pi+     PHSP;
+0.0001     psi(2S)  K0  pi0  pi0     PHSP;
+0.0001     psi(2S)  K+  pi-  pi0     PHSP;
+0.0004     psi(2S)  K_10              PHSP;
+#
+####
+0.000620000 psi(2S) K0                                      PHSP;  #[New mode added] #[Reconstructed PDG2011]
+
+Enddecay
+
+Decay anti-B0
+#
+0.000310000 psi(2S) K_S0                                    SVS; #[Reconstructed PDG2011]
+0.000310000 psi(2S) K_L0                                    SVS; #[Reconstructed PDG2011]
+#
+#
+0.000610000 psi(2S) anti-K*0                                SVV_HELAMP PKHminus PKphHminus PKHzero PKphHzero PKHplus PKphHplus; #[Reconstructed PDG2011]
+
+0.0004     psi(2S)  K-  pi+          PHSP;
+0.0002     psi(2S)  anti-K0   pi0          PHSP;
+0.0002     psi(2S)  anti-K0   pi+  pi-     PHSP;
+0.0001     psi(2S)  anti-K0   pi0  pi0     PHSP;
+0.0001     psi(2S)  K-  pi+  pi0     PHSP;
+0.0004     psi(2S)  anti-K_10              PHSP;
+#
+0.000620000 psi(2S) anti-K0                                 PHSP;  #[New mode added] #[Reconstructed PDG2011]
+####
+Enddecay
+
+Decay B-
+#
+#   B -> cc= s          sum = 1.92%
+#
+0.000646000 psi(2S) K-                                      SVS; #[Reconstructed PDG2011]
+0.000620000 psi(2S) K*-                                     SVV_HELAMP PKHminus PKphHminus PKHzero PKphHzero PKHplus PKphHplus; #[Reconstructed PDG2011]
+0.0004   psi(2S) anti-K0   pi-                   PHSP;
+0.0002   psi(2S) K-  pi0                   PHSP;
+0.001900000 psi(2S) K-      pi+     pi-                     PHSP; #[Reconstructed PDG2011]
+0.0001   psi(2S) K-  pi0  pi0              PHSP;
+0.0001   psi(2S) anti-K0   pi-  pi0              PHSP;
+0.0004     psi(2S)  K_1-              PHSP;
+#
+0.000025800 psi(2S) pi-                                     PHSP;  #[New mode added] #[Reconstructed PDG2011]
+
+Enddecay
+
+Decay B+
+#
+#   B -> cc= s          sum = 1.92%
+#
+0.000646000 psi(2S) K+                                      SVS; #[Reconstructed PDG2011]
+0.000620000 psi(2S) K*+                                     SVV_HELAMP PKHplus PKphHplus PKHzero PKphHzero PKHminus PKphHminus; #[Reconstructed PDG2011]
+0.0004   psi(2S) K0  pi+                   PHSP;
+0.0002   psi(2S) K+  pi0                   PHSP;
+0.001900000 psi(2S) K+      pi-     pi+                     PHSP; #[Reconstructed PDG2011]
+0.0001   psi(2S) K+  pi0  pi0              PHSP;
+0.0001   psi(2S) K0  pi+  pi0              PHSP;
+0.0004   psi(2S)  K_1+                     PHSP;
+#
+####
+0.000025800 psi(2S) pi+                                     PHSP;  #[New mode added] #[Reconstructed PDG2011]
+
+Enddecay
+
+Decay B_s0
+#					psi' = 0.34%  CLNS 94/1315
+0.000465    psi(2S)     eta'	                           SVS;
+0.000235    psi(2S)     eta                                SVS;
+0.000680000 psi(2S) phi                                     SVV_HELAMP  1.0 0.0 1.0 0.0 1.0 0.0; #[Reconstructed PDG2011]
+0.0003      psi(2S)     K-          K+                     PHSP;
+0.0003      psi(2S)     anti-K0     K0                     PHSP;
+0.0003      psi(2S)     K0          K-         pi+         PHSP;
+0.0003      psi(2S)     anti-K0     K0         pi0         PHSP;
+0.0003      psi(2S)     K-          K+         pi0         PHSP;
+0.00034   psi(2S) phi      pi+  pi-  PHSP;
+0.00034   psi(2S) phi      pi0  pi0  PHSP;
+0.0002    psi(2S) eta      pi+  pi-  PHSP;
+0.0002    psi(2S) eta      pi0  pi0  PHSP;
+0.0004    psi(2S) eta'     pi+  pi-  PHSP;
+0.0004    psi(2S) eta'     pi0  pi0  PHSP;
+0.0002    psi(2S) pi+      pi-       PHSP;
+0.0002    psi(2S) pi0      pi0       PHSP;
+#
+# PR LHCb 04/08/2004 : add Bs -> phi mu mu, phi e e
+0.0000023   phi    e+     e-                   BTOSLLALI;
+Enddecay
+
+Decay anti-B_s0
+#					psi' = 0.34%  CLNS 94/1315
+#
+0.000465  psi(2S) eta'               SVS;
+0.000235  psi(2S) eta                SVS;
+0.000680000 psi(2S) phi                                     SVV_HELAMP  1.0 0.0 1.0 0.0 1.0 0.0; #[Reconstructed PDG2011]
+0.0003    psi(2S) K-       K+        PHSP;
+0.0003    psi(2S) anti-K0  K0        PHSP;
+0.0003    psi(2S) anti-K0  K+    pi- PHSP;
+0.0003    psi(2S) anti-K0  K0    pi0 PHSP;
+0.0003    psi(2S) K-       K+    pi0 PHSP;
+0.00034   psi(2S) phi      pi+  pi-  PHSP;
+0.00034   psi(2S) phi      pi0  pi0  PHSP;
+0.0002    psi(2S) eta      pi+  pi-  PHSP;
+0.0002    psi(2S) eta      pi0  pi0  PHSP;
+0.0004    psi(2S) eta'     pi+  pi-  PHSP;
+0.0004    psi(2S) eta'     pi0  pi0  PHSP;
+0.0002    psi(2S) pi+      pi-       PHSP;
+0.0002    psi(2S) pi0      pi0       PHSP;
+#
+# PR LHCb 04/08/2004 : add Bs -> phi mu mu, phi e e
+0.0000023   phi    e-     e+                   BTOSLLALI;
+
+Enddecay
+
+Decay Lambda_b0
+  0.00038    Lambda0         psi(2S)                      PHSP;
+Enddecay
+
+Decay anti-Lambda_b0
+  0.00038    anti-Lambda0    psi(2S)                      PHSP;
+Enddecay
+
+Decay Xi_b-
+  0.00038    Xi-        psi(2S)                           PHSP;
+Enddecay
+
+Decay anti-Xi_b+
+  0.00038    anti-Xi+        psi(2S)                      PHSP;
+Enddecay
+
+Decay Omega_b-
+  0.00038    Omega-        psi(2S)                        PHSP;
+Enddecay
+
+Decay anti-Omega_b+
+  0.00038    anti-Omega+        psi(2S)                   PHSP;
+Enddecay
+
+Decay B_c-
+#                  SemiLeptonic Decays
+0.00094    psi(2S)    e-      anti-nu_e         PHOTOS  PHSP;
+Enddecay
+
+Decay B_c+
+#                  SemiLeptonic Decays
+0.00094    psi(2S)    e+          nu_e         PHOTOS  PHSP;
+Enddecay
+
+Decay psi(2S)
+1.000 J/psi   pi+     pi-         VVPIPI;
+Enddecay
+
+Decay J/psi
+### from DECAY.DEC
+1.000 e+   e-                PHOTOS VLL;
+Enddecay
+
+End

MC/config/PWGDQ/EvtGen/GeneratorEvtGen.C

@@ -42,6 +42,7 @@ enum DecayModeEvt { kEvtAll = 0,
                     kEvtElectronEM,
                     kEvtDiElectronEM,
                     kEvtGammaEM,
+                    kEvtBtoPsi2SToJpsiPiPi,
                     kEvtBeautyUpgrade };
 
 namespace o2
@@ -355,6 +356,9 @@ class GeneratorEvtGen : public T
       case kEvtBeautyUpgrade:
         SetDecayTable(Form("%s/BEAUTYUPGRADE.DEC", pathO2.Data()));
         break;
+      case kEvtBtoPsi2SToJpsiPiPi:
+        SetDecayTable(Form("%s/BTOPSITOJPSIPIPI.DEC", pathO2.Data()));
+        break;
     }
     return;
   };

MC/config/PWGDQ/external/generator/generator_pythia8_NonPromptSignals_gaptriggered_dq.C

@@ -351,3 +351,41 @@ FairGenerator*
   return gen;
 }
 
+FairGenerator*
+  GeneratorBeautyToPsiToJpsi_EvtGenMidY(double rapidityMin = -1.5, double rapidityMax = 1.5, bool verbose = false, TString pdgs = "511;521;531;541;5112;5122;5232;5132;5332")
+{
+  auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::GeneratorPythia8NonPromptInjectedGapTriggeredDQ>();
+  gen->setRapidity(rapidityMin, rapidityMax);
+  gen->setPDG(5);
+  gen->setRapidityHadron(rapidityMin,rapidityMax);
+  gen->setHadronMultiplicity(1);
+  TString pathO2table = gSystem->ExpandPathName("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/decayer/switchOffBhadrons.cfg");
+  gen->readFile(pathO2table.Data());
+  gen->setConfigMBdecays(pathO2table);
+  gen->setVerbose(verbose);
+
+  std::string spdg;
+  TObjArray* obj = pdgs.Tokenize(";");
+  gen->SetSizePdg(obj->GetEntriesFast());
+  for (int i = 0; i < obj->GetEntriesFast(); i++) {
+    spdg = obj->At(i)->GetName();
+    gen->AddPdg(std::stoi(spdg), i);
+    printf("PDG %d \n", std::stoi(spdg));
+    gen->addHadronPDGs(std::stoi(spdg));
+  }
+  gen->SetForceDecay(kEvtBtoPsi2SToJpsiPiPi);
+
+  // set random seed
+  gen->readString("Random:setSeed on");
+  uint random_seed;
+  unsigned long long int random_value = 0;
+  ifstream urandom("/dev/urandom", ios::in|ios::binary);
+  urandom.read(reinterpret_cast<char*>(&random_value), sizeof(random_seed));
+  gen->readString(Form("Random:seed = %d", random_value % 900000001));
+
+  // print debug
+  // gen->PrintDebug();
+
+  return gen;
+}
+

MC/config/PWGDQ/ini/GeneratorHF_bbbar_PsiToJpsi_midy_triggerGap.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/generator_pythia8_NonPromptSignals_gaptriggered_dq.C 
+funcName = GeneratorBeautyToPsiToJpsi_EvtGenMidY()
+
+### 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/common/pythia8/generator/pythia8_hf.cfg
+hooksFileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/hooks/pythia8_userhooks_qqbar.C
+hooksFuncName = pythia8_userhooks_bbbar(-1.5,1.5)

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

@@ -0,0 +1,104 @@
+int External()
+{
+    int checkPdgSignal[] = {100443};
+    int checkPdgDecay[] = {443, 211, -211}; 
+    int leptonPdg = 11;
+    Double_t rapidityWindow = 1.0;
+    std::string path{"o2sim_Kine.root"};
+    std::cout << "Check for\nsignal PDG " << checkPdgSignal[0] << "\n decay PDG " << checkPdgDecay[0] << ", " << checkPdgDecay[1] << ", " << checkPdgDecay[2] << "\n";
+    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 nAntileptons{};
+    int nLeptonPairs{};
+    int nLeptonPairsToBeDone{};
+    int nSignalJpsi{};
+    int nSignalPionsPos{};
+    int nSignalPionsNeg{};
+    int nSignalPsi2S{};
+    int nSignalJpsiWithinAcc{};
+    int nSignalPionsPosWithinAcc{};
+    int nSignalPionsNegWithinAcc{};
+    auto nEvents = tree->GetEntries();
+    o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine);
+    Int_t  bpdgs[] = {511, 521, 531, 5112, 5122, 5232, 5132};
+    Int_t sizePdg = sizeof(bpdgs)/sizeof(Int_t);
+    Bool_t hasBeautyMoth = kFALSE;
+
+    for (int i = 0; i < nEvents; i++) {
+        tree->GetEntry(i);
+        for (auto& track : *tracks) {
+            auto pdg = track.GetPdgCode();
+            auto rapidity =  track.GetRapidity();
+            auto idMoth = track.getMotherTrackId();
+            int idJpsi = -1; int IdChild0 = -1; int IdChild1 = -1;
+            if (pdg == leptonPdg) {
+                // count leptons
+                nLeptons++;
+            } else if(pdg == -leptonPdg) {
+                // count anti-leptons
+                nAntileptons++;
+            } else if (pdg == checkPdgSignal[0]) {
+                hasBeautyMoth  = kFALSE;
+                if(idMoth){ //  check beauty mother
+                    auto tdM = mcreader.getTrack(i, idMoth);
+                    for(int i=0; i<sizePdg; i++){ if (TMath::Abs(tdM->GetPdgCode()) == bpdgs[i] ) hasBeautyMoth = kTRUE; }
+                }
+                if(hasBeautyMoth){
+                    nSignalPsi2S++;
+                    for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
+                        auto pdgDau = tracks->at(j).GetPdgCode();
+                        if(TMath::Abs(pdgDau) == checkPdgDecay[0] ) { nSignalJpsi++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalJpsiWithinAcc++; idJpsi = j; }
+                        if(pdgDau == checkPdgDecay[1] ) { nSignalPionsPos++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalPionsPosWithinAcc++; }
+                        if(pdgDau == checkPdgDecay[2] ) { nSignalPionsNeg++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalPionsNegWithinAcc++; }
+                    }
+
+                    auto trackJpsi = tracks->at(idJpsi);
+                    for (int j{trackJpsi.getFirstDaughterTrackId()}; j <= trackJpsi.getLastDaughterTrackId(); ++j) {
+                        auto pdgDau = tracks->at(j).GetPdgCode();
+                        if(pdgDau == leptonPdg ) IdChild0 = j;
+                        if(pdgDau == -leptonPdg ) IdChild1 = j;
+                    }
+                    auto child0 = tracks->at(IdChild0);
+                    auto child1 = tracks->at(IdChild1);
+                    // check for parent-child relations
+                    auto pdg0 = child0.GetPdgCode();
+                    auto pdg1 = child1.GetPdgCode();
+                    if (std::abs(pdg0) == leptonPdg && std::abs(pdg1) == leptonPdg && pdg0 == -pdg1) {
+                        nLeptonPairs++;
+                        if (child0.getToBeDone() && child1.getToBeDone()) {
+                            nLeptonPairsToBeDone++;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    std::cout << "#events: " << nEvents << "\n"
+              << "#leptons: " << nLeptons << "\n"
+              << "#antileptons: " << nAntileptons << "\n"
+              << "#signal (jpsi <- psi2S): " << nSignalJpsi << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalJpsiWithinAcc << "\n"
+              << "#signal (pi+ <- psi2S): " << nSignalPionsPos << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalPionsPosWithinAcc << "\n"
+              << "#signal (pi- <- psi2S): " << nSignalPionsNeg << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalPionsNegWithinAcc << "\n"
+              << "#lepton pairs: " << nLeptonPairs << "\n"
+              << "#lepton pairs to be done: " << nLeptonPairs << "\n";
+
+
+    if (nLeptonPairs == 0 || nLeptons == 0 || nAntileptons == 0) {
+        std::cerr << "Number of leptons, number of anti-leptons as well as number of lepton pairs should all be greater than 1.\n";
+        return 1;
+    }
+    if (nLeptonPairs != nLeptonPairsToBeDone) {
+        std::cerr << "The number of lepton pairs should be the same as the number of lepton pairs which should be transported.\n";
+        return 1;
+    }
+
+    return 0;
+}