diff --git a/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/PSITOJPSIPIPI.DEC b/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/PSITOJPSIPIPI.DEC new file mode 100644 index 000000000..70bc72517 --- /dev/null +++ b/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/PSITOJPSIPIPI.DEC @@ -0,0 +1,8 @@ +Decay J/psi +### from DECAY.DEC +1.000 e+ e- PHOTOS VLL; +Enddecay +Decay psi(2S) +1.000 J/psi pi+ pi- VVPIPI; +Enddecay +End diff --git a/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/X3872TOJPSIPIPI.DEC b/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/X3872TOJPSIPIPI.DEC new file mode 100644 index 000000000..8451b0777 --- /dev/null +++ b/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/X3872TOJPSIPIPI.DEC @@ -0,0 +1,10 @@ +Decay X_1(3872) +#### Breit-Wigner function for the pion distribution (S-Wave approximation) +1.000 J/psi pi+ pi- XJPSIRHO0PIPI_S; +Enddecay +Decay J/psi +### from DECAY.DEC +1.000 e+ e- PHOTOS VLL; +Enddecay + +End diff --git a/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C b/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C index 2be7d7d8b..61bbf2c26 100644 --- a/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C +++ b/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C @@ -638,6 +638,75 @@ class O2_GeneratorParamPsiPbPb5TeV : public GeneratorTGenerator GeneratorParam* paramPsi = nullptr; }; + +class O2_GeneratorParamX3872MidY : public GeneratorTGenerator +{ + + public: + O2_GeneratorParamX3872MidY() : GeneratorTGenerator("ParamX3872MidY") + { + paramX3872 = new GeneratorParam(1, -1, PtX3872pp13TeV, YX3872pp13TeV, V2X3872pp13TeV, IpX3872pp13TeV); + paramX3872->SetMomentumRange(0., 1.e6); + paramX3872->SetPtRange(0., 1000.); + paramX3872->SetYRange(-1.0, 1.0); + paramX3872->SetPhiRange(0., 360.); + paramX3872->SetDecayer(new TPythia6Decayer()); // Pythia + paramX3872->SetForceDecay(kNoDecay); // particle left undecayed + setTGenerator(paramX3872); + }; + + ~O2_GeneratorParamX3872MidY() + { + delete paramX3872; + }; + + Bool_t Init() override + { + GeneratorTGenerator::Init(); + paramX3872->Init(); + return true; + } + + void SetNSignalPerEvent(Int_t nsig) { paramX3872->SetNumberParticles(nsig); } + + //-------------------------------------------------------------------------// + static Double_t PtX3872pp13TeV(const Double_t* px, const Double_t* /*dummy*/) + { + // prompt X3872 pT + // pp, 13TeV (tuned LHCb pp 13 TeV) + // + const Double_t kC = 7.64519e+00 ; + const Double_t kpt0 = 5.30628e+00; + const Double_t kn = 3.30887e+00; + Double_t pt = px[0]; + return kC * pt / TMath::Power((1. + (pt / kpt0) * (pt / kpt0)), kn); + } + + //-------------------------------------------------------------------------// + static Double_t YX3872pp13TeV(const Double_t* py, const Double_t* /*dummy*/) + { + // flat rapidity distribution assumed at midrapidity + return 1.; + } + + //-------------------------------------------------------------------------// + static Double_t V2X3872pp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/) + { + // X3872 v2 + return 0.; + } + + //-------------------------------------------------------------------------// + static Int_t IpX3872pp13TeV(TRandom*) + { + return 9920443; + } + + private: + GeneratorParam* paramX3872 = nullptr; +}; + + } // namespace eventgen } // namespace o2 @@ -741,6 +810,31 @@ FairGenerator* GeneratorCocktailPromptCharmoniaToMuonEvtGen_pp13TeV() return genCocktailEvtGen; } + +FairGenerator* + GeneratorParamPromptPsiToJpsiPiPiEvtGen_pp13TeV(TString pdgs = "100443") +{ + auto gen = new o2::eventgen::GeneratorEvtGen(); + gen->SetNSignalPerEvent(1); // number of jpsis per event + + 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)); + } + TString pathO2 = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen"); + gen->SetDecayTable(Form("%s/PSITOJPSIPIPI.DEC", pathO2.Data())); + + // print debug + gen->PrintDebug(); + + return gen; +} + + FairGenerator* GeneratorParamPromptJpsiToMuonEvtGen_pp13TeV(TString pdgs = "443") { @@ -844,6 +938,29 @@ FairGenerator* } +FairGenerator* + GeneratorParamX3872ToJpsiEvtGen_pp13TeV(TString pdgs = "9920443") +{ + auto gen = new o2::eventgen::GeneratorEvtGen(); + gen->SetNSignalPerEvent(1); // number of jpsis per event + + 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)); + } + TString pathO2 = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen"); + gen->SetDecayTable(Form("%s/X3872TOJPSIPIPI.DEC", pathO2.Data())); + + // print debug + gen->PrintDebug(); + + return gen; +} + diff --git a/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C b/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C index f0910dcad..ce6c7a7a5 100644 --- a/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C +++ b/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C @@ -45,7 +45,12 @@ public: case 7: // generate prompt charmonia cocktail at forward rapidity at 5TeV mGeneratorParam = (Generator*)GeneratorCocktailPromptCharmoniaToMuonEvtGen_PbPb5TeV(); break; - + case 8: // generate prompt X_1(3872) to Jpsi pi pi at midrapidity + mGeneratorParam = (Generator*)GeneratorParamX3872ToJpsiEvtGen_pp13TeV("9920443"); + break; + case 9: // generate prompt psi2S to Jpsi pi pi at midrapidity + mGeneratorParam = (Generator*)GeneratorParamPromptPsiToJpsiPiPiEvtGen_pp13TeV("100443"); + break; } mGeneratorParam->Init(); } diff --git a/MC/config/PWGDQ/ini/Generator_InjectedPromptPsi2SToJpsiPiPiMidy_TriggerGap.ini b/MC/config/PWGDQ/ini/Generator_InjectedPromptPsi2SToJpsiPiPiMidy_TriggerGap.ini new file mode 100644 index 000000000..30082ec0f --- /dev/null +++ b/MC/config/PWGDQ/ini/Generator_InjectedPromptPsi2SToJpsiPiPiMidy_TriggerGap.ini @@ -0,0 +1,6 @@ +### The external generator derives from GeneratorPythia8. +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C +funcName=GeneratorPythia8InjectedPromptCharmoniaGapTriggered(5,9) +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg diff --git a/MC/config/PWGDQ/ini/Generator_InjectedPromptX3782Midy_TriggerGap.ini b/MC/config/PWGDQ/ini/Generator_InjectedPromptX3782Midy_TriggerGap.ini new file mode 100644 index 000000000..097d4377d --- /dev/null +++ b/MC/config/PWGDQ/ini/Generator_InjectedPromptX3782Midy_TriggerGap.ini @@ -0,0 +1,6 @@ +### The external generator derives from GeneratorPythia8. +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C +funcName=GeneratorPythia8InjectedPromptCharmoniaGapTriggered(5,8) +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg diff --git a/MC/config/PWGDQ/ini/tests/Generator_InjectedPromptPsi2SToJpsiPiPiMidy_TriggerGap.C b/MC/config/PWGDQ/ini/tests/Generator_InjectedPromptPsi2SToJpsiPiPiMidy_TriggerGap.C new file mode 100644 index 000000000..9254fdab1 --- /dev/null +++ b/MC/config/PWGDQ/ini/tests/Generator_InjectedPromptPsi2SToJpsiPiPiMidy_TriggerGap.C @@ -0,0 +1,99 @@ +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); + Bool_t hasPsi2SMoth = 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]) { + // check daughters + std::cout << "Signal PDG: " << pdg << "\n"; + for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) { + auto pdgDau = tracks->at(j).GetPdgCode(); + std::cout << "Daughter " << j << " is: " << pdgDau << "\n"; + 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(); + std::cout << "Lepton daughter particles of mother " << trackJpsi.GetPdgCode() << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n"; + 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; +} diff --git a/MC/config/PWGDQ/ini/tests/Generator_InjectedPromptX3782Midy_TriggerGap.C b/MC/config/PWGDQ/ini/tests/Generator_InjectedPromptX3782Midy_TriggerGap.C new file mode 100644 index 000000000..836c8df30 --- /dev/null +++ b/MC/config/PWGDQ/ini/tests/Generator_InjectedPromptX3782Midy_TriggerGap.C @@ -0,0 +1,99 @@ +int External() +{ + int checkPdgSignal[] = {9920443}; // pdg code X3872 + 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 nSignalX3872{}; + int nSignalPionsPos{}; + int nSignalPionsNeg{}; + int nSignalPsi2S{}; + int nSignalX3872WithinAcc{}; + int nSignalPionsPosWithinAcc{}; + int nSignalPionsNegWithinAcc{}; + auto nEvents = tree->GetEntries(); + o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine); + Bool_t hasPsi2SMoth = 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 idX3872 = -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]) { + // check daughters + std::cout << "Signal PDG: " << pdg << "\n"; + for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) { + auto pdgDau = tracks->at(j).GetPdgCode(); + std::cout << "Daughter " << j << " is: " << pdgDau << "\n"; + if(TMath::Abs(pdgDau) == checkPdgDecay[0] ) { nSignalX3872++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalX3872WithinAcc++; idX3872 = 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 trackX3872 = tracks->at(idX3872); + for (int j{trackX3872.getFirstDaughterTrackId()}; j <= trackX3872.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(); + std::cout << "Lepton daughter particles of mother " << trackX3872.GetPdgCode() << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n"; + 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 <- X3872): " << nSignalX3872 << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalX3872WithinAcc << "\n" + << "#signal (pi+ <- X3872): " << nSignalPionsPos << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalPionsPosWithinAcc << "\n" + << "#signal (pi- <- X3872): " << 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; +}