Skip to content

Commit

Permalink
[DRAFT] PWGLF: Add continuously-increasing-with-Nch xi/om injection
Browse files Browse the repository at this point in the history
  • Loading branch information
@ddobrigk
ddobrigk committed Apr 23, 2024
1 parent 02d071a commit b3dcd18
Show file tree
Hide file tree
Showing 2 changed files with 356 additions and 0 deletions.
6 changes: 6 additions & 0 deletions MC/config/PWGLF/ini/GeneratorLF_Strangeness_PbPb5360_injection.ini
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
[GeneratorExternal]
fileName=${O2DPG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_extraStrangeness.C
funcName=generator_extraStrangeness()

[GeneratorPythia8]
config=${O2DPG_ROOT}/MC/config/common/pythia8/generator/pythia8_hi.cfg
350 changes: 350 additions & 0 deletions MC/config/PWGLF/pythia8/generator_pythia8_extraStrangeness.C
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,350 @@

#if !defined(__CLING__) || defined(__ROOTCLING__)
#include "Pythia8/Pythia.h"
#include "FairGenerator.h"
#include "FairPrimaryGenerator.h"
#include "Generators/GeneratorPythia8.h"
#include "TRandom3.h"
#include "TParticlePDG.h"
#include "TDatabasePDG.h"
#include "TMath.h"
#include <cmath>
using namespace Pythia8;
#endif

// Default pythia8 minimum bias generator
// Please do not change

class GeneratorPythia8ExtraStrangeness : public o2::eventgen::GeneratorPythia8
{
public:
/// Constructor
GeneratorPythia8ExtraStrangeness() {
genMinPt=0.0;
genMaxPt=20.0;
genminY=-1.5;
genmaxY=1.5;
genminEta=-1.5;
genmaxEta=1.5;

pdg=0;
E=0;
px=0;
py=0;
pz=0;
p=0;
y=0;
eta=0;
xProd=0;
yProd=0;
zProd=0;
xProd=0.; yProd=0.; zProd=0.;

fLVHelper = new TLorentzVector();

// Initialize the Xi and Omega spectra to be used at injection time
fSpectraXi = STAR_BlastWave("fSpectraXi", 1.32171, 20);
fSpectraXi ->SetNpx( 1000 );
fSpectraOm = STAR_BlastWave("fSpectraOm", 1.67245, 20);
fSpectraOm ->SetNpx( 1000 );

fSpectraXi->SetParameter(0,1.32171);
fSpectraXi->SetParameter(1,0.6615); //beta-max
fSpectraXi->SetParameter(2,0.0905); //T
fSpectraXi->SetParameter(3,0.7355); //n
fSpectraXi->SetParameter(4,1000); //norm (not relevant)

fSpectraOm->SetParameter(0,1.67245);
fSpectraOm->SetParameter(1,0.6615); //beta-max
fSpectraOm->SetParameter(2,0.0905); //T
fSpectraOm->SetParameter(3,0.7355); //n
fSpectraOm->SetParameter(4,1000); //norm (not relevant)

fLVHelper = new TLorentzVector();
}

/// Destructor
~GeneratorPythia8ExtraStrangeness() = default;

Double_t y2eta(Double_t pt, Double_t mass, Double_t y){
Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
return TMath::ASinH(mt / pt * TMath::SinH(y));
}

/// set 4-momentum
void set4momentum(double input_px, double input_py, double input_pz){
px = input_px;
py = input_py;
pz = input_pz;
E = sqrt( m*m+px*px+py*py+pz*pz );
fourMomentum.px(px);
fourMomentum.py(py);
fourMomentum.pz(pz);
fourMomentum.e(E);
p = sqrt( px*px+py*py+pz*pz );
y = 0.5*log( (E+pz)/(E-pz) );
eta = 0.5*log( (p+pz)/(p-pz) );
}

//__________________________________________________________________
Pythia8::Particle createParticle(){
//std::cout << "createParticle() mass " << m << " pdgCode " << pdg << std::endl;
Pythia8::Particle myparticle;
myparticle.id(pdg);
myparticle.status(11);
myparticle.px(px);
myparticle.py(py);
myparticle.pz(pz);
myparticle.e(E);
myparticle.m(m);
myparticle.xProd(xProd);
myparticle.yProd(yProd);
myparticle.zProd(zProd);

return myparticle;
}

//_________________________________________________________________________________
/// generate uniform eta and uniform momentum
void genSpectraMomentumEtaXi(double minP, double maxP, double minY, double maxY){
// random generator
std::unique_ptr<TRandom3> ranGenerator { new TRandom3() };
ranGenerator->SetSeed(0);

// generate transverse momentum
const double gen_pT = fSpectraXi->GetRandom(minP,maxP);

//Actually could be something else without loss of generality but okay
const double gen_phi = ranGenerator->Uniform(0,2*TMath::Pi());

// sample flat in rapidity, calculate eta
Double_t gen_Y=10, gen_eta=10;

while( gen_eta>genmaxEta || gen_eta<genminEta ){
gen_Y = ranGenerator->Uniform(minY,maxY);
//(Double_t pt, Double_t mass, Double_t y)
gen_eta = y2eta(gen_pT, m, gen_Y);
}

fLVHelper->SetPtEtaPhiM(gen_pT, gen_eta, gen_phi, m);
set4momentum(fLVHelper->Px(),fLVHelper->Py(),fLVHelper->Pz());
}

//_________________________________________________________________________________
/// generate uniform eta and uniform momentum
void genSpectraMomentumEtaOm(double minP, double maxP, double minY, double maxY){
// random generator
std::unique_ptr<TRandom3> ranGenerator { new TRandom3() };
ranGenerator->SetSeed(0);

// generate transverse momentum
const double gen_pT = fSpectraOm->GetRandom(minP,maxP);

//Actually could be something else without loss of generality but okay
const double gen_phi = ranGenerator->Uniform(0,2*TMath::Pi());

// sample flat in rapidity, calculate eta
Double_t gen_Y=10, gen_eta=10;

while( gen_eta>genmaxEta || gen_eta<genminEta ){
gen_Y = ranGenerator->Uniform(minY,maxY);
//(Double_t pt, Double_t mass, Double_t y)
gen_eta = y2eta(gen_pT, m, gen_Y);
}

fLVHelper->SetPtEtaPhiM(gen_pT, gen_eta, gen_phi, m);
set4momentum(fLVHelper->Px(),fLVHelper->Py(),fLVHelper->Pz());
}

Double_t STAR_BlastWave_Func(const Double_t *x, const Double_t *p) {
/* dN/dpt */

Double_t pt = x[0];
Double_t mass = p[0];
Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
Double_t beta_max = p[1];
Double_t temp = p[2];
Double_t n = p[3];
Double_t norm = p[4];

Double_t integral = 0;

// a bit manual but ok...
if(TMath::Abs(mass-1.32171)<0.002){
if (!fIntegrandXi)
fIntegrandXi = new TF1("fIntegrandXi", this, &GeneratorPythia8ExtraStrangeness::STAR_BlastWave_Integrand_Improved, 0., 1., 5, "GeneratorPythia8ExtraStrangeness", "STAR_BlastWave_Integrand_Improved");
fIntegrandXi->SetParameters(mt, pt, beta_max, temp, n);

integral = fIntegrandXi->Integral(0., 1.);
}
if(TMath::Abs(mass-1.67245)<0.002){
if (!fIntegrandOm)
fIntegrandOm = new TF1("fIntegrandOm", this, &GeneratorPythia8ExtraStrangeness::STAR_BlastWave_Integrand_Improved, 0., 1., 5, "GeneratorPythia8ExtraStrangeness", "STAR_BlastWave_Integrand_Improved");
fIntegrandOm->SetParameters(mt, pt, beta_max, temp, n);

integral = fIntegrandOm->Integral(0., 1.);
}
return norm * pt * integral;
}

//___________________________________________________________________
Double_t STAR_BlastWave_Integrand_Improved(const Double_t *x, const Double_t *p) {

/*
x[0] -> r (radius)
p[0] -> mT (transverse mass)
p[1] -> pT (transverse momentum)
p[2] -> beta_max (surface velocity)
p[3] -> T (freezout temperature)
p[4] -> n (velocity profile)
*/

Double_t r = x[0];
Double_t mt = p[0];
Double_t pt = p[1];
Double_t beta_max = p[2];
Double_t temp_1 = 1. / p[3];
Double_t n = p[4];

Double_t beta = beta_max * TMath::Power(r, n);
Double_t rho = TMath::ATanH(beta);
Double_t argI0 = pt * TMath::SinH(rho) * temp_1;
Double_t argK1 = mt * TMath::CosH(rho) * temp_1;
// if (argI0 > 100 || argI0 < -100)
// printf("r=%f, pt=%f, beta_max=%f, temp=%f, n=%f, mt=%f, beta=%f, rho=%f, argI0=%f, argK1=%f\n", r, pt, beta_max, 1. / temp_1, n, mt, beta, rho, argI0, argK1);
return r * mt * TMath::BesselI0(argI0) * TMath::BesselK1(argK1);

}

//___________________________________________________________________
TF1 *STAR_BlastWave(const Char_t *name, Double_t mass,Float_t upperlim, Double_t beta_max = 0.9, Double_t temp = 0.1, Double_t n = 1., Double_t norm = 1.e6) {

//new TF1("fSpectra",this ,&GeneratorPythia8GunPbPb::myLevyPt, 0.0,20,3, "GeneratorPythia8GunPbPb","myLevyPt");
TF1 *fBlastWave = new TF1(name, this, &GeneratorPythia8ExtraStrangeness::STAR_BlastWave_Func, 0., upperlim, 5, "GeneratorPythia8ExtraStrangeness", "STAR_BlastWave_Func");
fBlastWave->SetParameters(mass, beta_max, temp, n, norm);
fBlastWave->SetParNames("mass", "beta_max", "T", "n", "norm");
fBlastWave->FixParameter(0, mass);
fBlastWave->SetParLimits(1, 0.1, 0.9); // don't touch :) adding some 99 youu get floating point exception
fBlastWave->SetParLimits(2, 0.03,1.);//0.05, 1.); // no negative values!! for the following as well
fBlastWave->SetParLimits(3, 0.25,4.5); // was 2.5 // omega-->at limit even moving it to 4.5 but yield same
return fBlastWave;
}


//__________________________________________________________________
Bool_t generateEvent() override {

// Generate PYTHIA event
Bool_t lPythiaOK = kFALSE;
while (!lPythiaOK){
lPythiaOK = mPythia.next();
}

// characterise event
Long_t nParticles = mPythia.event.size();
Long_t nChargedParticlesAtMidRap = 0;
Long_t nPionsAtMidRap = 0;
for ( Long_t j=0; j < nParticles; j++ ) {
Int_t pypid = mPythia.event[j].id();
Float_t pyrap = mPythia.event[j].y();
Float_t pyeta = mPythia.event[j].eta();

// final only
if (!mPythia.event[j].isFinal()) continue;

if ( TMath::Abs(pyrap) < 0.5 && TMath::Abs(pypid)==211 ) nPionsAtMidRap++;
if ( TMath::Abs(pyeta) < 0.5 && TMath::Abs(mPythia.event[j].charge())>1e-5 ) nChargedParticlesAtMidRap++;
}

// now we have the multiplicity

//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// XI ABUNDANCE FIX
// FCN=0.35879 FROM MINOS STATUS=SUCCESSFUL 126 CALLS 634 TOTAL
// EDM=3.7456e-09 STRATEGY= 1 ERROR MATRIX ACCURATE
// EXT PARAMETER STEP FIRST
// NO. NAME VALUE ERROR SIZE DERIVATIVE
// 1 p0 4.74929e-03 3.29248e-04 -3.35914e-06 5.38225e+00
// 2 p1 -4.08255e-03 8.62587e-04 -2.02577e-05 2.45132e+00
// 3 p2 4.76660e+00 1.93593e+00 1.93593e+00 2.70369e-04
// Info in <TCanvas::MakeDefCanvas>: created default TCanvas with name c1
//Adjust relative abundance of multi-strange particles by injecting some
Double_t lExpectedXiToPion = TMath::Max(4.74929e-03 - 4.08255e-03*TMath::Exp(-nChargedParticlesAtMidRap/4.76660e+00) - 0.00211334,0.);
Double_t lExpectedXi = nPionsAtMidRap*lExpectedXiToPion;
Int_t lXiYield = gRandom->Poisson(3*lExpectedXi); //factor 3: fix the rapidity acceptance
m = 1.32171;
pdg = 3312;
cout<<"Adding extra xi: "<<lXiYield<<" (to reach average "<<Form("%.6f",lExpectedXi)<<" at this Nch = "<<nChargedParticlesAtMidRap<<", ratio: "<<Form("%.6f",lExpectedXiToPion)<<")"<<endl;
for(Int_t ii=0; ii<lXiYield; ii++){
pdg *= gRandom->Uniform()>0.5?+1:-1;
xProd=0.0;
yProd=0.0;
zProd=0.0;
genSpectraMomentumEtaXi(genMinPt,genMaxPt,genminY,genmaxY);
Pythia8::Particle lAddedParticle = createParticle();
mPythia.event.append(lAddedParticle);
//lAddedParticles++;
}
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// OMEGA ABUNDANCE FIX
//Adjust relative abundance of multi-strange particles by injecting some
Double_t lExpectedOmegaToPion = TMath::Max(8.55057e-04 - 7.38732e-04*TMath::Exp(-nChargedParticlesAtMidRap/2.40545e+01) - 6.56785e-05,0.);
Double_t lExpectedOmega = nPionsAtMidRap*lExpectedOmegaToPion;
Int_t lOmegaYield = gRandom->Poisson(3*lExpectedOmega); //factor 3: fix the rapidity acceptance
m = 1.67245;
pdg = 3334;
cout<<"Adding extra omegas: "<<lOmegaYield<<" (to reach average "<<Form("%.6f",lExpectedOmega)<<" at this Nch = "<<nChargedParticlesAtMidRap<<", ratio: "<<Form("%.6f",lExpectedOmegaToPion)<<")"<<endl;
for(Int_t ii=0; ii<lOmegaYield; ii++){
pdg *= gRandom->Uniform()>0.5?+1:-1;
xProd=0.0;
yProd=0.0;
zProd=0.0;
genSpectraMomentumEtaOm(genMinPt,genMaxPt,genminY,genmaxY);
Pythia8::Particle lAddedParticle = createParticle();
mPythia.event.append(lAddedParticle);
//lAddedParticles++;
}
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

return true;
}

private:

double genMinPt; /// minimum 3-momentum for generated particles
double genMaxPt; /// maximum 3-momentum for generated particles
double genminY; /// minimum pseudorapidity for generated particles
double genmaxY; /// maximum pseudorapidity for generated particles
double genminEta;
double genmaxEta;

Pythia8::Vec4 fourMomentum; /// four-momentum (px,py,pz,E)
double E; /// energy: sqrt( m*m+px*px+py*py+pz*pz ) [GeV/c]
double m; /// particle mass [GeV/c^2]
int pdg; /// particle pdg code
double px; /// x-component momentum [GeV/c]
double py; /// y-component momentum [GeV/c]
double pz; /// z-component momentum [GeV/c]
double p; /// momentum
double y; /// rapidity
double eta; /// pseudorapidity
double xProd; /// x-coordinate position production vertex [cm]
double yProd; /// y-coordinate position production vertex [cm]
double zProd; /// z-coordinate position production vertex [cm]
///
TF1 *fSpectraXi; /// TF1 to store more realistic shape of spectrum
TF1 *fSpectraOm; /// TF1 to store more realistic shape of spectrum

//BW integrand
TF1 *fIntegrandXi = NULL;
TF1 *fIntegrandOm = NULL;

TLorentzVector *fLVHelper;
};

FairGenerator *generator_extraStrangeness()
{
return new GeneratorPythia8ExtraStrangeness();
}

0 comments on commit b3dcd18

Please sign in to comment.