EnergyLossExtrapolationValidation.cxx

Go to the documentation of this file.

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

// EnergyLossValidation.cxx, (c) ATLAS Detector software
 
// Tracking
 
 
#include "TrkExAlgs/EnergyLossExtrapolationValidation.h"
#include "TrkExInterfaces/IExtrapolator.h"
#include "TrkSurfaces/CylinderSurface.h"
#include "TrkSurfaces/DiscSurface.h"
#include "TrkSurfaces/PlaneSurface.h"
#include "TrkEventPrimitives/ParticleHypothesis.h"
#include "TrkGeometry/TrackingGeometry.h"
#include "TrkGeometry/TrackingVolume.h"
#include "TrkVolumes/CylinderVolumeBounds.h"
#include "TrkTrack/TrackStateOnSurface.h"
// Validation mode - TTree includes
#include "TTree.h"
#include "GaudiKernel/ITHistSvc.h"
#include "GaudiKernel/MsgStream.h"
#include <cmath>
 
//================ Constructor =================================================
 
Trk::EnergyLossExtrapolationValidation::EnergyLossExtrapolationValidation(const std::string& name,
                                                                          ISvcLocator* pSvcLocator)
  : AthAlgorithm(name, pSvcLocator) {}
 
//================ Destructor =================================================
 
Trk::EnergyLossExtrapolationValidation::~EnergyLossExtrapolationValidation() {
  // clear random number generators
  delete m_gaussDist;
  delete m_flatDist;
  // clear data vectors
  delete m_theCylinders;
  delete m_theDiscs1;
  delete m_theDiscs2;
}
 
//================ Initialization =================================================
 
