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..9cc44c4d7
--- /dev/null
+++ b/MC/config/PWGDQ/external/generator/GeneratorCocktailStarlight_PbPb5TeV.C
@@ -0,0 +1,111 @@
+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;
+}
+o2::eventgen::GeneratorEvtGen* makeStarlightToEvtGenGenerator(std::string configuration,
+                                                               float energyCM,
+                                                               int beam1Z,
+                                                               int beam1A,
+                                                               int beam2Z,
+                                                               int beam2A,
+                                                               std::string extraParams = "",
+                                                               std::string dpmjetconf = "")
+{
+  auto generator = makeStarlightGenerator(configuration, energyCM, beam1Z, beam1A, beam2Z, beam2A, extraParams, dpmjetconf);
+  generator->SetPolarization(1); //Transversal
+  generator->SetSizePdg(2);
+  generator->AddPdg(443,0);
+  generator->AddPdg(100443,1);
+  TString pathO2 = gSystem->ExpandPathName("$O2DPG_MC_CONFIG_ROOT/MC/config/PWGUD/external/generator/DecayTablesEvtGen");
+  if      (configuration.find("Psi2sToMuPi") != std::string::npos) generator->SetDecayTable(Form("%s/PSI2S.MUMUPIPI.DEC",pathO2.Data()));
+  else if (configuration.find("Psi2sToElPi") != std::string::npos) generator->SetDecayTable(Form("%s/PSI2S.EEPIPI.DEC",pathO2.Data()));
+  else if (configuration.find("JpsiToElRad") != std::string::npos) generator->SetDecayTable(Form("%s/JPSI.EE.DEC",pathO2.Data()));
+  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 = makeStarlightToEvtGenGenerator("kCohJpsiToElRad", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genCohJpsi, 1);
+  auto genIncohJpsi = makeStarlightToEvtGenGenerator("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 = makeStarlightToEvtGenGenerator("kCohPsi2sToElPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genCohPsi2S, 1);
+  auto genIncohPsi2S = makeStarlightToEvtGenGenerator("kIncohPsi2sToElPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genIncohPsi2S, 1);
+
+  return genCocktailEvtGen;
+}
+
+FairGenerator* GeneratorCocktailStarlightCoherentMidy_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 = makeStarlightToEvtGenGenerator("kCohJpsiToElRad", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genCohJpsi, 1);
+  auto genTwoGammaToEl = makeStarlightGenerator("kTwoGammaToElLow", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genTwoGammaToEl, 1);
+  auto genCohPsi2S = makeStarlightToEvtGenGenerator("kCohPsi2sToElPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genCohPsi2S, 1);
+
+  return genCocktailEvtGen;
+}
+
+FairGenerator* GeneratorCocktailStarlightIncoherentMidy_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 genIncohJpsi = makeStarlightToEvtGenGenerator("kIncohJpsiToElRad", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genIncohJpsi, 1);
+  auto genIncohPsi2S = makeStarlightToEvtGenGenerator("kIncohPsi2sToElPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genIncohPsi2S, 1);
+
+  return genCocktailEvtGen;
+}
+
+FairGenerator* GeneratorCocktailStarlightCoherentFwdy_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("kCohJpsiToMu", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genCohJpsi, 1);
+  auto genTwoGammaToEl = makeStarlightGenerator("kTwoGammaToMuLow", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genTwoGammaToEl, 1);
+  auto genCohPsi2S = makeStarlightToEvtGenGenerator("kCohPsi2sToMuPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genCohPsi2S, 1);
+
+  return genCocktailEvtGen;
+}
+
+FairGenerator* GeneratorCocktailStarlightIncoherentFwdy_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 genIncohJpsi = makeStarlightGenerator("kIncohJpsiToMu", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genIncohJpsi, 1);
+  auto genIncohPsi2S = makeStarlightToEvtGenGenerator("kIncohPsi2sToMuPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genIncohPsi2S, 1);
+
+  return genCocktailEvtGen;
+}
\ No newline at end of file
diff --git a/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C b/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C
new file mode 100644
index 000000000..42aa445dd
--- /dev/null
+++ b/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C
@@ -0,0 +1,192 @@
+#include "FairGenerator.h"
+#include "Generators/GeneratorPythia8.h"
+#include "Pythia8/HeavyIons.h"
+#include "Pythia8/Pythia.h"
+#include "TRandom.h"
+#include "GeneratorCocktailStarlight_PbPb5TeV.C"
+#include <algorithm>
+#include <cmath>
+#include <string>
+#include <vector>
+
+using namespace o2::eventgen;
+using namespace Pythia8;
+
+namespace o2
+{
+namespace eventgen
+{
+
+
+class GeneratorPythia8StarlightTriggeredWithGap : public o2::eventgen::GeneratorPythia8 {
+public: 
+
+    /// default constructor
+    GeneratorPythia8StarlightTriggeredWithGap() = default;
+
+    /// constructor
+    GeneratorPythia8StarlightTriggeredWithGap(int inputTriggerRatio = 5, int gentype = 0) {
+        mGeneratedEvents = 0;
+        mInverseTriggerRatio = inputTriggerRatio; 
+        mNumSignals = 0;
+        switch (gentype) {
+            case 0: // generate photoproduced charmonia cocktail at mid rapidity at 5TeV
+                mGeneratorParam = (Generator*)GeneratorCocktailStarlightMidy_PbPb5TeV();
+                break;
+            case 1: // generate coherent charmonia cocktail at mid rapidity at 5TeV
+                mGeneratorParam = (Generator*)GeneratorCocktailStarlightCoherentMidy_PbPb5TeV();
+                break;
+            case 2: // generate incoherent charmonia cocktail at mid rapidity at 5TeV
+                mGeneratorParam = (Generator*)GeneratorCocktailStarlightIncoherentMidy_PbPb5TeV();
+                break;
+            case 3: // generate coherent charmonia cocktail at forward rapidity at 5TeV
+                mGeneratorParam = (Generator*)GeneratorCocktailStarlightCoherentFwdy_PbPb5TeV();
+                break;
+            case 4: // generate incoherent charmonia cocktail at forward rapidity at 5TeV
+                mGeneratorParam = (Generator*)GeneratorCocktailStarlightIncoherentFwdy_PbPb5TeV();
+                break;
+        }
+        mGeneratorParam->Init();  
+
+        addSubGenerator(0, "Minimum bias");
+        addSubGenerator(1, "event with injected signals");
+    }
+
+    /// Deconstructor
+    ~GeneratorPythia8StarlightTriggeredWithGap() = default;
+
+    void addSignalPDGs(int pdg) { mSignalsPDGs.push_back(pdg); };
+
+    void setRapidityRange(double valMin, double valMax)
+    {
+        mHadronRapidityMin = valMin;
+        mHadronRapidityMax = valMax;
+    };
+
+protected: 
+    Bool_t generateEvent() override
+    {
+        // GeneratorPythia8::generateEvent();
+        mNumSignals = 0;
+        bool genOk = false;
+        if (mGeneratedEvents % mInverseTriggerRatio == 0){ // add injected prompt signals to the stack
+            std::cout<<"generating event with injected signals"<<std::endl;
+            while (!genOk){ 
+                genOk = GeneratorPythia8::generateEvent();
+            }
+            double impactParameter = -1.;
+            if (mPythia.info.hiInfo) {
+                impactParameter = mPythia.info.hiInfo->b();
+                mNumSignals = getNumSignalsForImpactParameter(impactParameter);
+            } else {
+                std::cout<<"no heavy-ion info available; embedding no signal events"<<std::endl;
+            }
+            std::cout<<"MB event impact parameter: "<<impactParameter<<", embedding "<<mNumSignals<<" signal events"<<std::endl;
+            notifySubGenerator(mNumSignals > 0 ? 1 : 0);
+        } else { // gap event
+            while (!genOk) {
+                genOk = GeneratorPythia8::generateEvent();
+            }
+            notifySubGenerator(0);
+        }
+        mGeneratedEvents++;
+        std::cout<<"generated events: "<<mGeneratedEvents<<std::endl;
+        return true;
+    }
+
+    Bool_t importParticles() override
+    {
+        GeneratorPythia8::importParticles();
+
+        if ((mGeneratedEvents-1) % mInverseTriggerRatio == 0){ // add injected prompt signals to the stack
+            for (int isig = 0; isig < mNumSignals; isig++) {
+                bool found = false;
+                while (!found) {
+                    if (!mGeneratorParam->generateEvent()) {
+                        continue;
+                    }
+                    if (!mGeneratorParam->importParticles()) {
+                        mGeneratorParam->clearParticles();
+                        continue;
+                    }
+                    found = findSignalInAcceptance();
+                    if (found) {
+                        appendCurrentSignalParticles(isig);
+                    }
+                    mGeneratorParam->clearParticles();
+                }
+            }
+        }
+
+        return true;
+    }
+
+    bool findSignalInAcceptance() {
+        std::cout<<"loop over" << mGeneratorParam->getParticles().size()<<" particles"<<std::endl;
+        for (int pdg : mSignalsPDGs) {
+            std::cout<<"signal pdg: "<<pdg<<std::endl;
+        }
+        for (int ipart = 0; ipart < mGeneratorParam->getParticles().size(); ipart++) {
+            TParticle part = TParticle(mGeneratorParam->getParticles().at(ipart));
+            // make sure all signals are in the acceptance
+            for (int pdg : mSignalsPDGs) {
+                if (part.GetPdgCode() == pdg) {
+                    std::cout<<"found signal with pdg: "<<part.GetPdgCode()<<", mother: "<<part.GetFirstMother()<<std::endl;
+                    if (part.GetFirstMother() == -1) {
+                        if (part.Y() < mHadronRapidityMin || part.Y() > mHadronRapidityMax) {
+                            return false;
+                        }
+                    }
+                }
+            }
+        }
+        std::cout<<"generated signal in acceptance"<<std::endl;
+        return true;
+    }
+
+    int getNumSignalsForImpactParameter(double impactParameter) const
+    {
+        return std::max(0, static_cast<int>(std::lround(5.0 + 0.886202881 * std::pow(std::max(0.0, 17.5 - impactParameter), 1.7))));
+    }
+
+    void appendCurrentSignalParticles(int signalIndex)
+    {
+        int originalSize = mParticles.size();
+        std::cout<<"adding "<<mGeneratorParam->getParticles().size()<<" particles from signal event "<<signalIndex<<" to the stack"<<std::endl;
+        for(size_t ipart=0; ipart < mGeneratorParam->getParticles().size(); ipart++){
+            TParticle part = TParticle(mGeneratorParam->getParticles().at(ipart));
+            if(part.GetFirstMother() >= 0) part.SetFirstMother(part.GetFirstMother() + originalSize);
+            if(part.GetSecondMother() >= 0) part.SetLastMother(part.GetSecondMother() + originalSize);
+            if(part.GetFirstDaughter() >= 0) part.SetFirstDaughter(part.GetFirstDaughter() + originalSize);
+            if(part.GetLastDaughter() >= 0) part.SetLastDaughter(part.GetLastDaughter() + originalSize);
+            mParticles.push_back(part);
+            // encodeParticleStatusAndTracking method already called in GeneratorEvtGen.C
+        }
+    }
+
+private:
+    Generator* mGeneratorParam = nullptr;
+    unsigned long long mGeneratedEvents = 0;
+    int mInverseTriggerRatio = 1;
+    int mNumSignals = 0;
+    // Pythia8::Pythia pythiaMBgen; // minimum bias event
+    std::vector<int> mSignalsPDGs;
+    double mHadronRapidityMin;
+    double mHadronRapidityMax;
+};
+}
+}
+
+FairGenerator*
+  GeneratorPhotoproduction(int triggerGap, int gentype = 0, double rapidityMin = -1.5, double rapidityMax = 1.5)
+{
+  auto gen = new o2::eventgen::GeneratorPythia8StarlightTriggeredWithGap(triggerGap, gentype);
+  gen->setRapidityRange(rapidityMin, rapidityMax);
+  gen->addSignalPDGs(443); // J/Psi
+  gen->addSignalPDGs(100443); // Psi(2S)
+  gen->addSignalPDGs(11); // e
+//   auto seed = (gRandom->TRandom::GetSeed() % 900000000);
+//   gen->readString("Random:setSeed on");
+//   gen->readString("Random:seed " + std::to_string(seed));
+  return gen;
+}
\ No newline at end of file
diff --git a/MC/config/PWGDQ/ini/Generator_InjectedCoherentCharmoniaFwdy_PbPb5TeV.ini b/MC/config/PWGDQ/ini/Generator_InjectedCoherentCharmoniaFwdy_PbPb5TeV.ini
new file mode 100644
index 000000000..c69e9f661
--- /dev/null
+++ b/MC/config/PWGDQ/ini/Generator_InjectedCoherentCharmoniaFwdy_PbPb5TeV.ini
@@ -0,0 +1,5 @@
+[GeneratorExternal] 
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C 
+funcName = GeneratorPhotoproduction(2, 3, 2.0, 4.5)
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_pp_536tev.cfg
diff --git a/MC/config/PWGDQ/ini/Generator_InjectedCoherentCharmoniaMidy_PbPb5TeV.ini b/MC/config/PWGDQ/ini/Generator_InjectedCoherentCharmoniaMidy_PbPb5TeV.ini
new file mode 100644
index 000000000..939a7327e
--- /dev/null
+++ b/MC/config/PWGDQ/ini/Generator_InjectedCoherentCharmoniaMidy_PbPb5TeV.ini
@@ -0,0 +1,5 @@
+[GeneratorExternal] 
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C 
+funcName = GeneratorPhotoproduction(2, 1, -1.5, 1.5)
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_pp_536tev.cfg
diff --git a/MC/config/PWGDQ/ini/Generator_InjectedIncoherentCharmoniaFwdy_PbPb5TeV.ini b/MC/config/PWGDQ/ini/Generator_InjectedIncoherentCharmoniaFwdy_PbPb5TeV.ini
new file mode 100644
index 000000000..5024d65e9
--- /dev/null
+++ b/MC/config/PWGDQ/ini/Generator_InjectedIncoherentCharmoniaFwdy_PbPb5TeV.ini
@@ -0,0 +1,5 @@
+[GeneratorExternal] 
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C 
+funcName = GeneratorPhotoproduction(2, 4, 2.0, 4.5)
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_pp_536tev.cfg
diff --git a/MC/config/PWGDQ/ini/Generator_InjectedIncoherentCharmoniaMidy_PbPb5TeV.ini b/MC/config/PWGDQ/ini/Generator_InjectedIncoherentCharmoniaMidy_PbPb5TeV.ini
new file mode 100644
index 000000000..e7f4adb39
--- /dev/null
+++ b/MC/config/PWGDQ/ini/Generator_InjectedIncoherentCharmoniaMidy_PbPb5TeV.ini
@@ -0,0 +1,5 @@
+[GeneratorExternal] 
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C 
+funcName = GeneratorPhotoproduction(2, 2, -1.5, 1.5)
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_pp_536tev.cfg
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..5a6f58e5c
--- /dev/null
+++ b/MC/config/PWGDQ/ini/Generator_InjectedPhotoproucedCharmoniaMidy_PbPb5TeV.ini
@@ -0,0 +1,5 @@
+[GeneratorExternal] 
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C 
+funcName = GeneratorPhotoproduction(2, 0, -1.5, 1.5)
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_pp_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<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;
+}
\ No newline at end of file
diff --git a/MC/config/PWGDQ/pythia8/decayer/jpsiPsi2Sdielectron.cfg b/MC/config/PWGDQ/pythia8/decayer/jpsiPsi2Sdielectron.cfg
new file mode 100644
index 000000000..145dc8cec
--- /dev/null
+++ b/MC/config/PWGDQ/pythia8/decayer/jpsiPsi2Sdielectron.cfg
@@ -0,0 +1,9 @@
+### Force J/psi and psi(2S) -> e+ e-
+
+### J/psi -> e+ e-
+443:onMode = off		# switch off all decay channels
+443:onIfAll = 11 -11		# switch on e+ e- decay channel
+
+### psi(2S) -> e+ e-
+100443:onMode = off		# switch off all decay channels
+100443:onIfAll = 11 -11		# switch on e+ e- decay channel
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