Root/src/RootMeasurementIo.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 "ActsPlugins/Root/RootMeasurementIo.hpp"
0010
0011 #include "Acts/Geometry/GeometryIdentifier.hpp"
0012 #include "Acts/Utilities/Enumerate.hpp"
0013 #include "Acts/Utilities/VectorHelpers.hpp"
0014
0015 #include <TTree.h>
0016
0017 using namespace Acts;
0018
0019 ActsPlugins::RootMeasurementIo::RootMeasurementIo(const Config& config)
0020     : m_cfg(config) {
0021   clear();
0022 }
0023
0024 void ActsPlugins::RootMeasurementIo::connectForWrite(TTree& measurementTree) {
0025   measurementTree.Branch("event_nr", &m_measurementPayload.eventNr);
0026   measurementTree.Branch("volume_id", &m_measurementPayload.volumeID);
0027   measurementTree.Branch("layer_id", &m_measurementPayload.layerID);
0028   measurementTree.Branch("surface_id", &m_measurementPayload.surfaceID);
0029   measurementTree.Branch("extra_id", &m_measurementPayload.extraID);
0030
0031   for (auto ib : m_cfg.recoIndices) {
0032     measurementTree.Branch(("rec_" + bNames[ib]).c_str(),
0033                            &m_measurementPayload.recBound[ib]);
0034   }
0035   for (auto ib : m_cfg.recoIndices) {
0036     measurementTree.Branch(("var_" + bNames[ib]).c_str(),
0037                            &m_measurementPayload.varBound[ib]);
0038   }
0039
0040   measurementTree.Branch("rec_gx", &m_measurementPayload.recGx);
0041   measurementTree.Branch("rec_gy", &m_measurementPayload.recGy);
0042   measurementTree.Branch("rec_gz", &m_measurementPayload.recGz);
0043
0044   measurementTree.Branch("clus_size", &m_clusterPayload.nch);
0045   measurementTree.Branch("channel_value", &m_clusterPayload.chValue);
0046   // Both are allocated, but only relevant ones are set
0047   for (auto ib : m_cfg.clusterIndices) {
0048     if (static_cast<unsigned int>(ib) < 2) {
0049       measurementTree.Branch(("channel_" + bNames[ib]).c_str(),
0050                              &m_clusterPayload.chId[ib]);
0051       measurementTree.Branch(("clus_size_" + bNames[ib]).c_str(),
0052                              &m_clusterPayload.clusterSize[ib]);
0053     }
0054   }
0055
0056   for (unsigned int ib = 0; ib < eBoundSize; ++ib) {
0057     measurementTree.Branch(("true_" + bNames[ib]).c_str(),
0058                            &m_measurementPayload.trueBound[ib]);
0059   }
0060   measurementTree.Branch("true_x", &m_measurementPayload.trueGx);
0061   measurementTree.Branch("true_y", &m_measurementPayload.trueGy);
0062   measurementTree.Branch("true_z", &m_measurementPayload.trueGz);
0063   measurementTree.Branch("true_incident_phi",
0064                          &m_measurementPayload.incidentPhi);
0065   measurementTree.Branch("true_incident_theta",
0066                          &m_measurementPayload.incidentTheta);
0067
0068   for (auto ib : m_cfg.recoIndices) {
0069     measurementTree.Branch(("residual_" + bNames[ib]).c_str(),
0070                            &m_measurementPayload.residual[ib]);
0071   }
0072   for (auto ib : m_cfg.recoIndices) {
0073     measurementTree.Branch(("pull_" + bNames[ib]).c_str(),
0074                            &m_measurementPayload.pull[ib]);
0075   }
0076   clear();
0077 }
0078
0079 void ActsPlugins::RootMeasurementIo::fillIdentification(
0080     int evnt, const GeometryIdentifier& geoId) {
0081   m_measurementPayload.eventNr = evnt;
0082   m_measurementPayload.volumeID = static_cast<int>(geoId.volume());
0083   m_measurementPayload.layerID = static_cast<int>(geoId.layer());
0084   m_measurementPayload.surfaceID = static_cast<int>(geoId.sensitive());
0085   m_measurementPayload.extraID = static_cast<int>(geoId.extra());
0086 }
0087
0088 void ActsPlugins::RootMeasurementIo::fillTruthParameters(
0089     const Vector2& lp, const Vector4& xt, const Vector3& dir,
0090     const std::pair<double, double> angles) {
0091   m_measurementPayload.trueBound[eBoundLoc0] = lp[eBoundLoc0];
0092   m_measurementPayload.trueBound[eBoundLoc1] = lp[eBoundLoc1];
0093   m_measurementPayload.trueBound[eBoundPhi] = VectorHelpers::phi(dir);
0094   m_measurementPayload.trueBound[eBoundTheta] = VectorHelpers::theta(dir);
0095   m_measurementPayload.trueBound[eBoundTime] = xt[eTime];
0096
0097   m_measurementPayload.trueGx = xt[ePos0];
0098   m_measurementPayload.trueGy = xt[ePos1];
0099   m_measurementPayload.trueGz = xt[ePos2];
0100
0101   m_measurementPayload.incidentPhi = static_cast<float>(angles.first);
0102   m_measurementPayload.incidentTheta = static_cast<float>(angles.second);
0103 }
0104
0105 void ActsPlugins::RootMeasurementIo::fillBoundMeasurement(
0106     const std::vector<double>& measurements,
0107     const std::vector<double>& variances,
0108     const std::vector<unsigned int>& subspaceIndex) {
0109   if (measurements.size() != subspaceIndex.size() ||
0110       variances.size() != subspaceIndex.size()) {
0111     throw std::invalid_argument(
0112         "Measurement, covariance and access index size mismatch");
0113   }
0114
0115   for (auto [im, m] : enumerate(measurements)) {
0116     auto ib = subspaceIndex[im];
0117
0118     m_measurementPayload.recBound[ib] = static_cast<float>(m);
0119     m_measurementPayload.varBound[ib] = static_cast<float>(variances[im]);
0120     m_measurementPayload.residual[ib] =
0121         m_measurementPayload.recBound[ib] - m_measurementPayload.trueBound[ib];
0122     m_measurementPayload.pull[ib] =
0123         m_measurementPayload.residual[ib] /
0124         std::sqrt(m_measurementPayload.varBound[ib]);
0125   }
0126 }
0127
0128 void ActsPlugins::RootMeasurementIo::fillGlobalPosition(const Vector3& pos) {
0129   m_measurementPayload.recGx = pos.x();
0130   m_measurementPayload.recGy = pos.y();
0131   m_measurementPayload.recGz = pos.z();
0132 }
0133
0134 void ActsPlugins::RootMeasurementIo::fillCluster(
0135     const std::vector<std::tuple<int, int, float>>& channels) {
0136   m_clusterPayload.nch = static_cast<int>(channels.size());
0137   if (m_clusterPayload.nch == 0) {
0138     return;
0139   }
0140   for (auto [ch0, ch1, chv] : channels) {
0141     m_clusterPayload.chId[0].push_back(ch0);
0142     m_clusterPayload.chId[1].push_back(ch1);
0143     m_clusterPayload.chValue.push_back(chv);
0144   }
0145   // Calculate cluster size in 0 and 1 direction
0146   auto [min0, max0] = std::ranges::minmax_element(m_clusterPayload.chId[0]);
0147   auto [min1, max1] = std::ranges::minmax_element(m_clusterPayload.chId[1]);
0148   m_clusterPayload.clusterSize[0] = (*max0 - *min0 + 1);
0149   m_clusterPayload.clusterSize[1] = (*max1 - *min1 + 1);
0150 }
0151
0152 void ActsPlugins::RootMeasurementIo::clear() {
0153   for (unsigned int ib = 0; ib < eBoundSize; ++ib) {
0154     m_measurementPayload.trueBound[ib] =
0155         std::numeric_limits<float>::quiet_NaN();
0156     m_measurementPayload.recBound[ib] = std::numeric_limits<float>::quiet_NaN();
0157     m_measurementPayload.varBound[ib] = std::numeric_limits<float>::quiet_NaN();
0158     m_measurementPayload.residual[ib] = std::numeric_limits<float>::quiet_NaN();
0159     m_measurementPayload.pull[ib] = std::numeric_limits<float>::quiet_NaN();
0160   }
0161   m_measurementPayload.recGx = std::numeric_limits<float>::quiet_NaN();
0162   m_measurementPayload.recGy = std::numeric_limits<float>::quiet_NaN();
0163   m_measurementPayload.recGz = std::numeric_limits<float>::quiet_NaN();
0164   m_measurementPayload.trueGx = std::numeric_limits<float>::quiet_NaN();
0165   m_measurementPayload.trueGy = std::numeric_limits<float>::quiet_NaN();
0166   m_measurementPayload.trueGz = std::numeric_limits<float>::quiet_NaN();
0167   m_measurementPayload.incidentPhi = std::numeric_limits<float>::quiet_NaN();
0168   m_measurementPayload.incidentTheta = std::numeric_limits<float>::quiet_NaN();
0169
0170   m_clusterPayload.nch = 0;
0171   m_clusterPayload.clusterSize[0] = 0;
0172   m_clusterPayload.clusterSize[1] = 0;
0173   m_clusterPayload.chId[0].clear();
0174   m_clusterPayload.chId[1].clear();
0175   m_clusterPayload.chValue.clear();
0176 }