Root/src/RootMeasurementWriter.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/RootMeasurementWriter.hpp"
0010
0011 #include "Acts/Definitions/TrackParametrization.hpp"
0012 #include "ActsExamples/EventData/AverageSimHits.hpp"
0013 #include "ActsExamples/EventData/Index.hpp"
0014 #include "ActsExamples/EventData/Measurement.hpp"
0015 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0016 #include "ActsExamples/Utilities/Range.hpp"
0017
0018 #include <ios>
0019 #include <limits>
0020 #include <memory>
0021 #include <stdexcept>
0022 #include <utility>
0023
0024 #include <TFile.h>
0025 #include <TTree.h>
0026
0027 namespace ActsExamples {
0028
0029 struct RootMeasurementWriter::DigitizationTree {
0030   const std::array<std::string, Acts::eBoundSize> bNames = {
0031       "loc0", "loc1", "phi", "theta", "qop", "time"};
0032
0033   TTree* tree = nullptr;
0034
0035   // Identification parameters
0036   int eventNr = 0;
0037   int volumeID = 0;
0038   int layerID = 0;
0039   int surfaceID = 0;
0040   int extraID = 0;
0041
0042   // Reconstruction information
0043   float recBound[Acts::eBoundSize] = {};
0044   float varBound[Acts::eBoundSize] = {};
0045
0046   // Truth parameters
0047   float trueBound[Acts::eBoundSize] = {};
0048   float trueGx = 0.;
0049   float trueGy = 0.;
0050   float trueGz = 0.;
0051   float incidentPhi = 0.;
0052   float incidentTheta = 0.;
0053
0054   // Residuals and pulls
0055   float residual[Acts::eBoundSize] = {};
0056   float pull[Acts::eBoundSize] = {};
0057
0058   // Cluster information comprised of
0059   // nch :  number of channels
0060   // cSize : cluster size in loc0 and loc1
0061   // chId : channel identification
0062   // chValue: value/activation of the channel
0063   int nch = 0;
0064   int cSize[2] = {};
0065   std::array<std::vector<int>, 2> chId;
0066   std::vector<float> chValue;
0067
0068   /// Constructor from tree name
0069   DigitizationTree(const std::string& treeName,
0070                    const std::vector<Acts::BoundIndices>& recoIndices,
0071                    const std::vector<Acts::BoundIndices>& clusterIndices) {
0072     tree = new TTree(treeName.c_str(), treeName.c_str());
0073
0074     tree->Branch("event_nr", &eventNr);
0075     tree->Branch("volume_id", &volumeID);
0076     tree->Branch("layer_id", &layerID);
0077     tree->Branch("surface_id", &surfaceID);
0078     tree->Branch("extra_id", &extraID);
0079
0080     for (auto ib : recoIndices) {
0081       tree->Branch(("rec_" + bNames[ib]).c_str(), &recBound[ib]);
0082     }
0083     for (auto ib : recoIndices) {
0084       tree->Branch(("var_" + bNames[ib]).c_str(), &varBound[ib]);
0085     }
0086
0087     tree->Branch("clus_size", &nch);
0088     tree->Branch("channel_value", &chValue);
0089     // Both are allocated, but only relevant ones are set
0090     for (auto ib : clusterIndices) {
0091       if (static_cast<unsigned int>(ib) < 2) {
0092         tree->Branch(("channel_" + bNames[ib]).c_str(), &chId[ib]);
0093         tree->Branch(("clus_size_" + bNames[ib]).c_str(), &cSize[ib]);
0094       }
0095     }
0096
0097     for (unsigned int ib = 0; ib < Acts::eBoundSize; ++ib) {
0098       tree->Branch(("true_" + bNames[ib]).c_str(), &trueBound[ib]);
0099     }
0100     tree->Branch("true_x", &trueGx);
0101     tree->Branch("true_y", &trueGy);
0102     tree->Branch("true_z", &trueGz);
0103     tree->Branch("true_incident_phi", &incidentPhi);
0104     tree->Branch("true_incident_theta", &incidentTheta);
0105
0106     for (auto ib : recoIndices) {
0107       tree->Branch(("residual_" + bNames[ib]).c_str(), &residual[ib]);
0108     }
0109     for (auto ib : recoIndices) {
0110       tree->Branch(("pull_" + bNames[ib]).c_str(), &pull[ib]);
0111     }
0112
0113     clear();
0114   }
0115
0116   /// Convenience function to register idenfication
0117   ///
0118   /// @param eventNr The event number
0119   /// @param geoID The geometry identifier of the measurement
0120   void fillIdentification(int evnt, Acts::GeometryIdentifier geoId) {
0121     eventNr = evnt;
0122     volumeID = geoId.volume();
0123     layerID = geoId.layer();
0124     surfaceID = geoId.sensitive();
0125     extraID = geoId.extra();
0126   }
0127
0128   /// Convenience function to register the truth parameters
0129   ///
0130   /// @param lp The true local position
0131   /// @param xt The true 4D global position
0132   /// @param dir The true particle direction
0133   /// @param angles The incident angles
0134   /// @param qop The true q/p
0135   void fillTruthParameters(const Acts::Vector2& lp, const Acts::Vector4& xt,
0136                            const Acts::Vector3& dir,
0137                            const std::pair<double, double> angles) {
0138     trueBound[Acts::eBoundLoc0] = lp[Acts::eBoundLoc0];
0139     trueBound[Acts::eBoundLoc1] = lp[Acts::eBoundLoc1];
0140     trueBound[Acts::eBoundPhi] = Acts::VectorHelpers::phi(dir);
0141     trueBound[Acts::eBoundTheta] = Acts::VectorHelpers::theta(dir);
0142     trueBound[Acts::eBoundTime] = xt[Acts::eTime];
0143
0144     trueGx = xt[Acts::ePos0];
0145     trueGy = xt[Acts::ePos1];
0146     trueGz = xt[Acts::ePos2];
0147
0148     incidentPhi = angles.first;
0149     incidentTheta = angles.second;
0150   }
0151
0152   /// Convenience function to fill bound parameters
0153   ///
0154   /// @param m The measurement
0155   void fillBoundMeasurement(const ConstVariableBoundMeasurementProxy& m) {
0156     for (unsigned int i = 0; i < m.size(); ++i) {
0157       auto ib = m.subspaceIndexVector()[i];
0158
0159       recBound[ib] = m.parameters()[i];
0160       varBound[ib] = m.covariance()(i, i);
0161
0162       residual[ib] = recBound[ib] - trueBound[ib];
0163       pull[ib] = residual[ib] / std::sqrt(varBound[ib]);
0164     }
0165   }
0166
0167   /// Convenience function to fill the cluster information
0168   ///
0169   /// @param c The cluster
0170   void fillCluster(const Cluster& c) {
0171     nch = static_cast<int>(c.channels.size());
0172     cSize[0] = static_cast<int>(c.sizeLoc0);
0173     cSize[1] = static_cast<int>(c.sizeLoc1);
0174     for (auto ch : c.channels) {
0175       chId[0].push_back(static_cast<int>(ch.bin[0]));
0176       chId[1].push_back(static_cast<int>(ch.bin[1]));
0177       chValue.push_back(static_cast<float>(ch.activation));
0178     }
0179   }
0180
0181   /// Fill the tree
0182   void fill() { tree->Fill(); }
0183
0184   /// Clear the tree
0185   void clear() {
0186     for (unsigned int ib = 0; ib < Acts::eBoundSize; ++ib) {
0187       trueBound[ib] = std::numeric_limits<float>::quiet_NaN();
0188       recBound[ib] = std::numeric_limits<float>::quiet_NaN();
0189       varBound[ib] = std::numeric_limits<float>::quiet_NaN();
0190       residual[ib] = std::numeric_limits<float>::quiet_NaN();
0191       pull[ib] = std::numeric_limits<float>::quiet_NaN();
0192     }
0193     trueGx = std::numeric_limits<float>::quiet_NaN();
0194     trueGy = std::numeric_limits<float>::quiet_NaN();
0195     trueGz = std::numeric_limits<float>::quiet_NaN();
0196     incidentPhi = std::numeric_limits<float>::quiet_NaN();
0197     incidentTheta = std::numeric_limits<float>::quiet_NaN();
0198     nch = 0;
0199     cSize[0] = 0;
0200     cSize[1] = 0;
0201     chId[0].clear();
0202     chId[1].clear();
0203     chValue.clear();
0204   }
0205 };
0206
0207 RootMeasurementWriter::RootMeasurementWriter(
0208     const RootMeasurementWriter::Config& config, Acts::Logging::Level level)
0209     : WriterT(config.inputMeasurements, "RootMeasurementWriter", level),
0210       m_cfg(config) {
0211   // Input container for measurements is already checked by base constructor
0212   if (m_cfg.inputSimHits.empty()) {
0213     throw std::invalid_argument("Missing simulated hits input collection");
0214   }
0215   if (m_cfg.inputMeasurementSimHitsMap.empty()) {
0216     throw std::invalid_argument(
0217         "Missing hit-to-simulated-hits map input collection");
0218   }
0219
0220   m_inputClusters.maybeInitialize(m_cfg.inputClusters);
0221   m_inputSimHits.initialize(m_cfg.inputSimHits);
0222   m_inputMeasurementSimHitsMap.initialize(m_cfg.inputMeasurementSimHitsMap);
0223
0224   if (m_cfg.surfaceByIdentifier.empty()) {
0225     throw std::invalid_argument("Missing Surface-GeoID association map");
0226   }
0227   // Setup ROOT File
0228   m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0229   if (m_outputFile == nullptr) {
0230     throw std::ios_base::failure("Could not open '" + m_cfg.filePath + "'");
0231   }
0232
0233   m_outputFile->cd();
0234
0235   std::vector bIndices = {Acts::eBoundLoc0, Acts::eBoundLoc1, Acts::eBoundTime};
0236   m_outputTree =
0237       std::make_unique<DigitizationTree>("measurements", bIndices, bIndices);
0238 }
0239
0240 RootMeasurementWriter::~RootMeasurementWriter() {
0241   if (m_outputFile != nullptr) {
0242     m_outputFile->Close();
0243   }
0244 }
0245
0246 ProcessCode RootMeasurementWriter::finalize() {
0247   /// Close the file if it's yours
0248   m_outputFile->cd();
0249   m_outputTree->tree->Write();
0250   m_outputFile->Close();
0251
0252   return ProcessCode::SUCCESS;
0253 }
0254
0255 ProcessCode RootMeasurementWriter::writeT(
0256     const AlgorithmContext& ctx, const MeasurementContainer& measurements) {
0257   const auto& simHits = m_inputSimHits(ctx);
0258   const auto& hitSimHitsMap = m_inputMeasurementSimHitsMap(ctx);
0259
0260   const ClusterContainer* clusters = nullptr;
0261   if (!m_cfg.inputClusters.empty()) {
0262     clusters = &m_inputClusters(ctx);
0263   }
0264
0265   // Exclusive access to the tree while writing
0266   std::lock_guard<std::mutex> lock(m_writeMutex);
0267
0268   for (Index hitIdx = 0u; hitIdx < measurements.size(); ++hitIdx) {
0269     const ConstVariableBoundMeasurementProxy meas =
0270         measurements.getMeasurement(hitIdx);
0271
0272     Acts::GeometryIdentifier geoId = meas.geometryId();
0273     // find the corresponding surface
0274     auto surfaceItr = m_cfg.surfaceByIdentifier.find(geoId);
0275     if (surfaceItr == m_cfg.surfaceByIdentifier.end()) {
0276       continue;
0277     }
0278     const Acts::Surface& surface = *(surfaceItr->second);
0279
0280     // Fill the identification
0281     m_outputTree->fillIdentification(ctx.eventNumber, geoId);
0282
0283     // Find the contributing simulated hits
0284     auto indices = makeRange(hitSimHitsMap.equal_range(hitIdx));
0285     // Use average truth in the case of multiple contributing sim hits
0286     auto [local, pos4, dir] =
0287         averageSimHits(ctx.geoContext, surface, simHits, indices, logger());
0288     Acts::RotationMatrix3 rot =
0289         surface
0290             .referenceFrame(ctx.geoContext, pos4.segment<3>(Acts::ePos0), dir)
0291             .inverse();
0292     std::pair<double, double> angles =
0293         Acts::VectorHelpers::incidentAngles(dir, rot);
0294
0295     m_outputTree->fillTruthParameters(local, pos4, dir, angles);
0296     m_outputTree->fillBoundMeasurement(meas);
0297     if (clusters != nullptr) {
0298       const auto& c = (*clusters)[hitIdx];
0299       m_outputTree->fillCluster(c);
0300     }
0301
0302     m_outputTree->fill();
0303     m_outputTree->clear();
0304   }
0305
0306   return ProcessCode::SUCCESS;
0307 }
0308
0309 }  // namespace ActsExamples