Skip to content

Commit d853a91

Browse files
ddobrigkddobrigk
authored
LF: add multi-strange enriched generator (#1722)
1 parent 99e5525 commit d853a91

File tree

3 files changed

+277
-0
lines changed

3 files changed

+277
-0
lines changed

‎MC/config/PWGLF/ini/GeneratorLF_Strangeness_PbPb5360_injection.ini‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
[GeneratorExternal]
2+
fileName=${O2DPG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_extraStrangeness.C
3+
funcName=generator_extraStrangeness()
4+
5+
[GeneratorPythia8]
6+
config=${O2DPG_ROOT}/MC/config/common/pythia8/generator/pythia8_hi.cfg

‎MC/config/PWGLF/ini/tests/GeneratorLF_Strangeness_PbPb5360_injection.C‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
int External()
2+
{
3+
std::string path{"o2sim_Kine.root"};
4+
5+
TFile file(path.c_str(), "READ");
6+
if (file.IsZombie())
7+
{
8+
std::cerr << "Cannot open ROOT file " << path << "\n";
9+
return 1;
10+
}
11+
12+
auto tree = (TTree *)file.Get("o2sim");
13+
if (!tree)
14+
{
15+
std::cerr << "Cannot find tree o2sim in file " << path << "\n";
16+
return 1;
17+
}
18+
std::vector<o2::MCTrack> *tracks{};
19+
tree->SetBranchAddress("MCTrack", &tracks);
20+
21+
auto nEvents = tree->GetEntries();
22+
if (nEvents < 1)
23+
{
24+
std::cerr << "No events actually generated: not OK!";
25+
return 1;
26+
}
27+
return 0;
28+
}

‎MC/config/PWGLF/pythia8/generator_pythia8_extraStrangeness.C‎

Lines changed: 243 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,243 @@
1+
2+
#include "Pythia8/Pythia.h"
3+
#include "Pythia8/HeavyIons.h"
4+
#include "FairGenerator.h"
5+
#include "FairPrimaryGenerator.h"
6+
#include "Generators/GeneratorPythia8.h"
7+
#include "TRandom3.h"
8+
#include "TParticlePDG.h"
9+
#include "TDatabasePDG.h"
10+
11+
#include <map>
12+
#include <unordered_set>
13+
14+
class GeneratorPythia8ExtraStrangeness : public o2::eventgen::GeneratorPythia8
15+
{
16+
public:
17+
/// Constructor
18+
GeneratorPythia8ExtraStrangeness() {
19+
genMinPt=0.0;
20+
genMaxPt=20.0;
21+
genminY=-1.5;
22+
genmaxY=1.5;
23+
genminEta=-1.5;
24+
genmaxEta=1.5;
25+
26+
pdg=0;
27+
E=0;
28+
px=0;
29+
py=0;
30+
pz=0;
31+
p=0;
32+
y=0;
33+
eta=0;
34+
xProd=0;
35+
yProd=0;
36+
zProd=0;
37+
xProd=0.; yProd=0.; zProd=0.;
38+
39+
fLVHelper = std::make_unique<TLorentzVector>();
40+
}
41+
42+
Double_t y2eta(Double_t pt, Double_t mass, Double_t y){
43+
Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
44+
return TMath::ASinH(mt / pt * TMath::SinH(y));
45+
}
46+
47+
/// set 4-momentum
48+
void set4momentum(double input_px, double input_py, double input_pz){
49+
px = input_px;
50+
py = input_py;
51+
pz = input_pz;
52+
E = sqrt( m*m+px*px+py*py+pz*pz );
53+
fourMomentum.px(px);
54+
fourMomentum.py(py);
55+
fourMomentum.pz(pz);
56+
fourMomentum.e(E);
57+
p = sqrt( px*px+py*py+pz*pz );
58+
y = 0.5*log( (E+pz)/(E-pz) );
59+
eta = 0.5*log( (p+pz)/(p-pz) );
60+
}
61+
62+
//__________________________________________________________________
63+
Pythia8::Particle createParticle(){
64+
//std::cout << "createParticle() mass " << m << " pdgCode " << pdg << std::endl;
65+
Pythia8::Particle myparticle;
66+
myparticle.id(pdg);
67+
myparticle.status(11);
68+
myparticle.px(px);
69+
myparticle.py(py);
70+
myparticle.pz(pz);
71+
myparticle.e(E);
72+
myparticle.m(m);
73+
myparticle.xProd(xProd);
74+
myparticle.yProd(yProd);
75+
myparticle.zProd(zProd);
76+
77+
return myparticle;
78+
}
79+
80+
//_________________________________________________________________________________
81+
/// generate uniform eta and uniform momentum
82+
void genSpectraMomentumEtaXi(double minP, double maxP, double minY, double maxY){
83+
// random generator
84+
std::unique_ptr<TRandom3> ranGenerator { new TRandom3() };
85+
ranGenerator->SetSeed(0);
86+
87+
// generate transverse momentum
88+
const double gen_pT = ranGenerator->Uniform(0,5);
89+
90+
//Actually could be something else without loss of generality but okay
91+
const double gen_phi = ranGenerator->Uniform(0,2*TMath::Pi());
92+
93+
// sample flat in rapidity, calculate eta
94+
Double_t gen_Y=10, gen_eta=10;
95+
96+
while( gen_eta>genmaxEta || gen_eta<genminEta ){
97+
gen_Y = ranGenerator->Uniform(minY,maxY);
98+
//(Double_t pt, Double_t mass, Double_t y)
99+
gen_eta = y2eta(gen_pT, m, gen_Y);
100+
}
101+
102+
fLVHelper->SetPtEtaPhiM(gen_pT, gen_eta, gen_phi, m);
103+
set4momentum(fLVHelper->Px(),fLVHelper->Py(),fLVHelper->Pz());
104+
}
105+
106+
//_________________________________________________________________________________
107+
/// generate uniform eta and uniform momentum
108+
void genSpectraMomentumEtaOm(double minP, double maxP, double minY, double maxY){
109+
// random generator
110+
std::unique_ptr<TRandom3> ranGenerator { new TRandom3() };
111+
ranGenerator->SetSeed(0);
112+
113+
// generate transverse momentum
114+
const double gen_pT = ranGenerator->Uniform(0,5);
115+
116+
//Actually could be something else without loss of generality but okay
117+
const double gen_phi = ranGenerator->Uniform(0,2*TMath::Pi());
118+
119+
// sample flat in rapidity, calculate eta
120+
Double_t gen_Y=10, gen_eta=10;
121+
122+
while( gen_eta>genmaxEta || gen_eta<genminEta ){
123+
gen_Y = ranGenerator->Uniform(minY,maxY);
124+
//(Double_t pt, Double_t mass, Double_t y)
125+
gen_eta = y2eta(gen_pT, m, gen_Y);
126+
}
127+
128+
fLVHelper->SetPtEtaPhiM(gen_pT, gen_eta, gen_phi, m);
129+
set4momentum(fLVHelper->Px(),fLVHelper->Py(),fLVHelper->Pz());
130+
}
131+
132+
//__________________________________________________________________
133+
Bool_t generateEvent() override {
134+
135+
// Generate PYTHIA event
136+
Bool_t lPythiaOK = kFALSE;
137+
while (!lPythiaOK){
138+
lPythiaOK = mPythia.next();
139+
//if(TMath::Abs(mPythia.info.hiInfo->b()-8.5)>0.5) lPythiaOK = false; // select peripheral for flow
140+
}
141+
142+
// characterise event
143+
Long_t nParticles = mPythia.event.size();
144+
Long_t nChargedParticlesAtMidRap = 0;
145+
Long_t nPionsAtMidRap = 0;
146+
for ( Long_t j=0; j < nParticles; j++ ) {
147+
Int_t pypid = mPythia.event[j].id();
148+
Float_t pyrap = mPythia.event[j].y();
149+
Float_t pyeta = mPythia.event[j].eta();
150+
151+
// final only
152+
if (!mPythia.event[j].isFinal()) continue;
153+
154+
if ( TMath::Abs(pyrap) < 0.5 && TMath::Abs(pypid)==211 ) nPionsAtMidRap++;
155+
if ( TMath::Abs(pyeta) < 0.5 && TMath::Abs(mPythia.event[j].charge())>1e-5 ) nChargedParticlesAtMidRap++;
156+
}
157+
158+
// now we have the multiplicity
159+
160+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
161+
// XI ABUNDANCE FIX
162+
// FCN=0.35879 FROM MINOS STATUS=SUCCESSFUL 126 CALLS 634 TOTAL
163+
// EDM=3.7456e-09 STRATEGY= 1 ERROR MATRIX ACCURATE
164+
// EXT PARAMETER STEP FIRST
165+
// NO. NAME VALUE ERROR SIZE DERIVATIVE
166+
// 1 p0 4.74929e-03 3.29248e-04 -3.35914e-06 5.38225e+00
167+
// 2 p1 -4.08255e-03 8.62587e-04 -2.02577e-05 2.45132e+00
168+
// 3 p2 4.76660e+00 1.93593e+00 1.93593e+00 2.70369e-04
169+
// Info in <TCanvas::MakeDefCanvas>: created default TCanvas with name c1
170+
//Adjust relative abundance of multi-strange particles by injecting some
171+
Double_t lExpectedXiToPion = TMath::Max(4.74929e-03 - 4.08255e-03*TMath::Exp(-nChargedParticlesAtMidRap/4.76660e+00) - 0.00211334,0.);
172+
Double_t lExpectedXi = 5.0*nPionsAtMidRap*lExpectedXiToPion; // extra rich, factor 5
173+
Int_t lXiYield = gRandom->Poisson(3*lExpectedXi); //factor 3: fix the rapidity acceptance
174+
m = 1.32171;
175+
pdg = 3312;
176+
cout<<"Adding extra xi: "<<lXiYield<<" (to reach average "<<Form("%.6f",lExpectedXi)<<" at this Nch = "<<nChargedParticlesAtMidRap<<", ratio: "<<Form("%.6f",lExpectedXiToPion)<<")"<<endl;
177+
for(Int_t ii=0; ii<lXiYield; ii++){
178+
pdg *= gRandom->Uniform()>0.5?+1:-1;
179+
xProd=0.0;
180+
yProd=0.0;
181+
zProd=0.0;
182+
genSpectraMomentumEtaXi(genMinPt,genMaxPt,genminY,genmaxY);
183+
Pythia8::Particle lAddedParticle = createParticle();
184+
mPythia.event.append(lAddedParticle);
185+
//lAddedParticles++;
186+
}
187+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
188+
189+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
190+
// OMEGA ABUNDANCE FIX
191+
//Adjust relative abundance of multi-strange particles by injecting some
192+
Double_t lExpectedOmegaToPion = TMath::Max(8.55057e-04 - 7.38732e-04*TMath::Exp(-nChargedParticlesAtMidRap/2.40545e+01) - 6.56785e-05,0.);
193+
Double_t lExpectedOmega = 5.0*nPionsAtMidRap*lExpectedOmegaToPion; // extra rich, factor 5
194+
Int_t lOmegaYield = gRandom->Poisson(3*lExpectedOmega); //factor 3: fix the rapidity acceptance
195+
m = 1.67245;
196+
pdg = 3334;
197+
cout<<"Adding extra omegas: "<<lOmegaYield<<" (to reach average "<<Form("%.6f",lExpectedOmega)<<" at this Nch = "<<nChargedParticlesAtMidRap<<", ratio: "<<Form("%.6f",lExpectedOmegaToPion)<<")"<<endl;
198+
for(Int_t ii=0; ii<lOmegaYield; ii++){
199+
pdg *= gRandom->Uniform()>0.5?+1:-1;
200+
xProd=0.0;
201+
yProd=0.0;
202+
zProd=0.0;
203+
genSpectraMomentumEtaOm(genMinPt,genMaxPt,genminY,genmaxY);
204+
Pythia8::Particle lAddedParticle = createParticle();
205+
mPythia.event.append(lAddedParticle);
206+
//lAddedParticles++;
207+
}
208+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
209+
210+
return true;
211+
}
212+
213+
private:
214+
215+
double genMinPt; /// minimum 3-momentum for generated particles
216+
double genMaxPt; /// maximum 3-momentum for generated particles
217+
double genminY; /// minimum pseudorapidity for generated particles
218+
double genmaxY; /// maximum pseudorapidity for generated particles
219+
double genminEta;
220+
double genmaxEta;
221+
222+
Pythia8::Vec4 fourMomentum; /// four-momentum (px,py,pz,E)
223+
224+
double E; /// energy: sqrt( m*m+px*px+py*py+pz*pz ) [GeV/c]
225+
double m; /// particle mass [GeV/c^2]
226+
int pdg; /// particle pdg code
227+
double px; /// x-component momentum [GeV/c]
228+
double py; /// y-component momentum [GeV/c]
229+
double pz; /// z-component momentum [GeV/c]
230+
double p; /// momentum
231+
double y; /// rapidity
232+
double eta; /// pseudorapidity
233+
double xProd; /// x-coordinate position production vertex [cm]
234+
double yProd; /// y-coordinate position production vertex [cm]
235+
double zProd; /// z-coordinate position production vertex [cm]
236+
237+
std::unique_ptr<TLorentzVector> fLVHelper;
238+
};
239+
240+
FairGenerator *generator_extraStrangeness()
241+
{
242+
return new GeneratorPythia8ExtraStrangeness();
243+
}

0 commit comments

Comments
 (0)