StatusCode Trk::EnergyLossExtrapolationValidation::initialize() {
  // Code entered here will be executed once at program start.
  ATH_MSG_INFO("initialize()");
  ATH_MSG_DEBUG("initialize() m_materialCollectionValidation = " << m_materialCollectionValidation);
 
  // Get Extrapolator from ToolService
  if (m_extrapolator.retrieve().isFailure()) {
    ATH_MSG_FATAL("initialize() Could not retrieve Tool " << m_extrapolator << ". Exiting.");
    return StatusCode::FAILURE;
  }
 
  // create the new Trees
  m_validationTree = new TTree(m_validationTreeName.value().c_str(),
                               m_validationTreeDescription.value().c_str());
  m_validationRunTree = new TTree(m_validationRunTreeName.value().c_str(),
                                  m_validationRunTreeDescription.value().c_str());
 
  // the branches for the parameters
  m_validationTree->Branch("Entries", &m_entries, "entries/i");
  m_validationTree->Branch("Energy", m_energy, "energy[entries]/F");
  m_validationTree->Branch("EnergyLoss", m_energyLoss, "eLoss[entries]/F");
  m_validationTree->Branch("tInX0", m_parameterX0, "tinX0[entries]/F");
  m_validationTree->Branch("Radius", m_radius, "radius[entries]/F");
  m_validationTree->Branch("PosX", m_positionX, "posX[entries]/F");
  m_validationTree->Branch("PosY", m_positionY, "posY[entries]/F");
  m_validationTree->Branch("PosZ", m_positionZ, "posZ[entries]/F");
  m_validationTree->Branch("Eta", m_parameterEta, "eta[entries]/F");
  m_validationTree->Branch("Phi", m_parameterPhi, "phi[entries]/F");
  m_validationTree->Branch("Layer", m_layer, "layer[entries]/i");
 
  m_validationRunTree->Branch("Layers", &m_cylinders, "layers/i");
  m_validationRunTree->Branch("CylR", m_cylinderR, "cylR[layers]/F");
  m_validationRunTree->Branch("CylZ", m_cylinderZ, "cylZ[layers]/F");
  m_validationRunTree->Branch("Momentum", &m_momentum, "momentum/F");
  m_validationRunTree->Branch("UsePt", &m_usePt, "usePt/O");
  m_validationRunTree->Branch("MinEta", &m_minEta, "minEta/F");
  m_validationRunTree->Branch("MaxEta", &m_maxEta, "maxEta/F");
  m_validationRunTree->Branch("PDG", &m_pdg, "pdg/I");
  m_validationRunTree->Branch("Events", &m_events, "events/i");
  m_validationRunTree->Branch("AvgRecordedLayers", &m_avgRecordedLayers, "avgRecLayers/F");
 
  // now register the Trees
  SmartIF<ITHistSvc> tHistSvc {
    service("THistSvc")
  };
  if (!tHistSvc) {
    ATH_MSG_ERROR("initialize() Could not find Hist Service -> Switching ValidationMode Off !");
    delete m_validationTree;
    m_validationTree = nullptr;
    delete m_validationRunTree;
    m_validationRunTree = nullptr;
  }
  if ((tHistSvc->regTree(m_validationTreeFolder, m_validationTree)).isFailure()
      || (tHistSvc->regTree(m_validationRunTreeFolder, m_validationRunTree)).isFailure()) {
    ATH_MSG_ERROR("initialize() Could not register the validation Trees -> Switching ValidationMode Off !");
    delete m_validationTree;
    m_validationTree = nullptr;
    delete m_validationRunTree;
    m_validationRunTree = nullptr;
  }
 
  // initialize the random number generators
  ATH_MSG_INFO("initialize() RandomService = " << randSvc()->name());
  m_gaussDist = new Rndm::Numbers(randSvc(), Rndm::Gauss(0., 1.));
  m_flatDist = new Rndm::Numbers(randSvc(), Rndm::Flat(0., 1.));
 
  // initialize cylinders if they are not set in jobOptions
  double const s_cylInitR[TRKEXALGS_MAXPARAMETERS] = {
    0.5, 34.5, 250, 550, 1120, 4250, 13000, 0, 0, 0
  };
  double const s_cylInitZ[TRKEXALGS_MAXPARAMETERS] = {
    100, 10e6, 680, 2820, 3120, 6500, 22000, 0, 0, 0
  };
 
  // output of vector from jobOptions
  ATH_MSG_INFO("initialize() cylinder dimensions vector from jobOptions :");
  for (size_t lay = 0; lay < m_cylinders + 1; ++lay) {
    ATH_MSG_INFO(
      "initialize() m_cylinderVR[" << lay << "] = " << m_cylinderVR[lay] << "\t ... m_cylinderVZ[" << lay << "] = " <<
        m_cylinderVZ[lay]);
  }
  // transform vector (from jobOptions) into array (for ROOT tree)
  ATH_MSG_INFO("initialize() cylinder dimensions array for algorithm and ROOT tree :");
  for (size_t lay = 0; lay < m_cylinders + 1; ++lay) {
    m_cylinderR[lay] = m_cylinderVR[lay] > 0 ? m_cylinderVR[lay] : s_cylInitR[lay];
    m_cylinderZ[lay] = m_cylinderVZ[lay] > 0 ? m_cylinderVZ[lay] : s_cylInitZ[lay];
    // in "strict onion mode", constrain m_cylinders if the values don't make sense
    if (m_onion && lay > 0 && (m_cylinderR[lay] < m_cylinderR[lay - 1])) {
      ATH_MSG_WARNING(
        "initialize() layer " << lay << "dimensions are smaller than those of layer " << lay - 1 << " - constraining m_cylinders to " << lay -
          1);
      ATH_MSG_INFO("initialize() cutting off here :");
      ATH_MSG_INFO(
        "initialize() m_cylinderR[" << lay << "] = " << m_cylinderR[lay] << "\t ... m_cylinderZ[" << lay << "] = " <<
          m_cylinderZ[lay]);
      m_cylinders = lay - 1;
      break;
    }
    ATH_MSG_INFO(
      "initialize() m_cylinderR[" << lay << "] = " << m_cylinderR[lay] << "\t ... m_cylinderZ[" << lay << "] = " <<
        m_cylinderZ[lay]);
  }
 
  // fill data vector with cylinders once (in order not to create them every time)
  m_theCylinders = new DataVector<const Trk::CylinderSurface>();
  m_theDiscs1 = new DataVector<const Trk::DiscSurface>();
  m_theDiscs2 = new DataVector<const Trk::DiscSurface>();
  for (size_t lay = 0; lay < m_cylinders + 1; ++lay) {
    m_theCylinders->push_back(new Trk::CylinderSurface(Amg::Transform3D(), m_cylinderR[lay], m_cylinderZ[lay]));
    ATH_MSG_INFO("initialize() Cylinder " << lay << ": " << *m_theCylinders->at(lay));
    m_theDiscs1->push_back(new Trk::DiscSurface(*createTransform(0., 0., -m_cylinderZ[lay]), 0., m_cylinderR[lay]));
    ATH_MSG_INFO("initialize() Disc1 " << lay << ": " << *m_theDiscs1->at(lay));
    m_theDiscs2->push_back(new Trk::DiscSurface(*createTransform(0., 0., m_cylinderZ[lay]), 0., m_cylinderR[lay]));
    ATH_MSG_INFO("initialize() Disc2 " << lay << ": " << *m_theDiscs2->at(lay));
  }
 
  if (m_particleType == 0) m_pdg = 999;                                                // geantino
  else if (m_particleType == 1) m_pdg = 11; // electron
  else if (m_particleType == 2) m_pdg = 13; // muon-
  else if (m_particleType == 3) m_pdg = 211; // pion+
  else if (m_particleType == 4) m_pdg = 321; // kaon+
  ATH_MSG_INFO("initialize() ParticleType = " << m_particleType << " ... PDG = " << m_pdg);
 
  ATH_MSG_INFO("initialize() successful");
  return StatusCode::SUCCESS;
}
 
