Decay generator particles before transport

Generators/CMakeLists.txt

@@ -47,6 +47,8 @@ o2_add_library(Generators
                        $<$<BOOL:${pythia_FOUND}>:src/DecayerPythia8.cxx>
                        $<$<BOOL:${pythia_FOUND}>:src/GeneratorPythia8Param.cxx>
                        $<$<BOOL:${pythia_FOUND}>:src/DecayerPythia8Param.cxx>
+                       $<$<BOOL:${pythia_FOUND}>:src/EventDecayer.cxx>
+                       $<$<BOOL:${pythia_FOUND}>:src/EventDecayerParam.cxx>
                        $<$<BOOL:${HepMC3_FOUND}>:src/GeneratorHepMC.cxx>
                        $<$<BOOL:${HepMC3_FOUND}>:src/GeneratorHepMCParam.cxx>
                        $<$<BOOL:${HepMC3_FOUND}>:src/AODToHepMC.cxx>
@@ -103,6 +105,7 @@ if(pythia_FOUND)
               include/Generators/DecayerPythia8.h
               include/Generators/GeneratorPythia8Param.h
               include/Generators/DecayerPythia8Param.h
+              include/Generators/EventDecayerParam.h
               include/Generators/GeneratorFactory.h)
 endif()
 
@@ -147,6 +150,22 @@ if(doBuildSimulation)
                 PUBLIC_LINK_LIBRARIES O2::Generators)
   endif()
 
+  if(pythia_FOUND)
+    o2_add_test(EventDecayer NAME test_Generator_test_EventDecayer
+                SOURCES test/test_EventDecayer.cxx
+                COMPONENT_NAME Generator
+                ENVIRONMENT O2_ROOT=${CMAKE_BINARY_DIR}/stage
+                LABELS generator
+                PUBLIC_LINK_LIBRARIES O2::Generators)
+
+    o2_add_test(GeneratorEventDecayer NAME test_Generator_test_GeneratorEventDecayer
+                SOURCES test/test_GeneratorEventDecayer.cxx
+                COMPONENT_NAME Generator
+                ENVIRONMENT O2_ROOT=${CMAKE_BINARY_DIR}/stage
+                LABELS generator
+                PUBLIC_LINK_LIBRARIES O2::Generators)
+  endif()
+
   # o2_add_test(GeneratorPythia8Param NAME test_Generator_test_GeneratorPythia8Param
   #             SOURCES test/test_GeneratorPythia8Param.cxx
   #             COMPONENT_NAME Generator

Generators/include/Generators/EventDecayer.h

@@ -0,0 +1,90 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#ifndef ALICEO2_EVENTGEN_EVENTDECAYER_H_
+#define ALICEO2_EVENTGEN_EVENTDECAYER_H_
+
+#include "Pythia8/Pythia.h"
+#include <TParticle.h>
+#include <map>
+#include <string>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+namespace o2
+{
+namespace eventgen
+{
+
+/// Counts of what one or more EventDecayer::decay calls did, per PDG code
+struct EventDecayerReport {
+  std::map<int, long> decayed;                               // particles decayed by Pythia8
+  std::map<int, long> failed;                                // particles Pythia8 could not decay
+  std::map<int, long> unknownInO2;                           // particles to transport unknown to O2DatabasePDG
+  std::map<int, long> unknownInG4;                           // particles to transport unknown to Geant4
+  std::map<int, std::pair<double, double>> lifetimeMismatch; // pdg -> {pythia ctau mm, g4 ctau mm}
+  void merge(EventDecayerReport const& other);
+  void print() const;
+};
+
+/// One entry of the Geant4 particle table
+struct G4ParticleInfo {
+  std::string name;
+  double ctauMM;
+  bool shortLived;
+  bool hasDecayTable;
+};
+
+/// Reads the Geant4 particle table written by o2-sim-g4-particle-table
+std::unordered_map<int, G4ParticleInfo> loadG4ParticleTable(std::string const& path);
+
+/// Decays unstable generator particles with Pythia8 before transport, keeping the full chain
+class EventDecayer
+{
+ public:
+  bool init();
+  /// Sets the caller's units in cm, s, GeV/c and GeV, as the unit setters of Generator (default cm, s, GeV/c, GeV)
+  void setUnits(double positionUnit, double timeUnit, double momentumUnit = 1., double energyUnit = 1.)
+  {
+    mPositionUnit = positionUnit;
+    mTimeUnit = timeUnit;
+    mMomentumUnit = momentumUnit;
+    mEnergyUnit = energyUnit;
+  }
+  EventDecayerReport decay(std::vector<TParticle>& particles);
+  std::string configSummary() const;
+  bool vertexLimitsActive() const;
+  int seed() const { return mSeed; }
+
+ private:
+  bool decayable(int pdg) const;
+  bool decayOne(std::vector<TParticle>& particles, int index);
+  void checkParticles(std::vector<TParticle> const& particles, EventDecayerReport& report) const;
+  void checkLifetimes(EventDecayerReport& report) const;
+
+  mutable ::Pythia8::Pythia mPythia;                    // mutable since Pythia8 settings lookups are not const
+  std::vector<std::string> mConfigs;                    // expanded config paths
+  std::unordered_map<int, G4ParticleInfo> mG4Particles; // Geant4 particle table
+  bool mHasG4Table = false;
+  bool mVerbose = false;
+  double mPositionUnit = 1.; // caller's position unit in cm
+  double mTimeUnit = 1.;     // caller's time unit in s
+  double mMomentumUnit = 1.; // caller's momentum unit in GeV/c
+  double mEnergyUnit = 1.;   // caller's energy unit in GeV
+  double mTau0Max = 0.;      // tau0 limit in mm/c in use after init, 0 when off
+  int mSeed = 0;             // Pythia8 seed in use after init
+};
+
+} // namespace eventgen
+} // namespace o2
+
+#endif // ALICEO2_EVENTGEN_EVENTDECAYER_H_

