[PWGGAJE] Param Model: better multiplicity model

MC/config/PWGGAJE/external/generator/parametrisedJetModel/generator_pythia8_box_parametrisedModel_pythia6Fragmentation.C

@@ -26,18 +26,22 @@ using namespace Pythia8;
 // #include "SimulationDataFormat/MCEventHeader.h"
 
 // Input to simulation:
-// inputFilePathName file is expected to be a json file with the structure like so:
+// inputFilePathName file is expected to be a json file with the structure like
+// so:
 //  {
 //    "simLog": false,
 //    "sglGenRAA": 0.45,
 //    "sglGenTAA": 23.38,
-//    "sglCutoffSteepNess": 15,
+//    "sglCutoffSteepness": 15,
 //    "sglCutoffAbscissa": 10,
 //    "fallSpecterSlopeLog": 10.833274,
 //    "fallSpecterAffinePowerConstantTerm": -5.476804,
 //    "fallSpecterAffinePowerSlope": 0.001110,
 //    "bkgAveragePt": 0.670,
-//    "collTotalMultWithBkg": 2000
+//    "collMultPowerLawAmplitude": ?,
+//    "collMultPowerLawExponent": ?,
+//    "collMultMin": ?,
+//    "collMultMax": ?
 //  }
 // can be uploaded to grid using for example: alien.py cp
 // file:/local/path/parametrisedModel_PbPb_5p36TeV_cent0010.json
