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