HF Jet MC - fix weighting, adding new HF jet generator and its test

MC/config/PWGGAJE/ini/GeneratorHFJETrigger_ccbar.ini

@@ -0,0 +1,8 @@
+### The external generator derives from GeneratorPythia8.
+[GeneratorExternal]
+fileName=${O2DPG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
+funcName=GeneratorPythia8GapTriggeredCharm(3, -5, 5)
+
+[GeneratorPythia8]
+config=${O2DPG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_jet_charmtriggers_with_decays.cfg
+

MC/config/PWGGAJE/ini/tests/GeneratorHFJETrigger_ccbar.C

@@ -0,0 +1,125 @@
+int External() {
+  std::string path{"o2sim_Kine.root"};
+  int checkPdgQuark{4};
+  float ratioTrigger = 1. / 3; // one event triggered out of 3
+  std::vector<int> checkPdgHadron{411, 421, 431, 443, 4122, 4132, 4232, 4332};
+  std::map<int, std::vector<std::vector<int>>> checkHadronDecays{
+      // sorted pdg of daughters
+      {411, {{-321, 211, 211}, {-313, 211}, {211, 311}, {211, 333}}}, // D+
+      {421, {{-321, 211}, {-321, 111, 211}}},                         // D0
+      {431, {{211, 333}}},                                            // Ds+
+      {443, {{-11, 11}}},                                             // Jpsi
+      {4122,
+       {{-313, 2212}, {-321, 2224}, {211, 3124}, {-321, 211, 2212}}}, // Lc+
+      {4132, {{211, 3312}}},                                          // Xic0
+      {4232,
+       {{-313, 2212},
+        {-321, 3324},
+        {211, 211, 3312},
+        {-321, 211, 2212}}}, // Xic+
+      {4332, {{211, 3334}}}  // Omegac+
+  };
+
+  TFile file(path.c_str(), "READ");
+  if (file.IsZombie()) {
+    std::cerr << "Cannot open ROOT file " << path << "\n";
+    return 1;
+  }
+
+  auto tree = (TTree *)file.Get("o2sim");
+  std::vector<o2::MCTrack> *tracks{};
+  tree->SetBranchAddress("MCTrack", &tracks);
+
+  int nQuarks{}, nSignals{}, nSignalGoodDecay{};
+  auto nEvents = tree->GetEntries();
+
+  // Setting up event header to access event weight info
+  o2::dataformats::MCEventHeader *eventHeader = nullptr;
+  ;
+  tree->SetBranchAddress("MCEventHeader.", &eventHeader);
+  int weight1Counter = 0;
+  bool isvalid;
+
+  for (int i = 0; i < nEvents; i++) {
+    tree->GetEntry(i);
+
+    float eventWeight = eventHeader->getInfo<float>("weight", isvalid);
+    if (!isvalid) {
+      std::cerr << "Could not retrieve event weight from MCEventHeader\n";
+      return 1;
+    }
+    if (abs(eventWeight - 1) <
+        1E-5) { // checks that the event weights are not unitary
+      weight1Counter++;
+    }
+
+    for (auto &track : *tracks) {
+      auto pdg = track.GetPdgCode();
+      if (std::abs(pdg) == checkPdgQuark) {
+        nQuarks++;
+        continue;
+      }
+      if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(),
+                    std::abs(pdg)) != checkPdgHadron.end()) // found signal
+      {
+        // count signal PDG
+        nSignals++;
+
+        std::vector<int> pdgsDecay{};
+        std::vector<int> pdgsDecayAntiPart{};
+        for (int j{track.getFirstDaughterTrackId()};
+             j <= track.getLastDaughterTrackId(); ++j) {
+          auto pdgDau = tracks->at(j).GetPdgCode();
+          pdgsDecay.push_back(pdgDau);
+          if (pdgDau != 333) { // phi is antiparticle of itself
+            pdgsDecayAntiPart.push_back(-pdgDau);
+          } else {
+            pdgsDecayAntiPart.push_back(pdgDau);
+          }
+        }
+
+        std::sort(pdgsDecay.begin(), pdgsDecay.end());
+        std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end());
+
+        for (auto &decay : checkHadronDecays[std::abs(pdg)]) {
+          if (pdgsDecay == decay || pdgsDecayAntiPart == decay) {
+            nSignalGoodDecay++;
+            break;
+          }
+        }
+      }
+    }
+  }
+  std::cout << "--------------------------------\n";
+  std::cout << "# Events: " << nEvents << "\n";
+  std::cout << Form("# %d (anti)quarks: ", checkPdgQuark) << nQuarks << "\n";
+  std::cout << "# signal hadrons: " << nSignals << "\n";
+  std::cout << "# signal hadrons decaying in the correct channel: "
+            << nSignalGoodDecay << "\n";
+
+  if (nQuarks <
+      2 * nEvents *
+          ratioTrigger) // we expect anyway more because the same quark is
+                        // repeated several time, after each gluon radiation
+  {
+    std::cerr << "Number of generated (anti)quarks " << checkPdgQuark
+              << " lower than expected\n";
+    return 1;
+  }
+
+  float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
+  if (fracForcedDecays < 0.9) // we put some tolerance (e.g. due to oscillations
+                              // which might change the final state)
+  {
+    std::cerr << "Fraction of signals decaying into the correct channel "
+              << fracForcedDecays << " lower than expected\n";
+    return 1;
+  }
+
+  if (weight1Counter == nEvents) {
+    std::cerr << "All events have a unitary weight\n";
+    return 1;
+  }
+
+  return 0;
+}

