diff --git a/MC/config/PWGDQ/external/generator/GeneratorCocktailStarlight_PbPb5TeV.C b/MC/config/PWGDQ/external/generator/GeneratorCocktailStarlight_PbPb5TeV.C new file mode 100644 index 000000000..1463c9fb9 --- /dev/null +++ b/MC/config/PWGDQ/external/generator/GeneratorCocktailStarlight_PbPb5TeV.C @@ -0,0 +1,53 @@ +R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT/MC/config/PWGUD/external/generator) +R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT/MC/config/PWGDQ/EvtGen) +#include "GeneratorCocktail.C" +#include "GeneratorStarlight.C" +#include "GeneratorEvtGen.C" +#include + +namespace +{ +o2::eventgen::GeneratorStarlight_class* makeStarlightGenerator(std::string configuration, + float energyCM, + int beam1Z, + int beam1A, + int beam2Z, + int beam2A, + std::string extraParams = "", + std::string dpmjetconf = "") +{ + auto generator = new o2::eventgen::GeneratorStarlight_class(); + generator->selectConfiguration(configuration); + generator->setCollisionSystem(energyCM, beam1Z, beam1A, beam2Z, beam2A); + generator->setExtraParams(extraParams); + generator->setupDpmjet(dpmjetconf); + return generator; +} +} // namespace + +FairGenerator* GeneratorCocktailStarlightMidy_PbPb5TeV(float energyCM = 5360, int beam1Z = 82, int beam1A = 208, int beam2Z = 82, int beam2A = 208, std::string extrapars = "",std::string dpmjetconf = "") +{ + auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen(); + auto genCohJpsi = makeStarlightGenerator("kCohJpsiToElRad", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf); + genCocktailEvtGen->AddGenerator(genCohJpsi, 1); + auto genIncohJpsi = makeStarlightGenerator("kIncohJpsiToElRad", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf); + genCocktailEvtGen->AddGenerator(genIncohJpsi, 1); + auto genTwoGammaToEl = makeStarlightGenerator("kTwoGammaToElLow", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf); + genCocktailEvtGen->AddGenerator(genTwoGammaToEl, 1); + auto genCohPsi2S = makeStarlightGenerator("kCohPsi2sToElPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf); + genCocktailEvtGen->AddGenerator(genCohPsi2S, 1); + auto genIncohPsi2S = makeStarlightGenerator("kIncohPsi2sToElPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf); + genCocktailEvtGen->AddGenerator(genIncohPsi2S, 1); + + TString pdgs = "443;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)); + } + genCocktailEvtGen->SetForceDecay(kEvtPsiPrimeJpsiDiElectron); + return genCocktailEvtGen; +} \ No newline at end of file diff --git a/MC/config/PWGDQ/ini/Generator_InjectedPhotoproucedCharmoniaMidy_PbPb5TeV.ini b/MC/config/PWGDQ/ini/Generator_InjectedPhotoproucedCharmoniaMidy_PbPb5TeV.ini new file mode 100644 index 000000000..4d2931f03 --- /dev/null +++ b/MC/config/PWGDQ/ini/Generator_InjectedPhotoproucedCharmoniaMidy_PbPb5TeV.ini @@ -0,0 +1,8 @@ +[GeneratorExternal] +fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/GeneratorCocktailStarlight_PbPb5TeV.C +funcName = GeneratorCocktailStarlightMidy_PbPb5TeV() +[TriggerExternal] +fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C +funcName = selectMotherPartInAcc(-1.5,1.5) +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_PbPb_536tev.cfg diff --git a/MC/config/PWGDQ/ini/tests/Generator_InjectedPhotoproucedCharmoniaMidy_PbPb5TeV.C b/MC/config/PWGDQ/ini/tests/Generator_InjectedPhotoproucedCharmoniaMidy_PbPb5TeV.C new file mode 100644 index 000000000..fe55e024a --- /dev/null +++ b/MC/config/PWGDQ/ini/tests/Generator_InjectedPhotoproucedCharmoniaMidy_PbPb5TeV.C @@ -0,0 +1,104 @@ +int External() { + 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"); + if (!tree) + { + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; + return 1; + } + std::vector* tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + int checkPdgSignal[2] = {443, 100443}; + int checkPdgDecay = 11; + double rapiditymin = -1.5; double rapiditymax = 1.5; + + int nLeptons{}; + int nAntileptons{}; + int nLeptonPairs{}; + int nPrimaryLepton{}; + int nPrimaryAntilepton{}; + int nLeptonPairsToBeDone{}; + int nSignalJpsi{}; + int nSignalPsi2S{}; + int nPions{}; + int nPhotonsFromPsi{}; + int nPhotonsFromJPsi{}; + int nSignalJpsiWithinAcc{}; + int nSignalPsi2SWithinAcc{}; + + auto nEvents = tree->GetEntries(); + 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(); + if (pdg == checkPdgDecay) { + // count leptons + nLeptons++; + if (idMoth == -1) nPrimaryLepton++; + } else if (pdg == -checkPdgDecay) { + // count anti-leptons + nAntileptons++; + if (idMoth == -1) nPrimaryAntilepton++; + } else if (pdg == checkPdgSignal[0] || pdg == checkPdgSignal[1]) { + if (idMoth < 0) { + // count signal PDG + pdg == checkPdgSignal[0] ? nSignalJpsi++ : nSignalPsi2S++; + // count signal PDG within acceptance + if (rapidity > rapiditymin && rapidity < rapiditymax) { + pdg == checkPdgSignal[0] ? nSignalJpsiWithinAcc++ : nSignalPsi2SWithinAcc++; + } + } + int child0Id = track.getFirstDaughterTrackId(); + int child1Id = track.getLastDaughterTrackId(); + for (int iDaugh = child0Id; iDaugh <= child1Id; iDaugh++) { + auto daugh = tracks->at(iDaugh); + if (daugh.GetPdgCode() == 22) { + if (pdg == checkPdgSignal[0]) { + nPhotonsFromJPsi++; + } else if (pdg == checkPdgSignal[1]) { + nPhotonsFromPsi++; + } + } + if (TMath::Abs(daugh.GetPdgCode()) == 211) nPions++; + } + auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks); + auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks); + if (child0 != nullptr && child1 != nullptr) { + // check for parent-child relations + auto pdg0 = child0->GetPdgCode(); + auto pdg1 = child1->GetPdgCode(); + // std::cout << "First and last children of parent " << checkPdgSignal << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n"; + if (pdg == checkPdgSignal[0]) std::cout << "First and last children of parent Jpsi are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n"; + if (pdg == checkPdgSignal[1]) std::cout << "First and last children of parent Psi(2S) are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n"; + if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1) { + nLeptonPairs++; + if (child0->getToBeDone() && child1->getToBeDone()) { + nLeptonPairsToBeDone++; + } + } + } + } + } + } + + std::cout << "#events: " << nEvents << "\n" + << "#leptons: " << nLeptons << "\n" + << "#antileptons: " << nAntileptons << "\n" + << "#primary leptons: " << nPrimaryLepton << "\n" + << "#primary antileptons: " << nPrimaryAntilepton << "\n" + << "#signal (prompt Jpsi): " << nSignalJpsi << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalJpsiWithinAcc << "\n" + << "#signal (prompt Psi(2S)): " << nSignalPsi2S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalPsi2SWithinAcc << "\n" + << "#photons from Jpsi: " << nPhotonsFromJPsi <<" #photons from Psi(2S): "<< nPhotonsFromPsi<<"\n" + << "#pions from Psi(2S): "<< nPions<<"\n" + << "#lepton pairs: " << nLeptonPairs <<" #lepton pairs to be done: "<< nLeptonPairsToBeDone<<"\n"; + + return 0; +} \ No newline at end of file diff --git a/MC/run/PWGDQ/runPhotoproducedCharmonia_midy_PbPb5TeV.sh b/MC/run/PWGDQ/runPhotoproducedCharmonia_midy_PbPb5TeV.sh new file mode 100644 index 000000000..4347822d7 --- /dev/null +++ b/MC/run/PWGDQ/runPhotoproducedCharmonia_midy_PbPb5TeV.sh @@ -0,0 +1,34 @@ +#!/usr/bin/env bash + +RNDSEED=${RNDSEED:-0} +NWORKERS=${NWORKERS:-8} + + +export ALIEN_JDL_LPMPRODUCTIONTYPE=MC +export ALIEN_JDL_CPULIMIT=20 +export ALIEN_JDL_LPMANCHORPASSNAME=apass4 +export ALIEN_JDL_MCANCHOR=apass4 +export ALIEN_JDL_COLLISIONSYSTEM=PbPb +export ALIEN_JDL_LPMPASSNAME=apass4 +export ALIEN_JDL_LPMRUNNUMBER=544474 +export ALIEN_JDL_LPMANCHORRUN=544474 + +export ALIEN_JDL_LPMINTERACTIONTYPE=PbPb +export ALIEN_JDL_LPMPRODUCTIONTAG="TestProd" +export ALIEN_JDL_LPMANCHORPRODUCTION="TestProd" +export ALIEN_JDL_LPMANCHORYEAR=2023 + +export ALIEN_JDL_O2DPGWORKFLOWTARGET="aod" + +export PRODSPLIT=${ALIEN_O2DPG_GRIDSUBMIT_PRODSPLIT:-10} +export SPLITID=${ALIEN_O2DPG_GRIDSUBMIT_SUBJOBID:-10} +export CYCLE=0 +export NTIMEFRAMES=8 +export NSIGEVENTS=20 +export NBKGEVENTS=10 + +CONFIGNAME=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/ini/Generator_InjectedPhotoproucedCharmoniaMidy_PbPb5TeV.ini + +export ALIEN_JDL_ANCHOR_SIM_OPTIONS="-eCM 5360 -gen external -j ${NWORKERS} -tf ${NTIMEFRAMES} -e TGeant4 -seed 0 -ini ${CONFIGNAME} -genBkg pythia8 -procBkg \"heavy_ion\" -colBkg PbPb --embedding -nb $NBKGEVENTS --embeddPattern @0:e2" + +${O2DPG_ROOT}/MC/run/ANCHOR/anchorMC.sh \ No newline at end of file