Photoproduction injection in PbPb hadronic event

MC/config/PWGDQ/external/generator/GeneratorCocktailStarlight_PbPb5TeV.C

@@ -0,0 +1,112 @@
+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;
+}
+void configureEvtGen(o2::eventgen::GeneratorEvtGen<GeneratorCocktail>* generator, bool midrapidity)
+{
+  generator->SetPolarization(1); // Transverse
+  if (midrapidity) {
+    generator->SetSizePdg(2);
+    generator->AddPdg(443, 0);
+    generator->AddPdg(100443, 1);
+    TString pathO2 = gSystem->ExpandPathName("$O2DPG_MC_CONFIG_ROOT/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen");
+    generator->SetDecayTable(Form("%s/PSITOJPSIPIPI.DEC", pathO2.Data()));
+  } else {
+    generator->SetSizePdg(1);
+    generator->AddPdg(100443, 0);
+    TString pathO2 = gSystem->ExpandPathName("$O2DPG_MC_CONFIG_ROOT/MC/config/PWGUD/external/generator/DecayTablesEvtGen");
+    generator->SetDecayTable(Form("%s/PSI2S.MUMUPIPI.DEC", pathO2.Data()));
+  }
+}
+} // 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>();
+  configureEvtGen(genCocktailEvtGen, true);
+  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);
+
+  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>();
+  configureEvtGen(genCocktailEvtGen, true);
+  auto genCohJpsi = makeStarlightGenerator("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 = makeStarlightGenerator("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>();
+  configureEvtGen(genCocktailEvtGen, true);
+  auto genIncohJpsi = makeStarlightGenerator("kIncohJpsiToElRad", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genIncohJpsi, 1);
+  auto genIncohPsi2S = makeStarlightGenerator("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>();
+  configureEvtGen(genCocktailEvtGen, false);
+  auto genCohJpsi = makeStarlightGenerator("kCohJpsiToMu", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genCohJpsi, 1);
+  auto genTwoGammaToMu = makeStarlightGenerator("kTwoGammaToMuLow", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genTwoGammaToMu, 1);
+  auto genCohPsi2S = makeStarlightGenerator("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>();
+  configureEvtGen(genCocktailEvtGen, false);
+  auto genIncohJpsi = makeStarlightGenerator("kIncohJpsiToMu", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genIncohJpsi, 1);
+  auto genIncohPsi2S = makeStarlightGenerator("kIncohPsi2sToMuPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
+  genCocktailEvtGen->AddGenerator(genIncohPsi2S, 1);
+
+  return genCocktailEvtGen;
+}

MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C

@@ -0,0 +1,212 @@
+#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
+    {
+        const double photoproductionImpactParameterCut = 3.5; // tmp cut to embed signal > 10% centrality
+        if (impactParameter < photoproductionImpactParameterCut) {
+            std::cout<<"impact parameter: "<<impactParameter<<" fm, below photoproduction cut "
+                     <<photoproductionImpactParameterCut<<" fm; embedding no signal events"<<std::endl;
+            return 0;
+        }
+        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;
+}
+
+FairGenerator*
+  GeneratorPhotoproductionFwd(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(13); // mu
+//   auto seed = (gRandom->TRandom::GetSeed() % 900000000);
+//   gen->readString("Random:setSeed on");
+//   gen->readString("Random:seed " + std::to_string(seed));
+  return gen;
+}

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 = GeneratorPhotoproductionFwd(2, 3, -2.0, -4.5)
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_PbPb_536tev.cfg

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_PbPb_536tev.cfg

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 = GeneratorPhotoproductionFwd(2, 4, -2.0, -4.5)
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_PbPb_536tev.cfg

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_PbPb_536tev.cfg

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_PbPb_536tev.cfg