CETmaterial.cxx

Go to the documentation of this file.

/*
   Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
 */

// CETmaterial.cxx, (c) ATLAS Detector software
 
 
// Tracking
#include "TrkExAlgs/CETmaterial.h"
#include "TrkSurfaces/CylinderSurface.h"
#include "TrkSurfaces/DiscSurface.h"
#include "TrkSurfaces/PerigeeSurface.h"
#include "TrkGeometry/TrackingGeometry.h"
#include "TrkEventPrimitives/ParticleHypothesis.h"
#include "TrkTrack/TrackStateOnSurface.h"
#include "TrkVolumes/CylinderVolumeBounds.h"
#include "EventPrimitives/EventPrimitivesHelpers.h"
 
#include <fstream>
//================ Constructor =================================================
 
Trk::CETmaterial::CETmaterial(const std::string& name, ISvcLocator* pSvcLocator)
  :
  AthAlgorithm(name, pSvcLocator),
  m_chronoStatSvc("ChronoStatSvc", name)
{ }
 
//================ Destructor =================================================
 
Trk::CETmaterial::~CETmaterial() {
  delete m_err;
}
 
//================ Initialization =================================================
 
StatusCode Trk::CETmaterial::initialize() {
  // Code entered here will be executed once at program start.
 
  ATH_MSG_INFO("initialize()");
 
  // Get Extrapolator from ToolService
  if (m_extrapolator.retrieve().isFailure()) {
    ATH_MSG_FATAL("Could not retrieve Tool " << m_extrapolator << ". Exiting.");
    return StatusCode::FAILURE;
  }
  if (m_extraprec.retrieve().isFailure()) {
    ATH_MSG_FATAL("Could not retrieve Tool " << m_extraprec << ". Exiting.");
    return StatusCode::FAILURE;
  }
 
  ATH_MSG_INFO("initialize() successful");
  return StatusCode::SUCCESS;
}
 
//================ Finalization =================================================
 
StatusCode Trk::CETmaterial::finalize() {
  // Code entered here will be executed once at the end of the program run.
 
  if (m_chronoStatSvc) m_chronoStatSvc->chronoPrint("MS::scan");
 
  return StatusCode::SUCCESS;
}
 
//================ Execution ====================================================
 
