Merge pull request #39 from gconesab/master

Fix Pb PDG and automatic energy setting for Pb and HI pythia configuration

Author: gconesab

MC/bin/o2dpg_sim_workflow.py

@@ -5,10 +5,16 @@
 # It aims to handle the different MC possible configurations 
 # It just creates a workflow.json txt file, to execute the workflow one must execute right after
 #   ${O2DPG_ROOT}/MC/bin/o2_dpg_workflow_runner.py -f workflow.json 
+#
 # Execution examples:
-#   ./o2dpg_sim_workflow.py -e TGeant3 -nb 0 -ns 2 -j 8 -tf 1 -mod "-m TPC" -proc "jets" -ptTrigMin 3.5 -ptHatBin 3 -trigger "external" -ini "\$O2DPG_ROOT/MC/config/PWGGAJE/ini/trigger_decay_gamma.ini" --embedding False 
+#  - pp PYTHIA jets, 2 events, triggered on high pT decay photons on EMCal acceptance, eCMS 13 TeV
+#     ./o2dpg_sim_workflow.py -e TGeant3 -ns 2 -j 8 -tf 1 -mod "--skipModules ZDC" -col pp -eCM 13000 \
+#                             -proc "jets" -ptTrigMin 3.5 -acceptance 4 -ptHatBin 3 \
+#                             -trigger "external" -ini "\$O2DPG_ROOT/MC/config/PWGGAJE/ini/trigger_decay_gamma.ini"
 #
-#  ./o2dpg_sim_workflow.py -e TGeant3 -nb 0 -ns 2 -j 8 -tf 1 -mod "--skipModules ZDC" -proc "ccbar"  --embedding True 
+#  - pp PYTHIA ccbar events embedded into heavy-ion environment, 2 PYTHIA events into 1 bkg event, beams energy 2.510
+#     ./o2dpg_sim_workflow.py -e TGeant3 -nb 1 -ns 2 -j 8 -tf 1 -mod "--skipModules ZDC"  \
+#                             -col pp -eA 2.510 -proc "ccbar"  --embedding
 # 
 
 import argparse
@@ -25,7 +31,9 @@
 parser.add_argument('-ini',help='generator init parameters file, for example: ${O2DPG_ROOT}/MC/config/PWGHF/ini/GeneratorHF.ini', default='')
 parser.add_argument('-confKey',help='generator or trigger configuration key values, for example: GeneratorPythia8.config=pythia8.cfg', default='')
 
-parser.add_argument('-eCMS',help='CMS energy', default=5200.0)
+parser.add_argument('-eCM',help='CMS energy', default=-1)
+parser.add_argument('-eA',help='Beam A energy', default=6499.) #6369 PbPb, 2.510 pp 5 TeV, 4 pPb
+parser.add_argument('-eB',help='Beam B energy', default=-1)
 parser.add_argument('-col',help='collision sytem: pp, PbPb, pPb, Pbp, ...', default='pp')
 parser.add_argument('-ptHatBin',help='pT hard bin number', default=-1)
 parser.add_argument('-ptHatMin',help='pT hard minimum when no bin requested', default=0)
@@ -34,6 +42,7 @@
 
 parser.add_argument('-ptTrigMin',help='generated pT trigger minimum', default=0)
 parser.add_argument('-ptTrigMax',help='generated pT trigger maximum', default=-1)
+parser.add_argument('-acceptance',help='select particles within predefined acceptance in ${O2DPG_ROOT}/MC/run/common/detector_acceptance.C', default=0)
 
 parser.add_argument('--embedding',action='store_true', help='With embedding into background')
 parser.add_argument('-nb',help='number of background events / timeframe', default=20)
@@ -159,7 +168,9 @@ def getDPL_global_options(bigshm=False,nosmallrate=False):
    # function encapsulating the signal sim part
    # first argument is timeframe id
    RNDSEED=args.seed    # 0 means random seed !
-   ECMS=args.eCMS
+   ECMS=float(args.eCM)
+   EBEAMA=float(args.eA)
+   EBEAMB=float(args.eB)
    NSIGEVENTS=args.ns
    GENERATOR=args.gen
    INIFILE=''
