feat: Reading of nuclear interaction parametrisation (#695)

This PR adds a reader to the files written by #574.
Note: Since the reader, the internal storage and the application and closely linked to each other, all these files are part of this PR and should be discussed altogether. However, the Fatras process (NuclearInteraction.hpp/cpp) can be split into another PR.

Examples/Algorithms/Fatras/include/ActsExamples/Fatras/FatrasAlgorithm.hpp

@@ -14,6 +14,7 @@
 #include "ActsExamples/Framework/RandomNumbers.hpp"
 #include "ActsExamples/MagneticField/MagneticField.hpp"
 #include "ActsExamples/Utilities/OptionsFwd.hpp"
+#include "ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp"
 
 #include <memory>
 #include <string>
@@ -35,6 +36,9 @@ class FatrasAlgorithm final : public BareAlgorithm {
     std::string outputParticlesFinal;
     /// The simulated hits output collection.
     std::string outputSimHits;
+    /// Parametrisation of nuclear interaction
+    std::string imputParametrisationNuclearInteraction =
+        "nuclearInteractionParameters";
     /// Random number service.
     std::shared_ptr<const RandomNumbers> randomNumbers;
     /// The tracking geometry that should be used.

Examples/Io/NuclearInteractions/CMakeLists.txt

@@ -1,5 +1,6 @@
 add_library(
   ActsExamplesIoNuclearInteractions SHARED
+  src/NuclearInteractionOptions.cpp
   src/RootNuclearInteractionParametersWriter.cpp
   src/detail/NuclearInteractionParametrisation.cpp
 )
@@ -11,7 +12,7 @@ target_link_libraries(
   PUBLIC
     ActsCore ActsExamplesFramework
     Threads::Threads
-  PRIVATE ROOT::Core ROOT::Hist ROOT::Tree)
+  PRIVATE ROOT::Core ROOT::Hist ROOT::Tree Boost::program_options)
 
 install(
   TARGETS ActsExamplesIoNuclearInteractions

Examples/Io/NuclearInteractions/include/ActsExamples/Io/NuclearInteractions/NuclearInteractionOptions.hpp

@@ -0,0 +1,67 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2021 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/Utilities/Logger.hpp"
+#include "ActsExamples/Utilities/OptionsFwd.hpp"
+#include "ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp"
+#include "ActsFatras/Physics/NuclearInteraction/NuclearInteractionParameters.hpp"
+
+#include <utility>
+
+#include "boost/program_options.hpp"
+
+namespace ActsExamples {
+namespace Options {
+
+/// Add Fatras options.
+///
+/// @param desc The options description to add options to
+void addNuclearInteractionOptions(Description& desc);
+
+/// Reads the parametrisation and provides the parametrisation
+ActsFatras::detail::MultiParticleNuclearInteractionParametrisation
+readParametrisations(const std::string& fileName);
+
+/// Read Fatras options to create the algorithm config.
+///
+/// @param variables the variables to read from
+std::string readNuclearInteractionConfig(
+    const boost::program_options::variables_map& variables);
+
+/// Read the parametrisations.
+///
+/// @tparam simulator_t type of the simulation kernel
+/// @param [in] nuclearInteractionParametrisation File name of the
+/// parametrisations
+/// @param [in, out] simulator The simulation kernel
+template <typename simulator_t>
+void setNuclearInteractionParametrisations(
+    const std::string& nuclearInteractionParametrisation,
+    simulator_t& simulator) {
+  if (nuclearInteractionParametrisation.empty()) {
+    return;
+  }
+
+  auto& chargedNuclearInteraction =
+      simulator.charged.interactions
+          .template get<ActsFatras::NuclearInteraction>();
+  auto& neutralNuclearInteraction =
+      simulator.neutral.interactions
+          .template get<ActsFatras::NuclearInteraction>();
+
+  const auto mpp = readParametrisations(nuclearInteractionParametrisation);
+
+  chargedNuclearInteraction.multiParticleParameterisation = mpp;
+  neutralNuclearInteraction.multiParticleParameterisation = mpp;
+}
+
+}  // namespace Options
+}  // namespace ActsExamples

Examples/Io/NuclearInteractions/src/NuclearInteractionOptions.cpp

@@ -0,0 +1,234 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2021 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "ActsExamples/Io/NuclearInteractions/NuclearInteractionOptions.hpp"
+
+#include "ActsExamples/Utilities/Options.hpp"
+
+#include <string>
+
+#include <TFile.h>
+#include <TH1F.h>
+#include <TVectorF.h>
+
+namespace {
+std::pair<std::vector<float>, std::vector<uint32_t>> readHistogram(
+    const std::string&& binBordersName, const std::string&& binContentsName) {
+  std::vector<float>* binBorders;
+  std::vector<uint32_t>* binContents;
+  // Get the decomposed histogram
+  binBorders = (std::vector<float>*)gDirectory->Get(binBordersName.c_str());
+  binContents =
+      (std::vector<uint32_t>*)gDirectory->Get(binContentsName.c_str());
+  // Return the histogram if available
+  if (binBorders != nullptr && binContents != nullptr) {
+    return std::make_pair(*binBorders, *binContents);
+  }
+  return std::make_pair(std::vector<float>(), std::vector<uint32_t>());
+}
+
+Acts::ActsDynamicVector readVector(const std::string&& vectorName) {
+  std::vector<float>* vector;
+  // Get the vector
+  vector = (std::vector<float>*)gDirectory->Get(vectorName.c_str());
+  // Return the vector if available
+  if (vector != nullptr) {
+    Acts::ActsDynamicVector result;
+    const unsigned int sizeVec = vector->size();
+    result.resize(sizeVec);
+    for (unsigned int i = 0; i < sizeVec; i++)
+      result(i) = (*vector)[i];
+
+    return result;
+  }
+  return Acts::ActsDynamicVector();
+}
+
+void readKinematicParameters(
+    ActsFatras::detail::NuclearInteractionParameters& parameters,
+    TObject* folder, bool softInteractionParameters) {
+  if (folder->IsFolder()) {
+    // Find the momentum and invariant mass distributions
+    const char* distributionName = folder->GetName();
+    unsigned int mult = std::stoi(distributionName);
+    gDirectory->cd(distributionName);
+    std::vector<std::pair<std::vector<float>, std::vector<uint32_t>>>
+        momentumDistributions;
+    momentumDistributions.resize(mult + 1);
+    std::vector<std::pair<std::vector<float>, std::vector<uint32_t>>>
+        invariantMassDistributions;
+    invariantMassDistributions.resize(mult);
+    for (unsigned int i = 0; i < mult; i++) {
+      momentumDistributions[i] = readHistogram(
+          ("MomentumDistributionBinBorders_" + std::to_string(i)).c_str(),
+          ("MomentumDistributionBinContents_" + std::to_string(i)).c_str());
+
+      invariantMassDistributions[i] = readHistogram(
+          ("InvariantMassDistributionBinBorders_" + std::to_string(i)).c_str(),
+          ("InvariantMassDistributionBinContents_" + std::to_string(i))
+              .c_str());
+    }
+    momentumDistributions.back() = readHistogram(
+        ("MomentumDistributionBinBorders_" + std::to_string(mult)).c_str(),
+        ("MomentumDistributionBinContents_" + std::to_string(mult)).c_str());
+
+    // Get the eigenspace components for the kinematic parameters
+    Acts::ActsDynamicVector momentumEigenvalues =
+        readVector("MomentumEigenvalues");
+    Acts::ActsDynamicVector momentumEigenvectors =
+        readVector("MomentumEigenvectors");
+    Acts::ActsDynamicVector momentumMean = readVector("MomentumMean");
+    Acts::ActsDynamicVector invariantMassEigenvalues =
+        readVector("InvariantMassEigenvalues");
+    Acts::ActsDynamicVector invariantMassEigenvectors =
+        readVector("InvariantMassEigenvectors");
+    Acts::ActsDynamicVector invariantMassMean = readVector("InvariantMassMean");
+
+    // Test that a parametrisation is present
+    if (momentumEigenvalues.size() != 0 && momentumEigenvectors.size() != 0 &&
+        momentumMean.size() != 0 && invariantMassEigenvalues.size() != 0 &&
+        invariantMassEigenvectors.size() != 0 &&
+        invariantMassMean.size() != 0) {
+      // Prepare and store the kinematic parameters
+      ActsFatras::detail::NuclearInteractionParameters::
+          ParametersWithFixedMultiplicity kinematicParameters(
+              momentumDistributions, momentumEigenvalues, momentumEigenvectors,
+              momentumMean, invariantMassDistributions,
+              invariantMassEigenvalues, invariantMassEigenvectors,
+              invariantMassMean);
+      if (softInteractionParameters) {
+        if (mult >= parameters.softKinematicParameters.size())
+          parameters.softKinematicParameters.resize(mult + 1);
+        parameters.softKinematicParameters[mult] = kinematicParameters;
+      } else {
+        if (mult >= parameters.hardKinematicParameters.size())
+          parameters.hardKinematicParameters.resize(mult + 1);
+        parameters.hardKinematicParameters[mult] = kinematicParameters;
+      }
+    }
+    gDirectory->cd("..");
+  }
+}
+}  // namespace
+
+void ActsExamples::Options::addNuclearInteractionOptions(
+    ActsExamples::Options::Description& desc) {
+  using boost::program_options::value;
+
+  auto opt = desc.add_options();
+  opt("fatras-nuclear-interaction-parametrisation",
+      value<std::string>()->default_value({}),
+      "File containing parametrisations for nuclear interaction.");
+}
+
+std::string ActsExamples::Options::readNuclearInteractionConfig(
+    const boost::program_options::variables_map& variables) {
+  return variables["fatras-nuclear-interaction-parametrisation"]
+      .as<std::string>();
+}
+
+ActsFatras::detail::MultiParticleNuclearInteractionParametrisation
+ActsExamples::Options::readParametrisations(const std::string& fileName) {
+  // The collection
+  ActsFatras::detail::MultiParticleNuclearInteractionParametrisation mpp;
+
+  // Now read file
+  ActsFatras::detail::NuclearInteractionParametrisation parametrisation;
+  TFile tf(fileName.c_str(), "read");
+  gDirectory->cd();
+  auto listOfParticles = gDirectory->GetListOfKeys();
+  auto initialParticle = listOfParticles->First();
+  while (initialParticle) {
+    // Get the initial particle
+    char const* particleName = initialParticle->GetName();
+    gDirectory->cd(particleName);
+
+    // Walk over all initial momenta for a particle
+    auto listOfMomenta = gDirectory->GetListOfKeys();
+    auto initialMomentum = listOfMomenta->First();
+    while (initialMomentum) {
+      // Parameters for a fixed inital momentum
+      ActsFatras::detail::NuclearInteractionParameters parameters;
+      // Get the initial momentum
+      char const* nameMomentum = initialMomentum->GetName();
+      parameters.momentum = std::stof(nameMomentum);
+      gDirectory->cd(nameMomentum);
+
+      // Get the nuclear interaction probability
+      parameters.nuclearInteractionProbability = readHistogram(
+          "NuclearInteractionBinBorders", "NuclearInteractionBinContents");
+
+      // Get the soft interaction probability
+      TVectorF* softInteraction = (TVectorF*)gDirectory->Get("SoftInteraction");
+      parameters.softInteractionProbability = (*softInteraction)[0];
+
+      // Get the branching probabilities
+      std::vector<int> branchingPdgIds =
+          *((std::vector<int>*)gDirectory->Get("BranchingPdgIds"));
+      std::vector<int> targetPdgIds =
+          *((std::vector<int>*)gDirectory->Get("TargetPdgIds"));
+      std::vector<float> targetPdgProbability =
+          *((std::vector<float>*)gDirectory->Get("TargetPdgProbability"));
+      parameters.pdgMap.reserve(branchingPdgIds.size());
+      for (unsigned int i = 0; i < branchingPdgIds.size(); i++) {
+        auto it = parameters.pdgMap.begin();
+        while (it->first != branchingPdgIds[i] && it != parameters.pdgMap.end())
+          it++;
+
+        const auto target =
+            std::make_pair(targetPdgIds[i], targetPdgProbability[i]);
+        if (it != parameters.pdgMap.end()) {
+          it->second.push_back(target);
+        } else {
+          parameters.pdgMap.push_back(std::make_pair(
+              branchingPdgIds[i], std::vector<std::pair<int, float>>{target}));
+        }
+      }
+
+      // Get the soft distributions
+      gDirectory->cd("soft");
+      // Get the multiplicity distribution
+      parameters.softMultiplicity =
+          readHistogram("MultiplicityBinBorders", "MultiplicityBinContents");
+      // Get the distributions for each final state multiplicity
+      auto softList = gDirectory->GetListOfKeys();
+      auto softElement = softList->First();
+      while (softElement) {
+        readKinematicParameters(parameters, softElement, true);
+        softElement = softList->After(softElement);
+      }
+
+      // Get the hard distributions
+      gDirectory->cd("../hard");
+      // Get the multiplicity distribution
+      parameters.hardMultiplicity =
+          readHistogram("MultiplicityBinBorders", "MultiplicityBinContents");
+
+      // Get the distributions for each final state multiplicity
+      auto hardList = gDirectory->GetListOfKeys();
+      auto hardElement = hardList->First();
+      while (hardElement) {
+        readKinematicParameters(parameters, hardElement, false);
+        hardElement = hardList->After(hardElement);
+      }
+
+      initialMomentum = listOfMomenta->After(initialMomentum);
+      // Store the parametrisation
+      parametrisation.push_back(
+          std::make_pair(parameters.momentum, parameters));
+    }
+    tf.Close();
+
+    // Write to the collection to the EventStore
+    mpp.push_back(std::make_pair(std::stof(particleName), parametrisation));
+
+    initialParticle = listOfParticles->After(initialParticle);
+  }
+  // Return success flag
+  return mpp;
+}

Fatras/CMakeLists.txt

@@ -8,6 +8,7 @@ add_library(
   src/EventData/ProcessType.cpp
   src/Kernel/SimulationError.cpp
   src/Physics/BetheHeitler.cpp
+  src/Physics/NuclearInteraction/NuclearInteraction.cpp
   src/Physics/PhotonConversion.cpp
   src/Physics/StandardInteractions.cpp
   src/Utilities/LandauDistribution.cpp

Fatras/include/ActsFatras/EventData/ProcessType.hpp

@@ -21,6 +21,7 @@ enum class ProcessType : uint32_t {
   eDecay = 1,
   ePhotonConversion = 2,
   eBremsstrahlung = 3,
+  eNuclearInteraction = 4,
 };
 
 std::ostream &operator<<(std::ostream &os, ProcessType processType);