Store the monopole mass in the MCEventHeader and use it when reading kinematics

The monopole mass of a simulation (G4.monopoleMass) was only known to the
simulation itself. Consumers of the output (digitizer, reconstruction MC
labels, QC, analysis) took the compiled-in default of 100 GeV through
O2DatabasePDG, so MCTrack::GetMass(), GetEnergy(), GetRapidity() and
Get4Momentum() were wrong for monopoles of any other mass.

- O2HitMerger: apply the G4.* config key-values (they were not applied in
  the merger process) and write the mass to the MCEventHeader under
  o2::sim::MonopoleMassHeaderKey when G4.monopole is set.
- MCKinematicsReader: take the mass from the header when headers or
  tracks are loaded and hand it to O2DatabasePDG; warn if kinematics of
  different masses are read in the same process.
- O2DatabasePDG: setMonopoleMass()/getMonopoleMass(); Mass() returns the
  set mass for the monopole PDG codes. Unset by default, so the simulation
  itself keeps using TDatabasePDG.

Limitations: entries already registered in a TDatabasePDG keep their mass
(ROOT has no setter); G4 parameters given through --configFile are not
applied in the merger; files produced before this change have no key.

Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>

Author: tubagundem

DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h

@@ -34,6 +34,10 @@ constexpr int MonopolePdgAsymm = 4120000;
 /// Default monopole mass in GeV.
 constexpr double MonopoleMassDefaultGeV = 100.;
 
+/// Key under which the monopole mass (GeV) used in a simulation is stored in the MCEventHeader,
+/// so that readers of the kinematics can use the mass the events were produced with.
+constexpr const char* MonopoleMassHeaderKey = "monopoleMass";
+
 /// true for the monopole and antimonopole of both species
 constexpr bool isMonopole(int pdg) noexcept
 {

DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h

@@ -16,6 +16,7 @@
 #ifndef O2_O2DATABASEPDG_H
 #define O2_O2DATABASEPDG_H
 
+#include <atomic>
 #include <cmath>
 #include <string>
 #include "TDatabasePDG.h"
@@ -39,7 +40,8 @@ class O2DatabasePDG
     static bool initialized = false; // initialize this --> adds particles to TDatabasePDG;
     auto db = TDatabasePDG::Instance();
     if (!initialized) {
-      addALICEParticles(db);
+      // uses the monopole mass of the simulation if a reader has already set it
+      addALICEParticles(db, getMonopoleMass());
       initialized = true;
     }
     return db;
@@ -68,6 +70,26 @@ class O2DatabasePDG
     return true;
   }
 