//================ Finalisation =================================================
 
StatusCode Trk::EnergyLossExtrapolationValidation::finalize() {
  // Code entered here will be executed once at the end of the program run.
  ATH_MSG_INFO("finalize() ================== Output Statistics =========================");
  ATH_MSG_INFO("finalize() = Navigation : ");
  ATH_MSG_INFO("finalize() =  - breaks fwd : " << static_cast<double>(m_breaksForward) /
    static_cast<double>(m_triesForward)
                                               << " (" << m_breaksForward << "/" << m_triesForward << ")");
  ATH_MSG_INFO("finalize() =  - breaks bwd : " << static_cast<double>(m_breaksBack) / static_cast<double>(m_triesBack)
                                               << " (" << m_breaksBack << "/" << m_triesBack << ")");
  if (m_materialCollectionValidation) {
    ATH_MSG_INFO("finalize() = Material collection : ");
    ATH_MSG_INFO("finalize() =  - layer collected fwd : " << m_collectedLayerForward);
    ATH_MSG_INFO("finalize() =  - layer collected bwd : " << m_collectedLayerBack);
  }
  ATH_MSG_INFO("finalize() ==============================================================");
 
  m_avgRecordedLayers = m_events ? static_cast<float>(m_totalRecordedLayers) / static_cast<float>(m_events) : 0;
  ++m_cylinders;
  if (m_validationRunTree) m_validationRunTree->Fill();
  --m_cylinders;
 
  return StatusCode::SUCCESS;
}
 
//================ Execution ====================================================
 
