[PWGLF] Add hybrid PYTHIA8 generator for antideuteron coalescence in jets

Author: alcaliva

MC/config/PWGLF/hybrid/hybrid_coalescence.json

@@ -0,0 +1,24 @@
+
+{
+  "mode": "sequential",
+  "generators": [
+    {
+      "name": "pythia8",
+      "config": {
+        "config": "${O2_ROOT}/share/Generators/egconfig/pythia8_inel.cfg",
+        "hooksFileName": "",
+        "hooksFuncName": "",
+        "includePartonEvent": false,
+        "particleFilter": "",
+        "verbose": 0
+      }
+    },
+    {
+      "name": "external",
+      "config": {
+        "iniFile": "${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.ini"
+      }
+    }
+  ],
+  "fractions": [4, 1]
+}

MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.ini

@@ -0,0 +1,6 @@
+[GeneratorExternal]
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_coalescence.C
+funcName = generateCoalescence({-1000010020}, 1, 0.4)
+
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_jet.cfg

MC/config/PWGLF/ini/tests/GeneratorAntideuteronCoalescenceJet.C

@@ -0,0 +1,85 @@
+
+#if !defined(__CLING__) || defined(__ROOTCLING__)
+#include "SimulationDataFormat/MCTrack.h"
+#include "TFile.h"
+#include "TTree.h"
+
+#include <algorithm>
+#include <iostream>
+#include <string>
+#include <vector>
+#endif
+
+int External()
+{
+  std::string path{"o2sim_Kine.root"};
+
+  std::vector<int> injectedPDGs = {
+    -1000010020 // anti-deuteron
+  };
+
+  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");
+  if (!tree) {
+    std::cerr << "Cannot find tree o2sim in file " << path << "\n";
+    return 1;
+  }
+
+  std::vector<o2::MCTrack>* tracks = nullptr;
+  tree->SetBranchAddress("MCTrack", &tracks);
+
+  std::vector<int> nSignal(injectedPDGs.size(), 0);
+  std::vector<int> nEventsWithSignal(injectedPDGs.size(), 0);
+
+  const auto nEvents = tree->GetEntries();
+
+  for (Long64_t i = 0; i < nEvents; ++i) {
+    tree->GetEntry(i);
+
+    std::vector<bool> found(injectedPDGs.size(), false);
+
+    for (const auto& track : *tracks) {
+      const int pdg = track.GetPdgCode();
+
+      auto it = std::find(
+          injectedPDGs.begin(),
+          injectedPDGs.end(),
+          pdg);
+
+      if (it != injectedPDGs.end()) {
+        const auto index = std::distance(injectedPDGs.begin(), it);
+        ++nSignal[index];
+        found[index] = true;
+      }
+    }
+
+    for (size_t j = 0; j < injectedPDGs.size(); ++j) {
+      if (found[j]) {
+        ++nEventsWithSignal[j];
+      }
+    }
+  }
+
+  std::cout << "--------------------------------\n";
+  std::cout << "# Events: " << nEvents << "\n";
+
+  for (size_t i = 0; i < injectedPDGs.size(); ++i) {
+    std::cout << "PDG: " << injectedPDGs[i] << "\n";
+    std::cout << "# Generated particles: " << nSignal[i] << "\n";
+    std::cout << "# Events containing particle: "
+              << nEventsWithSignal[i] << "\n";
+
+    if (nSignal[i] == 0) {
+      std::cerr << "No generated particles with PDG "
+                << injectedPDGs[i] << "\n";
+      return 1;
+    }
+  }
+
+  return 0;
+}

MC/config/PWGLF/pythia8/generator_pythia8_coalescence.C

@@ -9,6 +9,7 @@
 #include "TDatabasePDG.h"
 #include "TSystem.h"
 #include "TMath.h"
+#include <algorithm>
 #include <cmath>
 #include <vector>
 #include <fstream>
@@ -18,20 +19,26 @@ using namespace Pythia8;
 
 R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT)
 #include "MC/config/common/external/generator/CoalescencePythia8.h"
-/// First version of the simple coalescence generator based PYTHIA8
 