Generators/include/Generators/EventDecayerParam.h

@@ -0,0 +1,40 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#ifndef ALICEO2_EVENTGEN_EVENTDECAYERPARAM_H_
+#define ALICEO2_EVENTGEN_EVENTDECAYERPARAM_H_
+
+#include "CommonUtils/ConfigurableParam.h"
+#include "CommonUtils/ConfigurableParamHelper.h"
+#include <string>
+
+namespace o2
+{
+namespace eventgen
+{
+
+/// Settings of the EventDecayer, which decays generator particles with Pythia8 before transport
+struct EventDecayerParam : public o2::conf::ConfigurableParamHelper<EventDecayerParam> {
+  bool enable = false;                                                               // decay step in Generator::ReadEvent
+  std::string config[8] = {"", "", "", "", "", "", "", ""};                          // Pythia8 config files, an empty slot takes DecayerPythia8.config
+  double tau0Max = 10.;                                                              // species with a larger tau0 in mm/c stay undecayed (<= 0 switches it off)
+  std::string skipPdg = "115 225 335 10221 10331 20223 20333 102132 102134 9010221"; // PDG codes never decayed, matched on |pdg|
+  std::string g4ParticleTable = "";                                                  // Geant4 particle table from o2-sim-g4-particle-table
+  int seed = 0;                                                                      // Pythia8 random seed (0 derives it from gRandom)
+  bool verbose = false;                                                              // list every Pythia8 decay
+  bool showChanged = false;                                                          // show changed particle data at init
+  O2ParamDef(EventDecayerParam, "EventDecayer");
+};
+
+} // end namespace eventgen
+} // end namespace o2
+
+#endif // ALICEO2_EVENTGEN_EVENTDECAYERPARAM_H_

Generators/include/Generators/Generator.h

@@ -35,6 +35,8 @@ namespace eventgen
 {
 
 class GenTPCLoopers; // Forward declaration
+class EventDecayer;
+struct EventDecayerReport;
 
 /*****************************************************************/
 /*****************************************************************/
@@ -181,6 +183,11 @@ class Generator : public FairGenerator
 
   // Loopers generator instance
   o2::eventgen::GenTPCLoopers* mTPCLoopersGen = nullptr;
+
+  // decay step after the trigger, created on first use when EventDecayer.enable is set
+  o2::eventgen::EventDecayer* mEventDecayer = nullptr;             //!
+  o2::eventgen::EventDecayerReport* mEventDecayerReport = nullptr; //!
+  Bool_t decayParticles(o2::dataformats::MCEventHeader* header);
 #ifdef GENERATORS_WITH_TPCLOOPERS
   bool initTPCLoopersGen();
 #endif

Generators/share/pythia8/decays/base.cfg

@@ -32,3 +32,14 @@ Init:showChangedParticleData off
 4232:mayDecay off	# Xi_c+
 4332:mayDecay off	# Omega_c0
 5122:mayDecay off	# Lambda_b0
+
+### lifetimes (PDG 2026), Pythia defaults are older measurements
+
+531:tau0  = 0.4542	# B_s0
+541:tau0  = 0.1529	# B_c+
+4132:tau0 = 0.04497	# Xi_c0
+4232:tau0 = 0.1358	# Xi_c+
+4332:tau0 = 0.08214	# Omega_c0
+5132:tau0 = 0.4707	# Xi_b-
+5232:tau0 = 0.4413	# Xi_b0
+5332:tau0 = 0.4917	# Omega_b-