Skip to content

Commit f139a34

Browse files
ddobrigkddobrigk
committed
Synthe-flow: move config histos to ccdb
1 parent 0fe6797 commit f139a34

File tree

3 files changed

+365
-7
lines changed

3 files changed

+365
-7
lines changed

‎MC/config/PWGLF/ini/GeneratorLF_SyntheFlowXi.ini‎

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

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

Lines changed: 20 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,9 @@
77
#include "TRandom3.h"
88
#include "TParticlePDG.h"
99
#include "TDatabasePDG.h"
10+
#include "CCDB/BasicCCDBManager.h"
11+
#include "TH1F.h"
12+
#include "TH1D.h"
1013

1114
#include <map>
1215
#include <unordered_set>
@@ -19,15 +22,25 @@ public:
1922
lutGen = new o2::eventgen::FlowMapper();
2023

2124
// -------- CONFIGURE SYNTHETIC FLOW ------------
22-
// specify a v2 vs pT here
23-
TFile *filehep = new TFile("/Users/daviddc/Downloads/HEPData-ins1116150-v1-Table_1.root", "READ");
24-
TH1D *hv = (TH1D*) filehep->Get("Table 1/Hist1D_y6");
25-
26-
TFile *fileEcc = new TFile("/Users/daviddc/Downloads/eccentricityvsb.root", "READ");
27-
TH1D *hEccentricities = (TH1D*) fileEcc->Get("hEccentricities");
25+
// establish connection to ccdb
26+
o2::ccdb::CcdbApi ccdb_api;
27+
ccdb_api.init("https://alice-ccdb.cern.ch");
28+
29+
// config was placed at midpoint of run 544122, retrieve that
30+
std::map<string, string> metadataRCT, headers;
31+
headers = ccdb_api.retrieveHeaders("RCT/Info/RunInformation/544122", metadataRCT, -1);
32+
int64_t tsSOR = atol(headers["SOR"].c_str());
33+
int64_t tsEOR = atol(headers["EOR"].c_str());
34+
int64_t midRun = 0.5*tsSOR+0.5*tsEOR;
35+
36+
map<string, string> metadata; // can be empty
37+
auto list = ccdb_api.retrieveFromTFileAny<TList>("Users/d/ddobrigk/syntheflow", metadata, midRun);
38+
39+
TH1D *hv2vspT = (TH1D*) list->FindObject("hFlowVsPt_ins1116150_v1_Table_1");
40+
TH1D *heccvsb = (TH1D*) list->FindObject("hEccentricityVsB");
2841

2942
cout<<"Generating LUT for flow test"<<endl;
30-
lutGen->CreateLUT(hv, hEccentricities);
43+
lutGen->CreateLUT(hv2vspT, heccvsb);
3144
cout<<"Finished creating LUT!"<<endl;
3245
// -------- END CONFIGURE SYNTHETIC FLOW ------------
3346
}

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