@@ -175,6 +186,7 @@ def getDPL_global_options(bigshm=False,nosmallrate=False):
 
    PTTRIGMIN=float(args.ptTrigMin)  
    PTTRIGMAX=float(args.ptTrigMax) 
+   PARTICLE_ACCEPTANCE=int(args.acceptance)
 
    ## Pt Hat productions
    WEIGHTPOW=int(args.weightPow)
@@ -219,24 +231,41 @@ def getDPL_global_options(bigshm=False,nosmallrate=False):
            PTHATMAX=hig_edge[PTHATBIN]
 
    # translate here collision type to PDG
-   # not sure this is what we want to do (GCB)
    COLTYPE=args.col
 
    if COLTYPE == 'pp':
       PDGA=2212 # proton
       PDGB=2212 # proton
 
    if COLTYPE == 'PbPb':
-      PDGA=2212 # Pb???? #---> to be checked (seems same as pp case)
-      PDGB=2212 # Pb????
+      PDGA=1000822080 # Pb
+      PDGB=1000822080 # Pb
+      if ECMS < 0:    # assign 5.02 TeV to Pb-Pb
+         print('o2dpg_sim_workflow: Set CM Energy to PbPb case 5.02 TeV')
+         ECMS=5020.0
 
    if COLTYPE == 'pPb':
-      PDGA=2212 # proton
-      PDGB=2212 # Pb????
+      PDGA=2212       # proton
+      PDGB=1000822080 # Pb
 
    if COLTYPE == 'Pbp':
-      PDGA=2212 # Pb????
-      PDGB=2212 # proton
+      PDGA=1000822080 # Pb
+      PDGB=2212       # proton
+
+   # If not set previously, set beam energy B equal to A
+   if EBEAMB < 0 and ECMS < 0:
+      EBEAMB=EBEAMA
+      print('o2dpg_sim_workflow: Set beam energy same in A and B beams')
+      if COLTYPE=="pPb" or COLTYPE=="Pbp":
+         print('o2dpg_sim_workflow: Careful! both beam energies are the same')
+
+   if ECMS > 0:
+      if COLTYPE=="pPb" or COLTYPE=="Pbp":
+         print('o2dpg_sim_workflow: Careful! ECM set for pPb/Pbp collisions!')
+
+   if ECMS < 0 and EBEAMA < 0 and EBEAMB < 0:
+      print('o2dpg_sim_workflow: Error! CM or Beam Energy not set!!!')
+      exit(1)
 
    # produce the signal configuration
    SGN_CONFIG_task=createTask(name='gensgnconf_'+str(tf), tf=tf, cwd=timeframeworkdir)

MC/config/PWGGAJE/hooks/jets_hook.C

@@ -7,8 +7,8 @@ R__ADD_INCLUDE_PATH($O2DPG_ROOT)
 
 /// \brief Select jet events within acceptance or associated parton flavor using Pythia Hooks.
 ///
-/// Select outoging jets on the 2->2 process, at least one in a selected acceptance and 
-/// with a  given PDG value.
+/// Select outoging parton/jets on the 2->2 process, at least one in a selected acceptance and
+/// optionally select the parton with a given PDG value.
 /// Only valid for PYTHIA8 and using Hooks
 ///
 /// \author Gustavo Conesa Balbastre (LPSC-IN2P3-CNRS)