@@ -46,36 +50,40 @@ using namespace Pythia8;
 class GeneratorParametrisedJetModel : public o2::eventgen::GeneratorPythia8 {
 public:
   /// constructor
-  GeneratorParametrisedJetModel(std::string inputSimParametersPath, std::string inputSimParametersFileName, bool generateUE, bool generateSignal = true) : mInputSimParametersPath{inputSimParametersPath}, mInputSimParametersFileName{inputSimParametersFileName}, mGenerateUE{generateUE}, mGenerateSignal{generateSignal} {
+  GeneratorParametrisedJetModel(std::string inputSimParametersFileName, bool generateUE, bool generateSignal = true) : mInputSimParametersFileName{inputSimParametersFileName}, mGenerateUE{generateUE}, mGenerateSignal{generateSignal} {
 
     if (!generateUE && !generateSignal) {
       cout << "Both the UE and the signal generation should not be off at the same time" << endl;
       exit(1);
     }
 
-    std::string inputFilePathName = "alien://" + inputSimParametersPath + inputSimParametersFileName;
-    if (!gGrid) {
-      TGrid::Connect("alien://");
-      if (!gGrid) {
-        LOG(fatal) << "AliEn connection failed, check token.";
-        exit(1);
-      }
-    }
-    // fetch and copy the .json file to the sim work directory
-    std::string outputPath = "./";
-    TString aliencp = Form("alien_cp alien://%s%s file:%s%s", inputSimParametersPath.c_str(),
-                            inputSimParametersFileName.c_str(), outputPath.c_str(), 
-                            inputSimParametersFileName.c_str()); // an internal operation in ROOT that discards the
-                                                                 // JSON file directly (smaller than 300 bytes and if
-                                                                 // a bit bigger it would discard the file because it
-                                                                 // understands it's not a ROOT file), thus one
-                                                                 // cannot use TFile::Cp()
-    if (gSystem->Exec(aliencp.Data()) != 0) {
-      cout << "Error: Sim parameters .json file " << inputFilePathName << " does not exist!" << endl;
-      exit(1);
-    }
+    // old: retrieving from GRID; very flexible but to be avoided because poor reproduceability if file is changed etc for CI tests and analyses
+    // std::string inputFilePathName = "alien://" + inputSimParametersPath + inputSimParametersFileName;
+    // if (!gGrid) {
+    //   TGrid::Connect("alien://");
+    //   if (!gGrid) {
+    //     LOG(fatal) << "AliEn connection failed, check token.";
+    //     exit(1);
+    //   }
+    // }
+    // // fetch and copy the .json file to the sim work directory
+    // std::string outputPath = "./";
+    // TString aliencp = Form("alien_cp alien://%s%s file:%s%s", inputSimParametersPath.c_str(),
+    //                         inputSimParametersFileName.c_str(), outputPath.c_str(), 
+    //                         inputSimParametersFileName.c_str()); // an internal operation in ROOT that discards the
+    //                                                              // JSON file directly (smaller than 300 bytes and if
+    //                                                              // a bit bigger it would discard the file because it
+    //                                                              // understands it's not a ROOT file), thus one
+    //                                                              // cannot use TFile::Cp()
+    // if (gSystem->Exec(aliencp.Data()) != 0) {
+    //   cout << "Error: Sim parameters .json file " << inputFilePathName << " does not exist!" << endl;
+    //   exit(1);
+    // }
+    std::string O2DPG_ROOT(getenv("O2DPG_MC_CONFIG_ROOT"));
+    std::string inputSimParametersPath = O2DPG_ROOT+"/MC/config/PWGGAJE/external/generator/parametrisedJetModel/paramConfigs/";
+    std::string inputFilePathName = inputSimParametersPath + inputSimParametersFileName;
     // open the file
-    std::FILE *fjson = std::fopen(inputSimParametersFileName.c_str(), "r");
+    std::FILE *fjson = std::fopen(inputFilePathName.c_str(), "r");
     if (!fjson) {
       cout << "Could not open sim parameters file " << inputFilePathName << endl;
       exit(1);
@@ -114,32 +122,50 @@ public:
     // get parameters for sim
     mSglGenRAA = jsonDocument[mConfigurableSimParameterNames.at(0).c_str()].GetDouble();
     mSglGenTAA = jsonDocument[mConfigurableSimParameterNames.at(1).c_str()].GetDouble();
-    mSglCutoffSteepNess = jsonDocument[mConfigurableSimParameterNames.at(2).c_str()].GetDouble();
+    msglCutoffSteepness = jsonDocument[mConfigurableSimParameterNames.at(2).c_str()].GetDouble();
     mSglCutoffAbscissa = jsonDocument[mConfigurableSimParameterNames.at(3).c_str()].GetDouble();
     mFallSpecterSlopeLog = jsonDocument[mConfigurableSimParameterNames.at(4).c_str()].GetDouble();
     mFallSpecterAffinePowerConstantTerm = jsonDocument[mConfigurableSimParameterNames.at(5).c_str()].GetDouble();
     mFallSpecterAffinePowerSlope = jsonDocument[mConfigurableSimParameterNames.at(6).c_str()].GetDouble();
     mBkgAveragePt = jsonDocument[mConfigurableSimParameterNames.at(7).c_str()].GetDouble();
-    mCollTotalMultWithBkg = jsonDocument[mConfigurableSimParameterNames.at(8).c_str()].GetDouble();
+    mCollMultPowerLawAmplitude =
+        jsonDocument[mConfigurableSimParameterNames.at(8).c_str()].GetDouble();
+    mCollMultPowerLawExponent =
+        jsonDocument[mConfigurableSimParameterNames.at(9).c_str()].GetDouble();
+    mCollMultMin =
+        jsonDocument[mConfigurableSimParameterNames.at(10).c_str()].GetInt();
+    mCollMultMax =
+        jsonDocument[mConfigurableSimParameterNames.at(11).c_str()].GetInt();
 
     // clean up
     std::fclose(fjson);
 
     cout << "param retrieved: mSglGenRAA = " << mSglGenRAA << endl;
     cout << "param retrieved: mSglGenTAA = " << mSglGenTAA << endl;
-    cout << "param retrieved: mSglCutoffSteepNess = " << mSglCutoffSteepNess << endl;
+    cout << "param retrieved: msglCutoffSteepness = " << msglCutoffSteepness << endl;
     cout << "param retrieved: mSglCutoffAbscissa = " << mSglCutoffAbscissa << endl;
     cout << "param retrieved: mFallSpecterSlopeLog = " << mFallSpecterSlopeLog << endl;
     cout << "param retrieved: mFallSpecterAffinePowerConstantTerm = " << mFallSpecterAffinePowerConstantTerm << endl;
     cout << "param retrieved: mFallSpecterAffinePowerSlope = " << mFallSpecterAffinePowerSlope << endl;
     cout << "param retrieved: mBkgAveragePt = " << mBkgAveragePt << endl;
-    cout << "param retrieved: mCollTotalMultWithBkg = " << mCollTotalMultWithBkg << endl;
-
-    // thermal background function
+    cout << "param retrieved: mCollMultPowerLawAmplitude = "
+         << mCollMultPowerLawAmplitude << endl;
+    cout << "param retrieved: mCollMultPowerLawExponent = "
+         << mCollMultPowerLawExponent << endl;
+    cout << "param retrieved: mcollMultMin = " << mCollMultMin << endl;
+    cout << "param retrieved: mcollMultMax = " << mCollMultMax << endl;
+
+    // thermal background pdf
     mBoltzmannPDF = new TF1("f1", "[0]*[0]*x*exp(-[0]*x)", mBkgGenPtMin, mPtInfinity);
     mBoltzmannPDF->SetParameter(0, 2. / mBkgAveragePt);
 
-    // jet signal function
+    // collision multiplicity pdf
+    mCollisionMultPDF =
+        new TF1("f1", "exp([0])*pow(x,[1])", mCollMultMin, mCollMultMax);
+    mCollisionMultPDF->SetParameter(0, mCollMultPowerLawAmplitude);
+    mCollisionMultPDF->SetParameter(1, mCollMultPowerLawExponent);
+
+    // jet signal pdf
     // this thesis says that the jet distrib used to sample parton pt is
     // actually full jet -> solves neutral particle fragments issue (better than
     // scaling) https://drupal.star.bnl.gov/STAR/files/phd_thesis_rusnak.pdf for
@@ -151,7 +177,7 @@ public:
     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
     mJetYieldFit->SetParameter(0, mSglGenRAA); // rAA (single value for all pt)
     mJetYieldFit->SetParameter(1, mSglGenTAA); // <tAA>
-    mJetYieldFit->SetParameter(2, mSglCutoffSteepNess); // steepness for the cutoff shape; higher is
+    mJetYieldFit->SetParameter(2, 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
     mJetYieldFit->SetParameter(3, mSglCutoffAbscissa); // cutoffpoint of hard population, in GeV/c
@@ -295,7 +321,24 @@ public:
       if (mDebug) {
         cout << "####################### Adding Thermal Background #######################" << endl;
       }
-      for (int iBkg{0}; iBkg < mCollTotalMultWithBkg - nHardParticles; ++iBkg) {
+
+      int mCollTotalMultWithBkg = 0;
+      if (std::abs(mCollMultMax - mCollMultMin) == 0) {
+        mCollTotalMultWithBkg = mCollMultMin;
+      } else {
+        mCollTotalMultWithBkg =
+            mCollisionMultPDF->GetRandom(mCollMultMin, mCollMultMax);
+      }
+
+      if (mDebug) {
+        cout << "mCollTotalMultWithBkg = " << mCollTotalMultWithBkg << endl;
+      }
+
+      int nBkgParticles = (mCollTotalMultWithBkg - nHardParticles) > 0
+                              ? mCollTotalMultWithBkg - nHardParticles
+                              : 0;
+
+      for (int iBkg{0}; iBkg < nBkgParticles; ++iBkg) {
         const double bkgPt = mBoltzmannPDF->GetRandom(mBkgGenPtMin, mPtInfinity);
         const double bkgEta = gRandom->Uniform(mGenMinEta, mGenMaxEta);
         const double bkgPhi = gRandom->Uniform(0, o2::constants::math::TwoPI);
@@ -339,28 +382,43 @@ public:
 
 private:
   bool mDebug = false; // setting to true will display particle lists
-  std::string mInputSimParametersPath, mInputSimParametersFileName; // input path and file name of .json used to read simulation parameters
+  std::string mInputSimParametersFileName; // input path and file name of .json used to read simulation parameters
 
   ////////////////////////////////////////////////
   ///////// Common signal and background parameters ////////
   ////////////////////////////////////////////////
 
   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
   const double mGenMinEta = -0.9; /// minimum pseudorapidity for generated particles
-  const double mGenMaxEta = +0.9; /// maximum pseudorapidity for generated particles
-  int mCollTotalMultWithBkg; /// total multiplicity of the collision
+  const double mGenMaxEta =
+      +0.9; /// maximum pseudorapidity for generated particles
+
+  TF1 *mCollisionMultPDF; /// TF1 to store pdf function from which collision
+                          /// multiplicity is drawn
+  double mCollMultPowerLawAmplitude; /// collision total multiplicity: power law
+                                     /// amplitude of the PDF
+  double mCollMultPowerLawExponent;  /// collision total multiplicity: power law
+                                     /// exponent of the PDF
+  int mCollMultMin; /// collision total multiplicity: minimum abscissa of the
+                    /// PDF
+  int mCollMultMax; /// collision total multiplicity: maximum abscissa of the
+                    /// PDF
+
   bool mGenerateSignal = true; /// boolean to request (or not) the generation of the jet signal
   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
   const std::vector<std::string> mConfigurableSimParameterNames = {
       "sglGenRAA",
       "sglGenTAA",
-      "sglCutoffSteepNess",
+      "sglCutoffSteepness",
       "sglCutoffAbscissa",
       "fallSpecterSlopeLog",
       "fallSpecterAffinePowerConstantTerm",
       "fallSpecterAffinePowerSlope",
       "bkgAveragePt",
-      "collTotalMultWithBkg"};
+      "collMultPowerLawAmplitude",
+      "collMultPowerLawExponent",
+      "collMultMin",
+      "collMultMax"};
 
   /////////////////////////////////////////////
   /////// Thermal background parameters ///////
@@ -383,7 +441,7 @@ private:
   // signal fit
   double mSglGenRAA;
   double mSglGenTAA;
-  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
+  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
   double mSglCutoffAbscissa;  /// minimum pt (in GeV/c) for jet signal distribution; it's a smooth cutoff
   double mFallSpecterSlopeLog;
   double mFallSpecterAffinePowerConstantTerm;
@@ -398,8 +456,7 @@ private:
 };
 
 ///___________________________________________________________
-FairGenerator *generateParametrisedJetModel(std::string inputSimParametersPath, 
-                                            std::string inputSimParametersFileName, 
+FairGenerator *generateParametrisedJetModel(std::string inputSimParametersFileName, 
                                             bool generateUE, bool generateSignal = true) {
-  return new GeneratorParametrisedJetModel(inputSimParametersPath, inputSimParametersFileName, generateUE, generateSignal);
+  return new GeneratorParametrisedJetModel(inputSimParametersFileName, generateUE, generateSignal);
 }

MC/config/PWGGAJE/external/generator/parametrisedJetModel/paramConfigs/parametrisedModel_PbPb_5p36TeV_cent0010.json

@@ -0,0 +1,15 @@
+{
+  "simLog": false,
+  "sglGenRAA": 0.45,
+  "sglGenTAA": 23.38,
+  "sglCutoffSteepness": 15,
+  "sglCutoffAbscissa": 4,
+  "fallSpecterSlopeLog": 6.619076,
+  "fallSpecterAffinePowerConstantTerm": -3.918165,
+  "fallSpecterAffinePowerSlope": -0.002032,
+  "bkgAveragePt": 0.764567,
+  "collMultPowerLawAmplitude": 20.66,
+  "collMultPowerLawExponent": -1.26971,
+  "collMultMin": 1681,
+  "collMultMax": 2400
+}

MC/config/PWGGAJE/external/generator/parametrisedJetModel/paramConfigs/parametrisedModel_PbPb_5p36TeV_cent5080.json

@@ -0,0 +1,15 @@
+{
+  "simLog": false,
+  "sglGenRAA": 0.8,
+  "sglGenTAA": 0.6933,
+  "sglCutoffSteepness": 15,
+  "sglCutoffAbscissa": 4,
+  "fallSpecterSlopeLog": 6.619076,
+  "fallSpecterAffinePowerConstantTerm": -3.918165,
+  "fallSpecterAffinePowerSlope": -0.002032,
+  "bkgAveragePt": 0.732677,
+  "collMultPowerLawAmplitude": 16.43,
+  "collMultPowerLawExponent": -0.683784,
+  "collMultMin": 47,
+  "collMultMax": 302
+}

MC/config/PWGGAJE/ini/GeneratorPythia8Box_parametrisedModel_pythia6Fragmentation_PbPb_cent0010_Bkg.ini

@@ -1,7 +1,7 @@
 ### bkg production using pythia8 box generator, to embed jet-jet production into
 [GeneratorExternal]
 fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/generator/parametrisedJetModel/generator_pythia8_box_parametrisedModel_pythia6Fragmentation.C
-funcName = generateParametrisedJetModel("/alice/cern.ch/user/a/alandou/Analysis/PWGJE/simSettings/combinatorialBkgClosure/", "parametrisedModel_PbPb_5p36TeV_cent0010.json", true, false)
+funcName = generateParametrisedJetModel("parametrisedModel_PbPb_5p36TeV_cent0010.json", true, false)
 
 [GeneratorPythia8]
 config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8box_parametrisedModel.cfg

MC/config/PWGGAJE/ini/GeneratorPythia8Box_parametrisedModel_pythia6Fragmentation_PbPb_cent0010_Sgl.ini

@@ -1,7 +1,7 @@
 ### bkg production using pythia8 box generator, to embed jet-jet production into
 [GeneratorExternal]
 fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/generator/parametrisedJetModel/generator_pythia8_box_parametrisedModel_pythia6Fragmentation.C
-funcName = generateParametrisedJetModel("/alice/cern.ch/user/a/alandou/Analysis/PWGJE/simSettings/combinatorialBkgClosure/", "parametrisedModel_PbPb_5p36TeV_cent0010.json", false, true)
+funcName = generateParametrisedJetModel("parametrisedModel_PbPb_5p36TeV_cent0010.json", false, true)
 
 [GeneratorPythia8]
 config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8box_parametrisedModel.cfg

MC/config/PWGGAJE/ini/GeneratorPythia8Box_parametrisedModel_pythia6Fragmentation_PbPb_cent0010_SglBkg.ini

@@ -1,7 +1,7 @@
 ### bkg production using pythia8 box generator, to embed jet-jet production into
 [GeneratorExternal]
 fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/generator/parametrisedJetModel/generator_pythia8_box_parametrisedModel_pythia6Fragmentation.C
-funcName = generateParametrisedJetModel("/alice/cern.ch/user/a/alandou/Analysis/PWGJE/simSettings/combinatorialBkgClosure/", "parametrisedModel_PbPb_5p36TeV_cent0010.json", true, true)
+funcName = generateParametrisedJetModel("parametrisedModel_PbPb_5p36TeV_cent0010.json", true, true)
 
 [GeneratorPythia8]
 config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8box_parametrisedModel.cfg

MC/config/PWGGAJE/ini/GeneratorPythia8Box_parametrisedModel_pythia6Fragmentation_PbPb_cent5080_Bkg.ini

@@ -1,7 +1,7 @@
 ### bkg production using pythia8 box generator, to embed jet-jet production into
 [GeneratorExternal]
 fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/generator/parametrisedJetModel/generator_pythia8_box_parametrisedModel_pythia6Fragmentation.C
-funcName = generateParametrisedJetModel("/alice/cern.ch/user/a/alandou/Analysis/PWGJE/simSettings/combinatorialBkgClosure/", "parametrisedModel_PbPb_5p36TeV_cent5080.json", true, false)
+funcName = generateParametrisedJetModel("parametrisedModel_PbPb_5p36TeV_cent5080.json", true, false)
 
 [GeneratorPythia8]
 config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8box_parametrisedModel.cfg

MC/config/PWGGAJE/ini/GeneratorPythia8Box_parametrisedModel_pythia6Fragmentation_PbPb_cent5080_Sgl.ini

@@ -1,7 +1,7 @@
 ### bkg production using pythia8 box generator, to embed jet-jet production into
 [GeneratorExternal]
 fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/generator/parametrisedJetModel/generator_pythia8_box_parametrisedModel_pythia6Fragmentation.C
-funcName = generateParametrisedJetModel("/alice/cern.ch/user/a/alandou/Analysis/PWGJE/simSettings/combinatorialBkgClosure/", "parametrisedModel_PbPb_5p36TeV_cent5080.json", false, true)
+funcName = generateParametrisedJetModel("parametrisedModel_PbPb_5p36TeV_cent5080.json", false, true)
 
 [GeneratorPythia8]
 config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8box_parametrisedModel.cfg

MC/config/PWGGAJE/ini/GeneratorPythia8Box_parametrisedModel_pythia6Fragmentation_PbPb_cent5080_SglBkg.ini

@@ -1,7 +1,7 @@
 ### bkg production using pythia8 box generator, to embed jet-jet production into
 [GeneratorExternal]
 fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/generator/parametrisedJetModel/generator_pythia8_box_parametrisedModel_pythia6Fragmentation.C
-funcName = generateParametrisedJetModel("/alice/cern.ch/user/a/alandou/Analysis/PWGJE/simSettings/combinatorialBkgClosure/", "parametrisedModel_PbPb_5p36TeV_cent5080.json", true, true)
+funcName = generateParametrisedJetModel("parametrisedModel_PbPb_5p36TeV_cent5080.json", true, true)
 
 [GeneratorPythia8]
 config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8box_parametrisedModel.cfg