Fix seeding and include X particle (#42)

* WIP: add files for X simulations

* Fix rpythia8

* Add B hadron productions

* Fix path of the KrKr database

* remove cfg which are migrated to O2DPG

* restore ccbar example config

* Revert "restore ccbar example config"

This reverts commit faabf5a.

* Revert "remove cfg which are migrated to O2DPG"

This reverts commit 73b92ab.

* revert to default LUTs

Author: ginnocen

examples/pythia8/decays/force_hadronic_B.cfg

@@ -0,0 +1,18 @@
+### Force golden B-hadron hadronic decay modes
+### Based on AliRoot AliDecayerPythia8::ForceBeautyUpgrade, latest commit 3154cf8 on Nov 30, 2020
+
+# Bs -> Ds- pi+
+531:onMode = off
+531:onIfMatch = 431 211
+# Lb: 50% to Lc any, 50% to Lc pion
+#FIXME: we should find a way to generate 50% of inclusive or Lc decays
+#FIXME: currently all the Lb -> Lc + pi
+5122:onMode = off
+#5122:onIfAll = 4122
+5122:onIfMatch = 4122 211
+
+511:onMode = off
+511:onIfMatch = 413 211
+
+521:onMode = off
+521:onIfMatch = 421 211

examples/pythia8/pythia8_bbbar.cfg

@@ -0,0 +1,17 @@
+### main
+Main:numberOfEvents 10000
+
+### service
+Next:numberShowEvent = 0
+
+### random
+Random:setSeed = on
+Random:seed = 123456789
+
+### beams
+Beams:idA 2212				# proton
+Beams:idB 2212 				# proton
+Beams:eCM 14000. 			# GeV
+
+### processes
+HardQCD:hardbbbar on			# g-g / q-qbar -> b-bbar 

examples/pythia8/pythia_onia_X.cfg

@@ -0,0 +1,64 @@
+### beams
+Beams:idA 2212                  # proton
+Beams:idB 2212                  # proton
+Beams:eCM 14000.                # GeV
+#663:addChannel = 1 0. 0 1 -1
+#9920443:addParticle = X(3872) 3 0 0 3.87196 0.00012
+Onia:all = on
+Random:setSeed = on
+
+Charmonium:states(3PJ)    = 10441,20443,445,9920443
+Charmonium:O(3PJ)[3P0(1)] = 0.05,0.05,0.05,0.05
+Charmonium:O(3PJ)[3S1(8)] = 0.0031,0.0031,0.0031,0.0031
+Charmonium:gg2ccbar(3PJ)[3PJ(1)]g    = off,off,off,on
+Charmonium:qg2ccbar(3PJ)[3PJ(1)]q    = off,off,off,on
+Charmonium:qqbar2ccbar(3PJ)[3PJ(1)]g = off,off,off,on
+Charmonium:gg2ccbar(3PJ)[3S1(8)]g    = off,off,off,on
+Charmonium:qg2ccbar(3PJ)[3S1(8)]q    = off,off,off,on
+Charmonium:qqbar2ccbar(3PJ)[3S1(8)]g = off,off,off,on
+
+### B0  -> J/psi X
+511:onMode = off
+511:onIfAny = 443
+
+###  B+/-  -> J/psi X
+521:onMode = off
+521:onIfAny = 443
+
+###B_s  -> J/psi X
+531:onMode = off
+531:onIfAny = 443
+
+##Lambda_b -> J/psi X
+5122:onMode = off
+5122:onIfAny = 443
+
+
+###J/psi -> mu+ mu-
+443:onMode = off
+443:onIfAll = 13 -13
+443:onIfAll = 11 -11 //e+e-
+
+###Psi(2S) -> mu+ mu-
+100443:onMode = off
+100443:onIfAll = 13 -13
+100443:onIfAll = 11 -11 //e+e-
+
+###Upsilon -> mu+ mu-
+553:onMode = off
+553:onIfAll = 13 -13
+553:onIfAll = 11 -11 //e+e-
+
+###Upsilon(2S) -> mu+ mu-
+100553:onMode = off
+100553:onIfAll = 13 -13
+100553:onIfAll = 11 -11 //e+e-
+
+###Upsilon(3S) -> mu+ mu-
+200553:onMode = off
+200553:onIfAll = 13 -13
+200553:onIfAll = 11 -11 //e+e-
+
+9920443:addChannel = 1 0 443 211 -211
+9920443:onMode = off
+9920443:onIfMatch = 443 211 -211

examples/scripts/createO2tables.py

@@ -266,7 +266,8 @@ def copy_and_link(file_name):
                 gen_log_file = f"gen.{run_number}.log"
                 hepmc_file = f"hepmcfile.{run_number}.hepmc"
                 custom_gen_option = f" --output {hepmc_file} --nevents {nevents}"
-                write_to_runner(custom_gen + custom_gen_option,
+                custom_seed = f" --seed {run_number}"
+                write_to_runner(custom_gen + custom_gen_option + custom_seed,
                                 log_file=gen_log_file)
                 write_to_runner(f"DelphesHepMC propagate.tcl {delphes_file} {hepmc_file}",
                                 log_file=delphes_log_file)

examples/scripts/default_configfile.ini

@@ -11,7 +11,7 @@ bfield = 5.
 sigmat = 0.020
 
 # radius [cm]
-radius = 100.
+radius = 20.
 
 # half length [cm]
 length = 200.
@@ -27,17 +27,22 @@ propagate_card = propagate.2kG.tcl
 
 # path to the LUTs
 lut_path = $DELPHESO2_ROOT/examples/smearing/luts/
+#lut_path = /data/Run5data/LUTs
 
 # tag of the LUTs to use in simulation
-lut_tag = 100cm.default
+#lut_tag = 100cm.default
+lut_tag = 20cm.scenario3
 
 # path of the DELPHES aod utilities
 aod_path = $DELPHESO2_ROOT/examples/aod/
 
 ## Pythia setup
-[CCBAR]
+[CCBAR_Dforced]
 generators = $DELPHESO2_ROOT/examples/pythia8/pythia8_ccbar.cfg $DELPHESO2_ROOT/examples/pythia8/decays/force_hadronic_D.cfg
 
+[BBBAR_BDforced]
+generators = $DELPHESO2_ROOT/examples/pythia8/pythia8_bbbar.cfg $DELPHESO2_ROOT/examples/pythia8/decays/force_hadronic_B.cfg $DELPHESO2_ROOT/examples/pythia8/decays/force_hadronic_D.cfg
+
 [CCBAR_Lc_pK0s] # (need O2DPG loaded)
 generators = $DELPHESO2_ROOT/examples/pythia8/pythia8_ccbar.cfg $O2DPG_ROOT/MC/config/PWGHF/pythia8/decayer/force_hadronic_D.cfg $O2DPG_ROOT/MC/config/PWGHF/pythia8/decayer/force_hadronic_D_forceLcChannel2.cfg
 
@@ -54,7 +59,7 @@ custom_gen = rpythia8 --config $O2DPG_ROOT/MC/config/common/pythia8/generator/py
 generators = $O2DPG_ROOT/MC/config/common/pythia8/generator/pythia8_hi.cfg
 
 [INEL_HI_KRKR] # (need O2DPG loaded)
-generators = $O2DPG_ROOT/MC/config/common/pythia8/generator/pythia8_hi_krkr.cfg
+generators = $O2DPG_ROOT/MC/config/common/pythia8/generator/pythia8_KrKr.cfg
 
 [GUN]
 custom_gen = rpythia8-gun --pdg 421 --px 1. --py 0. --pz 0. --xProd 0. --yProd 0. --zProd 0. --config $DELPHESO2_ROOT/examples/pythia8/decays/force_hadronic_D.cfg --decay
@@ -65,3 +70,5 @@ custom_gen = rpythia8-gun --pdg 4122 --px 1. --py 0. --pz 0. --xProd 1. --yProd
 [GUN_Lc_pK0s]
 custom_gen = rpythia8-gun --pdg 4122 --px 1. --py 0. --pz 0. --xProd 1. --yProd 0. --zProd 0. --config $O2DPG_ROOT/MC/config/PWGHF/pythia8/decayer/force_hadronic_D_forceLcChannel2.cfg --decay
 
+[ONIA_X] #(need O2DPG loaded)
+custom_gen = rpythia8 --config $DELPHESO2_ROOT/examples/pythia8/pythia_onia_X.cfg

rpythia8/rpythia8.cc

@@ -2,12 +2,13 @@
 #include <string>
 
 #include "Pythia8/Pythia.h"
+#include "Pythia8/ParticleData.h"
 #include "Pythia8Plugins/HepMC2.h"
 
 int main(int argc, char **argv)
 {
 
-  int nevents, inject_nevents;
+  int nevents, inject_nevents, seed;
   std::string config, output, inject_config;
 
   /** process arguments **/
@@ -21,6 +22,7 @@ int main(int argc, char **argv)
       ("output,o"       , po::value<std::string>(&output)->default_value("pythia8.hepmc"), "Output HepMC file")
       ("inject-config"  , po::value<std::string>(&inject_config), "Injected event configuration file")
       ("inject-nevents" , po::value<int>(&inject_nevents)->default_value(1), "Number of events to inject")
+      ("seed"           , po::value<int>(&seed)->default_value(1), "initial seed")
       ;
 
     po::variables_map vm;
@@ -40,11 +42,16 @@ int main(int argc, char **argv)
 
   // pythia, config and init
   Pythia8::Pythia pythia;
+  pythia.particleData.addParticle(9920443, "X(3872)", 3, 0, 0, 3.87196, 0.00012);
+
   Pythia8::Rndm   rndm;
   if (!config.empty() && !pythia.readFile(config)) {
     std::cout << "Error: could not read config file \"" << config << "\"" << std::endl;
     return 1;
   }
+  std::cout<<"Random:seed =" + std::to_string(seed)<<std::endl;
+  pythia.readString("Random:setSeed = on");
+  pythia.readString("Random:seed =" + std::to_string(seed));
   pythia.init();
   rndm.init();
 
@@ -61,6 +68,8 @@ int main(int argc, char **argv)
       std::cout << "Error: could not read config file \"" << inject_config << "\"" << std::endl;
       return 1;
     }
+    pythia_inj->readString("Random:setSeed = on");
+    pythia_inj->readString("Random:seed =" + std::to_string(seed));
     pythia_inj->init();
   }