+  // Sets the mass (GeV) of the monopole species for consumers of simulated data, which do not
+  // have access to G4Params.monopoleMass. It is taken from the MCEventHeader of the
+  // simulation (key o2::sim::MonopoleMassHeaderKey) by MCKinematicsReader.
+  // It only affects Mass() (hence MCTrack::GetMass() / GetEnergy() / GetRapidity()) and
+  // Instance() if called before the first Instance(): particles already registered in a
+  // TDatabasePDG keep their mass, since ROOT offers no way to change it afterwards.
+  // Not used (and not needed) in the simulation itself, where TDatabasePDG is registered
+  // with G4Params.monopoleMass.
+  static void setMonopoleMass(double massGeV) { monopoleMassOverride().store(massGeV); }
+
+  // mass set by setMonopoleMass(), or the default if it was never set
+  static double getMonopoleMass()
+  {
+    const double mass = monopoleMassOverride().load();
+    return mass > 0. ? mass : o2::sim::MonopoleMassDefaultGeV;
+  }
+
+  // true if a monopole mass was set by setMonopoleMass()
+  static bool hasMonopoleMassOverride() { return monopoleMassOverride().load() > 0.; }
+
   // get particle's (if any) mass
   static Double_t MassImpl(TParticlePDG* particle, bool& success)
   {
@@ -82,6 +104,11 @@ class O2DatabasePDG
   // determine particle to get mass for based on PDG
   static Double_t Mass(int pdg, bool& success, TDatabasePDG* db = O2DatabasePDG::Instance())
   {
+    if (hasMonopoleMassOverride() && o2::sim::isMonopole(pdg)) {
+      // the table may have been filled with the default mass before the simulation mass was known
+      success = true;
+      return getMonopoleMass();
+    }
     if (pdg < IONBASELOW || pdg > IONBASEHIGH) {
       // not an ion, return immediately
       return MassImpl(db->GetParticle(pdg), success);
@@ -99,6 +126,13 @@ class O2DatabasePDG
   O2DatabasePDG() = delete;
 
  private:
+  // negative: not set (the TDatabasePDG content is used)
+  static std::atomic<double>& monopoleMassOverride()
+  {
+    static std::atomic<double> mass{-1.};
+    return mass;
+  }
+
   static constexpr int IONBASELOW{1000000000};
   static constexpr int IONBASEHIGH{1099999999};
 };
@@ -161,7 +195,7 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa
   db->AddParticle("A_0+    ", " ", 0.9960, kFALSE, 0.0, +3, "meson", 9000211);
   db->AddParticle("A_0-    ", " ", 0.9960, kFALSE, 0.0, -3, "meson", -9000211);
 
-  //db->AddParticle("F0P0    ", " ", 0.9960, kFALSE, 0.0, 0, "meson",  9010221);
+  // db->AddParticle("F0P0    ", " ", 0.9960, kFALSE, 0.0, 0, "meson",  9010221);
 
   db->AddParticle("KDL_2+  ", " ", 1.773, kFALSE, 0.0, +3, "meson", 10325);
   db->AddParticle("KDL_2-  ", " ", 1.773, kFALSE, 0.0, -3, "meson", -10325);
@@ -261,7 +295,7 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa
   db->AddParticle("FeedbackPhoton", "FeedbackPhoton", 0, kFALSE,
                   0, 0, "Special", kspe + 51);
 
-  //Hyper nuclei and exotica
+  // Hyper nuclei and exotica
   ionCode = 1010010030;
   if (!db->GetParticle(ionCode)) {
     db->AddParticle("HyperTriton", "HyperTriton", 2.991134, kFALSE,
@@ -274,25 +308,25 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa
                     2.5e-15, 3, "Ion", ionCode);
   }
 
-  //hyper hydrogen 4 ground state
+  // hyper hydrogen 4 ground state
   ionCode = 1010010040;
   if (!db->GetParticle(ionCode)) {
     db->AddParticle("Hyperhydrog4", "Hyperhydrog4", 3.922434, kFALSE,
                     2.5e-15, 3, "Ion", ionCode);
   }
-  //anti hyper hydrogen 4 ground state
+  // anti hyper hydrogen 4 ground state
   ionCode = -1010010040;
   if (!db->GetParticle(ionCode)) {
     db->AddParticle("AntiHyperhydrog4", "AntiHyperhydrog4", 3.922434, kFALSE,
                     2.5e-15, 3, "Ion", ionCode);
   }
-  //hyper helium 4 ground state
+  // hyper helium 4 ground state
   ionCode = 1010020040;
   if (!db->GetParticle(ionCode)) {
     db->AddParticle("Hyperhelium4", "Hyperhelium4", 3.921728, kFALSE,
                     2.5e-15, 6, "Ion", ionCode);
   }
-  //anti hyper helium 4 ground state
+  // anti hyper helium 4 ground state
   ionCode = -1010020040;
   if (!db->GetParticle(ionCode)) {
     db->AddParticle("AntiHyperhelium4", "AntiHyperhelium4", 3.921728, kFALSE,
@@ -642,9 +676,9 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa
                     0.040, -3, "Resonance", -ionCode);
   }
 
-  //Additional Hidden Strangeness Pentaquarks
+  // Additional Hidden Strangeness Pentaquarks
 
-  //Ps +/-
+  // Ps +/-
   ionCode = 9322132;
   if (!db->GetParticle(ionCode)) {
     db->AddParticle("Ps_Plus_1870", "Ps_Plus_1870", 1.870, kFALSE,
@@ -682,7 +716,7 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa
                     0.10, -3, "Resonance", -ionCode);
   }
 
-  //Ps0
+  // Ps0
   ionCode = 9322131;
   if (!db->GetParticle(ionCode)) {
     db->AddParticle("Ps_0_1870", "Ps_0_1870", 1.870, kFALSE,
@@ -786,14 +820,14 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa
 
   // d*(2380) - dibaryon resonance
 
-  ionCode = 900010020; //Arbitrary choice - as deuteron but with leading 9 instead of 10
+  ionCode = 900010020; // Arbitrary choice - as deuteron but with leading 9 instead of 10
   if (!db->GetParticle(ionCode)) {
     db->AddParticle("d*_2380", "d*_2380", 2.38, kFALSE,
                     0.070, 3, "Resonance", ionCode);
   }
   db->AddAntiParticle("d*_2380_bar", -ionCode);
 
-  //Sexaquark (uuddss): compact, neutral and stable hypothetical bound state (arxiv.org/abs/1708.08951)
+  // Sexaquark (uuddss): compact, neutral and stable hypothetical bound state (arxiv.org/abs/1708.08951)
   ionCode = 900000020;
   if (!db->GetParticle(ionCode)) {
     db->AddParticle("Sexaquark", "Sexaquark", 2.0, kTRUE, 0.0, 0, "Special", ionCode);
@@ -826,4 +860,4 @@ inline void O2DatabasePDG::addParticlesFromExternalFile(TDatabasePDG* db)
 
 } // namespace o2
 
-#endif //O2_O2DATABASEPDG_H
+#endif // O2_O2DATABASEPDG_H

Steer/src/MCKinematicsReader.cxx

@@ -13,6 +13,8 @@
 #include "Steer/MCKinematicsReader.h"
 #include "SimulationDataFormat/MCEventHeader.h"
 #include "SimulationDataFormat/TrackReference.h"
+#include "SimulationDataFormat/MonopoleParticles.h"
+#include "SimulationDataFormat/O2DatabasePDG.h"
 #include <TChain.h>
 #include <stdexcept>
 #include <string>
@@ -82,6 +84,29 @@ void MCKinematicsReader::initTracksForSource(int source) const
   }
 }
 
+namespace
+{
+// Makes the monopole mass used by the simulation, which is stored in the event header, the
+// one of O2DatabasePDG; otherwise MCTrack::GetMass() (and GetEnergy(), GetRapidity(), ...) would
+// give the default mass for monopoles of a run with another G4.monopoleMass.
+void applyMonopoleMass(o2::dataformats::MCEventHeader const& header)
+{
+  if (!header.hasInfo(o2::sim::MonopoleMassHeaderKey)) {
+    return;
+  }
+  bool valid = false;
+  const double mass = header.getInfo<double>(o2::sim::MonopoleMassHeaderKey, valid);
+  if (!valid || !(mass > 0.)) {
+    return;
+  }
+  if (o2::O2DatabasePDG::hasMonopoleMassOverride() && std::abs(o2::O2DatabasePDG::getMonopoleMass() - mass) > 1.e-6 * mass) {
+    LOG(warn) << "Kinematics with different monopole masses (" << o2::O2DatabasePDG::getMonopoleMass() << " and " << mass
+              << " GeV) are read in the same process; monopole masses will be wrong for one of them";
+  }
+  o2::O2DatabasePDG::setMonopoleMass(mass);
+}
+} // namespace
+
 void MCKinematicsReader::loadTracksForSourceAndEvent(int source, int event) const
 {
   auto chain = mInputChains[source];
@@ -95,6 +120,17 @@ void MCKinematicsReader::loadTracksForSourceAndEvent(int source, int event) cons
       mTracks[source][event] = new std::vector<o2::MCTrack>;
       *mTracks[source][event] = *loadtracks;
       delete loadtracks;
+      // the monopole mass has to be known before the tracks are used
+      if (auto hbr = chain->GetBranch("MCEventHeader.")) {
+        o2::dataformats::MCEventHeader* header = nullptr;
+        hbr->SetAddress(&header);
+        hbr->GetEntry(event);
+        if (header) {
+          applyMonopoleMass(*header);
+        }
+        delete header;
+        hbr->ResetAddress();
+      }
     }
   }
 }