StatusCode Trk::EnergyLossExtrapolationValidation::execute(const EventContext& ctx) {
  // get the overall dimensions
  if (!m_highestVolume) {
    // get TrackingGeometry and highest volume
    const Trk::TrackingGeometry* trackingGeometry = m_extrapolator->trackingGeometry();
    m_highestVolume = trackingGeometry ? trackingGeometry->highestTrackingVolume() : nullptr;
    const Trk::CylinderVolumeBounds* cylBounds = m_highestVolume ?
                                                 dynamic_cast<const Trk::CylinderVolumeBounds*>(&(m_highestVolume->
                                                                                                   volumeBounds())) :
                                                 nullptr;
    // bail out
    if (!cylBounds) {
      ATH_MSG_WARNING("execute() No highest TrackingVolume / no VolumeBounds ... pretty useless!");
      return StatusCode::SUCCESS;
    }
  }
 
 
  m_entries = 0;
  for (size_t par = 0; par < TRKEXALGS_MAXPARAMETERS; ++par) {
    // -----------> start parameters
 
    m_parameterP[par] = 0.;
    m_energy[par] = 0.;
    m_energyLoss[par] = 0.;
    m_parameterX0[par] = 0.;
    m_radius[par] = 0.;
    m_positionX[par] = 0.;
    m_positionY[par] = 0.;
    m_positionZ[par] = 0.;
    m_parameterEta[par] = 0.;
    m_parameterPhi[par] = 0.;
    m_parameterTheta[par] = 0.;
    m_layer[par] = 0;
  }
 
  // the local start parameters
  // are adopted for planar and straight line surfaces
  m_parameterPhi[0] = M_PI * (2 * m_flatDist->shoot() - 1);
  m_parameterEta[0] = m_minEta + m_flatDist->shoot() * (m_maxEta - m_minEta);
  m_parameterTheta[0] = 2. * std::atan(std::exp(-m_parameterEta[0]));
 
  double charge = -1.;
  m_parameterP[0] = m_momentum;
  // convert transverse momentum (pt) to momentum (p) if flag is set: p = pt/sin(theta)
  if (m_usePt) m_parameterP[0] /= std::sin(m_parameterTheta[0]);
  m_parameterQoverP[0] = charge / m_parameterP[0];
 
  double mass = Trk::ParticleMasses::mass[m_particleType];
 
  double energy1 = sqrt(m_parameterP[0] * m_parameterP[0] + mass * mass);
 
  double newX0 = 0;
 
  const Trk::PerigeeSurface perigeeSurface;
  // the initial perigee with random numbers
  Trk::Perigee startParameters(0, // m_parameterLoc1[0],
                               0, // m_parameterLoc2[0],
                               m_parameterPhi[0],
                               m_parameterTheta[0],
                               m_parameterQoverP[0],
                               perigeeSurface);
 
  ATH_MSG_VERBOSE("execute() Start Parameters : " << startParameters);
  ATH_MSG_DEBUG(
    "execute() Start Parameters : [phi,eta] = [ " << startParameters.momentum().phi() << ", " << startParameters.eta() <<
      " ]");
 
  // --------------- propagate to find an intersection ---------------------
 
  // fill the TrackParameters vector with extrapolation from startParameters to dummy cylinder surface
  std::unique_ptr<Trk::TrackParameters> lastParameters {};
  std::unique_ptr<Trk::TrackParameters> newParameters {};
 
  if (!m_materialCollectionValidation) {
    lastParameters = m_extrapolator->extrapolate(
      ctx,
      startParameters,
      *(m_theCylinders->at(0)),
      Trk::alongMomentum,
      true,
      static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
  } else { // material collection validation
    // get the vector of TrackStateOnSurfaces back
    const std::vector<const Trk::TrackStateOnSurface*>* collectedMaterial =
      m_extrapolator->extrapolateM(
        ctx,
        startParameters,
        *(m_theCylinders->at(0)),
        Trk::alongMomentum,
        true,
        static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
 
    // get the last one and clone it
    if (collectedMaterial && !collectedMaterial->empty()) {
      // get the last track state on surface & clone the destination parameters
      const Trk::TrackStateOnSurface* destinationState = collectedMaterial->back();
      lastParameters.reset(destinationState->trackParameters() ? destinationState->trackParameters()->clone() : nullptr);
      m_collectedLayerForward += collectedMaterial->size();
      // delete the layers / cleanup
      for (const auto* tsos : *collectedMaterial) {
        newX0 += tsos->materialEffectsOnTrack() ? tsos->materialEffectsOnTrack()->thicknessInX0() : 0;
        delete tsos;
      }
      ATH_MSG_VERBOSE("execute() newX0 = " << newX0);
    }
  }
 
 
  for (size_t lay = 1; lay < m_cylinders + 1; ++lay) {
    if (!m_onion) newX0 = 0;
 
    if (m_onion) {
      // safety check
      if (!lastParameters) {
        ATH_MSG_WARNING("execute() Layer " << lay << ": start parameters for cylinder NOT found - skip event !");
        break;
      }
      ATH_MSG_VERBOSE("execute() Layer " << lay << ": start parameters for cylinder found: " << *lastParameters);
    }
 
    // trying to extrapolate to cylinder barrel
    newParameters.reset();
    if (!m_materialCollectionValidation) {
      newParameters = m_extrapolator->extrapolate(
        ctx,
        m_onion ? *lastParameters : startParameters,
        *(m_theCylinders->at(lay)),
        Trk::alongMomentum,
        true,
        static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
    } else { // material collection validation
      // get the vector of TrackStateOnSurfaces back
      const std::vector<const Trk::TrackStateOnSurface*>* collectedMaterial =
        m_extrapolator->extrapolateM(ctx,
                                     m_onion ? *lastParameters : startParameters,
                                     *(m_theCylinders->at(lay)),
                                     Trk::alongMomentum,
                                     true,
                                     static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
 
      // get the last one and clone it
      if (collectedMaterial && !collectedMaterial->empty()) {
        // get the last track state on surface & clone the destination parameters
        const Trk::TrackStateOnSurface* destinationState = collectedMaterial->back();
        newParameters.reset(destinationState->trackParameters() ? destinationState->trackParameters()->clone() : nullptr);
        if (m_onion) m_collectedLayerForward += collectedMaterial->size();
        else m_collectedLayerForward = collectedMaterial->size(); // TODO: shouldn't there be something else here?
        // delete the layers / cleanup
        for (const auto* tsos : *collectedMaterial) {
          newX0 += tsos->materialEffectsOnTrack() ? tsos->materialEffectsOnTrack()->thicknessInX0() : 0;
          delete tsos;
        }
        ATH_MSG_VERBOSE("execute() newX0 = " << newX0);
      }
    }
 
    // no intersection with cylinder barrel, now trying disc endcaps
    if (!newParameters) {
      if (!m_materialCollectionValidation) {
        newParameters = m_extrapolator->extrapolate(
          ctx,
          m_onion ? *lastParameters : startParameters,
          (m_parameterEta[0] < 0) ? *(m_theDiscs1->at(lay))
          : *(m_theDiscs2->at(lay)),
          Trk::alongMomentum,
          true,
          static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
      } else { // material collection validation
        // get the vector of TrackStateOnSurfaces back
        const std::vector<const Trk::TrackStateOnSurface*>* collectedMaterial =
          m_extrapolator->extrapolateM(ctx,
                                       m_onion ? *lastParameters : startParameters,
                                       (m_parameterEta[0] < 0) ? *(m_theDiscs1->at(lay)) : *(m_theDiscs2->at(lay)),
                                       Trk::alongMomentum,
                                       true,
                                       static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
 
        // get the last one and clone it
        if (collectedMaterial && !collectedMaterial->empty()) {
          // get the last track state on surface & clone the destination parameters
          const Trk::TrackStateOnSurface* destinationState = collectedMaterial->back();
          newParameters.reset(
            destinationState->trackParameters() ? destinationState->trackParameters()->clone() : nullptr);
          if (m_onion) m_collectedLayerForward += collectedMaterial->size();
          else m_collectedLayerForward = collectedMaterial->size(); // TODO: shouldn't there be something else here?
          // delete the layers / cleanup
          for (const auto* tsos : *collectedMaterial) {
            newX0 += tsos->materialEffectsOnTrack() ? tsos->materialEffectsOnTrack()->thicknessInX0() : 0;
            delete tsos;
          }
          ATH_MSG_VERBOSE("execute() newX0 = " << newX0);
        }
      }
    }
 
    // still no intersection
    if (!newParameters) {
      ATH_MSG_WARNING("execute() Layer " << lay << " intersection did not work !");
    } else if (m_highestVolume && newParameters && !(m_highestVolume->inside(newParameters->position()))) {
      ATH_MSG_WARNING("execute() Layer " << lay << " intersection is outside the known world !");
    } else {
      // get the current surface intersection position
      const Amg::Vector3D& newPosition = newParameters->position();
      ATH_MSG_VERBOSE("execute() Track Parameters at layer " << lay << ": " << *newParameters);
      ATH_MSG_DEBUG(
        "execute() Track Parameters at layer " << lay << ": [r,z] = [ " << newPosition.perp() << ", " <<
          newPosition.z());
 
      // record the surface parameters
      ++m_triesForward;
      ++m_entries;
      m_parameterPhi[m_entries] = newParameters->parameters()[Trk::phi];
      m_parameterEta[m_entries] = newParameters->momentum().eta();
      m_parameterTheta[m_entries] = newParameters->parameters()[Trk::theta];
      m_parameterP[m_entries] = newParameters->momentum().mag();
      m_parameterX0[m_entries] = (float) newX0;
      ATH_MSG_DEBUG("execute() Layer " << lay << ": cumulated X0 = " << m_parameterX0[m_entries]);
 
      // get the current energy and calculate energy loss
      m_energy[m_entries] = sqrt(m_parameterP[m_entries] * m_parameterP[m_entries] + mass * mass);
      m_energyLoss[m_entries] = energy1 - m_energy[m_entries];
      ATH_MSG_DEBUG("execute() Layer " << lay << ": cumulated Energy Loss = " << m_energyLoss[m_entries]);
 
      // record the current layer ID
      m_layer[m_entries] = lay;
      // record the current position
      m_radius[m_entries] = newPosition.perp();
      m_positionX[m_entries] = newPosition.x();
      m_positionY[m_entries] = newPosition.y();
      m_positionZ[m_entries] = newPosition.z();
    }
    lastParameters = std::move(newParameters);
  }
 
 
  m_totalRecordedLayers += m_entries;
  ++m_events;
  // increase m_entries once more before the fill (to account for the "start layer" at index 0 with initial track
  // parameters)
  ++m_entries;
 
  // fill the event tree
  if (m_validationTree) m_validationTree->Fill();
 
  // memory cleanup
  ATH_MSG_DEBUG("execute() deleting DataVector parameters ... ");
 
  return StatusCode::SUCCESS;
}
 
//============================================================================================
std::unique_ptr<Amg::Transform3D> Trk::EnergyLossExtrapolationValidation::createTransform(
  double x, double y, double z, double phi, double theta, double alphaZ) {
  if (phi != 0. && theta != 0.) {
    // create the Start Surface
    Amg::Vector3D surfacePosition(x, y, z);
    // z direction
    Amg::Vector3D surfaceZdirection(cos(phi) * sin(theta),
                                    sin(phi) * sin(theta),
                                    cos(theta));
    // the global z axis
    Amg::Vector3D zAxis(0., 0., 1.);
    // the y direction
    Amg::Vector3D surfaceYdirection(zAxis.cross(surfaceZdirection));
    // the x direction
    Amg::Vector3D surfaceXdirection(surfaceYdirection.cross(surfaceZdirection));
    // the rotation
    Amg::RotationMatrix3D surfaceRotation;
    surfaceRotation.col(0) = surfaceXdirection;
    surfaceRotation.col(1) = surfaceYdirection;
    surfaceRotation.col(2) = surfaceZdirection;
    // return it
    if (alphaZ == 0.) {
      return std::make_unique<Amg::Transform3D>(surfaceRotation, surfacePosition);
    }
    Amg::Transform3D nominalTransform(surfaceRotation, surfacePosition);
    return std::make_unique<Amg::Transform3D>(nominalTransform * Amg::AngleAxis3D(alphaZ, zAxis));
  }
 
  return std::make_unique<Amg::Transform3D>(Amg::Translation3D(x, y, z));
}