From 65e6c224eb4294b0a8a2cd393855d55a6bdcc500 Mon Sep 17 00:00:00 2001 From: ffionda Date: Sat, 21 Dec 2024 12:12:21 +0100 Subject: [PATCH] update for generating a cocktail of X3872 + Psi2S to Jpsi pi pi --- .../X3872ANDPSI2STOJPSIPIPI.DEC | 15 +++ .../generator/GeneratorPromptCharmonia.C | 25 +++++ ...ithInjectedPromptSignals_gaptriggered_dq.C | 3 + ...ctedPromptX3782andPsi2SMidy_TriggerGap.ini | 6 + ...jectedPromptX3782andPsi2SMidy_TriggerGap.C | 104 ++++++++++++++++++ 5 files changed, 153 insertions(+) create mode 100644 MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/X3872ANDPSI2STOJPSIPIPI.DEC create mode 100644 MC/config/PWGDQ/ini/Generator_InjectedPromptX3782andPsi2SMidy_TriggerGap.ini create mode 100644 MC/config/PWGDQ/ini/tests/Generator_InjectedPromptX3782andPsi2SMidy_TriggerGap.C diff --git a/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/X3872ANDPSI2STOJPSIPIPI.DEC b/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/X3872ANDPSI2STOJPSIPIPI.DEC new file mode 100644 index 000000000..53f19edc6 --- /dev/null +++ b/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/X3872ANDPSI2STOJPSIPIPI.DEC @@ -0,0 +1,15 @@ +Decay X_1(3872) +#### Breit-Wigner function for the pion distribution (S-Wave approximation) +1.000 J/psi pi+ pi- XJPSIRHO0PIPI_S; +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 diff --git a/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C b/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C index c7a62623c..840d450cf 100644 --- a/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C +++ b/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C @@ -962,5 +962,30 @@ FairGenerator* } +FairGenerator* GeneratorCocktailX3872AndPsi2StoJpsi_pp13TeV() +{ + auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen(); + + auto genX3872 = new o2::eventgen::O2_GeneratorParamX3872MidY; + genX3872->SetNSignalPerEvent(1); // number of jpsis per event + auto genPsi2S = new o2::eventgen::O2_GeneratorParamPsiMidY; + genPsi2S->SetNSignalPerEvent(1); // number of jpsis per event + genCocktailEvtGen->AddGenerator(genX3872, 1); // add J/psi generator + genCocktailEvtGen->AddGenerator(genPsi2S, 1); // add Psi(2S) generator + + TString pdgs = "9920443;100443"; + std::string spdg; + TObjArray* obj = pdgs.Tokenize(";"); + genCocktailEvtGen->SetSizePdg(obj->GetEntriesFast()); + for (int i = 0; i < obj->GetEntriesFast(); i++) { + spdg = obj->At(i)->GetName(); + genCocktailEvtGen->AddPdg(std::stoi(spdg), i); + printf("PDG %d \n", std::stoi(spdg)); + } + TString pathO2 = gSystem->ExpandPathName("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen"); + genCocktailEvtGen->SetDecayTable(Form("%s/X3872ANDPSI2STOJPSIPIPI.DEC", pathO2.Data())); + + return genCocktailEvtGen; +} 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 ce6c7a7a5..c0d6bfd85 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 @@ -51,6 +51,9 @@ public: case 9: // generate prompt psi2S to Jpsi pi pi at midrapidity mGeneratorParam = (Generator*)GeneratorParamPromptPsiToJpsiPiPiEvtGen_pp13TeV("100443"); break; + case 10: // generate cocktail of prompt X_1(3872) and psi2S to Jpsi pi pi at midrapidity + mGeneratorParam = (Generator*)GeneratorCocktailX3872AndPsi2StoJpsi_pp13TeV(); + break; } mGeneratorParam->Init(); } diff --git a/MC/config/PWGDQ/ini/Generator_InjectedPromptX3782andPsi2SMidy_TriggerGap.ini b/MC/config/PWGDQ/ini/Generator_InjectedPromptX3782andPsi2SMidy_TriggerGap.ini new file mode 100644 index 000000000..566002ac8 --- /dev/null +++ b/MC/config/PWGDQ/ini/Generator_InjectedPromptX3782andPsi2SMidy_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,10) +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg diff --git a/MC/config/PWGDQ/ini/tests/Generator_InjectedPromptX3782andPsi2SMidy_TriggerGap.C b/MC/config/PWGDQ/ini/tests/Generator_InjectedPromptX3782andPsi2SMidy_TriggerGap.C new file mode 100644 index 000000000..87674d2f1 --- /dev/null +++ b/MC/config/PWGDQ/ini/tests/Generator_InjectedPromptX3782andPsi2SMidy_TriggerGap.C @@ -0,0 +1,104 @@ +int External() +{ + int checkPdgSignal[] = {9920443,100443}; // pdg code X3872 + TString PdgSignalName[] = {"X(3872)", "Psi2S"}; + int checkPdgDecay[] = {443, 211, -211}; + int leptonPdg = 11; + Double_t rapidityWindow = 1.0; + std::string path{"o2sim_Kine.root"}; + for(int iSig =0; iSig < 2; iSig++) std::cout << "Check for\nsignal PDG " << checkPdgSignal[iSig] << "\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[]={0,0}; + int nAntileptons[]={0,0}; + int nLeptonPairs[]={0,0}; + int nLeptonPairsToBeDone[]={0,0}; + int nSignalX3872[]={0,0}; + int nSignalPionsPos[]={0,0}; + int nSignalPionsNeg[]={0,0}; + int nSignalPsi2S{}; + int nSignalX3872WithinAcc[]={0,0}; + int nSignalPionsPosWithinAcc[]={0,0}; + int nSignalPionsNegWithinAcc[]={0,0}; + 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; + for(int iSig=0; iSig<2; iSig++) { + if (pdg == leptonPdg) { + // count leptons + nLeptons[iSig]++; + } else if(pdg == -leptonPdg) { + // count anti-leptons + nAntileptons[iSig]++; + } else if (pdg == checkPdgSignal[iSig]) { + // 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[iSig]++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalX3872WithinAcc[iSig]++; idX3872 = j; } + if(pdgDau == checkPdgDecay[1] ) { nSignalPionsPos[iSig]++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalPionsPosWithinAcc[iSig]++; } + if(pdgDau == checkPdgDecay[2] ) { nSignalPionsNeg[iSig]++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalPionsNegWithinAcc[iSig]++; } + } + + 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[iSig]++; + if (child0.getToBeDone() && child1.getToBeDone()) { + nLeptonPairsToBeDone[iSig]++; + } + } + } + } + } + } + + std::cout << "#events: " << nEvents << "\n"; + for(int iSig=0; iSig < 2; iSig++){ + std::cout << "#leptons from " << PdgSignalName[iSig] << ": " << nLeptons[iSig] << "\n" + << "#antileptons from " << PdgSignalName[iSig] << ": " << nAntileptons[iSig] << "\n" + << "#signal (jpsi <-" << PdgSignalName[iSig] <<"): " << nSignalX3872[iSig] << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalX3872WithinAcc[iSig] << "\n" + << "#signal (pi+ <-" << PdgSignalName[iSig] <<"): " << nSignalPionsPos[iSig] << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalPionsPosWithinAcc[iSig] << "\n" + << "#signal (pi- <-" << PdgSignalName[iSig] <<"): " << nSignalPionsNeg[iSig] << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalPionsNegWithinAcc[iSig] << "\n" + << "#lepton pairs from " << PdgSignalName[iSig] <<": " << nLeptonPairs[iSig] << "\n" + << "#lepton pairs to be done from " << PdgSignalName[iSig] <<": " << nLeptonPairs[iSig] << "\n"; + + + if (nLeptonPairs[iSig] == 0 || nLeptons[iSig] == 0 || nAntileptons[iSig] == 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[iSig] != nLeptonPairsToBeDone[iSig]) { + 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; +}