NuclearInteractions/src/RootNuclearInteractionParametersWriter.cpp

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