EnergyLossExtrapolationValidation.cxx

Go to the documentation of this file.

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

// EnergyLossValidation.cxx, (c) ATLAS Detector software
 
// Tracking
#include <cmath>
 
#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"
 
 
//================ 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 = Gaudi::Hive::currentContext();
    // 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
    const Trk::TrackParameters* lastParameters = nullptr;
    const Trk::TrackParameters* newParameters = nullptr;
 
    if (!m_materialCollectionValidation) {
 
        lastParameters = m_extrapolator->extrapolate(
          ctx,
          startParameters,
          *(m_theCylinders->at(0)),
          Trk::alongMomentum,
          true,
          static_cast<Trk::ParticleHypothesis>(m_particleType.value())).release();
 
    } 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 = 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 = nullptr;
        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())).release();
 
        } 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 = 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())).release();
 
            } 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 = 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 = 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));
}