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 "TF1.h"
8+ #include "TRandom3.h"
9+ //#include "TParticlePDG.h"
10+ //#include "TDatabasePDG.h"
11+
12+ #include <map>
13+ #include <unordered_set>
14+
15+ class GeneratorPythia8ExtraStrangeness : public o2 ::eventgen ::GeneratorPythia8
16+ {
17+ public :
18+ /// default constructor
19+ GeneratorPythia8ExtraStrangeness () = default ;
20+
21+ /// Constructor
22+ GeneratorPythia8ExtraStrangeness (int input_pdg )
23+ {
24+ genMinPt = 0.0 ;
25+ genMaxPt = 12.0 ;
26+ genminY = -1.5 ;
27+ genmaxY = 1.5 ;
28+ genminEta = -1.5 ;
29+ genmaxEta = 1.5 ;
30+
31+ pdg = input_pdg ;
32+ m = 0 ;
33+ E = 0 ;
34+ px = 0 ;
35+ py = 0 ;
36+ pz = 0 ;
37+ p = 0 ;
38+ y = 0 ;
39+ eta = 0 ;
40+ xProd = 0. ; yProd = 0. ; zProd = 0. ;
41+
42+ fLVHelper = std ::make_unique < TLorentzVector > (
43+
44+ randomizePDGsign = false;
45+
46+ fSpectra = makeLevySpectrum ("fSpectra" , genMinPt , genMaxPt , mPythia .particleData .m0 (input_pdg ), 0.30 , 7. );
47+ }
48+
49+ /// randomize the PDG code sign of core particle
50+ void setRandomizePDGsign () { randomizePDGsign = true; }
51+
52+ Double_t y2eta (Double_t pt , Double_t mass , Double_t y ){
53+ Double_t mt = TMath ::Sqrt (mass * mass + pt * pt );
54+ return TMath ::ASinH (mt / pt * TMath ::SinH (y ));
55+ }
56+
57+
58+ /// set mass
59+ double sampleMass (int input_pdg )
60+ {
61+ auto& pd = mPythia .particleData ;
62+ const double mass = pd .mWidth (input_pdg ) > 0. ? pd .mSel (input_pdg ) : pd .m0 (input_pdg );
63+ return mass ;
64+ }
65+
66+ static double myLevyPt (double * pt , double * par )
67+ {
68+ const double lMass = par [0 ];
69+ const double ldNdy = par [1 ];
70+ const double lTemp = par [2 ];
71+ const double lPower = par [3 ];
72+
73+ const double lBigCoef = ((lPower - 1 )* (lPower - 2 )) / (lPower * lTemp * (lPower * lTemp + lMass * (lPower - 2 )));
74+ const double lInPower = 1 + (TMath ::Sqrt (pt [0 ]* pt [0 ]+ lMass * lMass )- lMass ) / (lPower * lTemp );
75+
76+ return ldNdy * pt [0 ] * lBigCoef * TMath ::Power (lInPower , - lPower );
77+ }
78+
79+ /// build a Levy-Tsallis pT spectrum for a particle of given mass
80+ TF1 * makeLevySpectrum (const char * name , double ptMin , double ptMax , double mass , double T , double n , double norm = 1. )
81+ {
82+ TF1 * f = new TF1 (name , myLevyPt , ptMin , ptMax , 4 );
83+ f -> SetParNames ("mass" , "norm" , "T" , "n" );
84+ f -> FixParameter (0 , mass ); // mass [GeV/c^2]
85+ f -> SetParameter (1 , norm ); // normalization (irrelevant for GetRandom)
86+ f -> SetParameter (2 , T ); // slope parameter [GeV]
87+ f -> SetParameter (3 , n ); // power-law exponent
88+ f -> SetNpx (1000 );
89+ return f ;
90+ }
91+
92+ /// set 4-momentum
93+ void set4momentum (double input_px , double input_py , double input_pz ){
94+ px = input_px ;
95+ py = input_py ;
96+ pz = input_pz ;
97+ E = sqrt ( m * m + px * px + py * py + pz * pz );
98+ fourMomentum .px (px );
99+ fourMomentum .py (py );
100+ fourMomentum .pz (pz );
101+ fourMomentum .e (E );
102+ p = sqrt ( px * px + py * py + pz * pz );
103+ y = 0.5 * log ( (E + pz )/(E - pz ) );
104+ eta = 0.5 * log ( (p + pz )/(p - pz ) );
105+ }
106+
107+
108+ //_________________________________________________________________________________
109+ /// generate uniform eta and uniform momentum
110+ void genSpectraMomentumEta (double minPt , double maxPt , double minY , double maxY ){
111+ // random generator
112+ std ::unique_ptr < TRandom3 > ranGenerator { new TRandom3 () };
113+ ranGenerator -> SetSeed (0 );
114+
115+ // generate transverse momentum
116+ //const double gen_pT = ranGenerator->Uniform(minPt, maxPt);
117+ const double gen_pT = fSpectra -> GetRandom (minPt ,maxPt );
118+
119+ //Actually could be something else without loss of generality but okay
120+ const double gen_phi = ranGenerator -> Uniform (0 ,2 * TMath ::Pi ());
121+
122+ // sample flat in rapidity, calculate eta
123+ Double_t gen_Y = 10 , gen_eta = 10 ;
124+
125+ while ( gen_eta > genmaxEta || gen_eta < genminEta ){
126+ gen_Y = ranGenerator -> Uniform (minY ,maxY );
127+ gen_eta = y2eta (gen_pT , m , gen_Y );
128+ }
129+
130+ fLVHelper -> SetPtEtaPhiM (gen_pT , gen_eta , gen_phi , m );
131+ set4momentum (fLVHelper -> Px (),fLVHelper -> Py (),fLVHelper -> Pz ());
132+ }
133+
134+ //__________________________________________________________________
135+ Pythia8 ::Particle createParticle (){
136+ //std::cout << "createParticle() mass " << m << " pdgCode " << pdg << std::endl;
137+ Pythia8 ::Particle myparticle ;
138+ myparticle .id (pdg );
139+ myparticle .status (11 );
140+ myparticle .px (px );
141+ myparticle .py (py );
142+ myparticle .pz (pz );
143+ myparticle .e (E );
144+ myparticle .m (m );
145+ myparticle .xProd (xProd );
146+ myparticle .yProd (yProd );
147+ myparticle .zProd (zProd );
148+
149+ return myparticle ;
150+ }
151+
152+ //__________________________________________________________________
153+ int randomizeSign ()
154+ {
155+ std ::unique_ptr < TRandom3 > gen_random {new TRandom3 (0 )};
156+ const float n = gen_random -> Uniform (-1 , 1 );
157+
158+ return n / abs (n );
159+ }
160+
161+ //__________________________________________________________________
162+ Bool_t generateEvent () override {
163+ // Generate PYTHIA event
164+
165+ mPythia .event .reset ();
166+
167+ /// go to next Pythia event
168+ Bool_t lPythiaOK = kFALSE ;
169+ while (!lPythiaOK ){
170+ lPythiaOK = mPythia .next ();
171+ }
172+
173+ /// reset event
174+ //mPythia.event.reset();
175+
176+ /// create and append the desired particle
177+ m = sampleMass (pdg );
178+ genSpectraMomentumEta (genMinPt , genMaxPt , genminY , genmaxY );
179+
180+ if (randomizePDGsign )
181+ pdg *= randomizeSign ();
182+
183+ Pythia8 ::Particle particle = createParticle ();
184+ mPythia .event .append (particle );
185+ mPythia .moreDecays ();
186+
187+ return true;
188+ }
189+
190+ private :
191+
192+ double genMinPt ; /// minimum 3-momentum for generated particles
193+ double genMaxPt ; /// maximum 3-momentum for generated particles
194+ double genminY ; /// minimum pseudorapidity for generated particles
195+ double genmaxY ; /// maximum pseudorapidity for generated particles
196+ double genminEta ;
197+ double genmaxEta ;
198+
199+ Pythia8 ::Vec4 fourMomentum ; /// four-momentum (px,py,pz,E)
200+ //std::unique_ptr<o2::eventgen::FlowMapper> lutGen;
201+
202+ double E ; /// energy: sqrt( m*m+px*px+py*py+pz*pz ) [GeV/c]
203+ double m ; /// particle mass [GeV/c^2]
204+ int pdg ; /// particle pdg code
205+ double px ; /// x-component momentum [GeV/c]
206+ double py ; /// y-component momentum [GeV/c]
207+ double pz ; /// z-component momentum [GeV/c]
208+ double p ; /// momentum
209+ double y ; /// rapidity
210+ double eta ; /// pseudorapidity
211+ double xProd ; /// x-coordinate position production vertex [cm]
212+ double yProd ; /// y-coordinate position production vertex [cm]
213+ double zProd ; /// z-coordinate position production vertex [cm]
214+
215+ bool randomizePDGsign ; /// bool to randomize the PDG code of the core particle
216+
217+ TF1 * fSpectra = nullptr ; /// TF1 to store more realistic shape of spectrum
218+ std ::unique_ptr < TLorentzVector > fLVHelper ;
219+ };
220+
221+ FairGenerator * generator_Lstar ()
222+ {
223+ auto myGen = new GeneratorPythia8ExtraStrangeness (3124 );
224+ myGen -> setRandomizePDGsign (); // randomization of PDG switched on
225+ return myGen ;
226+ }
0 commit comments