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