Root/src/RootNuclearInteractionParametersWriter.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Io/Root/RootNuclearInteractionParametersWriter.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Common.hpp"
0013 #include "ActsExamples/EventData/SimParticle.hpp"
0014 #include "ActsFatras/EventData/Particle.hpp"
0015
0016 #include <algorithm>
0017 #include <cstdint>
0018 #include <iterator>
0019 #include <limits>
0020 #include <memory>
0021 #include <stdexcept>
0022 #include <tuple>
0023 #include <unordered_map>
0024 #include <utility>
0025
0026 #include <TAxis.h>
0027 #include <TDirectory.h>
0028 #include <TFile.h>
0029 #include <TH1.h>
0030 #include <TVectorT.h>
0031
0032 namespace ActsExamples {
0033 struct AlgorithmContext;
0034 }  // namespace ActsExamples
0035
0036 namespace {
0037
0038 /// @brief This method labels events as either soft or hard and sets the
0039 /// multiplicity
0040 ///
0041 /// @param [in] events The events that will be labeled
0042 void labelEvents(
0043     std::vector<
0044         ActsExamples::detail::NuclearInteractionParametrisation::EventFraction>&
0045         events) {
0046   namespace Parametrisation =
0047       ActsExamples::detail::NuclearInteractionParametrisation;
0048   // Search for the highest momentum particles per event
0049   for (Parametrisation::EventFraction& event : events) {
0050     double maxMom = 0.;
0051     double maxMomOthers = 0.;
0052     // Walk over all final state particles
0053     for (const ActsExamples::SimParticle& p : event.finalParticles) {
0054       // Search for the maximum in particles with the same PDG ID as the
0055       // interacting one
0056       if (p.pdg() == event.initialParticle.pdg()) {
0057         maxMom = std::max(p.absoluteMomentum(), maxMom);
0058         // And the maximum among the others
0059       } else {
0060         maxMomOthers = std::max(p.absoluteMomentum(), maxMomOthers);
0061       }
0062     }
0063
0064     // Label the initial momentum
0065     event.initialMomentum = event.initialParticle.absoluteMomentum();
0066
0067     // Label the indication that the interacting particle carries most of the
0068     // momentum
0069     event.soft = (maxMom > maxMomOthers);
0070
0071     // Get the final state p_T
0072     double pt = 0.;
0073     Acts::Vector2 ptVec(0., 0.);
0074     for (const ActsExamples::SimParticle& p : event.finalParticles) {
0075       Acts::Vector2 particlePt =
0076           p.fourMomentum().template segment<2>(Acts::eMom0);
0077       ptVec[0] += particlePt[0];
0078       ptVec[1] += particlePt[1];
0079     }
0080     pt = ptVec.norm();
0081
0082     // Use the final state p_T as veto for the soft label
0083     if (event.soft && pt <= event.interactingParticle.transverseMomentum()) {
0084       event.soft = false;
0085     }
0086
0087     // Store the multiplicity
0088     event.multiplicity = event.finalParticles.size();
0089   }
0090 }
0091
0092 /// @brief This method builds components for transforming a probability
0093 /// distribution into components that represent the cumulative probability
0094 /// distribution
0095 ///
0096 /// @param [in] hist The probability distribution
0097 ///
0098 /// @return Tuple containing the bin borders, the non-normalised cumulative
0099 /// distribution and the sum over all entries
0100 std::tuple<std::vector<float>, std::vector<double>, double>
0101 buildNotNormalisedMap(TH1F const* hist) {
0102   // Retrieve the number of bins & borders
0103   const int nBins = hist->GetNbinsX();
0104   std::vector<float> histoBorders(nBins + 1);
0105
0106   // Fill the cumulative histogram
0107   double integral = 0.;
0108   std::vector<double> temp_HistoContents(nBins);
0109   for (int iBin = 0; iBin < nBins; iBin++) {
0110     float binval = hist->GetBinContent(iBin + 1);
0111     // Avoid negative bin values
0112     if (binval < 0) {
0113       binval = 0.;
0114     }
0115     // Store the value
0116     integral += binval;
0117     temp_HistoContents[iBin] = integral;
0118   }
0119
0120   // Ensure that content is available
0121   if (integral == 0.) {
0122     histoBorders.clear();
0123     temp_HistoContents.clear();
0124     return {histoBorders, temp_HistoContents, integral};
0125   }
0126
0127   // Set the bin borders
0128   for (int iBin = 1; iBin <= nBins; iBin++) {
0129     histoBorders[iBin - 1] = hist->GetXaxis()->GetBinLowEdge(iBin);
0130   }
0131   histoBorders[nBins] = hist->GetXaxis()->GetXmax();
0132
0133   return {histoBorders, temp_HistoContents, integral};
0134 }
0135
0136 /// @brief This function combines neighbouring bins with the same value
0137 ///
0138 /// @param [in, out] histoBorders The borders of the bins
0139 /// @param [in, out] histoContents The content of each bin
0140 void reduceMap(std::vector<float>& histoBorders,
0141                std::vector<std::uint32_t>& histoContents) {
0142   for (auto cit = histoContents.cbegin(); cit != histoContents.cend(); cit++) {
0143     while (std::next(cit, 1) != histoContents.end() &&
0144            *cit == *std::next(cit, 1)) {
0145       const auto distance =
0146           std::distance(histoContents.cbegin(), std::next(cit, 1));
0147       // Remove the bin
0148       histoBorders.erase(histoBorders.begin() + distance);
0149       histoContents.erase(histoContents.begin() + distance);
0150     }
0151   }
0152 }
0153
0154 /// @brief This method transforms a probability distribution into components
0155 /// that represent the cumulative probability distribution
0156 ///
0157 /// @note This method is used to avoid Root dependencies afterwards by
0158 /// decomposing the histogram
0159 /// @param [in] hist The probability distribution
0160 ///
0161 /// @return Pair containing the bin borders and the bin content
0162 std::pair<std::vector<float>, std::vector<std::uint32_t>> buildMap(
0163     TH1F const* hist) {
0164   // Build the components
0165   std::tuple<std::vector<float>, std::vector<double>, double> map =
0166       buildNotNormalisedMap(hist);
0167   const int nBins = hist->GetNbinsX();
0168   std::vector<double>& histoContents = std::get<1>(map);
0169
0170   // Fast exit if the histogram is empty
0171   if (histoContents.empty()) {
0172     return {std::get<0>(map), std::vector<std::uint32_t>()};
0173   }
0174
0175   // Set the bin content
0176   std::vector<std::uint32_t> normalisedHistoContents(nBins);
0177   const double invIntegral = 1. / std::get<2>(map);
0178   for (int iBin = 0; iBin < nBins; ++iBin) {
0179     normalisedHistoContents[iBin] =
0180         static_cast<unsigned int>(std::numeric_limits<std::uint32_t>::max() *
0181                                   (histoContents[iBin] * invIntegral));
0182   }
0183
0184   auto histoBorders = std::get<0>(map);
0185   reduceMap(histoBorders, normalisedHistoContents);
0186   return {histoBorders, normalisedHistoContents};
0187 }
0188
0189 /// @brief This method transforms a probability distribution into components
0190 /// that represent the cumulative probability distribution
0191 ///
0192 /// @note This method is used to avoid Root dependencies afterwards by
0193 /// decomposing the histogram
0194 /// @param [in] hist The probability distribution
0195 /// @param [in] integral Scaling factor of the distribution
0196 /// @note If @p integral is less than the actual integral of the histogram then
0197 /// the latter is used.
0198 ///
0199 /// @return Pair containing the bin borders and the bin content
0200 std::pair<std::vector<float>, std::vector<std::uint32_t>> buildMap(
0201     TH1F const* hist, double integral) {
0202   // Build the components
0203   std::tuple<std::vector<float>, std::vector<double>, double> map =
0204       buildNotNormalisedMap(hist);
0205
0206   const int nBins = hist->GetNbinsX();
0207   std::vector<double>& histoContents = std::get<1>(map);
0208
0209   // Fast exit if the histogram is empty
0210   if (histoContents.empty()) {
0211     return {std::get<0>(map), std::vector<std::uint32_t>()};
0212   }
0213
0214   // Set the bin content
0215   std::vector<std::uint32_t> normalisedHistoContents(nBins);
0216   const double invIntegral = 1. / std::max(integral, std::get<2>(map));
0217   for (int iBin = 0; iBin < nBins; ++iBin) {
0218     normalisedHistoContents[iBin] =
0219         static_cast<unsigned int>(std::numeric_limits<std::uint32_t>::max() *
0220                                   (histoContents[iBin] * invIntegral));
0221   }
0222
0223   std::vector<float> histoBorders = std::get<0>(map);
0224   reduceMap(histoBorders, normalisedHistoContents);
0225
0226   return {histoBorders, normalisedHistoContents};
0227 }
0228
0229 /// @brief This method builds decomposed cumulative probability distributions
0230 /// out of a vector of probability distributions
0231 ///
0232 /// @param [in] histos Vector of probability distributions
0233 ///
0234 /// @return Vector containing the decomposed cumulative probability
0235 /// distributions
0236 std::vector<std::pair<std::vector<float>, std::vector<std::uint32_t>>>
0237 buildMaps(const std::vector<TH1F*>& histos) {
0238   std::vector<std::pair<std::vector<float>, std::vector<std::uint32_t>>> maps;
0239   for (auto& h : histos) {
0240     maps.push_back(buildMap(h));
0241   }
0242   return maps;
0243 }
0244
0245 /// @brief This method parametrises and stores recursively a parametrisation of
0246 /// the final state kinematics in a nuclear interaction
0247 ///
0248 /// @param [in] eventFractionCollection The event storage
0249 /// @param [in] interactionType The interaction type that will be parametrised
0250 /// @param [in] cfg Configuration that steers the binning of histograms
0251 inline void recordKinematicParametrisation(
0252     const std::vector<
0253         ActsExamples::detail::NuclearInteractionParametrisation::EventFraction>&
0254         eventFractionCollection,
0255     bool interactionType, unsigned int multiplicity,
0256     const ActsExamples::RootNuclearInteractionParametersWriter::Config& cfg) {
0257   namespace Parametrisation =
0258       ActsExamples::detail::NuclearInteractionParametrisation;
0259   gDirectory->mkdir(std::to_string(multiplicity).c_str());
0260   gDirectory->cd(std::to_string(multiplicity).c_str());
0261
0262   // Parametrise the momentum und invarian mass distributions
0263   const auto momentumParameters = Parametrisation::buildMomentumParameters(
0264       eventFractionCollection, multiplicity, interactionType, cfg.momentumBins);
0265   std::vector<Parametrisation::CumulativeDistribution> distributionsMom =
0266       momentumParameters.second;
0267
0268   const auto invariantMassParameters =
0269       Parametrisation::buildInvariantMassParameters(
0270           eventFractionCollection, multiplicity, interactionType,
0271           cfg.invariantMassBins);
0272   std::vector<Parametrisation::CumulativeDistribution> distributionsInvMass =
0273       invariantMassParameters.second;
0274
0275   // Fast exit in case of no events
0276   if (!distributionsMom.empty() && !distributionsInvMass.empty()) {
0277     if (multiplicity > 1) {
0278       // Write the eigenspace components for the momenta
0279       Parametrisation::EigenspaceComponents esComponentsMom =
0280           momentumParameters.first;
0281
0282       auto momEigenVal = std::get<0>(esComponentsMom);
0283       auto momEigenVec = std::get<1>(esComponentsMom);
0284       auto momMean = std::get<2>(esComponentsMom);
0285       std::vector<float> momVecVal(momEigenVal.data(),
0286                                    momEigenVal.data() + momEigenVal.size());
0287       std::vector<float> momVecVec(momEigenVec.data(),
0288                                    momEigenVec.data() + momEigenVec.size());
0289       std::vector<float> momVecMean(momMean.data(),
0290                                     momMean.data() + momMean.size());
0291       gDirectory->WriteObject(&momVecVal, "MomentumEigenvalues");
0292       gDirectory->WriteObject(&momVecVec, "MomentumEigenvectors");
0293       gDirectory->WriteObject(&momVecMean, "MomentumMean");
0294
0295       // Write the eigenspace components for the invariant masses
0296       Parametrisation::EigenspaceComponents esComponentsInvMass =
0297           invariantMassParameters.first;
0298
0299       auto invMassEigenVal = std::get<0>(esComponentsInvMass);
0300       auto invMassEigenVec = std::get<1>(esComponentsInvMass);
0301       auto invMassMean = std::get<2>(esComponentsInvMass);
0302       std::vector<float> invMassVecVal(
0303           invMassEigenVal.data(),
0304           invMassEigenVal.data() + invMassEigenVal.size());
0305       std::vector<float> invMassVecVec(
0306           invMassEigenVec.data(),
0307           invMassEigenVec.data() + invMassEigenVec.size());
0308       std::vector<float> invMassVecMean(
0309           invMassMean.data(), invMassMean.data() + invMassMean.size());
0310       gDirectory->WriteObject(&invMassVecVal, "InvariantMassEigenvalues");
0311       gDirectory->WriteObject(&invMassVecVec, "InvariantMassEigenvectors");
0312       gDirectory->WriteObject(&invMassVecMean, "InvariantMassMean");
0313     }
0314
0315     const auto momDistributions = buildMaps(distributionsMom);
0316     const auto invMassDistributions = buildMaps(distributionsInvMass);
0317
0318     // Write the distributions
0319     for (unsigned int i = 0; i <= multiplicity; i++) {
0320       if (cfg.writeOptionalHistograms) {
0321         gDirectory->WriteObject(
0322             distributionsMom[i],
0323             ("MomentumDistributionHistogram_" + std::to_string(i)).c_str());
0324       }
0325       gDirectory->WriteObject(
0326           &momDistributions[i].first,
0327           ("MomentumDistributionBinBorders_" + std::to_string(i)).c_str());
0328       gDirectory->WriteObject(
0329           &momDistributions[i].second,
0330           ("MomentumDistributionBinContents_" + std::to_string(i)).c_str());
0331     }
0332     for (unsigned int i = 0; i < multiplicity; i++) {
0333       if (cfg.writeOptionalHistograms) {
0334         gDirectory->WriteObject(
0335             distributionsInvMass[i],
0336             ("InvariantMassDistributionHistogram_" + std::to_string(i))
0337                 .c_str());
0338       }
0339       gDirectory->WriteObject(
0340           &invMassDistributions[i].first,
0341           ("InvariantMassDistributionBinBorders_" + std::to_string(i)).c_str());
0342       gDirectory->WriteObject(
0343           &invMassDistributions[i].second,
0344           ("InvariantMassDistributionBinContents_" + std::to_string(i))
0345               .c_str());
0346     }
0347   }
0348   gDirectory->cd("..");
0349 }
0350 }  // namespace
0351
0352 ActsExamples::RootNuclearInteractionParametersWriter::
0353     RootNuclearInteractionParametersWriter(
0354         const ActsExamples::RootNuclearInteractionParametersWriter::Config&
0355             config,
0356         Acts::Logging::Level level)
0357     : WriterT(config.inputSimulationProcesses,
0358               "RootNuclearInteractionParametersWriter", level),
0359       m_cfg(config) {
0360   if (m_cfg.inputSimulationProcesses.empty()) {
0361     throw std::invalid_argument("Missing input collection");
0362   }
0363   if (m_cfg.filePath.empty()) {
0364     throw std::invalid_argument("Missing output filename for parameters");
0365   }
0366 }
0367
0368 ActsExamples::RootNuclearInteractionParametersWriter::
0369     ~RootNuclearInteractionParametersWriter() = default;
0370
0371 ActsExamples::ProcessCode
0372 ActsExamples::RootNuclearInteractionParametersWriter::finalize() {
0373   namespace Parametrisation = detail::NuclearInteractionParametrisation;
0374   if (m_eventFractionCollection.empty()) {
0375     return ProcessCode::ABORT;
0376   }
0377
0378   // Exclusive access to the file while writing
0379   std::lock_guard<std::mutex> lock(m_writeMutex);
0380
0381   // The file
0382   TFile* tf = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0383   gDirectory->cd();
0384   gDirectory->mkdir(
0385       std::to_string(m_eventFractionCollection[0].initialParticle.pdg())
0386           .c_str());
0387   gDirectory->cd(
0388       std::to_string(m_eventFractionCollection[0].initialParticle.pdg())
0389           .c_str());
0390   gDirectory->mkdir(
0391       std::to_string(m_eventFractionCollection[0].initialMomentum).c_str());
0392   gDirectory->cd(
0393       std::to_string(m_eventFractionCollection[0].initialMomentum).c_str());
0394   gDirectory->mkdir("soft");
0395   gDirectory->mkdir("hard");
0396
0397   // Write the nuclear interaction probability
0398   ACTS_DEBUG("Starting parametrisation of nuclear interaction probability");
0399   const auto nuclearInteractionProbability =
0400       Parametrisation::cumulativeNuclearInteractionProbability(
0401           m_eventFractionCollection, m_cfg.interactionProbabilityBins);
0402
0403   if (m_cfg.writeOptionalHistograms) {
0404     gDirectory->WriteObject(nuclearInteractionProbability,
0405                             "NuclearInteractionHistogram");
0406   }
0407   const auto mapNIprob =
0408       buildMap(nuclearInteractionProbability, m_cfg.nSimulatedEvents);
0409   gDirectory->WriteObject(&mapNIprob.first, "NuclearInteractionBinBorders");
0410   gDirectory->WriteObject(&mapNIprob.second, "NuclearInteractionBinContents");
0411   ACTS_DEBUG("Nuclear interaction probability parametrised");
0412
0413   ACTS_DEBUG("Starting calculation of probability of interaction type");
0414   // Write the interaction type probability
0415   const auto softProbability =
0416       Parametrisation::softProbability(m_eventFractionCollection);
0417
0418   gDirectory->WriteObject(&softProbability, "SoftInteraction");
0419   ACTS_DEBUG("Calculation of probability of interaction type finished");
0420
0421   // Write the PDG id production distribution
0422   ACTS_DEBUG(
0423       "Starting calculation of transition probabilities between PDG IDs");
0424   const auto pdgIdMap =
0425       Parametrisation::cumulativePDGprobability(m_eventFractionCollection);
0426   std::vector<int> branchingPdgIds;
0427   std::vector<int> targetPdgIds;
0428   std::vector<float> targetPdgProbability;
0429   for (const auto& [targetKey, targetValue] : pdgIdMap) {
0430     for (const auto& [producedKey, producedValue] : targetValue) {
0431       branchingPdgIds.push_back(targetKey);
0432       targetPdgIds.push_back(producedKey);
0433       targetPdgProbability.push_back(producedValue);
0434     }
0435   }
0436
0437   gDirectory->WriteObject(&branchingPdgIds, "BranchingPdgIds");
0438   gDirectory->WriteObject(&targetPdgIds, "TargetPdgIds");
0439   gDirectory->WriteObject(&targetPdgProbability, "TargetPdgProbability");
0440   ACTS_DEBUG(
0441       "Calculation of transition probabilities between PDG IDs finished");
0442
0443   // Write the multiplicity and kinematics distribution
0444   ACTS_DEBUG("Starting parametrisation of multiplicity probabilities");
0445   const auto multiplicity = Parametrisation::cumulativeMultiplicityProbability(
0446       m_eventFractionCollection, m_cfg.multiplicityMax);
0447   ACTS_DEBUG("Parametrisation of multiplicity probabilities finished");
0448
0449   gDirectory->cd("soft");
0450   if (m_cfg.writeOptionalHistograms) {
0451     gDirectory->WriteObject(multiplicity.first, "MultiplicityHistogram");
0452   }
0453   const auto multProbSoft = buildMap(multiplicity.first);
0454   gDirectory->WriteObject(&multProbSoft.first, "MultiplicityBinBorders");
0455   gDirectory->WriteObject(&multProbSoft.second, "MultiplicityBinContents");
0456   for (unsigned int i = 1; i <= m_cfg.multiplicityMax; i++) {
0457     ACTS_DEBUG("Starting parametrisation of final state kinematics for soft " +
0458                std::to_string(i) + " particle(s) final state");
0459     recordKinematicParametrisation(m_eventFractionCollection, true, i, m_cfg);
0460     ACTS_DEBUG("Parametrisation of final state kinematics for soft " +
0461                std::to_string(i) + " particle(s) final state finished");
0462   }
0463   gDirectory->cd("../hard");
0464   if (m_cfg.writeOptionalHistograms) {
0465     gDirectory->WriteObject(multiplicity.second, "MultiplicityHistogram");
0466   }
0467   const auto multProbHard = buildMap(multiplicity.second);
0468   gDirectory->WriteObject(&multProbHard.first, "MultiplicityBinBorders");
0469   gDirectory->WriteObject(&multProbHard.second, "MultiplicityBinContents");
0470
0471   for (unsigned int i = 1; i <= m_cfg.multiplicityMax; i++) {
0472     ACTS_DEBUG("Starting parametrisation of final state kinematics for hard " +
0473                std::to_string(i) + " particle(s) final state");
0474     recordKinematicParametrisation(m_eventFractionCollection, false, i, m_cfg);
0475     ACTS_DEBUG("Parametrisation of final state kinematics for hard " +
0476                std::to_string(i) + " particle(s) final state finished");
0477   }
0478   gDirectory->cd();
0479   tf->Write();
0480   tf->Close();
0481   return ProcessCode::SUCCESS;
0482 }
0483
0484 ActsExamples::ProcessCode
0485 ActsExamples::RootNuclearInteractionParametersWriter::writeT(
0486     const AlgorithmContext& /*ctx*/,
0487     const ExtractedSimulationProcessContainer& event) {
0488   // Convert the tuple to use additional categorisation variables
0489   std::vector<detail::NuclearInteractionParametrisation::EventFraction>
0490       eventFractions;
0491   eventFractions.reserve(event.size());
0492   for (const auto& e : event) {
0493     eventFractions.emplace_back(e);
0494   }
0495   // Categorise the event
0496   labelEvents(eventFractions);
0497   // Store the event
0498   m_eventFractionCollection.insert(m_eventFractionCollection.end(),
0499                                    eventFractions.begin(),
0500                                    eventFractions.end());
0501
0502   return ProcessCode::SUCCESS;
0503 }