Skip to content

Commit 9132c20

Browse files
ypwanggypwangg
authored
Photoproduction injection in PbPb hadronic event
1 parent 6393923 commit 9132c20

15 files changed

Lines changed: 1039 additions & 0 deletions

‎MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/PSITOJPSITOMU.DEC‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
Decay psi(2S)
2+
1.000 J/psi pi+ pi- VVPIPI;
3+
Enddecay
4+
5+
Decay J/psi
6+
1.000 mu+ mu- PHOTOS VLL;
7+
Enddecay
8+
9+
End

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

Lines changed: 121 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,121 @@
1+
R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT/MC/config/PWGUD/external/generator)
2+
R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT/MC/config/PWGDQ/EvtGen)
3+
#include "GeneratorCocktail.C"
4+
#include "GeneratorStarlight.C"
5+
#include "GeneratorEvtGen.C"
6+
#include <string>
7+
8+
namespace
9+
{
10+
o2::eventgen::GeneratorStarlight_class* makeStarlightGenerator(std::string configuration,
11+
float energyCM,
12+
int beam1Z,
13+
int beam1A,
14+
int beam2Z,
15+
int beam2A,
16+
std::string extraParams = "",
17+
std::string dpmjetconf = "")
18+
{
19+
auto generator = new o2::eventgen::GeneratorStarlight_class();
20+
generator->selectConfiguration(configuration);
21+
generator->setCollisionSystem(energyCM, beam1Z, beam1A, beam2Z, beam2A);
22+
generator->setExtraParams(extraParams);
23+
// generator->setupDpmjet(dpmjetconf);
24+
return generator;
25+
}
26+
void configureEvtGen(o2::eventgen::GeneratorEvtGen<GeneratorCocktail>* generator, bool midrapidity)
27+
{
28+
generator->SetPolarization(1); // Transverse
29+
generator->SetSizePdg(2);
30+
generator->AddPdg(443, 0);
31+
generator->AddPdg(100443, 1);
32+
TString pathO2 = gSystem->ExpandPathName("$O2DPG_MC_CONFIG_ROOT/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen");
33+
if (midrapidity) {
34+
// generator->SetSizePdg(2);
35+
// generator->AddPdg(443, 0);
36+
// generator->AddPdg(100443, 1);
37+
// TString pathO2 = gSystem->ExpandPathName("$O2DPG_MC_CONFIG_ROOT/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen");
38+
generator->SetDecayTable(Form("%s/PSITOJPSIPIPI.DEC", pathO2.Data()));
39+
} else {
40+
// generator->SetSizePdg(1);
41+
// generator->AddPdg(100443, 0);
42+
// generator->SetSizePdg(2);
43+
// generator->AddPdg(443, 0);
44+
// generator->AddPdg(100443, 1);
45+
// TString pathO2 = gSystem->ExpandPathName("$O2DPG_MC_CONFIG_ROOT/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen");
46+
generator->SetDecayTable(Form("%s/PSITOJPSITOMU.DEC", pathO2.Data()));
47+
// TString pathO2 = gSystem->ExpandPathName("$O2DPG_MC_CONFIG_ROOT/MC/config/PWGUD/external/generator/DecayTablesEvtGen");
48+
// generator->SetDecayTable(Form("%s/PSI2S.MUMUPIPI.DEC", pathO2.Data()));
49+
}
50+
}
51+
} // namespace
52+
53+
FairGenerator* GeneratorCocktailStarlightMidy_PbPb5TeV(float energyCM = 5360, int beam1Z = 82, int beam1A = 208, int beam2Z = 82, int beam2A = 208, std::string extrapars = "",std::string dpmjetconf = "")
54+
{
55+
auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen<GeneratorCocktail>();
56+
configureEvtGen(genCocktailEvtGen, true);
57+
auto genCohJpsi = makeStarlightGenerator("kCohJpsiToElRad", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
58+
genCocktailEvtGen->AddGenerator(genCohJpsi, 1);
59+
auto genIncohJpsi = makeStarlightGenerator("kIncohJpsiToElRad", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
60+
genCocktailEvtGen->AddGenerator(genIncohJpsi, 1);
61+
auto genTwoGammaToEl = makeStarlightGenerator("kTwoGammaToElLow", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
62+
genCocktailEvtGen->AddGenerator(genTwoGammaToEl, 1);
63+
auto genCohPsi2S = makeStarlightGenerator("kCohPsi2sToElPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
64+
genCocktailEvtGen->AddGenerator(genCohPsi2S, 1);
65+
auto genIncohPsi2S = makeStarlightGenerator("kIncohPsi2sToElPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
66+
genCocktailEvtGen->AddGenerator(genIncohPsi2S, 1);
67+
68+
return genCocktailEvtGen;
69+
}
70+
71+
FairGenerator* GeneratorCocktailStarlightCoherentMidy_PbPb5TeV(float energyCM = 5360, int beam1Z = 82, int beam1A = 208, int beam2Z = 82, int beam2A = 208, std::string extrapars = "",std::string dpmjetconf = "")
72+
{
73+
auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen<GeneratorCocktail>();
74+
configureEvtGen(genCocktailEvtGen, true);
75+
auto genCohJpsi = makeStarlightGenerator("kCohJpsiToElRad", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
76+
genCocktailEvtGen->AddGenerator(genCohJpsi, 1);
77+
auto genTwoGammaToEl = makeStarlightGenerator("kTwoGammaToElLow", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
78+
genCocktailEvtGen->AddGenerator(genTwoGammaToEl, 1);
79+
auto genCohPsi2S = makeStarlightGenerator("kCohPsi2sToElPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
80+
genCocktailEvtGen->AddGenerator(genCohPsi2S, 1);
81+
82+
return genCocktailEvtGen;
83+
}
84+
85+
FairGenerator* GeneratorCocktailStarlightIncoherentMidy_PbPb5TeV(float energyCM = 5360, int beam1Z = 82, int beam1A = 208, int beam2Z = 82, int beam2A = 208, std::string extrapars = "",std::string dpmjetconf = "")
86+
{
87+
auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen<GeneratorCocktail>();
88+
configureEvtGen(genCocktailEvtGen, true);
89+
auto genIncohJpsi = makeStarlightGenerator("kIncohJpsiToElRad", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
90+
genCocktailEvtGen->AddGenerator(genIncohJpsi, 1);
91+
auto genIncohPsi2S = makeStarlightGenerator("kIncohPsi2sToElPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
92+
genCocktailEvtGen->AddGenerator(genIncohPsi2S, 1);
93+
94+
return genCocktailEvtGen;
95+
}
96+
97+
FairGenerator* GeneratorCocktailStarlightCoherentFwdy_PbPb5TeV(float energyCM = 5360, int beam1Z = 82, int beam1A = 208, int beam2Z = 82, int beam2A = 208, std::string extrapars = "",std::string dpmjetconf = "")
98+
{
99+
auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen<GeneratorCocktail>();
100+
configureEvtGen(genCocktailEvtGen, false);
101+
auto genCohJpsi = makeStarlightGenerator("kCohJpsiToElRad", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
102+
genCocktailEvtGen->AddGenerator(genCohJpsi, 1);
103+
auto genTwoGammaToMu = makeStarlightGenerator("kTwoGammaToMuLow", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
104+
genCocktailEvtGen->AddGenerator(genTwoGammaToMu, 1);
105+
auto genCohPsi2S = makeStarlightGenerator("kCohPsi2sToMuPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
106+
genCocktailEvtGen->AddGenerator(genCohPsi2S, 1);
107+
108+
return genCocktailEvtGen;
109+
}
110+
111+
FairGenerator* GeneratorCocktailStarlightIncoherentFwdy_PbPb5TeV(float energyCM = 5360, int beam1Z = 82, int beam1A = 208, int beam2Z = 82, int beam2A = 208, std::string extrapars = "",std::string dpmjetconf = "")
112+
{
113+
auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen<GeneratorCocktail>();
114+
configureEvtGen(genCocktailEvtGen, false);
115+
auto genIncohJpsi = makeStarlightGenerator("kIncohJpsiToElRad", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
116+
genCocktailEvtGen->AddGenerator(genIncohJpsi, 1);
117+
auto genIncohPsi2S = makeStarlightGenerator("kIncohPsi2sToMuPi", energyCM, beam1Z, beam1A, beam2Z, beam2A, extrapars, dpmjetconf);
118+
genCocktailEvtGen->AddGenerator(genIncohPsi2S, 1);
119+
120+
return genCocktailEvtGen;
121+
}

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

Lines changed: 212 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,212 @@
1+
#include "FairGenerator.h"
2+
#include "Generators/GeneratorPythia8.h"
3+
#include "Pythia8/HeavyIons.h"
4+
#include "Pythia8/Pythia.h"
5+
#include "TRandom.h"
6+
#include "GeneratorCocktailStarlight_PbPb5TeV.C"
7+
#include <algorithm>
8+
#include <cmath>
9+
#include <string>
10+
#include <vector>
11+
12+
using namespace o2::eventgen;
13+
using namespace Pythia8;
14+
15+
namespace o2
16+
{
17+
namespace eventgen
18+
{
19+
20+
21+
class GeneratorPythia8StarlightTriggeredWithGap : public o2::eventgen::GeneratorPythia8 {
22+
public:
23+
24+
/// default constructor
25+
GeneratorPythia8StarlightTriggeredWithGap() = default;
26+
27+
/// constructor
28+
GeneratorPythia8StarlightTriggeredWithGap(int inputTriggerRatio = 5, int gentype = 0) {
29+
mGeneratedEvents = 0;
30+
mInverseTriggerRatio = inputTriggerRatio;
31+
mNumSignals = 0;
32+
switch (gentype) {
33+
case 0: // generate photoproduced charmonia cocktail at mid rapidity at 5TeV
34+
mGeneratorParam = (Generator*)GeneratorCocktailStarlightMidy_PbPb5TeV();
35+
break;
36+
case 1: // generate coherent charmonia cocktail at mid rapidity at 5TeV
37+
mGeneratorParam = (Generator*)GeneratorCocktailStarlightCoherentMidy_PbPb5TeV();
38+
break;
39+
case 2: // generate incoherent charmonia cocktail at mid rapidity at 5TeV
40+
mGeneratorParam = (Generator*)GeneratorCocktailStarlightIncoherentMidy_PbPb5TeV();
41+
break;
42+
case 3: // generate coherent charmonia cocktail at forward rapidity at 5TeV
43+
mGeneratorParam = (Generator*)GeneratorCocktailStarlightCoherentFwdy_PbPb5TeV();
44+
break;
45+
case 4: // generate incoherent charmonia cocktail at forward rapidity at 5TeV
46+
mGeneratorParam = (Generator*)GeneratorCocktailStarlightIncoherentFwdy_PbPb5TeV();
47+
break;
48+
}
49+
mGeneratorParam->Init();
50+
51+
addSubGenerator(0, "Minimum bias");
52+
addSubGenerator(1, "event with injected signals");
53+
}
54+
55+
/// Deconstructor
56+
~GeneratorPythia8StarlightTriggeredWithGap() = default;
57+
58+
void addSignalPDGs(int pdg) { mSignalsPDGs.push_back(pdg); };
59+
60+
void setRapidityRange(double valMin, double valMax)
61+
{
62+
mHadronRapidityMin = valMin;
63+
mHadronRapidityMax = valMax;
64+
};
65+
66+
protected:
67+
Bool_t generateEvent() override
68+
{
69+
// GeneratorPythia8::generateEvent();
70+
mNumSignals = 0;
71+
bool genOk = false;
72+
if (mGeneratedEvents % mInverseTriggerRatio == 0){ // add injected prompt signals to the stack
73+
std::cout<<"generating event with injected signals"<<std::endl;
74+
while (!genOk){
75+
genOk = GeneratorPythia8::generateEvent();
76+
}
77+
double impactParameter = -1.;
78+
if (mPythia.info.hiInfo) {
79+
impactParameter = mPythia.info.hiInfo->b();
80+
mNumSignals = getNumSignalsForImpactParameter(impactParameter);
81+
} else {
82+
std::cout<<"no heavy-ion info available; embedding no signal events"<<std::endl;
83+
}
84+
std::cout<<"MB event impact parameter: "<<impactParameter<<", embedding "<<mNumSignals<<" signal events"<<std::endl;
85+
notifySubGenerator(mNumSignals > 0 ? 1 : 0);
86+
} else { // gap event
87+
while (!genOk) {
88+
genOk = GeneratorPythia8::generateEvent();
89+
}
90+
notifySubGenerator(0);
91+
}
92+
mGeneratedEvents++;
93+
std::cout<<"generated events: "<<mGeneratedEvents<<std::endl;
94+
return true;
95+
}
96+
97+
Bool_t importParticles() override
98+
{
99+
GeneratorPythia8::importParticles();
100+
101+
if ((mGeneratedEvents-1) % mInverseTriggerRatio == 0){ // add injected prompt signals to the stack
102+
for (int isig = 0; isig < mNumSignals; isig++) {
103+
bool found = false;
104+
while (!found) {
105+
if (!mGeneratorParam->generateEvent()) {
106+
continue;
107+
}
108+
if (!mGeneratorParam->importParticles()) {
109+
mGeneratorParam->clearParticles();
110+
continue;
111+
}
112+
found = findSignalInAcceptance();
113+
if (found) {
114+
appendCurrentSignalParticles(isig);
115+
}
116+
mGeneratorParam->clearParticles();
117+
}
118+
}
119+
}
120+
121+
return true;
122+
}
123+
124+
bool findSignalInAcceptance() {
125+
// std::cout<<"loop over" << mGeneratorParam->getParticles().size()<<" particles"<<std::endl;
126+
for (int pdg : mSignalsPDGs) {
127+
// std::cout<<"signal pdg: "<<pdg<<std::endl;
128+
}
129+
for (int ipart = 0; ipart < mGeneratorParam->getParticles().size(); ipart++) {
130+
TParticle part = TParticle(mGeneratorParam->getParticles().at(ipart));
131+
// make sure all signals are in the acceptance
132+
for (int pdg : mSignalsPDGs) {
133+
if (part.GetPdgCode() == pdg) {
134+
// std::cout<<"found signal with pdg: "<<part.GetPdgCode()<<", mother: "<<part.GetFirstMother()<<std::endl;
135+
if (part.GetFirstMother() == -1) {
136+
if (part.Y() < mHadronRapidityMin || part.Y() > mHadronRapidityMax) {
137+
return false;
138+
}
139+
}
140+
}
141+
}
142+
}
143+
std::cout<<"generated signal in acceptance"<<std::endl;
144+
return true;
145+
}
146+
147+
int getNumSignalsForImpactParameter(double impactParameter) const
148+
{
149+
// const double photoproductionImpactParameterCut = 3.5; // tmp cut to embed signal > 10% centrality
150+
// if (impactParameter < photoproductionImpactParameterCut) {
151+
// std::cout<<"impact parameter: "<<impactParameter<<" fm, below photoproduction cut "
152+
// <<photoproductionImpactParameterCut<<" fm; embedding no signal events"<<std::endl;
153+
// return 0;
154+
// }
155+
return std::max(0, static_cast<int>(std::lround(5.0 + 0.886202881 * std::pow(std::max(0.0, 17.5 - impactParameter), 1.7))));
156+
}
157+
158+
void appendCurrentSignalParticles(int signalIndex)
159+
{
160+
int originalSize = mParticles.size();
161+
// std::cout<<"adding "<<mGeneratorParam->getParticles().size()<<" particles from signal event "<<signalIndex<<" to the stack"<<std::endl;
162+
for(size_t ipart=0; ipart < mGeneratorParam->getParticles().size(); ipart++){
163+
TParticle part = TParticle(mGeneratorParam->getParticles().at(ipart));
164+
if(part.GetFirstMother() >= 0) part.SetFirstMother(part.GetFirstMother() + originalSize);
165+
if(part.GetSecondMother() >= 0) part.SetLastMother(part.GetSecondMother() + originalSize);
166+
if(part.GetFirstDaughter() >= 0) part.SetFirstDaughter(part.GetFirstDaughter() + originalSize);
167+
if(part.GetLastDaughter() >= 0) part.SetLastDaughter(part.GetLastDaughter() + originalSize);
168+
mParticles.push_back(part);
169+
// encodeParticleStatusAndTracking method already called in GeneratorEvtGen.C
170+
}
171+
}
172+
173+
private:
174+
Generator* mGeneratorParam = nullptr;
175+
unsigned long long mGeneratedEvents = 0;
176+
int mInverseTriggerRatio = 1;
177+
int mNumSignals = 0;
178+
// Pythia8::Pythia pythiaMBgen; // minimum bias event
179+
std::vector<int> mSignalsPDGs;
180+
double mHadronRapidityMin;
181+
double mHadronRapidityMax;
182+
};
183+
}
184+
}
185+
186+
FairGenerator*
187+
GeneratorPhotoproduction(int triggerGap, int gentype = 0, double rapidityMin = -1.5, double rapidityMax = 1.5)
188+
{
189+
auto gen = new o2::eventgen::GeneratorPythia8StarlightTriggeredWithGap(triggerGap, gentype);
190+
gen->setRapidityRange(rapidityMin, rapidityMax);
191+
gen->addSignalPDGs(443); // J/Psi
192+
gen->addSignalPDGs(100443); // Psi(2S)
193+
gen->addSignalPDGs(11); // e
194+
// auto seed = (gRandom->TRandom::GetSeed() % 900000000);
195+
// gen->readString("Random:setSeed on");
196+
// gen->readString("Random:seed " + std::to_string(seed));
197+
return gen;
198+
}
199+
200+
FairGenerator*
201+
GeneratorPhotoproductionFwd(int triggerGap, int gentype = 0, double rapidityMin = -1.5, double rapidityMax = 1.5)
202+
{
203+
auto gen = new o2::eventgen::GeneratorPythia8StarlightTriggeredWithGap(triggerGap, gentype);
204+
gen->setRapidityRange(rapidityMin, rapidityMax);
205+
gen->addSignalPDGs(443); // J/Psi
206+
gen->addSignalPDGs(100443); // Psi(2S)
207+
gen->addSignalPDGs(13); // mu
208+
// auto seed = (gRandom->TRandom::GetSeed() % 900000000);
209+
// gen->readString("Random:setSeed on");
210+
// gen->readString("Random:seed " + std::to_string(seed));
211+
return gen;
212+
}

‎MC/config/PWGDQ/ini/Generator_InjectedCoherentCharmoniaFwdy_PbPb5TeV.ini‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
[GeneratorExternal]
2+
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C
3+
funcName = GeneratorPhotoproductionFwd(2, 3, -4.5, -2.0)
4+
[GeneratorPythia8]
5+
config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_PbPb_536tev.cfg

‎MC/config/PWGDQ/ini/Generator_InjectedCoherentCharmoniaMidy_PbPb5TeV.ini‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
[GeneratorExternal]
2+
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C
3+
funcName = GeneratorPhotoproduction(2, 1, -1.5, 1.5)
4+
[GeneratorPythia8]
5+
config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_PbPb_536tev.cfg

‎MC/config/PWGDQ/ini/Generator_InjectedIncoherentCharmoniaFwdy_PbPb5TeV.ini‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
[GeneratorExternal]
2+
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C
3+
funcName = GeneratorPhotoproductionFwd(2, 4, -4.5, -2.0)
4+
[GeneratorPythia8]
5+
config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_PbPb_536tev.cfg

‎MC/config/PWGDQ/ini/Generator_InjectedIncoherentCharmoniaMidy_PbPb5TeV.ini‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
[GeneratorExternal]
2+
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C
3+
funcName = GeneratorPhotoproduction(2, 2, -1.5, 1.5)
4+
[GeneratorPythia8]
5+
config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_PbPb_536tev.cfg

‎MC/config/PWGDQ/ini/Generator_InjectedPhotoproucedCharmoniaMidy_PbPb5TeV.ini‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
[GeneratorExternal]
2+
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C
3+
funcName = GeneratorPhotoproduction(2, 0, -1.5, 1.5)
4+
[GeneratorPythia8]
5+
config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator/pythia8_PbPb_536tev.cfg

0 commit comments

Comments
 (0)