StatusCode Trk::CETmaterial::execute(const EventContext& ctx) {
  ATH_MSG_INFO("execute()");
  // retrieve outer boundary cylinder surface
  if (!m_outerBoundary) {
    m_trackingGeometry = m_extrapolator->trackingGeometry();
    m_outerBoundary = &(m_trackingGeometry->highestTrackingVolume()->boundarySurfaces()[2]->surfaceRepresentation());
    if (!m_outerBoundary) {
      ATH_MSG_FATAL("Could not retrieve cylinder boundary  from " << m_extrapolator << ". Exiting.");
      return StatusCode::FAILURE;
    }
    ATH_MSG_INFO("boundary retrieved ");
  }
 
  if (m_chronoStatSvc) m_chronoStatSvc->chronoStart("MS::scan");
 
  // scan
  std::vector<const TrackStateOnSurface*> material;
  std::vector<const TrackStateOnSurface*> matPrec;
  double phi = -M_PI;
  for (unsigned int it = 0; it < m_numScan + 1; it++) {
    // the initial perigee
    double z0 = m_minZ0 + (m_maxZ0 - m_minZ0) / m_numScan * it;
    phi += 1 * Gaudi::Units::deg;
    if (phi > M_PI) phi -= 2 * M_PI;
 
    double theta = m_minTheta + (m_maxTheta - m_minTheta) / m_numScan * it;
    double p = m_minP + (m_maxP - m_minP) / m_numScan * it;
    Trk::PerigeeSurface surface(Amg::Vector3D(0., 0., 0.));
    Trk::Perigee initialPerigee(0., z0, phi, theta, m_charge / p, surface);
 
    const Trk::TrackParameters* seed = initialPerigee.clone();
    const Trk::PerigeeSurface& pSf = initialPerigee.associatedSurface();
    material.clear();
    matPrec.clear();
    const Trk::TrackParameters* currPar = seed;
    const Trk::TrackParameters* precPar = seed;
    if (m_domsentry) {
      if (!m_msentry) {
        m_msentry = m_trackingGeometry->trackingVolume("Calo::Containers::Calorimeter");
      }
      if (m_msentry) {
        auto msEntry = m_extrapolator->extrapolateToVolume(
          ctx,
          *currPar,
          *m_msentry,
          Trk::alongMomentum,
          static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
        if (msEntry) {
          printMat(theta, phi,currPar->momentum().mag() - msEntry->momentum().mag(),
            Amg::error(msEntry->covariance()->inverse().eval(), Trk::theta),
            Amg::error(msEntry->covariance()->inverse().eval(), Trk::phi));
          const std::vector<const Trk::TrackStateOnSurface*>* mmsentry = m_extrapolator->extrapolateM(ctx,
            *currPar, msEntry->associatedSurface(), Trk::alongMomentum, false, 
            static_cast<Trk:: ParticleHypothesis>( m_particleType.value()));
          if (mmsentry) {
            for (const auto& entry : *mmsentry) {
              if (entry) {
                ATH_MSG_DEBUG("position:eloss:"
                              << entry->trackParameters()->position() << ":"
                              << entry->trackParameters()->momentum().mag() - currPar->momentum().mag());
              }
            }
            currPar = (mmsentry->back()) ?  mmsentry->back()->trackParameters() : msEntry.release();
            const std::vector<const Trk::TrackStateOnSurface*>* peri = m_extrapolator->extrapolateM(ctx,
              *currPar, pSf, Trk::oppositeMomentum, false, 
              static_cast<Trk::ParticleHypothesis>( m_particleType.value()));
            ATH_MSG_INFO("material scan:backward:");
            if (peri) {
              ATH_MSG_DEBUG("trPar vector size:" << peri->size());
            } else {
              ATH_MSG_ERROR("Perigee pointer is null in CETmaterial.cxx");
              return StatusCode::FAILURE;
            }
            for (const auto& entry : *peri) {
              if (entry && entry->trackParameters()) {
                ATH_MSG_DEBUG("position:eloss:"
                  << entry->trackParameters()->position() << ":"
                  << entry->trackParameters()->momentum().mag() - msEntry->momentum().mag());
              }
            }
 
            if (peri->back() && peri->back()->trackParameters()) {
              ATH_MSG_INFO("extrapolation to perigee:input: "
                   << initialPerigee.parameters()[0] << ","
                   << initialPerigee.parameters()[1] << ","
                   << initialPerigee.parameters()[2] << ","
                   << initialPerigee.parameters()[3] << ","
                   << initialPerigee.momentum().mag());
              ATH_MSG_INFO("extrapolation to perigee:output: "
                   << peri->back()->trackParameters()->parameters()[0] << ","
                   << peri->back()->trackParameters()->parameters()[1] << ","
                   << peri->back()->trackParameters()->parameters()[2] << ","
                   << peri->back()->trackParameters()->parameters()[3] << ","
                   << peri->back()->trackParameters()->momentum().mag());
            } else {
              ATH_MSG_ERROR("extrapolation to perigee failed for input parameters: " << msEntry->parameters());
            }
            delete peri;
            peri = nullptr;
          } else {
            ATH_MSG_ERROR("extrapolation to MSentry failed for input parameters: " << currPar->parameters());
            printMat(theta, phi, 0.);
          }
        }
      }
      delete currPar;
      currPar = nullptr;
      continue;
    }
    if (m_checkStepWise) {
      double matApp = 0.;
      while (currPar) {
        std::pair<std::unique_ptr<Trk::TrackParameters>,const Trk::Layer*> next = 
          m_extrapolator->extrapolateToNextActiveLayerM(
          ctx, *currPar, Trk::alongMomentum, true, material,
          static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
        const Trk::TrackParameters* nextPar = next.first.release();
        const Trk::Layer* lay = next.second;
        currPar = nextPar;
        if (m_doprecision && precPar && currPar) {
          // try to extrapolate to the same surface
          const std::vector<const Trk::TrackStateOnSurface*>* nextPrec = m_extraprec->extrapolateM(
            ctx, *precPar, currPar->associatedSurface(), Trk::alongMomentum, false,
            static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
          delete precPar;
          precPar = nullptr;
          // collect material
          if (nextPrec) {
            for (const auto* i : *nextPrec) {
              const Trk::MaterialEffectsBase* mEff = i->materialEffectsOnTrack();
              const Trk::TrackParameters* trPar = i->trackParameters();
              if (mEff && trPar) {
                matApp += mEff->thicknessInX0();
              }
            }
          }
          // stop of extrapolation failed
          if (!lay || !nextPrec || nextPrec->empty() || !nextPrec->back()) break;
          precPar = nextPrec->back()->trackParameters();
          double mat = 0.;
          if (!material.empty()) for (auto& i : material) {
              if (i->materialEffectsOnTrack()) mat += i->materialEffectsOnTrack()->thicknessInX0();
            }
          if (precPar) printMatComp(theta, phi, currPar,
            lay->enclosingDetachedTrackingVolume()->name(), mat, matApp,
            currPar->parameters()[0] - precPar->parameters()[0],
            currPar->parameters()[1] - precPar->parameters()[1]);
          else if (currPar) {
            //precPar is nullptr here
            ATH_MSG_INFO("expected layer not reached:" << currPar->position());
          }
        }
        if (nextPar && m_printActive) {
          int id = 0;
          if (lay) id = lay->layerType();
          double matc = 0.;
          if (!material.empty()) for (auto& i : material) {
              if (i->materialEffectsOnTrack()) matc += i->materialEffectsOnTrack()->thicknessInX0();
            }
          else ATH_MSG_INFO("mat & error:" << theta << "," << phi << "," << matc << ","
                                           << Amg::error(nextPar->covariance()->inverse().eval(), Trk::theta) << ","
                                           << Amg::error(nextPar->covariance()->inverse().eval(), Trk::phi));
 
          printMatPrec(theta, phi, nextPar, nextPar, matc, id, "unknown");
        }
        if (!lay) break;
      }
      if (m_printMaterial) {
        double mat = 0.;
        if (!material.empty()) for (auto& i : material) {
            if (i->materialEffectsOnTrack()) {
              mat += i->materialEffectsOnTrack()->thicknessInX0();
            }
          }
        printMat(theta, phi, mat);
      }
    } else {
      const std::vector<const Trk::TrackStateOnSurface*>* destParameters = m_extrapolator->extrapolateM(
        ctx, *currPar, *m_outerBoundary, Trk::alongMomentum, false,
        static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
 
      if (m_printMaterial) {
        double mat = 0.;
        if (destParameters) for (const auto* destParameter : *destParameters) {
            const Trk::MaterialEffectsBase* mEff = destParameter->materialEffectsOnTrack();
            const Trk::TrackParameters* trPar = destParameter->trackParameters();
            if (mEff && trPar) {
              mat += mEff->thicknessInX0();
              // find volume
              std::vector<const Trk::DetachedTrackingVolume*> detVols =
                m_extrapolator->trackingGeometry()->lowestDetachedTrackingVolumes(trPar->position());
              if (!detVols.empty()){ 
                printMatScan(theta, phi, trPar->position().perp(), trPar->position().z(), mEff->thicknessInX0(), (detVols)[0]->name());
              } else {
                printMatScan(theta, phi, trPar->position().perp(), trPar->position().z(),
                   mEff->thicknessInX0(), m_extrapolator->trackingGeometry()->lowestStaticTrackingVolume(
                  trPar->position())->volumeName());
              }
            }
          }
        
        printMat(theta, phi, mat);
      }
 
      if (!destParameters || destParameters->empty()) {
        ATH_MSG_ERROR("extrapolation to outer boundary failed for input parameters: " << initialPerigee.parameters());
      } else if (destParameters->back()->trackParameters()) {
        // forward extrapolation ok
        if (m_backward) {
          material.clear();
          const std::vector<const Trk::TrackStateOnSurface*>* peri = m_extrapolator->extrapolateM(
            ctx, *(destParameters->back()->trackParameters()),
            pSf,
            Trk::oppositeMomentum,
            false,
            static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
 
          if (peri) {
            ATH_MSG_INFO("trPar vector size:" << peri->size());
            for (unsigned int i = 0; i < peri->size(); i++){
              ATH_MSG_INFO("position:" << i << "," << (*peri)[i]->trackParameters()->position());
            }
            ATH_MSG_INFO("extrapolation to perigee:input: " << initialPerigee.parameters());
            ATH_MSG_INFO("extrapolation to perigee:output: " << peri->back()->trackParameters()->parameters());
          } else {
            ATH_MSG_ERROR(
              "extrapolation to perigee failed for input parameters: " <<
            destParameters->back()->trackParameters()->parameters());
          }
          delete peri;
          peri = nullptr;
        }
      }
 
      delete destParameters;
      destParameters = nullptr;
    }
  }
 
  if (m_chronoStatSvc) m_chronoStatSvc->chronoStop("MS::scan");
 
  return StatusCode::SUCCESS;
}
 
//============================================================================================
 
void Trk::CETmaterial::printMat(double theta, double phi, double mat, double dtheta, double dphi) const {
  std::ofstream myfilemat;
  myfilemat.open(m_matTotFile, std::ios::app);
  myfilemat << theta << " " << phi << " " << mat << " " << dtheta << " " << dphi << std::endl;
}
 
void Trk::CETmaterial::printMatScan(double theta, double phi, double r, double z, double mat,
                                    const std::string& name) const {
  std::ofstream myfilemat;
  myfilemat.open(m_matScanFile, std::ios::app);
  myfilemat << theta << " " << phi << " " << r << " " << z << " " << mat << " " << name << std::endl;
}
 
void Trk::CETmaterial::printMatPrec(double theta, double phi, const Trk::TrackParameters* nextPar,
                                    const Trk::TrackParameters* mdest, double mat, int id, const std::string& /*name*/) {
  ; // dummy to get rid of warning message (unused variable name)
  std::ofstream myfilemat;
  myfilemat.open(m_matActiveFile, std::ios::app);
  if (!m_th && !m_ph) {
    m_th = theta;
    m_ph = phi;
    m_id = id;
    m_matSaved = mat;
    delete m_next;
    m_next = nextPar->clone();
    delete m_err;
    m_err = nullptr;
    if (mdest) {
      m_err = new Amg::MatrixX;
      *m_err = mdest->covariance()->inverse().eval();
    }
    return;
  }
 
  if (theta != m_th || phi != m_ph) {
    if (m_err && m_id > 0) {
      myfilemat << m_th << " " << m_ph << " " << 1 << " " << m_id << " " << m_matSaved << std::endl;
      myfilemat << m_next->parameters()[Trk::locX] << " " << m_next->parameters()[Trk::locY] << " " <<
      m_next->parameters()[Trk::phi]
                << " " << m_next->parameters()[Trk::theta] << " " << m_next->parameters()[Trk::qOverP] << std::endl;
      myfilemat << Amg::error(*m_err, Trk::locX) << " " << Amg::error(*m_err, Trk::locY)
                << " " << Amg::error(*m_err, Trk::phi) << " " << Amg::error(*m_err, Trk::theta)
                << " " << Amg::error(*m_err, Trk::qOverP) << std::endl;
    } else {
      myfilemat << m_th << " " << m_ph << " " << 0 << " " << m_id << std::endl;
      myfilemat << m_next->parameters()[Trk::locX] << " " << m_next->parameters()[Trk::locY] << " " <<
      m_next->parameters()[Trk::phi]
                << " " << m_next->parameters()[Trk::theta] << " " << m_next->parameters()[Trk::qOverP] << std::endl;
    }
    m_th = theta;
    m_ph = phi;
    m_id = id;
    m_matSaved = mat;
    delete m_next;
    m_next = nextPar->clone();
    delete m_err;
    m_err = nullptr;
 
    if (mdest) {
      m_err = new Amg::MatrixX;
      *m_err = mdest->covariance()->inverse().eval();
    }
    return;
  }
 
  // update data
  if (id > 1) {
    m_th = theta;
    m_ph = phi;
    m_id = id;
    m_matSaved = mat;
    delete m_next;
    m_next = nextPar->clone();
    delete m_err;
    m_err = nullptr;
    if (mdest) {
      m_err = new Amg::MatrixX;
      *m_err = mdest->covariance()->inverse().eval();
    }
  }
}
 
void Trk::CETmaterial::printMatComp(double theta, double phi, const Trk::TrackParameters* currPar,
   const std::string& name, double mat, double matApp, double dx, double dy) const {
  std::ofstream myfilemat;
  myfilemat.open(m_matCompFile, std::ios::app);
  myfilemat << theta << " " << phi << " " << currPar->position().perp() << " " << currPar->position().z() << " " <<
  name.substr(0, 2)
            << " " << mat << " " << matApp << " " << dx << " " << dy << std::endl;
}