Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:10:01 

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017-2018 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/Root/RootMaterialTrackWriter.hpp"
0010 
0011 #include "Acts/Geometry/GeometryIdentifier.hpp"
0012 #include "Acts/Geometry/TrackingVolume.hpp"
0013 #include "Acts/Geometry/Volume.hpp"
0014 #include "Acts/Material/Material.hpp"
0015 #include "Acts/Material/MaterialInteraction.hpp"
0016 #include "Acts/Material/MaterialSlab.hpp"
0017 #include "Acts/Surfaces/CylinderBounds.hpp"
0018 #include "Acts/Surfaces/RadialBounds.hpp"
0019 #include "Acts/Surfaces/Surface.hpp"
0020 #include "Acts/Surfaces/SurfaceBounds.hpp"
0021 #include "Acts/Utilities/Intersection.hpp"
0022 #include "Acts/Utilities/Logger.hpp"
0023 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0024 
0025 #include <algorithm>
0026 #include <cstddef>
0027 #include <ios>
0028 #include <stdexcept>
0029 #include <type_traits>
0030 
0031 #include <TFile.h>
0032 #include <TTree.h>
0033 
0034 using Acts::VectorHelpers::eta;
0035 using Acts::VectorHelpers::perp;
0036 using Acts::VectorHelpers::phi;
0037 
0038 namespace ActsExamples {
0039 
0040 RootMaterialTrackWriter::RootMaterialTrackWriter(
0041     const RootMaterialTrackWriter::Config& config, Acts::Logging::Level level)
0042     : WriterT(config.inputMaterialTracks, "RootMaterialTrackWriter", level),
0043       m_cfg(config) {
0044   // An input collection name and tree name must be specified
0045   if (m_cfg.inputMaterialTracks.empty()) {
0046     throw std::invalid_argument("Missing input collection");
0047   } else if (m_cfg.treeName.empty()) {
0048     throw std::invalid_argument("Missing tree name");
0049   }
0050 
0051   // Setup ROOT I/O
0052   m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0053   if (m_outputFile == nullptr) {
0054     throw std::ios_base::failure("Could not open '" + m_cfg.filePath + "'");
0055   }
0056 
0057   m_outputFile->cd();
0058   m_outputTree =
0059       new TTree(m_cfg.treeName.c_str(), "TTree from RootMaterialTrackWriter");
0060   if (m_outputTree == nullptr) {
0061     throw std::bad_alloc();
0062   }
0063 
0064   // Set the branches
0065   m_outputTree->Branch("event_id", &m_eventId);
0066   m_outputTree->Branch("v_x", &m_v_x);
0067   m_outputTree->Branch("v_y", &m_v_y);
0068   m_outputTree->Branch("v_z", &m_v_z);
0069   m_outputTree->Branch("v_px", &m_v_px);
0070   m_outputTree->Branch("v_py", &m_v_py);
0071   m_outputTree->Branch("v_pz", &m_v_pz);
0072   m_outputTree->Branch("v_phi", &m_v_phi);
0073   m_outputTree->Branch("v_eta", &m_v_eta);
0074   m_outputTree->Branch("t_X0", &m_tX0);
0075   m_outputTree->Branch("t_L0", &m_tL0);
0076   m_outputTree->Branch("mat_x", &m_step_x);
0077   m_outputTree->Branch("mat_y", &m_step_y);
0078   m_outputTree->Branch("mat_z", &m_step_z);
0079   m_outputTree->Branch("mat_dx", &m_step_dx);
0080   m_outputTree->Branch("mat_dy", &m_step_dy);
0081   m_outputTree->Branch("mat_dz", &m_step_dz);
0082   m_outputTree->Branch("mat_step_length", &m_step_length);
0083   m_outputTree->Branch("mat_X0", &m_step_X0);
0084   m_outputTree->Branch("mat_L0", &m_step_L0);
0085   m_outputTree->Branch("mat_A", &m_step_A);
0086   m_outputTree->Branch("mat_Z", &m_step_Z);
0087   m_outputTree->Branch("mat_rho", &m_step_rho);
0088 
0089   if (m_cfg.prePostStep) {
0090     m_outputTree->Branch("mat_sx", &m_step_sx);
0091     m_outputTree->Branch("mat_sy", &m_step_sy);
0092     m_outputTree->Branch("mat_sz", &m_step_sz);
0093     m_outputTree->Branch("mat_ex", &m_step_ex);
0094     m_outputTree->Branch("mat_ey", &m_step_ey);
0095     m_outputTree->Branch("mat_ez", &m_step_ez);
0096   }
0097   if (m_cfg.storeSurface) {
0098     m_outputTree->Branch("sur_id", &m_sur_id);
0099     m_outputTree->Branch("sur_type", &m_sur_type);
0100     m_outputTree->Branch("sur_x", &m_sur_x);
0101     m_outputTree->Branch("sur_y", &m_sur_y);
0102     m_outputTree->Branch("sur_z", &m_sur_z);
0103     m_outputTree->Branch("sur_pathCorrection", &m_sur_pathCorrection);
0104     m_outputTree->Branch("sur_range_min", &m_sur_range_min);
0105     m_outputTree->Branch("sur_range_max", &m_sur_range_max);
0106   }
0107   if (m_cfg.storeVolume) {
0108     m_outputTree->Branch("vol_id", &m_vol_id);
0109   }
0110 }
0111 
0112 RootMaterialTrackWriter::~RootMaterialTrackWriter() {
0113   if (m_outputFile != nullptr) {
0114     m_outputFile->Close();
0115   }
0116 }
0117 
0118 ProcessCode RootMaterialTrackWriter::finalize() {
0119   // write the tree and close the file
0120   ACTS_INFO("Writing ROOT output File : " << m_cfg.filePath);
0121 
0122   m_outputFile->cd();
0123   m_outputTree->Write();
0124   m_outputFile->Close();
0125 
0126   return ProcessCode::SUCCESS;
0127 }
0128 
0129 ProcessCode RootMaterialTrackWriter::writeT(
0130     const AlgorithmContext& ctx,
0131     const std::unordered_map<std::size_t, Acts::RecordedMaterialTrack>&
0132         materialTracks) {
0133   // Exclusive access to the tree while writing
0134   std::lock_guard<std::mutex> lock(m_writeMutex);
0135 
0136   m_eventId = ctx.eventNumber;
0137   // Loop over the material tracks and write them out
0138   for (auto& [idTrack, mtrack] : materialTracks) {
0139     // Clearing the vector first
0140     m_step_sx.clear();
0141     m_step_sy.clear();
0142     m_step_sz.clear();
0143     m_step_x.clear();
0144     m_step_y.clear();
0145     m_step_z.clear();
0146     m_step_ex.clear();
0147     m_step_ey.clear();
0148     m_step_ez.clear();
0149     m_step_dx.clear();
0150     m_step_dy.clear();
0151     m_step_dz.clear();
0152     m_step_length.clear();
0153     m_step_X0.clear();
0154     m_step_L0.clear();
0155     m_step_A.clear();
0156     m_step_Z.clear();
0157     m_step_rho.clear();
0158 
0159     m_sur_id.clear();
0160     m_sur_type.clear();
0161     m_sur_x.clear();
0162     m_sur_y.clear();
0163     m_sur_z.clear();
0164     m_sur_pathCorrection.clear();
0165     m_sur_range_min.clear();
0166     m_sur_range_max.clear();
0167 
0168     m_vol_id.clear();
0169 
0170     auto materialInteractions = mtrack.second.materialInteractions;
0171     if (m_cfg.collapseInteractions) {
0172       std::vector<Acts::MaterialInteraction> collapsed;
0173 
0174       Acts::Vector3 positionSum = Acts::Vector3::Zero();
0175       double pathCorrectionSum = 0;
0176 
0177       for (std::size_t start = 0, end = 0; end < materialInteractions.size();
0178            ++end) {
0179         const auto& mintStart = materialInteractions[start];
0180         const auto& mintEnd = materialInteractions[end];
0181 
0182         positionSum += mintEnd.position;
0183         pathCorrectionSum += mintEnd.pathCorrection;
0184 
0185         const bool same = mintStart.materialSlab.material() ==
0186                           mintEnd.materialSlab.material();
0187         const bool last = end == materialInteractions.size() - 1;
0188 
0189         if (!same || last) {
0190           auto mint = mintStart;
0191           mint.position = positionSum / (end - start);
0192           mint.pathCorrection = pathCorrectionSum;
0193 
0194           collapsed.push_back(mint);
0195 
0196           start = end;
0197           positionSum = Acts::Vector3::Zero();
0198           pathCorrectionSum = 0;
0199         }
0200       }
0201 
0202       materialInteractions = std::move(collapsed);
0203     }
0204 
0205     // Reserve the vector then
0206     std::size_t mints = materialInteractions.size();
0207     m_step_sx.reserve(mints);
0208     m_step_sy.reserve(mints);
0209     m_step_sz.reserve(mints);
0210     m_step_x.reserve(mints);
0211     m_step_y.reserve(mints);
0212     m_step_z.reserve(mints);
0213     m_step_ex.reserve(mints);
0214     m_step_ey.reserve(mints);
0215     m_step_ez.reserve(mints);
0216     m_step_dx.reserve(mints);
0217     m_step_dy.reserve(mints);
0218     m_step_dz.reserve(mints);
0219     m_step_length.reserve(mints);
0220     m_step_X0.reserve(mints);
0221     m_step_L0.reserve(mints);
0222     m_step_A.reserve(mints);
0223     m_step_Z.reserve(mints);
0224     m_step_rho.reserve(mints);
0225 
0226     m_sur_id.reserve(mints);
0227     m_sur_type.reserve(mints);
0228     m_sur_x.reserve(mints);
0229     m_sur_y.reserve(mints);
0230     m_sur_z.reserve(mints);
0231     m_sur_pathCorrection.reserve(mints);
0232     m_sur_range_min.reserve(mints);
0233     m_sur_range_max.reserve(mints);
0234 
0235     m_vol_id.reserve(mints);
0236 
0237     // reset the global counter
0238     if (m_cfg.recalculateTotals) {
0239       m_tX0 = 0.;
0240       m_tL0 = 0.;
0241     } else {
0242       m_tX0 = mtrack.second.materialInX0;
0243       m_tL0 = mtrack.second.materialInL0;
0244     }
0245 
0246     // set the track information at vertex
0247     m_v_x = mtrack.first.first.x();
0248     m_v_y = mtrack.first.first.y();
0249     m_v_z = mtrack.first.first.z();
0250     m_v_px = mtrack.first.second.x();
0251     m_v_py = mtrack.first.second.y();
0252     m_v_pz = mtrack.first.second.z();
0253     m_v_phi = phi(mtrack.first.second);
0254     m_v_eta = eta(mtrack.first.second);
0255 
0256     // and now loop over the material
0257     for (const auto& mint : materialInteractions) {
0258       auto direction = mint.direction.normalized();
0259 
0260       // The material step position information
0261       m_step_x.push_back(mint.position.x());
0262       m_step_y.push_back(mint.position.y());
0263       m_step_z.push_back(mint.position.z());
0264       m_step_dx.push_back(direction.x());
0265       m_step_dy.push_back(direction.y());
0266       m_step_dz.push_back(direction.z());
0267 
0268       if (m_cfg.prePostStep) {
0269         Acts::Vector3 prePos =
0270             mint.position - 0.5 * mint.pathCorrection * direction;
0271         Acts::Vector3 posPos =
0272             mint.position + 0.5 * mint.pathCorrection * direction;
0273 
0274         m_step_sx.push_back(prePos.x());
0275         m_step_sy.push_back(prePos.y());
0276         m_step_sz.push_back(prePos.z());
0277         m_step_ex.push_back(posPos.x());
0278         m_step_ey.push_back(posPos.y());
0279         m_step_ez.push_back(posPos.z());
0280       }
0281 
0282       // Store surface information
0283       if (m_cfg.storeSurface) {
0284         const Acts::Surface* surface = mint.surface;
0285         if (mint.intersectionID.value() != 0) {
0286           m_sur_id.push_back(mint.intersectionID.value());
0287           m_sur_pathCorrection.push_back(mint.pathCorrection);
0288           m_sur_x.push_back(mint.intersection.x());
0289           m_sur_y.push_back(mint.intersection.y());
0290           m_sur_z.push_back(mint.intersection.z());
0291         } else if (surface != nullptr) {
0292           auto sfIntersection =
0293               surface
0294                   ->intersect(ctx.geoContext, mint.position, mint.direction,
0295                               Acts::BoundaryCheck(true))
0296                   .closest();
0297           m_sur_id.push_back(surface->geometryId().value());
0298           m_sur_pathCorrection.push_back(1.0);
0299           m_sur_x.push_back(sfIntersection.position().x());
0300           m_sur_y.push_back(sfIntersection.position().y());
0301           m_sur_z.push_back(sfIntersection.position().z());
0302         } else {
0303           m_sur_id.push_back(Acts::GeometryIdentifier().value());
0304           m_sur_x.push_back(0);
0305           m_sur_y.push_back(0);
0306           m_sur_z.push_back(0);
0307           m_sur_pathCorrection.push_back(1.0);
0308         }
0309         if (surface != nullptr) {
0310           m_sur_type.push_back(surface->type());
0311           const Acts::SurfaceBounds& surfaceBounds = surface->bounds();
0312           const Acts::RadialBounds* radialBounds =
0313               dynamic_cast<const Acts::RadialBounds*>(&surfaceBounds);
0314           const Acts::CylinderBounds* cylinderBounds =
0315               dynamic_cast<const Acts::CylinderBounds*>(&surfaceBounds);
0316           if (radialBounds != nullptr) {
0317             m_sur_range_min.push_back(radialBounds->rMin());
0318             m_sur_range_max.push_back(radialBounds->rMax());
0319           } else if (cylinderBounds != nullptr) {
0320             m_sur_range_min.push_back(
0321                 -cylinderBounds->get(Acts::CylinderBounds::eHalfLengthZ));
0322             m_sur_range_max.push_back(
0323                 cylinderBounds->get(Acts::CylinderBounds::eHalfLengthZ));
0324           } else {
0325             m_sur_range_min.push_back(0);
0326             m_sur_range_max.push_back(0);
0327           }
0328         } else {
0329           m_sur_type.push_back(-1);
0330           m_sur_range_min.push_back(0);
0331           m_sur_range_max.push_back(0);
0332         }
0333       }
0334 
0335       // store volume information
0336       if (m_cfg.storeVolume) {
0337         Acts::GeometryIdentifier vlayerID;
0338         if (!mint.volume.empty()) {
0339           vlayerID = mint.volume.geometryId();
0340           m_vol_id.push_back(vlayerID.value());
0341         } else {
0342           vlayerID.setVolume(0);
0343           vlayerID.setBoundary(0);
0344           vlayerID.setLayer(0);
0345           vlayerID.setApproach(0);
0346           vlayerID.setSensitive(0);
0347           m_vol_id.push_back(vlayerID.value());
0348         }
0349       }
0350 
0351       // the material information
0352       const auto& mprops = mint.materialSlab;
0353       m_step_length.push_back(mprops.thickness());
0354       m_step_X0.push_back(mprops.material().X0());
0355       m_step_L0.push_back(mprops.material().L0());
0356       m_step_A.push_back(mprops.material().Ar());
0357       m_step_Z.push_back(mprops.material().Z());
0358       m_step_rho.push_back(mprops.material().massDensity());
0359       // re-calculate if defined to do so
0360       if (m_cfg.recalculateTotals) {
0361         m_tX0 += mprops.thicknessInX0();
0362         m_tL0 += mprops.thicknessInL0();
0363       }
0364     }
0365     // write to
0366     m_outputTree->Fill();
0367   }
0368 
0369   // return success
0370   return ProcessCode::SUCCESS;
0371 }
0372 
0373 }  // namespace ActsExamples
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!