@@ -26,18 +26,22 @@ using namespace Pythia8;
2626// #include "SimulationDataFormat/MCEventHeader.h"
2727
2828// Input to simulation:
29- // inputFilePathName file is expected to be a json file with the structure like so:
29+ // inputFilePathName file is expected to be a json file with the structure like
30+ // so:
3031// {
3132// "simLog": false,
3233// "sglGenRAA": 0.45,
3334// "sglGenTAA": 23.38,
34- // "sglCutoffSteepNess ": 15,
35+ // "sglCutoffSteepness ": 15,
3536// "sglCutoffAbscissa": 10,
3637// "fallSpecterSlopeLog": 10.833274,
3738// "fallSpecterAffinePowerConstantTerm": -5.476804,
3839// "fallSpecterAffinePowerSlope": 0.001110,
3940// "bkgAveragePt": 0.670,
40- // "collTotalMultWithBkg": 2000
41+ // "collMultPowerLawAmplitude": ?,
42+ // "collMultPowerLawExponent": ?,
43+ // "collMultMin": ?,
44+ // "collMultMax": ?
4145// }
4246// can be uploaded to grid using for example: alien.py cp
4347// file:/local/path/parametrisedModel_PbPb_5p36TeV_cent0010.json
@@ -46,36 +50,40 @@ using namespace Pythia8;
4650class GeneratorParametrisedJetModel : public o2 ::eventgen ::GeneratorPythia8 {
4751public :
4852 /// constructor
49- GeneratorParametrisedJetModel (std ::string inputSimParametersPath , std :: string inputSimParametersFileName , bool generateUE , bool generateSignal = true) : mInputSimParametersPath { inputSimParametersPath }, mInputSimParametersFileName {inputSimParametersFileName }, mGenerateUE {generateUE }, mGenerateSignal {generateSignal } {
53+ GeneratorParametrisedJetModel (std ::string inputSimParametersFileName , bool generateUE , bool generateSignal = true) : mInputSimParametersFileName {inputSimParametersFileName }, mGenerateUE {generateUE }, mGenerateSignal {generateSignal } {
5054
5155 if (!generateUE && !generateSignal ) {
5256 cout << "Both the UE and the signal generation should not be off at the same time" << endl ;
5357 exit (1 );
5458 }
5559
56- std ::string inputFilePathName = "alien://" + inputSimParametersPath + inputSimParametersFileName ;
57- if (!gGrid ) {
58- TGrid ::Connect ("alien://" );
59- if (!gGrid ) {
60- LOG (fatal ) << "AliEn connection failed, check token." ;
61- exit (1 );
62- }
63- }
64- // fetch and copy the .json file to the sim work directory
65- std ::string outputPath = "./" ;
66- TString aliencp = Form ("alien_cp alien://%s%s file:%s%s" , inputSimParametersPath .c_str (),
67- inputSimParametersFileName .c_str (), outputPath .c_str (),
68- inputSimParametersFileName .c_str ()); // an internal operation in ROOT that discards the
69- // JSON file directly (smaller than 300 bytes and if
70- // a bit bigger it would discard the file because it
71- // understands it's not a ROOT file), thus one
72- // cannot use TFile::Cp()
73- if (gSystem -> Exec (aliencp .Data ()) != 0 ) {
74- cout << "Error: Sim parameters .json file " << inputFilePathName << " does not exist!" << endl ;
75- exit (1 );
76- }
60+ // old: retrieving from GRID; very flexible but to be avoided because poor reproduceability if file is changed etc for CI tests and analyses
61+ // std::string inputFilePathName = "alien://" + inputSimParametersPath + inputSimParametersFileName;
62+ // if (!gGrid) {
63+ // TGrid::Connect("alien://");
64+ // if (!gGrid) {
65+ // LOG(fatal) << "AliEn connection failed, check token.";
66+ // exit(1);
67+ // }
68+ // }
69+ // // fetch and copy the .json file to the sim work directory
70+ // std::string outputPath = "./";
71+ // TString aliencp = Form("alien_cp alien://%s%s file:%s%s", inputSimParametersPath.c_str(),
72+ // inputSimParametersFileName.c_str(), outputPath.c_str(),
73+ // inputSimParametersFileName.c_str()); // an internal operation in ROOT that discards the
74+ // // JSON file directly (smaller than 300 bytes and if
75+ // // a bit bigger it would discard the file because it
76+ // // understands it's not a ROOT file), thus one
77+ // // cannot use TFile::Cp()
78+ // if (gSystem->Exec(aliencp.Data()) != 0) {
79+ // cout << "Error: Sim parameters .json file " << inputFilePathName << " does not exist!" << endl;
80+ // exit(1);
81+ // }
82+ std ::string O2DPG_ROOT (getenv ("O2DPG_MC_CONFIG_ROOT" ));
83+ std ::string inputSimParametersPath = O2DPG_ROOT + "/MC/config/PWGGAJE/external/generator/parametrisedJetModel/paramConfigs/" ;
84+ std ::string inputFilePathName = inputSimParametersPath + inputSimParametersFileName ;
7785 // open the file
78- std ::FILE * fjson = std ::fopen (inputSimParametersFileName .c_str (), "r" );
86+ std ::FILE * fjson = std ::fopen (inputFilePathName .c_str (), "r" );
7987 if (!fjson ) {
8088 cout << "Could not open sim parameters file " << inputFilePathName << endl ;
8189 exit (1 );
@@ -114,32 +122,50 @@ public:
114122 // get parameters for sim
115123 mSglGenRAA = jsonDocument [mConfigurableSimParameterNames .at (0 ).c_str ()].GetDouble ();
116124 mSglGenTAA = jsonDocument [mConfigurableSimParameterNames .at (1 ).c_str ()].GetDouble ();
117- mSglCutoffSteepNess = jsonDocument [mConfigurableSimParameterNames .at (2 ).c_str ()].GetDouble ();
125+ msglCutoffSteepness = jsonDocument [mConfigurableSimParameterNames .at (2 ).c_str ()].GetDouble ();
118126 mSglCutoffAbscissa = jsonDocument [mConfigurableSimParameterNames .at (3 ).c_str ()].GetDouble ();
119127 mFallSpecterSlopeLog = jsonDocument [mConfigurableSimParameterNames .at (4 ).c_str ()].GetDouble ();
120128 mFallSpecterAffinePowerConstantTerm = jsonDocument [mConfigurableSimParameterNames .at (5 ).c_str ()].GetDouble ();
121129 mFallSpecterAffinePowerSlope = jsonDocument [mConfigurableSimParameterNames .at (6 ).c_str ()].GetDouble ();
122130 mBkgAveragePt = jsonDocument [mConfigurableSimParameterNames .at (7 ).c_str ()].GetDouble ();
123- mCollTotalMultWithBkg = jsonDocument [mConfigurableSimParameterNames .at (8 ).c_str ()].GetDouble ();
131+ mCollMultPowerLawAmplitude =
132+ jsonDocument [mConfigurableSimParameterNames .at (8 ).c_str ()].GetDouble ();
133+ mCollMultPowerLawExponent =
134+ jsonDocument [mConfigurableSimParameterNames .at (9 ).c_str ()].GetDouble ();
135+ mCollMultMin =
136+ jsonDocument [mConfigurableSimParameterNames .at (10 ).c_str ()].GetInt ();
137+ mCollMultMax =
138+ jsonDocument [mConfigurableSimParameterNames .at (11 ).c_str ()].GetInt ();
124139
125140 // clean up
126141 std ::fclose (fjson );
127142
128143 cout << "param retrieved: mSglGenRAA = " << mSglGenRAA << endl ;
129144 cout << "param retrieved: mSglGenTAA = " << mSglGenTAA << endl ;
130- cout << "param retrieved: mSglCutoffSteepNess = " << mSglCutoffSteepNess << endl ;
145+ cout << "param retrieved: msglCutoffSteepness = " << msglCutoffSteepness << endl ;
131146 cout << "param retrieved: mSglCutoffAbscissa = " << mSglCutoffAbscissa << endl ;
132147 cout << "param retrieved: mFallSpecterSlopeLog = " << mFallSpecterSlopeLog << endl ;
133148 cout << "param retrieved: mFallSpecterAffinePowerConstantTerm = " << mFallSpecterAffinePowerConstantTerm << endl ;
134149 cout << "param retrieved: mFallSpecterAffinePowerSlope = " << mFallSpecterAffinePowerSlope << endl ;
135150 cout << "param retrieved: mBkgAveragePt = " << mBkgAveragePt << endl ;
136- cout << "param retrieved: mCollTotalMultWithBkg = " << mCollTotalMultWithBkg << endl ;
137-
138- // thermal background function
151+ cout << "param retrieved: mCollMultPowerLawAmplitude = "
152+ << mCollMultPowerLawAmplitude << endl ;
153+ cout << "param retrieved: mCollMultPowerLawExponent = "
154+ << mCollMultPowerLawExponent << endl ;
155+ cout << "param retrieved: mcollMultMin = " << mCollMultMin << endl ;
156+ cout << "param retrieved: mcollMultMax = " << mCollMultMax << endl ;
157+
158+ // thermal background pdf
139159 mBoltzmannPDF = new TF1 ("f1" , "[0]*[0]*x*exp(-[0]*x)" , mBkgGenPtMin , mPtInfinity );
140160 mBoltzmannPDF -> SetParameter (0 , 2. / mBkgAveragePt );
141161
142- // jet signal function
162+ // collision multiplicity pdf
163+ mCollisionMultPDF =
164+ new TF1 ("f1" , "exp([0])*pow(x,[1])" , mCollMultMin , mCollMultMax );
165+ mCollisionMultPDF -> SetParameter (0 , mCollMultPowerLawAmplitude );
166+ mCollisionMultPDF -> SetParameter (1 , mCollMultPowerLawExponent );
167+
168+ // jet signal pdf
143169 // this thesis says that the jet distrib used to sample parton pt is
144170 // actually full jet -> solves neutral particle fragments issue (better than
145171 // scaling) https://drupal.star.bnl.gov/STAR/files/phd_thesis_rusnak.pdf for
@@ -151,7 +177,7 @@ public:
151177 mJetYieldFit = new TF1 ("f2" , "[0]*[1]*exp(-exp(-[2]*(x-[3]))) * exp([4])*pow(x, [5]+[6]*x)" , 0 , mPtInfinity ); // RAA * TAA * sigmoid(cutoff at [2]) * fit to fulljetSpectrum in pp PYTHIA
152178 mJetYieldFit -> SetParameter (0 , mSglGenRAA ); // rAA (single value for all pt)
153179 mJetYieldFit -> SetParameter (1 , mSglGenTAA ); // <tAA>
154- mJetYieldFit -> SetParameter (2 , mSglCutoffSteepNess ); // steepness for the cutoff shape; higher is
180+ mJetYieldFit -> SetParameter (2 , msglCutoffSteepness ); // steepness for the cutoff shape; higher is
155181 // steeper, but looks more and more like a
156182 // Heavyside step function as it gets to 50 or 100
157183 mJetYieldFit -> SetParameter (3 , mSglCutoffAbscissa ); // cutoffpoint of hard population, in GeV/c
@@ -295,7 +321,24 @@ public:
295321 if (mDebug ) {
296322 cout << "####################### Adding Thermal Background #######################" << endl ;
297323 }
298- for (int iBkg {0 }; iBkg < mCollTotalMultWithBkg - nHardParticles ; ++ iBkg ) {
324+
325+ int mCollTotalMultWithBkg = 0 ;
326+ if (std ::abs (mCollMultMax - mCollMultMin ) == 0 ) {
327+ mCollTotalMultWithBkg = mCollMultMin ;
328+ } else {
329+ mCollTotalMultWithBkg =
330+ mCollisionMultPDF -> GetRandom (mCollMultMin , mCollMultMax );
331+ }
332+
333+ if (mDebug ) {
334+ cout << "mCollTotalMultWithBkg = " << mCollTotalMultWithBkg << endl ;
335+ }
336+
337+ int nBkgParticles = (mCollTotalMultWithBkg - nHardParticles ) > 0
338+ ? mCollTotalMultWithBkg - nHardParticles
339+ : 0 ;
340+
341+ for (int iBkg {0 }; iBkg < nBkgParticles ; ++ iBkg ) {
299342 const double bkgPt = mBoltzmannPDF -> GetRandom (mBkgGenPtMin , mPtInfinity );
300343 const double bkgEta = gRandom -> Uniform (mGenMinEta , mGenMaxEta );
301344 const double bkgPhi = gRandom -> Uniform (0 , o2 ::constants ::math ::TwoPI );
@@ -339,28 +382,43 @@ public:
339382
340383private :
341384 bool mDebug = false; // setting to true will display particle lists
342- std ::string mInputSimParametersPath , mInputSimParametersFileName ; // input path and file name of .json used to read simulation parameters
385+ std ::string mInputSimParametersFileName ; // input path and file name of .json used to read simulation parameters
343386
344387 ////////////////////////////////////////////////
345388 ///////// Common signal and background parameters ////////
346389 ////////////////////////////////////////////////
347390
348391 const double mPtInfinity = 300 ; // maximum pt (in GeV/c) for generated particles, and upper pT limit for integral and TF1 purposes; too high and GetRandom struggles
349392 const double mGenMinEta = -0.9 ; /// minimum pseudorapidity for generated particles
350- const double mGenMaxEta = +0.9 ; /// maximum pseudorapidity for generated particles
351- int mCollTotalMultWithBkg ; /// total multiplicity of the collision
393+ const double mGenMaxEta =
394+ +0.9 ; /// maximum pseudorapidity for generated particles
395+
396+ TF1 * mCollisionMultPDF ; /// TF1 to store pdf function from which collision
397+ /// multiplicity is drawn
398+ double mCollMultPowerLawAmplitude ; /// collision total multiplicity: power law
399+ /// amplitude of the PDF
400+ double mCollMultPowerLawExponent ; /// collision total multiplicity: power law
401+ /// exponent of the PDF
402+ int mCollMultMin ; /// collision total multiplicity: minimum abscissa of the
403+ /// PDF
404+ int mCollMultMax ; /// collision total multiplicity: maximum abscissa of the
405+ /// PDF
406+
352407 bool mGenerateSignal = true; /// boolean to request (or not) the generation of the jet signal
353408 bool mGenerateUE = false; /// boolean to request (or not) embedding of the jet signal inside underlying event modelled by a thermal background; if mGenerateSignal = false, only the UE is generated
354409 const std ::vector < std ::string > mConfigurableSimParameterNames = {
355410 "sglGenRAA" ,
356411 "sglGenTAA" ,
357- "sglCutoffSteepNess " ,
412+ "sglCutoffSteepness " ,
358413 "sglCutoffAbscissa" ,
359414 "fallSpecterSlopeLog" ,
360415 "fallSpecterAffinePowerConstantTerm" ,
361416 "fallSpecterAffinePowerSlope" ,
362417 "bkgAveragePt" ,
363- "collTotalMultWithBkg" };
418+ "collMultPowerLawAmplitude" ,
419+ "collMultPowerLawExponent" ,
420+ "collMultMin" ,
421+ "collMultMax" };
364422
365423 /////////////////////////////////////////////
366424 /////// Thermal background parameters ///////
@@ -383,7 +441,7 @@ private:
383441 // signal fit
384442 double mSglGenRAA ;
385443 double mSglGenTAA ;
386- double mSglCutoffSteepNess ; /// steepness for the cutoff shape; higher is steeper, but looks more and more like a Heavyside step function as it gets to 50 or 100
444+ double msglCutoffSteepness ; /// steepness for the cutoff shape; higher is steeper, but looks more and more like a Heavyside step function as it gets to 50 or 100
387445 double mSglCutoffAbscissa ; /// minimum pt (in GeV/c) for jet signal distribution; it's a smooth cutoff
388446 double mFallSpecterSlopeLog ;
389447 double mFallSpecterAffinePowerConstantTerm ;
@@ -398,8 +456,7 @@ private:
398456};
399457
400458///___________________________________________________________
401- FairGenerator * generateParametrisedJetModel (std ::string inputSimParametersPath ,
402- std ::string inputSimParametersFileName ,
459+ FairGenerator * generateParametrisedJetModel (std ::string inputSimParametersFileName ,
403460 bool generateUE , bool generateSignal = true) {
404- return new GeneratorParametrisedJetModel (inputSimParametersPath , inputSimParametersFileName , generateUE , generateSignal );
461+ return new GeneratorParametrisedJetModel (inputSimParametersFileName , generateUE , generateSignal );
405462}
0 commit comments