+/// First version of the simple coalescence generator based on PYTHIA8
 class GeneratorPythia8Coalescence : public o2::eventgen::GeneratorPythia8
 {
 public:
   /// Constructor
-  GeneratorPythia8Coalescence(std::vector<unsigned int> pdgList, int input_trigger_ratio = 1, double coal_momentum = 0.4)
+  GeneratorPythia8Coalescence(std::vector<int> pdgList, int input_trigger_ratio = 1, double coal_momentum = 0.4)
       : o2::eventgen::GeneratorPythia8()
   {
     fmt::printf(">> Coalescence generator %d\n", input_trigger_ratio);
     mInverseTriggerRatio = input_trigger_ratio;
     mCoalMomentum = coal_momentum;
     mPdgList = pdgList;
+
+    // CoalescencePythia8 accepts unsigned PDG codes
+    for (auto pdg : mPdgList) {
+      mAbsPdgList.push_back(static_cast<unsigned int>(std::abs(pdg)));
+    }
   }
+
   /// Destructor
   ~GeneratorPythia8Coalescence() = default;
 
@@ -45,55 +52,86 @@ public:
 protected:
   bool generateEvent() override
   {
-    fmt::printf(">> Generating event %d\n", mGeneratedEvents);
+    fmt::printf(">> Generating event %llu\n",
+                static_cast<unsigned long long>(mGeneratedEvents));
+
     // Simple straightforward check to alternate generators
     if (mGeneratedEvents % mInverseTriggerRatio == 0)
     {
-      fmt::printf(">> Generating coalescence event %d\n", mGeneratedEvents);
+      fmt::printf(">> Generating coalescence event %llu\n",
+                  static_cast<unsigned long long>(mGeneratedEvents));
+
       bool genOk = false;
       int localCounter{0};
+
       while (!genOk)
       {
         if (GeneratorPythia8::generateEvent())
         {
-          genOk = CoalescencePythia8(mPythia.event, mPdgList, mCoalMomentum);
+          // Store the event size before coalescence
+          const int originalSize = mPythia.event.size();
+
+          // Run coalescence with absolute PDG codes
+          bool coalOk = CoalescencePythia8(
+              mPythia.event, mAbsPdgList, false, mCoalMomentum);
+
+          if (coalOk && mPythia.event.size() > originalSize) {
+            for (int i = originalSize; i < mPythia.event.size(); ++i) {
+              const int producedPdg = mPythia.event[i].id();
+              if (std::find(mPdgList.begin(), mPdgList.end(), producedPdg) != mPdgList.end()) {
+                genOk = true;
+                break;
+              }
+            }
+          }
         }
+
         localCounter++;
       }
-      fmt::printf(">> Coalescence successful after %i generations\n", localCounter);
+
+      fmt::printf(">> Coalescence successful after %i generations\n",
+                  localCounter);
       std::cout << std::endl << std::endl;
       notifySubGenerator(1);
     }
     else
     {
-      fmt::printf(">> Generating minimum-bias event %d\n", mGeneratedEvents);
+      fmt::printf(">> Generating minimum-bias event %llu\n",
+                  static_cast<unsigned long long>(mGeneratedEvents));
+
       // Generate minimum-bias event
       bool genOk = false;
       while (!genOk)
       {
         genOk = GeneratorPythia8::generateEvent();
       }
+
       notifySubGenerator(0);
     }
+
     mGeneratedEvents++;
     return true;
   }
 
-
 private:
-  std::vector<unsigned int> mPdgList; /// list of pdg codes to be generated
-  float mCoalMomentum = 0.4; /// coalescence momentum
-  uint64_t mGeneratedEvents = 0; /// number of events generated so far
-  int mInverseTriggerRatio = 1;  /// injection gap
+  std::vector<int> mPdgList;                /// Signed PDG codes
+  std::vector<unsigned int> mAbsPdgList;    /// Absolute PDG codes for coalescence
+  float mCoalMomentum = 0.4;                /// Coalescence momentum
+  uint64_t mGeneratedEvents = 0;           /// Number of generated events
+  int mInverseTriggerRatio = 1;            /// Injection gap
 };
 
-///___________________________________________________________
-FairGenerator *generateCoalescence(std::vector<unsigned int> pdgList, int input_trigger_ratio, double coal_momentum = 0.4)
+/// Factory function
+FairGenerator* generateCoalescence(std::vector<int> pdgList,
+                                  int input_trigger_ratio,
+                                  double coal_momentum = 0.4)
 {
-  auto myGen = new GeneratorPythia8Coalescence(pdgList, input_trigger_ratio, coal_momentum);
+  auto myGen = new GeneratorPythia8Coalescence(
+      pdgList, input_trigger_ratio, coal_momentum);
+
   auto seed = (gRandom->TRandom::GetSeed() % 900000000);
   myGen->readString("Random:setSeed on");
   myGen->readString("Random:seed " + std::to_string(seed));
+
   return myGen;
 }
-///___________________________________________________________