Root/src/RootMuonSpacePointWriter.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 // This file is part of the ACTS project.
0010 //
0011 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0012 //
0013 // This Source Code Form is subject to the terms of the Mozilla Public
0014 // License, v. 2.0. If a copy of the MPL was not distributed with this
0015 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0016
0017 #include "ActsExamples/Io/Root/RootMuonSpacePointWriter.hpp"
0018
0019 #include "Acts/Definitions/Units.hpp"
0020 #include "Acts/Seeding/detail/CompSpacePointAuxiliaries.hpp"
0021 #include "Acts/Surfaces/LineBounds.hpp"
0022 #include "Acts/Surfaces/RectangleBounds.hpp"
0023 #include "Acts/Surfaces/Surface.hpp"
0024 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0025 #include "Acts/Utilities/Enumerate.hpp"
0026 #include "Acts/Utilities/VectorHelpers.hpp"
0027
0028 #include "TFile.h"
0029 #include "TTree.h"
0030
0031 using namespace Acts;
0032 using namespace Acts::UnitLiterals;
0033 using namespace Acts::VectorHelpers;
0034 /// @brief Converts a surface Identifier to the one of the surrounding volume
0035 /// @param id: Surface identifier to convert
0036 constexpr GeometryIdentifier toChamberId(const GeometryIdentifier& id) {
0037   return GeometryIdentifier{}.withVolume(id.volume()).withLayer(id.layer());
0038 }
0039 /// @brief Pushes value back to a vector and applies a static cast at the same
0040 /// @param vec: Vector into which the value is pushed back
0041 /// @param val: Value to push
0042 template <typename T, typename T1>
0043 void castPush(std::vector<T>& vec, const T1& val) {
0044   // MARK: fpeMaskBegin(FLTUND, 1, #-1)
0045   const T castedVal = static_cast<T>(val);
0046   vec.push_back(castedVal);
0047   // MARK: fpeMaskEnd(FLTUND)
0048 }
0049 /// @brief Converts an angle in radians into an angle in degree
0050 constexpr double inDeg(const double radians) {
0051   if (Acts::abs(radians) < std::numeric_limits<float>::epsilon()) {
0052     return 0.;
0053   }
0054   using namespace Acts::UnitLiterals;
0055   return radians / 1._degree;
0056 }
0057 namespace ActsExamples {
0058 RootMuonSpacePointWriter::RootMuonSpacePointWriter(const Config& config,
0059                                                    Logging::Level level)
0060     : WriterT(config.inputSpacePoints, "RootMuonSpacePointWriter", level),
0061       m_cfg{config} {
0062   // inputParticles is already checked by base constructor
0063   if (m_cfg.filePath.empty()) {
0064     throw std::invalid_argument("Missing file path");
0065   }
0066   if (m_cfg.treeName.empty()) {
0067     throw std::invalid_argument("Missing tree name");
0068   }
0069
0070   // open root file and create the tree
0071   m_file.reset(TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str()));
0072   if (m_file == nullptr) {
0073     throw std::ios_base::failure("Could not open '" + m_cfg.filePath + "'");
0074   }
0075   m_file->cd();
0076   m_tree = new TTree(m_cfg.treeName.c_str(), m_cfg.treeName.c_str());
0077
0078   m_tree->Branch("event_id", &m_eventId);
0079   m_tree->Branch("spacePoint_bucketId", &m_bucketId);
0080   m_tree->Branch("spacePoint_geometryId", &m_geometryId);
0081   m_tree->Branch("spacePoint_muonId", &m_muonId);
0082   m_tree->Branch("spacePoint_localPosX", &m_localPositionX);
0083   m_tree->Branch("spacePoint_localPosY", &m_localPositionY);
0084   m_tree->Branch("spacePoint_localPosZ", &m_localPositionZ);
0085
0086   m_tree->Branch("spacePoint_sensorDirTheta", &m_sensorDirectionTheta);
0087   m_tree->Branch("spacePoint_sensorDirPhi", &m_sensorDirectionPhi);
0088
0089   m_tree->Branch("spacePoint_toNextDirTheta", &m_toNextSensorTheta);
0090   m_tree->Branch("spacePoint_toNextDirPhi", &m_toNextSensorPhi);
0091
0092   m_tree->Branch("spacePoint_covLoc0", &m_covLoc0);
0093   m_tree->Branch("spacePoint_covLoc1", &m_covLoc1);
0094   m_tree->Branch("spacePoint_covT", &m_covLocT);
0095
0096   m_tree->Branch("spacePoint_driftRadius", &m_driftR);
0097   m_tree->Branch("spacePoint_time", &m_time);
0098
0099   if (m_cfg.writeGlobal && m_cfg.trackingGeometry == nullptr) {
0100     throw std::runtime_error(
0101         "RootMuonSpacePointWriter() - Global coordinates can only be written "
0102         "with a tracking geometry");
0103   }
0104   if (m_cfg.writeGlobal) {
0105     m_tree->Branch("global_positionX", &m_globalPosX);
0106     m_tree->Branch("global_positionY", &m_globalPosY);
0107     m_tree->Branch("global_positionZ", &m_globalPosZ);
0108
0109     m_tree->Branch("global_lowEdgeX", &m_lowEdgeX);
0110     m_tree->Branch("global_lowEdgeY", &m_lowEdgeY);
0111     m_tree->Branch("global_lowEdgeZ", &m_lowEdgeZ);
0112
0113     m_tree->Branch("global_highEdgeX", &m_highEdgeX);
0114     m_tree->Branch("global_highEdgeY", &m_highEdgeY);
0115     m_tree->Branch("global_highEdgeZ", &m_highEdgeZ);
0116   }
0117 }
0118
0119 RootMuonSpacePointWriter::~RootMuonSpacePointWriter() = default;
0120 ProcessCode RootMuonSpacePointWriter::finalize() {
0121   m_file->cd();
0122   m_file->Write();
0123   m_file.reset();
0124   ACTS_INFO("Wrote muon spacepoints to tree '" << m_cfg.treeName << "' in '"
0125                                                << m_cfg.filePath << "'");
0126
0127   return ProcessCode::SUCCESS;
0128 }
0129 ProcessCode RootMuonSpacePointWriter::writeT(
0130     const AlgorithmContext& ctx, const MuonSpacePointContainer& hits) {
0131   std::lock_guard lock{m_mutex};
0132   m_eventId = ctx.eventNumber;
0133   const Acts::GeometryContext gctx{};
0134   for (const auto& [counter, bucket] : enumerate(hits)) {
0135     for (const MuonSpacePoint& writeMe : bucket) {
0136       ACTS_VERBOSE("Dump space point " << writeMe);
0137
0138       castPush(m_bucketId, counter);
0139       castPush(m_geometryId, writeMe.geometryId().value());
0140       castPush(m_muonId, writeMe.id().toInt());
0141       castPush(m_localPositionX, writeMe.localPosition().x());
0142       castPush(m_localPositionY, writeMe.localPosition().y());
0143       castPush(m_localPositionZ, writeMe.localPosition().z());
0144
0145       castPush(m_sensorDirectionTheta, inDeg(theta(writeMe.sensorDirection())));
0146       castPush(m_sensorDirectionPhi, inDeg(phi(writeMe.sensorDirection())));
0147
0148       castPush(m_toNextSensorTheta, inDeg(theta(writeMe.toNextSensor())));
0149       castPush(m_toNextSensorPhi, inDeg(phi(writeMe.toNextSensor())));
0150
0151       const auto& cov = writeMe.covariance();
0152       {
0153         using namespace Experimental::detail;
0154         using enum CompSpacePointAuxiliaries::ResidualIdx;
0155         castPush(m_covLoc0, cov[toUnderlying(nonBending)]);
0156         castPush(m_covLoc1, cov[toUnderlying(bending)]);
0157         castPush(m_covLocT, cov[toUnderlying(time)]);
0158       }
0159       castPush(m_driftR, writeMe.driftRadius());
0160       castPush(m_time, writeMe.time());
0161       if (!m_cfg.writeGlobal || writeMe.geometryId() == GeometryIdentifier{}) {
0162         continue;
0163       }
0164       const Surface* surface =
0165           m_cfg.trackingGeometry->findSurface(writeMe.geometryId());
0166       assert(surface != nullptr);
0167       const TrackingVolume* chambVol =
0168           m_cfg.trackingGeometry->findVolume(toChamberId(writeMe.geometryId()));
0169       assert(chambVol != nullptr);
0170
0171       const Vector3 globPos = chambVol->transform() *
0172                               AngleAxis3{-90._degree, Vector3::UnitZ()} *
0173                               writeMe.localPosition();
0174       castPush(m_globalPosX, globPos.x());
0175       castPush(m_globalPosY, globPos.y());
0176       castPush(m_globalPosZ, globPos.z());
0177
0178       const auto& bounds{surface->bounds()};
0179       const auto& trf{surface->transform(gctx)};
0180       Acts::Vector3 lowEdge{Vector3::Zero()};
0181       Acts::Vector3 highEdge{Vector3::Zero()};
0182       switch (bounds.type()) {
0183         using enum Acts::SurfaceBounds::BoundsType;
0184         case eLine: {
0185           const auto& lBounds = static_cast<const LineBounds&>(bounds);
0186           const double l = lBounds.get(LineBounds::eHalfLengthZ);
0187           lowEdge = trf * (-l * Vector3::UnitZ());
0188           highEdge = trf * (l * Vector3::UnitZ());
0189           break;
0190         }
0191         case eRectangle: {
0192           const auto& rBounds = static_cast<const RectangleBounds&>(bounds);
0193           const double l =
0194               rBounds.get(writeMe.measuresLoc1() ? RectangleBounds::eMaxX
0195                                                  : RectangleBounds::eMaxY);
0196           const auto dimIdx =
0197               static_cast<std::uint32_t>(writeMe.measuresLoc1() == false);
0198           lowEdge = trf * (-l * Vector3::Unit(dimIdx));
0199           highEdge = trf * (l * Vector3::Unit(dimIdx));
0200           break;
0201         }
0202         case eTrapezoid: {
0203           ACTS_WARNING(__FILE__ << ":" << __LINE__ << " Implement me");
0204           break;
0205         }
0206         /// Other surface bound types are not supported
0207         default:
0208           ACTS_ERROR("Unsupported bounds " << surface->toString(gctx));
0209           return ProcessCode::ABORT;
0210       }
0211       castPush(m_lowEdgeX, lowEdge.x());
0212       castPush(m_lowEdgeY, lowEdge.y());
0213       castPush(m_lowEdgeZ, lowEdge.z());
0214       castPush(m_highEdgeX, highEdge.x());
0215       castPush(m_highEdgeY, highEdge.y());
0216       castPush(m_highEdgeZ, highEdge.z());
0217     }
0218   }
0219   m_tree->Fill();
0220
0221   m_geometryId.clear();
0222   m_bucketId.clear();
0223   m_muonId.clear();
0224   m_localPositionX.clear();
0225   m_localPositionY.clear();
0226   m_localPositionZ.clear();
0227   m_sensorDirectionTheta.clear();
0228   m_sensorDirectionPhi.clear();
0229   m_toNextSensorTheta.clear();
0230   m_toNextSensorPhi.clear();
0231   m_covLoc0.clear();
0232   m_covLoc1.clear();
0233   m_covLocT.clear();
0234   m_driftR.clear();
0235   m_time.clear();
0236
0237   m_globalPosX.clear();
0238   m_globalPosY.clear();
0239   m_globalPosZ.clear();
0240
0241   m_lowEdgeX.clear();
0242   m_lowEdgeY.clear();
0243   m_lowEdgeZ.clear();
0244
0245   m_highEdgeX.clear();
0246   m_highEdgeY.clear();
0247   m_highEdgeZ.clear();
0248   return ProcessCode::SUCCESS;
0249 }
0250
0251 }  // namespace ActsExamples