Root/src/RootMaterialTrackIo.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 "Acts/Plugins/Root/RootMaterialTrackIo.hpp"
0010
0011 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0012 #include "Acts/Surfaces/Surface.hpp"
0013 #include "Acts/Utilities/VectorHelpers.hpp"
0014
0015 #include <TChain.h>
0016 #include <TTree.h>
0017
0018 void Acts::RootMaterialTrackIo::connectForRead(TChain& materialChain) {
0019   materialChain.SetBranchAddress("event_id", &m_eventId);
0020   materialChain.SetBranchAddress("v_x", &m_summaryPayload.vX);
0021   materialChain.SetBranchAddress("v_y", &m_summaryPayload.vY);
0022   materialChain.SetBranchAddress("v_z", &m_summaryPayload.vZ);
0023   materialChain.SetBranchAddress("v_px", &m_summaryPayload.vPx);
0024   materialChain.SetBranchAddress("v_py", &m_summaryPayload.vPy);
0025   materialChain.SetBranchAddress("v_pz", &m_summaryPayload.vPz);
0026   materialChain.SetBranchAddress("v_phi", &m_summaryPayload.vPhi);
0027   materialChain.SetBranchAddress("v_eta", &m_summaryPayload.vEta);
0028   materialChain.SetBranchAddress("t_X0", &m_summaryPayload.tX0);
0029   materialChain.SetBranchAddress("t_L0", &m_summaryPayload.tL0);
0030   materialChain.SetBranchAddress("mat_x", &m_stepPayload.stepXPtr);
0031   materialChain.SetBranchAddress("mat_y", &m_stepPayload.stepYPtr);
0032   materialChain.SetBranchAddress("mat_z", &m_stepPayload.stepZPtr);
0033   materialChain.SetBranchAddress("mat_dx", &m_stepPayload.stepDxPtr);
0034   materialChain.SetBranchAddress("mat_dy", &m_stepPayload.stepDyPtr);
0035   materialChain.SetBranchAddress("mat_dz", &m_stepPayload.stepDzPtr);
0036   materialChain.SetBranchAddress("mat_step_length",
0037                                  &m_stepPayload.stepLengthPtr);
0038   materialChain.SetBranchAddress("mat_X0", &m_stepPayload.stepMatX0Ptr);
0039   materialChain.SetBranchAddress("mat_L0", &m_stepPayload.stepMatL0Ptr);
0040   materialChain.SetBranchAddress("mat_A", &m_stepPayload.stepMatAPtr);
0041   materialChain.SetBranchAddress("mat_Z", &m_stepPayload.stepMatZPtr);
0042   materialChain.SetBranchAddress("mat_rho", &m_stepPayload.stepMatRhoPtr);
0043   if (m_cfg.surfaceInfo) {
0044     materialChain.SetBranchAddress("sur_id", &m_surfacePayload.surfaceIdPtr);
0045     materialChain.SetBranchAddress("sur_x", &m_surfacePayload.surfaceXPtr);
0046     materialChain.SetBranchAddress("sur_y", &m_surfacePayload.surfaceYPtr);
0047     materialChain.SetBranchAddress("sur_z", &m_surfacePayload.surfaceZPtr);
0048     materialChain.SetBranchAddress("sur_pathCorrection",
0049                                    &m_surfacePayload.surfacePathCorrectionPtr);
0050   }
0051 }
0052
0053 void Acts::RootMaterialTrackIo::connectForWrite(TTree& materialTree) {
0054   // This sets the branch addresses for the material track
0055   // Set the branches
0056   materialTree.Branch("event_id", &m_eventId);
0057   materialTree.Branch("v_x", &m_summaryPayload.vX);
0058   materialTree.Branch("v_y", &m_summaryPayload.vY);
0059   materialTree.Branch("v_z", &m_summaryPayload.vZ);
0060   materialTree.Branch("v_px", &m_summaryPayload.vPx);
0061   materialTree.Branch("v_py", &m_summaryPayload.vPy);
0062   materialTree.Branch("v_pz", &m_summaryPayload.vPz);
0063   materialTree.Branch("v_phi", &m_summaryPayload.vPhi);
0064   materialTree.Branch("v_eta", &m_summaryPayload.vEta);
0065   materialTree.Branch("t_X0", &m_summaryPayload.tX0);
0066   materialTree.Branch("t_L0", &m_summaryPayload.tL0);
0067   materialTree.Branch("mat_x", &m_stepPayload.stepX);
0068   materialTree.Branch("mat_y", &m_stepPayload.stepY);
0069   materialTree.Branch("mat_z", &m_stepPayload.stepZ);
0070   materialTree.Branch("mat_r", &m_stepPayload.stepR);
0071   materialTree.Branch("mat_dx", &m_stepPayload.stepDx);
0072   materialTree.Branch("mat_dy", &m_stepPayload.stepDy);
0073   materialTree.Branch("mat_dz", &m_stepPayload.stepDz);
0074   materialTree.Branch("mat_step_length", &m_stepPayload.stepLength);
0075   materialTree.Branch("mat_X0", &m_stepPayload.stepMatX0);
0076   materialTree.Branch("mat_L0", &m_stepPayload.stepMatL0);
0077   materialTree.Branch("mat_A", &m_stepPayload.stepMatA);
0078   materialTree.Branch("mat_Z", &m_stepPayload.stepMatZ);
0079   materialTree.Branch("mat_rho", &m_stepPayload.stepMatRho);
0080
0081   if (m_cfg.prePostStepInfo) {
0082     materialTree.Branch("mat_sx", &m_stepPayload.stepXs);
0083     materialTree.Branch("mat_sy", &m_stepPayload.stepYs);
0084     materialTree.Branch("mat_sz", &m_stepPayload.stepZs);
0085     materialTree.Branch("mat_ex", &m_stepPayload.stepXe);
0086     materialTree.Branch("mat_ey", &m_stepPayload.stepYe);
0087     materialTree.Branch("mat_ez", &m_stepPayload.stepZe);
0088   }
0089   if (m_cfg.surfaceInfo) {
0090     materialTree.Branch("sur_id", &m_surfacePayload.surfaceId);
0091     materialTree.Branch("sur_x", &m_surfacePayload.surfaceX);
0092     materialTree.Branch("sur_y", &m_surfacePayload.surfaceY);
0093     materialTree.Branch("sur_z", &m_surfacePayload.surfaceZ);
0094     materialTree.Branch("sur_r", &m_surfacePayload.surfaceR);
0095     materialTree.Branch("sur_distance", &m_surfacePayload.surfaceDistance);
0096     materialTree.Branch("sur_pathCorrection",
0097                         &m_surfacePayload.surfacePathCorrection);
0098   }
0099   if (m_cfg.volumeInfo) {
0100     materialTree.Branch("vol_id", &m_volumePayload.volumeId);
0101   }
0102 }
0103
0104 void Acts::RootMaterialTrackIo::write(
0105     const GeometryContext& gctx, std::uint32_t eventNum,
0106     const RecordedMaterialTrack& materialTrack) {
0107   m_eventId = eventNum;
0108   // Clearing the vector first
0109   m_stepPayload.stepXs.clear();
0110   m_stepPayload.stepYs.clear();
0111   m_stepPayload.stepZs.clear();
0112   m_stepPayload.stepX.clear();
0113   m_stepPayload.stepY.clear();
0114   m_stepPayload.stepZ.clear();
0115   m_stepPayload.stepR.clear();
0116   m_stepPayload.stepXe.clear();
0117   m_stepPayload.stepYe.clear();
0118   m_stepPayload.stepZe.clear();
0119   m_stepPayload.stepDx.clear();
0120   m_stepPayload.stepDy.clear();
0121   m_stepPayload.stepDz.clear();
0122   m_stepPayload.stepLength.clear();
0123   m_stepPayload.stepMatX0.clear();
0124   m_stepPayload.stepMatL0.clear();
0125   m_stepPayload.stepMatA.clear();
0126   m_stepPayload.stepMatZ.clear();
0127   m_stepPayload.stepMatRho.clear();
0128
0129   m_surfacePayload.surfaceId.clear();
0130   m_surfacePayload.surfaceX.clear();
0131   m_surfacePayload.surfaceY.clear();
0132   m_surfacePayload.surfaceZ.clear();
0133   m_surfacePayload.surfaceR.clear();
0134   m_surfacePayload.surfaceDistance.clear();
0135   m_surfacePayload.surfacePathCorrection.clear();
0136
0137   m_volumePayload.volumeId.clear();
0138
0139   auto materialInteractions = materialTrack.second.materialInteractions;
0140
0141   // Reserve the vector then
0142   std::size_t mints = materialInteractions.size();
0143   m_stepPayload.stepXs.reserve(mints);
0144   m_stepPayload.stepYs.reserve(mints);
0145   m_stepPayload.stepZs.reserve(mints);
0146   m_stepPayload.stepX.reserve(mints);
0147   m_stepPayload.stepY.reserve(mints);
0148   m_stepPayload.stepZ.reserve(mints);
0149   m_stepPayload.stepR.reserve(mints);
0150   m_stepPayload.stepXe.reserve(mints);
0151   m_stepPayload.stepYe.reserve(mints);
0152   m_stepPayload.stepZe.reserve(mints);
0153   m_stepPayload.stepDx.reserve(mints);
0154   m_stepPayload.stepDy.reserve(mints);
0155   m_stepPayload.stepDz.reserve(mints);
0156   m_stepPayload.stepLength.reserve(mints);
0157   m_stepPayload.stepMatX0.reserve(mints);
0158   m_stepPayload.stepMatL0.reserve(mints);
0159   m_stepPayload.stepMatA.reserve(mints);
0160   m_stepPayload.stepMatZ.reserve(mints);
0161   m_stepPayload.stepMatRho.reserve(mints);
0162
0163   m_surfacePayload.surfaceId.reserve(mints);
0164   m_surfacePayload.surfaceX.reserve(mints);
0165   m_surfacePayload.surfaceY.reserve(mints);
0166   m_surfacePayload.surfaceZ.reserve(mints);
0167   m_surfacePayload.surfaceR.reserve(mints);
0168   m_surfacePayload.surfaceDistance.reserve(mints);
0169   m_surfacePayload.surfacePathCorrection.reserve(mints);
0170
0171   m_volumePayload.volumeId.reserve(mints);
0172
0173   // reset the global counter
0174   if (m_cfg.recalculateTotals) {
0175     m_summaryPayload.tX0 = 0.;
0176     m_summaryPayload.tL0 = 0.;
0177   } else {
0178     m_summaryPayload.tX0 = materialTrack.second.materialInX0;
0179     m_summaryPayload.tL0 = materialTrack.second.materialInL0;
0180   }
0181
0182   // set the track information at vertex
0183   m_summaryPayload.vX = materialTrack.first.first.x();
0184   m_summaryPayload.vY = materialTrack.first.first.y();
0185   m_summaryPayload.vZ = materialTrack.first.first.z();
0186   m_summaryPayload.vPx = materialTrack.first.second.x();
0187   m_summaryPayload.vPy = materialTrack.first.second.y();
0188   m_summaryPayload.vPz = materialTrack.first.second.z();
0189   m_summaryPayload.vPhi = VectorHelpers::phi(materialTrack.first.second);
0190   m_summaryPayload.vEta = VectorHelpers::eta(materialTrack.first.second);
0191
0192   // and now loop over the material
0193   for (const auto& mint : materialInteractions) {
0194     auto direction = mint.direction.normalized();
0195
0196     // The material step position information
0197     m_stepPayload.stepX.push_back(mint.position.x());
0198     m_stepPayload.stepY.push_back(mint.position.y());
0199     m_stepPayload.stepZ.push_back(mint.position.z());
0200     m_stepPayload.stepR.push_back(VectorHelpers::perp(mint.position));
0201     m_stepPayload.stepDx.push_back(direction.x());
0202     m_stepPayload.stepDy.push_back(direction.y());
0203     m_stepPayload.stepDz.push_back(direction.z());
0204
0205     if (m_cfg.prePostStepInfo) {
0206       Acts::Vector3 prePos =
0207           mint.position - 0.5 * mint.pathCorrection * direction;
0208       Acts::Vector3 posPos =
0209           mint.position + 0.5 * mint.pathCorrection * direction;
0210
0211       m_stepPayload.stepXs.push_back(prePos.x());
0212       m_stepPayload.stepYs.push_back(prePos.y());
0213       m_stepPayload.stepZs.push_back(prePos.z());
0214       m_stepPayload.stepXe.push_back(posPos.x());
0215       m_stepPayload.stepYe.push_back(posPos.y());
0216       m_stepPayload.stepZe.push_back(posPos.z());
0217     }
0218
0219     // Store surface information
0220     if (m_cfg.surfaceInfo) {
0221       const Acts::Surface* surface = mint.surface;
0222       if (mint.intersectionID.value() != 0) {
0223         m_surfacePayload.surfaceId.push_back(mint.intersectionID.value());
0224         m_surfacePayload.surfacePathCorrection.push_back(mint.pathCorrection);
0225         m_surfacePayload.surfaceX.push_back(mint.intersection.x());
0226         m_surfacePayload.surfaceY.push_back(mint.intersection.y());
0227         m_surfacePayload.surfaceZ.push_back(mint.intersection.z());
0228         m_surfacePayload.surfaceR.push_back(
0229             VectorHelpers::perp(mint.intersection));
0230         m_surfacePayload.surfaceDistance.push_back(
0231             (mint.position - mint.intersection).norm());
0232       } else if (surface != nullptr) {
0233         auto sfIntersection =
0234             surface
0235                 ->intersect(gctx, mint.position, mint.direction,
0236                             Acts::BoundaryTolerance::None())
0237                 .closest();
0238         m_surfacePayload.surfaceId.push_back(surface->geometryId().value());
0239         m_surfacePayload.surfacePathCorrection.push_back(1.0);
0240         m_surfacePayload.surfaceX.push_back(sfIntersection.position().x());
0241         m_surfacePayload.surfaceY.push_back(sfIntersection.position().y());
0242         m_surfacePayload.surfaceZ.push_back(sfIntersection.position().z());
0243       } else {
0244         m_surfacePayload.surfaceId.push_back(
0245             Acts::GeometryIdentifier().value());
0246         m_surfacePayload.surfaceX.push_back(0);
0247         m_surfacePayload.surfaceY.push_back(0);
0248         m_surfacePayload.surfaceZ.push_back(0);
0249         m_surfacePayload.surfacePathCorrection.push_back(1.0);
0250       }
0251     }
0252
0253     // store volume information
0254     if (m_cfg.volumeInfo) {
0255       Acts::GeometryIdentifier vlayerID;
0256       if (!mint.volume.empty()) {
0257         vlayerID = mint.volume.geometryId();
0258         m_volumePayload.volumeId.push_back(vlayerID.value());
0259       } else {
0260         vlayerID = vlayerID.withVolume(0)
0261                        .withBoundary(0)
0262                        .withLayer(0)
0263                        .withApproach(0)
0264                        .withSensitive(0);
0265         m_volumePayload.volumeId.push_back(vlayerID.value());
0266       }
0267     }
0268
0269     // the material information
0270     const auto& mprops = mint.materialSlab;
0271     m_stepPayload.stepLength.push_back(mprops.thickness());
0272     m_stepPayload.stepMatX0.push_back(mprops.material().X0());
0273     m_stepPayload.stepMatL0.push_back(mprops.material().L0());
0274     m_stepPayload.stepMatA.push_back(mprops.material().Ar());
0275     m_stepPayload.stepMatZ.push_back(mprops.material().Z());
0276     m_stepPayload.stepMatRho.push_back(mprops.material().massDensity());
0277     // re-calculate if defined to do so
0278     if (m_cfg.recalculateTotals) {
0279       m_summaryPayload.tX0 += mprops.thicknessInX0();
0280       m_summaryPayload.tL0 += mprops.thicknessInL0();
0281     }
0282   }
0283 }
0284
0285 Acts::RecordedMaterialTrack Acts::RootMaterialTrackIo::read() const {
0286   Acts::RecordedMaterialTrack rmTrack;
0287   // Fill the position and momentum
0288   rmTrack.first.first = Acts::Vector3(m_summaryPayload.vX, m_summaryPayload.vY,
0289                                       m_summaryPayload.vZ);
0290   rmTrack.first.second = Acts::Vector3(
0291       m_summaryPayload.vPx, m_summaryPayload.vPy, m_summaryPayload.vPz);
0292
0293   // Fill the individual steps
0294   std::size_t msteps = m_stepPayload.stepLength.size();
0295   rmTrack.second.materialInteractions.reserve(msteps);
0296   rmTrack.second.materialInX0 = 0.;
0297   rmTrack.second.materialInL0 = 0.;
0298
0299   for (std::size_t is = 0; is < msteps; ++is) {
0300     float s = m_stepPayload.stepLength[is];
0301     if (s == 0.) {
0302       continue;
0303     }
0304
0305     float mX0 = m_stepPayload.stepMatX0[is];
0306     float mL0 = m_stepPayload.stepMatL0[is];
0307
0308     rmTrack.second.materialInX0 += s / mX0;
0309     rmTrack.second.materialInL0 += s / mL0;
0310     /// Fill the position & the material
0311     Acts::MaterialInteraction mInteraction;
0312     mInteraction.position =
0313         Acts::Vector3(m_stepPayload.stepX[is], m_stepPayload.stepY[is],
0314                       m_stepPayload.stepZ[is]);
0315     mInteraction.direction =
0316         Acts::Vector3(m_stepPayload.stepDx[is], m_stepPayload.stepDy[is],
0317                       m_stepPayload.stepDz[is]);
0318     mInteraction.materialSlab = Acts::MaterialSlab(
0319         Acts::Material::fromMassDensity(mX0, mL0, m_stepPayload.stepMatA[is],
0320                                         m_stepPayload.stepMatZ[is],
0321                                         m_stepPayload.stepMatRho[is]),
0322         s);
0323     if (m_cfg.surfaceInfo) {
0324       // add the surface information to the interaction this allows the
0325       // mapping to be speed up
0326       mInteraction.intersectionID =
0327           Acts::GeometryIdentifier(m_surfacePayload.surfaceId[is]);
0328       mInteraction.intersection = Acts::Vector3(m_surfacePayload.surfaceX[is],
0329                                                 m_surfacePayload.surfaceY[is],
0330                                                 m_surfacePayload.surfaceZ[is]);
0331       mInteraction.pathCorrection = m_surfacePayload.surfacePathCorrection[is];
0332     } else {
0333       mInteraction.intersectionID = Acts::GeometryIdentifier();
0334       mInteraction.intersection = Acts::Vector3(0, 0, 0);
0335     }
0336     rmTrack.second.materialInteractions.push_back(std::move(mInteraction));
0337   }
0338   return rmTrack;
0339 }