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