Lines changed: 339 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,339 @@
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+
#include "CCDB/BasicCCDBManager.h"
11+
#include "TH1F.h"
12+
#include "TH1D.h"
13+
14+
#include <map>
15+
#include <unordered_set>
16+
17+
class GeneratorPythia8SyntheFlowXi : public o2::eventgen::GeneratorPythia8
18+
{
19+
public:
20+
/// Constructor
21+
GeneratorPythia8SyntheFlowXi() {
22+
lutGen = new o2::eventgen::FlowMapper();
23+
24+
// -------- CONFIGURE SYNTHETIC FLOW ------------
25+
// establish connection to ccdb
26+
o2::ccdb::CcdbApi ccdb_api;
27+
ccdb_api.init("https://alice-ccdb.cern.ch");
28+
29+
// config was placed at midpoint of run 544122, retrieve that
30+
std::map<string, string> metadataRCT, headers;
31+
headers = ccdb_api.retrieveHeaders("RCT/Info/RunInformation/544122", metadataRCT, -1);
32+
int64_t tsSOR = atol(headers["SOR"].c_str());
33+
int64_t tsEOR = atol(headers["EOR"].c_str());
34+
int64_t midRun = 0.5*tsSOR+0.5*tsEOR;
35+
36+
map<string, string> metadata; // can be empty
37+
auto list = ccdb_api.retrieveFromTFileAny<TList>("Users/d/ddobrigk/syntheflow", metadata, midRun);
38+
39+
TH1D *hv2vspT = (TH1D*) list->FindObject("hFlowVsPt_ins1116150_v1_Table_1");
40+
TH1D *heccvsb = (TH1D*) list->FindObject("hEccentricityVsB");
41+
42+
cout<<"Generating LUT for flow test"<<endl;
43+
lutGen->CreateLUT(hv2vspT, heccvsb);
44+
cout<<"Finished creating LUT!"<<endl;
45+
// -------- END CONFIGURE SYNTHETIC FLOW ------------
46+
47+
genMinPt=0.0;
48+
genMaxPt=20.0;
49+
genminY=-1.5;
50+
genmaxY=1.5;
51+
genminEta=-1.5;
52+
genmaxEta=1.5;
53+
54+
pdg=0;
55+
E=0;
56+
px=0;
57+
py=0;
58+
pz=0;
59+
p=0;
60+
y=0;
61+
eta=0;
62+
xProd=0;
63+
yProd=0;
64+
zProd=0;
65+
xProd=0.; yProd=0.; zProd=0.;
66+
67+
fLVHelper = std::make_unique<TLorentzVector>();
68+
69+
fSpectrumXi = std::make_unique<TF1>("fSpectrumXi", this, &GeneratorPythia8ExtraStrangeness::boltzPlusPower, 0., genMaxPt, 5, "GeneratorPythia8ExtraStrangeness", "boltzPlusPower");
70+
71+
fSpectrumXi->FixParameter(0, 1.32171);
72+
fSpectrumXi->FixParameter(1, 4.84e-1);
73+
fSpectrumXi->FixParameter(2, 111.9);
74+
fSpectrumXi->FixParameter(3, -2.56511e+00);
75+
fSpectrumXi->FixParameter(4, 1.14011e-04);
76+
77+
fSpectrumOm = std::make_unique<TF1>("fSpectrumOm", this, &GeneratorPythia8ExtraStrangeness::boltzPlusPower, 0., genMaxPt, 5, "GeneratorPythia8ExtraStrangeness", "boltzPlusPower");
78+
79+
fSpectrumOm->FixParameter(0, 1.67245e+00);
80+
fSpectrumOm->FixParameter(1, 5.18174e-01);
81+
fSpectrumOm->FixParameter(2, 1.73747e+01);
82+
fSpectrumOm->FixParameter(3, -2.56681e+00);
83+
fSpectrumOm->FixParameter(4, 1.87513e-04);
84+
}
85+
86+
/// Destructor
87+
~GeneratorPythia8SyntheFlowXi() = default;
88+
89+
Double_t y2eta(Double_t pt, Double_t mass, Double_t y){
90+
Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
91+
return TMath::ASinH(mt / pt * TMath::SinH(y));
92+
}
93+
94+
/// set 4-momentum
95+
void set4momentum(double input_px, double input_py, double input_pz){
96+
px = input_px;
97+
py = input_py;
98+
pz = input_pz;
99+
E = sqrt( m*m+px*px+py*py+pz*pz );
100+
fourMomentum.px(px);
101+
fourMomentum.py(py);
102+
fourMomentum.pz(pz);
103+
fourMomentum.e(E);
104+
p = sqrt( px*px+py*py+pz*pz );
105+
y = 0.5*log( (E+pz)/(E-pz) );
106+
eta = 0.5*log( (p+pz)/(p-pz) );
107+
}
108+
109+
//__________________________________________________________________
110+
Pythia8::Particle createParticle(){
111+
//std::cout << "createParticle() mass " << m << " pdgCode " << pdg << std::endl;
112+
Pythia8::Particle myparticle;
113+
myparticle.id(pdg);
114+
myparticle.status(11);
115+
myparticle.px(px);
116+
myparticle.py(py);
117+
myparticle.pz(pz);
118+
myparticle.e(E);
119+
myparticle.m(m);
120+
myparticle.xProd(xProd);
121+
myparticle.yProd(yProd);
122+
myparticle.zProd(zProd);
123+
124+
return myparticle;
125+
}
126+
127+
//_________________________________________________________________________________
128+
/// generate uniform eta and uniform momentum
129+
void genSpectraMomentumEtaXi(double minP, double maxP, double minY, double maxY){
130+
// random generator
131+
std::unique_ptr<TRandom3> ranGenerator { new TRandom3() };
132+
ranGenerator->SetSeed(0);
133+
134+
// generate transverse momentum
135+
const double gen_pT = fSpectrumXi->GetRandom(genMinPt,genMaxPt);
136+
137+
//Actually could be something else without loss of generality but okay
138+
const double gen_phi = ranGenerator->Uniform(0,2*TMath::Pi());
139+
140+
// sample flat in rapidity, calculate eta
141+
Double_t gen_Y=10, gen_eta=10;
142+
143+
while( gen_eta>genmaxEta || gen_eta<genminEta ){
144+
gen_Y = ranGenerator->Uniform(minY,maxY);
145+
//(Double_t pt, Double_t mass, Double_t y)
146+
gen_eta = y2eta(gen_pT, m, gen_Y);
147+
}
148+
149+
fLVHelper->SetPtEtaPhiM(gen_pT, gen_eta, gen_phi, m);
150+
set4momentum(fLVHelper->Px(),fLVHelper->Py(),fLVHelper->Pz());
151+
}
152+
153+
//_________________________________________________________________________________
154+
/// generate uniform eta and uniform momentum
155+
void genSpectraMomentumEtaOm(double minP, double maxP, double minY, double maxY){
156+
// random generator
157+
std::unique_ptr<TRandom3> ranGenerator { new TRandom3() };
158+
ranGenerator->SetSeed(0);
159+
160+
// generate transverse momentum
161+
const double gen_pT = fSpectrumOm->GetRandom(genMinPt,genMaxPt);
162+
163+
//Actually could be something else without loss of generality but okay
164+
const double gen_phi = ranGenerator->Uniform(0,2*TMath::Pi());
165+
166+
// sample flat in rapidity, calculate eta
167+
Double_t gen_Y=10, gen_eta=10;
168+
169+
while( gen_eta>genmaxEta || gen_eta<genminEta ){
170+
gen_Y = ranGenerator->Uniform(minY,maxY);
171+
//(Double_t pt, Double_t mass, Double_t y)
172+
gen_eta = y2eta(gen_pT, m, gen_Y);
173+
}
174+
175+
fLVHelper->SetPtEtaPhiM(gen_pT, gen_eta, gen_phi, m);
176+
set4momentum(fLVHelper->Px(),fLVHelper->Py(),fLVHelper->Pz());
177+
}
178+
179+
//_________________________________________________________________________________
180+
/// shape function
181+
Double_t boltzPlusPower(const Double_t *x, const Double_t *p)
182+
{
183+
// a plain parametrization. not meant to be physics worthy.
184+
// adjusted to match preliminary 5 TeV shape.
185+
186+
Double_t pt = x[0];
187+
Double_t mass = p[0];
188+
Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
189+
Double_t T = p[1];
190+
Double_t norm = p[2];
191+
192+
Double_t lowptpart = mt * TMath::Exp(-mt / T);
193+
Double_t highptpart = p[4]*TMath::Power(x[0], p[3]);
194+
195+
Double_t mixup = 1./(1.+TMath::Exp((x[0]-4.5)/.1));
196+
197+
//return pt * norm * (mixup * mt * TMath::Exp(-mt / T) + (1.-mixup)*highptpart) ;
198+
return pt * norm * (mt * TMath::Exp(-mt / T) + (1.-mixup)*highptpart) ;
199+
}
200+
201+
//__________________________________________________________________
202+
Bool_t generateEvent() override {
203+
204+
// Generate PYTHIA event
205+
Bool_t lPythiaOK = kFALSE;
206+
while (!lPythiaOK){
207+
lPythiaOK = mPythia.next();
208+
}
209+
210+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
211+
// add extra xi/omega content
212+
// characterise event
213+
Long_t nParticles = mPythia.event.size();
214+
Long_t nChargedParticlesAtMidRap = 0;
215+
Long_t nPionsAtMidRap = 0;
216+
for ( Long_t j=0; j < nParticles; j++ ) {
217+
Int_t pypid = mPythia.event[j].id();
218+
Float_t pyrap = mPythia.event[j].y();
219+
Float_t pyeta = mPythia.event[j].eta();
220+
221+
// final only
222+
if (!mPythia.event[j].isFinal()) continue;
223+
224+
if ( TMath::Abs(pyrap) < 0.5 && TMath::Abs(pypid)==211 ) nPionsAtMidRap++;
225+
if ( TMath::Abs(pyeta) < 0.5 && TMath::Abs(mPythia.event[j].charge())>1e-5 ) nChargedParticlesAtMidRap++;
226+
}
227+
228+
// now we have the multiplicity
229+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
230+
// XI ABUNDANCE FIX
231+
//Adjust relative abundance of multi-strange particles by injecting some
232+
Double_t lExpectedXiToPion = TMath::Max(4.74929e-03 - 4.08255e-03*TMath::Exp(-nChargedParticlesAtMidRap/4.76660e+00) - 0.00211334,0.);
233+
Double_t lExpectedXi = 5.0*nPionsAtMidRap*lExpectedXiToPion; // extra rich, factor 5
234+
Int_t lXiYield = gRandom->Poisson(3*lExpectedXi); //factor 3: fix the rapidity acceptance
235+
m = 1.32171;
236+
pdg = 3312;
237+
cout<<"Adding extra xi: "<<lXiYield<<" (to reach average "<<Form("%.6f",lExpectedXi)<<" at this Nch = "<<nChargedParticlesAtMidRap<<", ratio: "<<Form("%.6f",lExpectedXiToPion)<<")"<<endl;
238+
for(Int_t ii=0; ii<lXiYield; ii++){
239+
pdg *= gRandom->Uniform()>0.5?+1:-1;
240+
xProd=0.0;
241+
yProd=0.0;
242+
zProd=0.0;
243+
genSpectraMomentumEtaXi(genMinPt,genMaxPt,genminY,genmaxY);
244+
Pythia8::Particle lAddedParticle = createParticle();
245+
mPythia.event.append(lAddedParticle);
246+
//lAddedParticles++;
247+
}
248+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
249+
250+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
251+
// OMEGA ABUNDANCE FIX
252+
//Adjust relative abundance of multi-strange particles by injecting some
253+
Double_t lExpectedOmegaToPion = TMath::Max(8.55057e-04 - 7.38732e-04*TMath::Exp(-nChargedParticlesAtMidRap/2.40545e+01) - 6.56785e-05,0.);
254+
Double_t lExpectedOmega = 5.0*nPionsAtMidRap*lExpectedOmegaToPion; // extra rich, factor 5
255+
Int_t lOmegaYield = gRandom->Poisson(3*lExpectedOmega); //factor 3: fix the rapidity acceptance
256+
m = 1.67245;
257+
pdg = 3334;
258+
cout<<"Adding extra omegas: "<<lOmegaYield<<" (to reach average "<<Form("%.6f",lExpectedOmega)<<" at this Nch = "<<nChargedParticlesAtMidRap<<", ratio: "<<Form("%.6f",lExpectedOmegaToPion)<<")"<<endl;
259+
for(Int_t ii=0; ii<lOmegaYield; ii++){
260+
pdg *= gRandom->Uniform()>0.5?+1:-1;
261+
xProd=0.0;
262+
yProd=0.0;
263+
zProd=0.0;
264+
genSpectraMomentumEtaOm(genMinPt,genMaxPt,genminY,genmaxY);
265+
Pythia8::Particle lAddedParticle = createParticle();
266+
mPythia.event.append(lAddedParticle);
267+
//lAddedParticles++;
268+
}
269+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
270+
271+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
272+
// loop over the entire event record and rotate all particles
273+
// synthetic flow exercise
274+
// first: get event plane
275+
float eventPlaneAngle = mPythia.info.hiInfo->phi();
276+
float impactParameter = mPythia.info.hiInfo->b();
277+
278+
for ( Long_t j=0; j < mPythia.event.size(); j++ ) {
279+
float pyphi = mPythia.event[j].phi();
280+
float pypT = mPythia.event[j].pT();
281+
282+
// calculate delta with EP
283+
float deltaPhiEP = pyphi - eventPlaneAngle;
284+
float shift = 0.0;
285+
while(deltaPhiEP<0.0){
286+
deltaPhiEP += 2*TMath::Pi();
287+
shift += 2*TMath::Pi();
288+
}
289+
while(deltaPhiEP>2*TMath::Pi()){
290+
deltaPhiEP -= 2*TMath::Pi();
291+
shift -= 2*TMath::Pi();
292+
}
293+
float newDeltaPhiEP = lutGen->MapPhi(deltaPhiEP, impactParameter, pypT);
294+
float pyphiNew = newDeltaPhiEP - shift + eventPlaneAngle;
295+
296+
if(pyphiNew>TMath::Pi())
297+
pyphiNew -= 2.0*TMath::Pi();
298+
if(pyphiNew<-TMath::Pi())
299+
pyphiNew += 2.0*TMath::Pi();
300+
mPythia.event[j].rot(0.0, pyphiNew-pyphi);
301+
}
302+
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
303+
304+
return true;
305+
}
306+
307+
private:
308+
double genMinPt; /// minimum 3-momentum for generated particles
309+
double genMaxPt; /// maximum 3-momentum for generated particles
310+
double genminY; /// minimum pseudorapidity for generated particles
311+
double genmaxY; /// maximum pseudorapidity for generated particles
312+
double genminEta;
313+
double genmaxEta;
314+
315+
Pythia8::Vec4 fourMomentum; /// four-momentum (px,py,pz,E)
316+
317+
double E; /// energy: sqrt( m*m+px*px+py*py+pz*pz ) [GeV/c]
318+
double m; /// particle mass [GeV/c^2]
319+
int pdg; /// particle pdg code
320+
double px; /// x-component momentum [GeV/c]
321+
double py; /// y-component momentum [GeV/c]
322+
double pz; /// z-component momentum [GeV/c]
323+
double p; /// momentum
324+
double y; /// rapidity
325+
double eta; /// pseudorapidity
326+
double xProd; /// x-coordinate position production vertex [cm]
327+
double yProd; /// y-coordinate position production vertex [cm]
328+
double zProd; /// z-coordinate position production vertex [cm]
329+
330+
std::unique_ptr<TLorentzVector> fLVHelper;
331+
std::unique_ptr<TF1> fSpectrumXi;
332+
std::unique_ptr<TF1> fSpectrumOm;
333+
o2::eventgen::FlowMapper *lutGen; // for mapping phi angles
334+
};
335+
336+
FairGenerator *generator_syntheFlowXi()
337+
{
338+
return new GeneratorPythia8SyntheFlowXi();
339+
}

0 commit comments

Comments
 (0)