diff --git a/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/BTOPSITOJPSIPIPI.DEC b/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/BTOPSITOJPSIPIPI.DEC new file mode 100644 index 000000000..82102bf60 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/MC/config/PWGDQ/EvtGen/GeneratorEvtGen.C b/MC/config/PWGDQ/EvtGen/GeneratorEvtGen.C index da0b1aa6a..18abd5ddf 100644 --- a/MC/config/PWGDQ/EvtGen/GeneratorEvtGen.C +++ b/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; }; diff --git a/MC/config/PWGDQ/external/generator/generator_pythia8_NonPromptSignals_gaptriggered_dq.C b/MC/config/PWGDQ/external/generator/generator_pythia8_NonPromptSignals_gaptriggered_dq.C index eb26e89f3..b2050175b 100644 --- a/MC/config/PWGDQ/external/generator/generator_pythia8_NonPromptSignals_gaptriggered_dq.C +++ b/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(); + 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(&random_value), sizeof(random_seed)); + gen->readString(Form("Random:seed = %d", random_value % 900000001)); + + // print debug + // gen->PrintDebug(); + + return gen; +} + diff --git a/MC/config/PWGDQ/ini/GeneratorHF_bbbar_PsiToJpsi_midy_triggerGap.ini b/MC/config/PWGDQ/ini/GeneratorHF_bbbar_PsiToJpsi_midy_triggerGap.ini new file mode 100644 index 000000000..a891cc8d1 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/MC/config/PWGDQ/ini/tests/GeneratorHF_bbbar_PsiToJpsi_midy_triggerGap.C b/MC/config/PWGDQ/ini/tests/GeneratorHF_bbbar_PsiToJpsi_midy_triggerGap.C new file mode 100644 index 000000000..5c50cd1af --- /dev/null +++ b/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* 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; iGetPdgCode()) == 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; +} \ No newline at end of file