Add pythia8 box gun generator

MC/config/PWGHF/ini/GeneratorHFOmegaCEmb.ini

@@ -0,0 +1,6 @@
+[GeneratorExternal]
+fileName=${O2DPG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_box.C
+funcName=generatePythia8Box(4332, 3)
+
+[GeneratorPythia8]
+config=${O2DPG_ROOT}/MC/config/PWGHF/pythia8/generator/omega_c.cfg

MC/config/PWGHF/pythia8/generator/omega_c.cfg

@@ -0,0 +1,18 @@
+### Omega_c
+#4332:all = Omega_c0 Omega_cbar0 2 0 0 2.69520 0. 0. 0. 10
+
+ProcessLevel:all = off
+
+### changing the ctau value in mm/c
+4332:tau0=0.08000000000
+
+### add OmegaC decay absent in PYTHIA8 decay table
+4332:addChannel = 1 0.0001 0 3334 211
+
+### force the OmegaC to decay in the Omega_c -> Omega pi channel
+4332:onMode = off
+4332:onIfMatch = 3334 211
+
+### switch off Omega and Lambda decay channel (treated in GEANT)
+3334:onMode = off
+3122:onMode = off

MC/config/PWGLF/ini/GeneratorLFOmegaEmb.ini

@@ -0,0 +1,6 @@
+[GeneratorExternal]
+fileName=${O2DPG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_box.C
+funcName=generatePythia8Box(3334, 3)
+
+[GeneratorPythia8]
+config=${O2DPG_ROOT}/MC/config/PWGLF/pythia8/generator/omega.cfg

MC/config/PWGLF/pythia8/generator/omega.cfg

@@ -0,0 +1,3 @@
+### switch off Omega decay channel (treated in GEANT)
+ProcessLevel:all = off
+3334:onMode = off

MC/config/PWGLF/pythia8/generator_pythia8_box.C

@@ -0,0 +1,99 @@
+#include "Pythia8/Pythia.h"
+#include "FairGenerator.h"
+#include "FairPrimaryGenerator.h"
+#include "Generators/GeneratorPythia8.h"
+#include "TRandom3.h"
+#include "TParticlePDG.h"
+#include "TDatabasePDG.h"
+#include "TMath.h"
+#include <cmath>
+
+using namespace Pythia8;
+
+class GeneratorPythia8Box : public o2::eventgen::GeneratorPythia8
+{
+public:
+  /// constructor
+  GeneratorPythia8Box(int input_pdg, int nInject = 1, float ptMin = 1, float ptMax = 10) : pdg{input_pdg}, nParticles{nInject}, genMinPt{ptMin}, genMaxPt{ptMax}, m{getMass(input_pdg)}
+  {
+  }
+
+  ///  Destructor
+  ~GeneratorPythia8Box() = default;
+
+  /// randomize the PDG code sign of core particle
+  void setRandomizePDGsign(bool val) { randomizePDGsign = val; }
+
+  /// get mass from TParticlePDG
+  static double getMass(int input_pdg)
+  {
+    double mass = 0;
+    if (TDatabasePDG::Instance())
+    {
+      TParticlePDG *particle = TDatabasePDG::Instance()->GetParticle(input_pdg);
+      if (particle)
+      {
+        mass = particle->Mass();
+      }
+      else
+      {
+        std::cout << "===> Unknown particle requested with PDG " << input_pdg << ", mass set to 0" << std::endl;
+      }
+    }
+    return mass;
+  }
+
+  Bool_t generateEvent() override
+  {
+    mPythia.event.reset();
+
+    static int sign{1};
+    for (int i{0}; i < nParticles; ++i)
+    {
+      const double pt = gRandom->Uniform(genMinPt, genMaxPt);
+      const double eta = gRandom->Uniform(genMinEta, genMaxEta);
+      const double phi = gRandom->Uniform(0, TMath::TwoPi());
+      const double px{pt * std::cos(phi)};
+      const double py{pt * std::sin(phi)};
+      const double pz{pt * std::sinh(eta)};
+      const double et{std::hypot(std::hypot(pt, pz), m)};
+      sign *= randomizePDGsign ? -1 : 1;
+
+      Particle myparticle;
+      myparticle.id(pdg);
+      myparticle.status(11);
+      myparticle.px(px);
+      myparticle.py(py);
+      myparticle.pz(pz);
+      myparticle.e(et);
+      myparticle.m(m);
+      myparticle.xProd(0);
+      myparticle.yProd(0);
+      myparticle.zProd(0);
+
+      mPythia.event.append(myparticle);
+    }
+    mPythia.next();
+    return true;
+  }
+
+  //__________________________________________________________________
+
+private:
+  double genMinPt = 0.5;  /// minimum 3-momentum for generated particles
+  double genMaxPt = 12.;  /// maximum 3-momentum for generated particles
+  double genMinEta = -1.; /// minimum pseudorapidity for generated particles
+  double genMaxEta = +1.; /// maximum pseudorapidity for generated particles
+
+  double m = 0;       /// particle mass [GeV/c^2]
+  int pdg = 0;        /// particle pdg code
+  int nParticles = 1; /// Number of injected particles
+
+  bool randomizePDGsign = true; /// bool to randomize the PDG code of the core particle
+};
+
+///___________________________________________________________
+FairGenerator *generatePythia8Box(int pdg, int nInject, float ptMin = 1, float ptMax = 10)
+{
+  return new GeneratorPythia8Box(pdg, nInject, ptMin, ptMax);
+}

MC/run/PWGHF/run_OmegaCInjected.sh

@@ -0,0 +1,32 @@
+#!/bin/bash
+
+#
+# Inject OmegaC signal into background with a pythia8 box gun generator
+
+# make sure O2DPG + O2 is loaded
+
+[ ! "${O2DPG_ROOT}" ] && echo "Error: This needs O2DPG loaded" && exit 1
+[ ! "${O2_ROOT}" ] && echo "Error: This needs O2 loaded" && exit 1
+
+# ----------- LOAD UTILITY FUNCTIONS --------------------------
+. ${O2_ROOT}/share/scripts/jobutils.sh
+
+# ----------- START ACTUAL JOB  -----------------------------
+
+NWORKERS=${NWORKERS:-8}
+MODULES="--skipModules ZDC"
+SIMENGINE=${SIMENGINE:-TGeant4}
+NSIGEVENTS=${NSIGEVENTS:-1}
+NBKGEVENTS=${NBKGEVENTS:-1}
+
+NTIMEFRAMES=${NTIMEFRAMES:-1}
+SYSTEM=${SYSTEM:-pp}
+ENERGY=${ENERGY:-13600}
+[[ ${SPLITID} != "" ]] && SEED="-seed ${SPLITID}" || SEED=""
+
+${O2DPG_ROOT}/MC/bin/o2dpg_sim_workflow.py -e ${SIMENGINE} ${SEED} -eCM 13600 -col pp -colBkg pp -gen external -genBkg pythia8 -procBkg inel -j ${NWORKERS} -ns ${NSIGEVENTS} -nb ${NBKGEVENTS} -tf ${NTIMEFRAMES} -interactionRate 500000 -confKey "Diamond.width[2]=6." -mod "--skipModules ZDC" \
+        --embedding -ini $O2DPG_ROOT/MC/config/PWGHF/ini/GeneratorHFOmegaCEmb.ini 
+
+
+# run workflow
+${O2DPG_ROOT}/MC/bin/o2_dpg_workflow_runner.py -f workflow.json -tt aod --cpu-limit 64

MC/run/PWGLF/run_OmegaInjected.sh

@@ -0,0 +1,32 @@
+#!/bin/bash
+
+#
+# Inject Omega signal into background with a pythia8 box gun generator
+
+# make sure O2DPG + O2 is loaded
+
+[ ! "${O2DPG_ROOT}" ] && echo "Error: This needs O2DPG loaded" && exit 1
+[ ! "${O2_ROOT}" ] && echo "Error: This needs O2 loaded" && exit 1
+
+# ----------- LOAD UTILITY FUNCTIONS --------------------------
+. ${O2_ROOT}/share/scripts/jobutils.sh
+
+# ----------- START ACTUAL JOB  -----------------------------
+
+NWORKERS=${NWORKERS:-8}
+MODULES="--skipModules ZDC"
+SIMENGINE=${SIMENGINE:-TGeant4}
+NSIGEVENTS=${NSIGEVENTS:-1}
+NBKGEVENTS=${NBKGEVENTS:-1}
+
+NTIMEFRAMES=${NTIMEFRAMES:-1}
+SYSTEM=${SYSTEM:-pp}
+ENERGY=${ENERGY:-13600}
+[[ ${SPLITID} != "" ]] && SEED="-seed ${SPLITID}" || SEED=""
+
+${O2DPG_ROOT}/MC/bin/o2dpg_sim_workflow.py -e ${SIMENGINE} ${SEED} -eCM 13600 -col pp -colBkg pp -gen external -genBkg pythia8 -procBkg inel -j ${NWORKERS} -ns ${NSIGEVENTS} -nb ${NBKGEVENTS} -tf ${NTIMEFRAMES} -interactionRate 500000 -confKey "Diamond.width[2]=6." -mod "--skipModules ZDC" \
+        --embedding -ini $O2DPG_ROOT/MC/config/PWGLF/ini/GeneratorLFOmegaEmb.ini 
+
+
+# run workflow
+${O2DPG_ROOT}/MC/bin/o2_dpg_workflow_runner.py -f workflow.json -tt aod --cpu-limit 64

test/PWGHF/checkKine.C

@@ -0,0 +1,76 @@
+int checkKine(std::string const &path, int checkPdgSignal, std::vector<int> checkPdgDecays)
+{
+    std::cout << "Check for\nsignal PDG " << checkPdgSignal;
+    for (auto pdg : checkPdgDecays)
+    {
+        std::cout << "\ndecay PDG " << pdg;
+    }
+    std::cout << "\n";
+
+    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 nSignal{};
+    std::vector<int> nDecays(checkPdgDecays.size(), 0);
+    auto nEvents = tree->GetEntries();
+
+    for (int i = 0; i < nEvents; i++)
+    {
+        tree->GetEntry(i);
+        for (auto &track : *tracks)
+        {
+            auto pdg = track.GetPdgCode();
+            if (pdg == checkPdgSignal)
+            {
+                // count signal PDG
+                nSignal++;
+
+                for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j)
+                {
+                    auto pdgDau = tracks->at(j).GetPdgCode();
+                    // count decay PDGs
+                    for (int i = 0, n = checkPdgDecays.size(); i < n; ++i)
+                    {
+                        if (pdgDau == checkPdgDecays[i])
+                        {
+                            nDecays[i]++;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    std::cout << "--------------------------------\n";
+    std::cout << "# Events: " << nEvents << "\n";
+    std::cout << "# Mother " << checkPdgSignal << ": " << nSignal << "\n";
+
+    for (int i = 0; i < checkPdgDecays.size(); i++)
+    {
+        std::cout << "# Daughter " << checkPdgDecays[i] << ": " << nDecays[i] << "\n";
+    }
+
+    if (nSignal != nEvents * 3)
+    {
+        std::cerr << "Number of generated" << checkPdgSignal << "lower than expected\n";
+        return 1;
+    }
+
+    for (int i = 0; i < checkPdgDecays.size(); i++)
+    {
+        if (nDecays[i] != nSignal)
+        {
+            std::cerr << "Number of generated" << checkPdgDecays[i] << "lower than expected\n";
+            return 1;
+        }
+    }
+
+    return 0;
+}

test/test-pythiabox_generator.sh

@@ -0,0 +1,26 @@
+
+#!/bin/bash
+
+########################################################################
+# A very basic test to test a custom generator implementation of PWGDQ #
+########################################################################
+
+[[ -z ${O2_ROOT+x} ]] && { echo "O2_ROOT not loaded." ; exit 1 ; }
+[[ -z ${O2DPG_ROOT+x} ]] && { echo "O2DPG_ROOT not loaded." ; exit 1 ; }
+
+
+# run the workflow up to signal simulation (included)
+# so far that is what we need for the subsequent test on the kinematics file
+NSIGEVENTS=1 NBKGEVENTS=1 NTIMEFRAMES=1 TARGETTASK=sgnsim ${O2DPG_ROOT}/MC/run/PWGHF/run_OmegaCInjected.sh > sim_workflow.log 2>&1
+
+# quit aleady if this didn't work
+[ "$?" != "0" ] && { cat sim_workflow.log ; exit 1 ; }
+
+# now, do some checks, for now only on kinematics
+root -l -b -q ${O2DPG_ROOT}/test/PWGHF/checkKine.C'("tf1/sgn_1_Kine.root", 4332, {3334, 211})' > checkKine.log 2>&1
+[ "$?" != "0" ] && { cat checkKine.log ; exit 1 ; }
+
+# Other potential tests...
+
+
+exit 0