MC/config/PWGGAJE/pythia8/generator/pythia8_jet_charmtriggers_with_decays.cfg

@@ -0,0 +1,117 @@
+### author: Aimeric Landou (aimeric.landou@cern.ch)
+### based on Fabrizio Grosa's HF template (fabrizio.grosa@cern.ch)
+### since:  January 2023
+
+### beams
+Beams:idA 2212			# proton
+Beams:idB 2212 			# proton
+Beams:eCM 13600. 		# GeV
+
+# ### processes
+# SoftQCD:inelastic on		# all inelastic processes
+
+### decays
+ParticleDecays:limitTau0 on	
+ParticleDecays:tau0Max 10.
+
+### Force golden charm hadrons decay modes for trigger studies
+
+### add D+ decays absent in PYTHIA8 decay table and set BRs from PDG for other
+411:oneChannel = 1 0.0752 0 -321 211 211
+411:addChannel = 1 0.0104 0 -313 211
+411:addChannel = 1 0.0156 0 311 211
+411:addChannel = 1 0.00276 0 333 211
+## add Lc decays absent in PYTHIA8 decay table and set BRs from PDG for other
+4122:oneChannel = 1 0.0196 100 2212 -313
+4122:addChannel = 1 0.0108 100 2224 -321
+4122:addChannel = 1 0.022 100 3124 211
+4122:addChannel = 1 0.035 0 2212 -321 211
+### add Xic+ decays absent in PYTHIA8 decay table
+4232:addChannel = 1 0.2 0 2212 -313
+4232:addChannel = 1 0.2 0 2212 -321 211
+4232:addChannel = 1 0.2 0 3324 211
+4232:addChannel = 1 0.2 0 3312 211 211
+### add Xic0 decays absent in PYTHIA8 decay table
+4132:addChannel = 1 0.2 0 3312 211
+
+### K* -> K pi
+313:onMode = off
+313:onIfAll = 321 211
+### for Ds -> Phi pi+
+333:onMode = off
+333:onIfAll = 321 321
+### for D0 -> rho0 pi+ k-
+113:onMode = off
+113:onIfAll = 211 211
+### for Lambda_c -> Delta++ K-
+2224:onMode = off
+2224:onIfAll = 2212 211
+### for Lambda_c -> Lambda(1520) K-
+102134:onMode = off
+102134:onIfAll = 2212 321
+
+### switch off all decay channels
+411:onMode = off
+421:onMode = off
+431:onMode = off
+4112:onMode = off
+4122:onMode = off
+4232:onMode = off
+4132:onMode = off
+443:onMode = off
+4332:onMode = off
+
+### D0 -> K pi
+421:onIfMatch = 321 211
+
+### D+/- -> K pi pi 
+411:onIfMatch = 321 211 211
+### D+/- -> K* pi
+411:onIfMatch = 313 211
+### D+/- -> phi pi
+411:onIfMatch = 333 211
+
+### D_s -> Phi pi
+431:onIfMatch = 333 211
+
+### Lambda_c -> p K*
+4122:onIfMatch = 2212 313
+### Lambda_c -> Delta K
+4122:onIfMatch = 2224 321
+### Lambda_c -> Lambda(1520) pi
+4122:onIfMatch = 3124 211
+### Lambda_c -> p K pi
+4122:onIfMatch = 2212 321 211
+
+### Xic+ -> pK*0
+4232:onIfMatch = 2212 313
+### Xic+ -> p K- pi+
+4232:onIfMatch = 2212 321 211
+### Xic+ -> Xi*0 pi+, Xi*->Xi- pi+
+4232:onIfMatch = 3324 211
+### Xic+ -> Xi- pi+ pi+
+4232:onIfMatch = 3312 211 211
+
+### Xic0 -> Xi- pi+
+4132:onIfMatch = 3312 211
+
+### Omega_c -> Omega pi
+4332:onIfMatch = 3334 211
+
+### Jpsi -> ee
+443:onIfMatch = 11 11
+
+# Correct Lb decay length (wrong in PYTHIA8 decay table)
+5122:tau0 = 4.41000e-01
+
+# Correct OmegaC decay length (wrong in PYTHIA8 decay table)
+4332:tau0 = 0.08000000000
+
+
+### Jets options
+PhaseSpace:pTHatMin = 5.0
+PhaseSpace:pTHatMax = 300.0
+PhaseSpace:bias2Selection = on
+PhaseSpace:bias2SelectionPow = 6.0
+SoftQCD:inelastic off		# all inelastic processes
+HardQCD:all on