@@ -120,6 +156,7 @@ void MCKinematicsReader::loadHeadersForSource(int source) const
       for (int event = 0; event < br->GetEntries(); ++event) {
         br->GetEntry(event);
         mHeaders[source][event] = *header;
+        applyMonopoleMass(mHeaders[source][event]);
       }
       delete header;
       header = nullptr;

run/O2HitMerger.cxx

@@ -41,6 +41,9 @@
 
 #include "PrimaryServerState.h"
 #include <SimConfig/SimConfig.h>
+#include <SimConfig/G4Params.h>
+#include <CommonUtils/ConfigurableParam.h>
+#include <SimulationDataFormat/MonopoleParticles.h>
 #include <TPCSimulation/Detector.h>
 #include <ITSSimulation/Detector.h>
 #include <MFTSimulation/Detector.h>
@@ -147,6 +150,23 @@ void O2HitMerger::InitTask()
     // we didn't manage to get a configuration --> better to fail
     LOG(fatal) << "No configuration received. Aborting";
   }
+  // The config key-values are not applied to the ConfigurableParams in this process. The merger needs the G4 ones
+  // to record the monopole mass in the event headers (G4.monopole, G4.monopoleMass), so apply those.
+  // Note: G4 parameters given via a config file (--configFile) are not applied here.
+  {
+    std::string g4Tokens;
+    std::stringstream tokens(o2::conf::SimConfig::Instance().getKeyValueString());
+    std::string token;
+    while (std::getline(tokens, token, ';')) {
+      const auto first = token.find_first_not_of(" \t");
+      if (first != std::string::npos && token.compare(first, 3, "G4.") == 0) {
+        g4Tokens += token.substr(first) + ";";
+      }
+    }
+    if (!g4Tokens.empty()) {
+      o2::conf::ConfigurableParam::updateFromString(g4Tokens);
+    }
+  }
   mAsService = o2::conf::SimConfig::Instance().asService();
   mForwardKine = o2::conf::SimConfig::Instance().forwardKine();
   mWriteToDisc = o2::conf::SimConfig::Instance().writeToDisc();
@@ -694,6 +714,12 @@ bool O2HitMerger::mergeAndFlushData()
       }
     }
 
+    // record the mass used for the monopoles, such that readers of the kinematics
+    // (MCTrack::GetMass() and derived quantities) do not fall back to the default
+    if (eventheader && o2::conf::G4Params::Instance().monopole) {
+      eventheader->putInfo<double>(o2::sim::MonopoleMassHeaderKey, o2::conf::G4Params::Instance().monopoleMass);
+    }
+
     // now see which events can be discarded in any case due to no hits
     if (confref.isFilterOutNoHitEvents()) {
       if (eventheader && eventheader->getMCEventStats().getNHits() == 0) {