@@ -22,6 +22,14 @@ class UserHooks_jets : public Pythia8::UserHooks
   ~UserHooks_jets() = default;
   bool canVetoPartonLevel() override { return true; };
   bool doVetoPartonLevel(const Pythia8::Event& event) override {
+
+//    for (Int_t id = 0; id < 10; id++)
+//    {
+//      printf("parton %d PDG %d, status %d, mother %d, E %2.2f, pT %2.2f, eta %2.2f, phi %2.2f\n",
+//             id, event[id].id() , event[id].status(), event[id].mother1(),
+//             event[id].e()  , event[id].pT(),
+//             event[id].eta(), event[id].phi()*TMath::RadToDeg());
+//    }
 
     // Get the outgoing 2->2 partons. 
     // The jets are in position 5 or 6.
@@ -37,39 +45,37 @@ class UserHooks_jets : public Pythia8::UserHooks
     if ( mOutPartonPDG > 0 && TMath::Abs(event[id1].id()) != mOutPartonPDG ) 
       okpdg1 = false;
 
-    if ( acc1 && okpdg1 )
-    {
-      printf("--- Accepted event, jet 1 ---\n");
-      printf("PDG %d, status %d, mother %d, E %2.2f, pT %2.2f, eta %2.2f, phi %2.2f\n", 
-             event[id1].id() , event[id1].status(), event[id1].mother1(),
-             event[id1].e()  , event[id1].pT(),
-             event[id1].eta(), event[id1].phi()*TMath::RadToDeg());
-      
-      return false;
-    }
-    
     // Check the second jet
     //
     bool acc2 = detector_acceptance(mAcceptance, event[id2].phi(), event[id2].eta())  ; 
     bool okpdg2 = true;
 
     if ( mOutPartonPDG > 0 && TMath::Abs(event[id2].id()) != mOutPartonPDG ) 
       okpdg2 = false;
-    
-    if ( acc2 && okpdg2 )
+
+    //printf("acc1 %d, acc2 %d, okpdg1 %d, okpdf2 %d\n",acc1,acc2,okpdg1,okpdg2);
+
+    if ( (acc1 || acc2) && (okpdg1 || okpdg2) )
     {
-      printf("--- Accepted event, jet 2 ---\n");
-      printf("PDG %d, status %d, mother %d, E %2.2f, pT %2.2f, eta %2.2f, phi %2.2f\n", 
+      printf("--- Accepted event ---\n");
+      printf("\t --- jet 1 ---\n");
+      printf("\t PDG %d, status %d, mother %d, E %2.2f, pT %2.2f, eta %2.2f, phi %2.2f\n",
+             event[id1].id() , event[id1].status(), event[id1].mother1(),
+             event[id1].e()  , event[id1].pT(),
+             event[id1].eta(), event[id1].phi()*TMath::RadToDeg());
+
+      printf("\t --- jet 2 ---\n");
+      printf("\t PDG %d, status %d, mother %d, E %2.2f, pT %2.2f, eta %2.2f, phi %2.2f\n",
              event[id2].id() , event[id2].status(), event[id2].mother1(),
              event[id2].e()  , event[id2].pT(),
              event[id2].eta(), event[id2].phi()*TMath::RadToDeg());
-      
+
       return false;
     }
 
-    // Jets not found
+    // Jets rejected
     //
-    printf("--- Rejected event ---\n");
+    printf("\t --- Rejected event ---\n");
 
     return true;
 
@@ -97,9 +103,9 @@ pythia8_userhooks_jets(int acc = 0, int pdgPar = 0)
     printf("Select outgoing partons with pdg = %d\n",pdgPar);
   }
 
-  if ( !acc && gSystem->Getenv("CONFIG_DETECTOR_ACCEPTANCE") )
+  if ( !acc && gSystem->Getenv("PARTICLE_ACCEPTANCE") )
   {
-    acc = atoi(gSystem->Getenv("CONFIG_DETECTOR_ACCEPTANCE"));
+    acc = atoi(gSystem->Getenv("PARTICLE_ACCEPTANCE"));
     printf("Requested acceptance %d\n",acc);
   }
 

MC/config/PWGGAJE/hooks/prompt_gamma_hook.C

@@ -121,9 +121,9 @@ Pythia8::UserHooks*
     printf("Select outgoing partons with pdg = %d\n",pdgPar);
   }
 
-  if ( !acc && gSystem->Getenv("CONFIG_DETECTOR_ACCEPTANCE") )
+  if ( !acc && gSystem->Getenv("PARTICLE_ACCEPTANCE") )
   {
-    acc = atoi(gSystem->Getenv("CONFIG_DETECTOR_ACCEPTANCE"));
+    acc = atoi(gSystem->Getenv("PARTICLE_ACCEPTANCE"));
     printf("Requested acceptance %d\n",acc);
   }
 