MC/run/PWGGAJE/run_jets_HF_ccbar.sh

@@ -15,42 +15,21 @@
 
 RNDSEED=${RNDSEED:-0}   # [default = 0] time-based random seed
 
-NSIGEVENTS=${NSIGEVENTS:-10}
+NSIGEVENTS=${NSIGEVENTS:-5}
 NTIMEFRAMES=${NTIMEFRAMES:-1}
 NWORKERS=${NWORKERS:-8}
 MODULES="--skipModules ZDC" #"PIPE ITS TPC EMCAL"
 CONFIG_ENERGY=${CONFIG_ENERGY:-13600.0}
 SIMENGINE=${SIMENGINE:-TGeant4}
-WEIGHTPOW=${WEIGHTPOW:-6.0}
 [[ ${SPLITID} != "" ]] && SEED="-seed ${SPLITID}" || SEED=""
 
-# Default for weighted productions
-PTHATMIN=${PTHATMIN:-5.0}
-PTHATMAX=${PTHATMAX:-300.0}
-
-# Define the pt hat bin arrays
-pthatbin_loweredges=(0 5 7 9 12 16 21 28 36 45 57 70 85 99 115 132 150 169 190 212 235)
-pthatbin_higheredges=( 5 7 9 12 16 21 28 36 45 57 70 85 99 115 132 150 169 190 212 235 -1)
-
-# Recover environmental vars for pt binning
-#PTHATBIN=${PTHATBIN:-1}
-
-if [ -z "$PTHATBIN" ]; then
-    echo "Open Pt-hat range set"
-else
-  PTHATMIN=${pthatbin_loweredges[$PTHATBIN]}
-  PTHATMAX=${pthatbin_higheredges[$PTHATBIN]}
-fi
-
-
-#ccbar filter
-${O2DPG_ROOT}/MC/bin/o2dpg_sim_workflow.py -eCM ${CONFIG_ENERGY} -col pp -gen external -proc "jets"           \
-                                            -ptHatMin ${PTHATMIN} -ptHatMax ${PTHATMAX}                       \
-                                            -tf ${NTIMEFRAMES} -ns ${NSIGEVENTS} -e ${SIMENGINE}              \
-                                            -j ${NWORKERS} -mod "--skipModules ZDC"                           \
-                                            -interactionRate 500000 -confKey "Diamond.width[2]=6." ${SEED}    \
-                                            -ini $O2DPG_ROOT/MC/config/PWGHF/ini/GeneratorHFTrigger_ccbar.ini \
-                                            -weightPow ${WEIGHTPOW}
+#ccbar filter and bias2SelectionPow and PtHat settings are in the ini file given below
+${O2DPG_ROOT}/MC/bin/o2dpg_sim_workflow.py -eCM ${CONFIG_ENERGY} -col pp -gen external -proc "jets"                 \
+                                            -tf ${NTIMEFRAMES} -ns ${NSIGEVENTS} -e ${SIMENGINE}                    \
+                                            -j ${NWORKERS} -mod "--skipModules ZDC"                                 \
+                                            -interactionRate 500000 -confKey "Diamond.width[2]=6." ${SEED}          \
+                                            -ini $O2DPG_ROOT/MC/config/PWGGAJE/ini/GeneratorHFJETrigger_ccbar.ini
+
 
 # run workflow
 # allow increased timeframe parallelism with --cpu-limit 32