Code for photoproduction at midrapidity

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 <string>
+
+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<GeneratorCocktail>();
+  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;
+}

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

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<o2::MCTrack>* 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;
+}

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