MC/config/PWGGAJE/trigger/decay_gamma_jets.C

@@ -27,8 +27,8 @@ o2::eventgen::Trigger decay_gamma_jets( )
     // Select photons within acceptance
     //
     Int_t acceptance = 0;
-    if ( gSystem->Getenv("CONFIG_DETECTOR_ACCEPTANCE") )
-      acceptance = atoi(gSystem->Getenv("CONFIG_DETECTOR_ACCEPTANCE"));
+    if ( gSystem->Getenv("PARTICLE_ACCEPTANCE") )
+      acceptance = atoi(gSystem->Getenv("PARTICLE_ACCEPTANCE"));
     //printf("Requested acceptance %d\n",acceptance);
 
     // Particle loop

MC/config/PWGGAJE/trigger/prompt_gamma.C

@@ -71,8 +71,8 @@ o2::eventgen::Trigger prompt_gamma( )
     // Select photons within acceptance
     //
     Int_t acceptance = 0;
-    if ( gSystem->Getenv("CONFIG_DETECTOR_ACCEPTANCE") )
-      acceptance = atoi(gSystem->Getenv("CONFIG_DETECTOR_ACCEPTANCE"));
+    if ( gSystem->Getenv("PARTICLE_ACCEPTANCE") )
+      acceptance = atoi(gSystem->Getenv("PARTICLE_ACCEPTANCE"));
     //printf("Requested acceptance %d\n",acceptance);
 
     if ( detector_acceptance(acceptance, gamma.Phi(),gamma.Eta()) ) 

MC/config/common/pythia8/utils/mkpy8cfg.py

@@ -12,13 +12,19 @@
                    help='The random seed')
 
 parser.add_argument('--idA', type=int, default='2212',
-                   help='PDG code of projectile beam')
+                   help='PDG code of projectile beam A')
 
 parser.add_argument('--idB', type=int, default='2212',
-                   help='PDG code of target beam')
+                   help='PDG code of target beam B')
 
-parser.add_argument('--eCM', type=float, default='13000.',
-                    help='Centre-of-mass energy')
+parser.add_argument('--eA', type=float, default='6499.',
+                    help='Energy of beam A')
+
+parser.add_argument('--eB', type=float, default='6499.',
+                    help='Energy of beam B')
+
+parser.add_argument('--eCM', type=float, default='-1',
+                    help='Centre-of-mass energy (careful!, better use beam energy)')
 
 parser.add_argument('--process', default='inel', choices=['none', 'inel', 'ccbar', 'bbbar', 'heavy', 'jets', 'dirgamma'],
                     help='Process to switch on')
@@ -76,7 +82,14 @@
 fout.write('### beams \n')
 fout.write('Beams:idA = %d \n' % (args.idA))
 fout.write('Beams:idB = %d \n' % (args.idB))
-fout.write('Beams:eCM = %f \n' % (args.eCM))
+if args.eCM > 0:
+   fout.write('Beams:eCM = %f \n' % (args.eCM))
+elif args.eA > 0 and args.eB > 0:
+   fout.write('Beams:eA = %f \n' % (args.eA))
+   fout.write('Beams:eB = %f \n' % (args.eB))
+else:
+   print('mkpy8cfg.py: Error, CM or Beam Energy not set!!!')
+   exit(1)
 fout.write('\n')
 
 ### processes
@@ -94,6 +107,14 @@
     fout.write('PromptPhoton:all = on \n')
 fout.write('\n')
 
+### heavy ion  settings (valid for Pb-Pb 5520 only)
+if args.idA==1000822080 and args.idB==1000822080:
+    fout.write('### heavy-ion settings (valid for Pb-Pb 5520 only) \n')
+    fout.write('HeavyIon:SigFitNGen = 0 \n')
+    fout.write('HeavyIon:SigFitDefPar = 13.88,1.84,0.22,0.0,0.0,0.0,0.0,0.0 \n')
+    fout.write('HeavyIon:bWidth = 14.48 \n')
+fout.write('\n')
+
 ### decays
 fout.write('### decays \n')
 fout.write('ParticleDecays:limitTau0 = on \n') ### we will need to put some parameters for these settings