ExtrapolationCell.h

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
 
// ExtrapolationCell.h, (c) ATLAS Detector software
 
#ifndef TRKEXUTILS_EXTRAPOLATIONCELL_H
#define TRKEXUTILS_EXTRAPOLATIONCELL_H
 
#include "TrkDetDescrUtils/GeometrySignature.h"
#include "TrkEventPrimitives/ParticleHypothesis.h"
#include "TrkEventPrimitives/PropDirection.h"
#include "TrkExUtils/MaterialUpdateMode.h"
#include "TrkEventPrimitives/TransportJacobian.h"
#include "TrkGeometry/MagneticFieldProperties.h"
#include "TrkGeometry/MaterialProperties.h"
#include "TrkMaterialOnTrack/EnergyLoss.h"
#include "TrkSurfaces/Surface.h"
 
#ifndef TRKEXUTILS_CHECKPATHMACRO
#define TRKEXUTILS_CHECKPATHMACRO
#define reachedLimit(current, limit, tolerance)                                \
  ((limit) > 0 &&                                                              \
   (((current) < (limit))                                                      \
      ? ((current) - (limit)) * ((current) - (limit)) / ((limit) * (limit)) <  \
          (tolerance) * (tolerance)                                            \
      : true))
#endif
 
namespace Trk {
 
class MaterialProperties;
class TrackingVolume;
class Layer;
// class Surface;
 
class ExtrapolationMode
{
 
public:
  enum eMode
  {
    Direct = 1,      // call propagator directly, no navigation
    Destination = 2, // try to hit the destination, if not other means to stop
    StopWithPathLimit = 3, // stop when the path limit is reached
    StopWithMaterialLimitX0 =
      4, // stop when the material limit is reached in X0
    StopWithMaterialLimitL0 =
      5,                   // stop when the material limit is reached in L0
    StopAtBoundary = 6,    // stop at the next ID / Calo / MS boundary
    CollectSensitive = 7,  // collect parameters on sensitive elements
    CollectPassive = 8,    // collect parameters on passive layers
    CollectBoundary = 9,   // collect parameters on boundary parameters
    CollectMaterial = 10,  // collect all material on the way
    CollectJacobians = 11, // collect the transport jacobians
    CollectPathSteps = 12, // collect the single path steps
    AvoidFallback = 13,    // don't fallback to propagation
    FATRAS = 14            // force initial radialDirection to be outward
  };
};
 
class ExtrapolationConfig
{
public:
  ExtrapolationConfig(unsigned int evalue = 0)
    : m_value(evalue)
  {}
 
  ExtrapolationConfig(const ExtrapolationConfig& eConfig)
    : m_value(eConfig.m_value)
  {}
 
  void addMode(ExtrapolationMode::eMode em)
  {
    // set the bit corresponding to this mode
    m_value |= (1 << int(em));
  }
 
  bool checkMode(ExtrapolationMode::eMode em) const
  {
    // check if the bit is set or not
    return (m_value & (1 << int(em)));
  }
 
private:
  unsigned int m_value;
};
 
class ExtrapolationCode
{
 
public:
  enum eCode
  {
    Unset = 0,              // no code set yet
    InProgress = 1,         // successful : extrapolation in process
    SuccessDestination = 2, // successful : destination reached
    SuccessBoundaryReached =
      3, // successful : boundary reached & configured to do so
    SuccessPathLimit =
      4, // successful : path limit reached & configured to do so
    SuccessMaterialLimit =
      5,           // successful : material limit reached & configured to do so
    Recovered = 6, // successful : recovered & configured to do so
    FailureDestination = 7, // failure    : could not reach destination
    FailureLoop = 8,        // failure    : loop or oscillation between volumes
    FailureNavigation = 9,  // failure    : general navigation failure
    FailureUpdateKill = 10, // failure    : updated track under threshold
    FailureConfiguration = 11, // failure    : general configuration failure
    LeftKnownWorld = 12        // successful ? failure ? if we just knew ...
  };
 
  eCode code;
 
  /* create a simple extrapolation code */
  ExtrapolationCode(eCode c)
    : code(c)
  {}
 
  ExtrapolationCode& operator=(const eCode& ec)
  {
    code = ec;
    return (*this);
  }
 
  bool operator==(const eCode& ec) const { return (ec == code); }
 
  bool operator!=(const eCode& ec) const { return (ec != code); }
 
  bool inProgress() const { return (code == InProgress); }
 
  bool isSuccess() const { return (code > InProgress && code < Recovered); }
 
  bool isSuccessBeforeDestination() const
  {
    return (code > SuccessDestination && code < Recovered);
  }
 
  bool isSuccessOrRecovered() const
  {
    return (code > InProgress && code <= FailureDestination);
  }
 
  bool isFailure() const { return (code > Recovered); }
 
  bool isFailureOrRecovered() const { return (code > SuccessMaterialLimit); };
 
  const std::string& toString() const { return s_ecodeNames[code]; }
 
private:
  static const std::vector<std::string> s_ecodeNames;
};
 
template<class T>
class ExtrapolationStep
{
public:
 const T* parameters;     
 const Surface* surface;  
 const Layer*
     layer;  
 ExtrapolationConfig stepConfiguration;  
 const MaterialProperties* material;     
 Amg::Vector3D materialPosition;  
 double materialScaling;  
 const TransportJacobian*
     transportJacobian;  
 double pathLength;      
 float time;             
 
 ExtrapolationStep(const T* pars = 0, const Surface* sf = nullptr,
                   const ExtrapolationConfig& eConfig = ExtrapolationConfig(),
                   const MaterialProperties* mprop = nullptr,
                   const TransportJacobian* tjac = nullptr,
                   double pLength = -1.)
     : parameters(pars),
       surface(sf),
       layer(nullptr),
       stepConfiguration(eConfig),
       material(mprop),
       materialPosition(Amg::Vector3D(0., 0., 0.)),
       materialScaling(1.),
       transportJacobian(tjac),
       pathLength(pLength),
       time(0) {}
};
 
template<class T>
class ExtrapolationCell
{
public:
  T* startParameters; 
  const TrackingVolume*
    startVolume; 
  const Layer* startLayer; 
 
  T* endParameters; 
  const TrackingVolume* endVolume; 
  const Layer* endLayer; 
  const Surface*
    endSurface; 
 
  T* leadParameters; 
  const TrackingVolume*
    leadVolume;           
  const Layer* leadLayer; 
  const Surface* leadLayerSurface; 
 
  T* lastBoundaryParameters; 
  const Surface*
    lastBoundarySurface; 
 
  T* lastLeadParameters; 
  PropDirection propDirection; 
  int radialDirection; 
 
  GeometrySignature
    nextGeometrySignature; 
 
  int navigationStep; 
  double pathLength;  
  double pathLimit;   
 
  double materialX0;      
  double materialLimitX0; 
  double materialL0;      
  double materialLimitL0; 
  int materialProcess;    
 
  ParticleHypothesis
    pHypothesis; 
  MagneticFieldProperties
    mFieldMode; 
  MaterialUpdateMode materialUpdateMode; 
  bool navigationCurvilinear;  
  bool sensitiveCurvilinear;   
  bool destinationCurvilinear; 
 
  std::vector<ExtrapolationStep<T>>
    extrapolationSteps; 
 
  ExtrapolationConfig
    extrapolationConfiguration; 
 
  float time;        
  EnergyLoss* eLoss; 
  float zOaTrX;      
  float zX;          
 
  ExtrapolationCell(T& sParameters,
                    PropDirection pDir = alongMomentum,
                    unsigned int econfig = 1)
    : startParameters(&sParameters)
    , startVolume(nullptr)
    , startLayer(nullptr)
    , endParameters(0)
    , endVolume(nullptr)
    , endLayer(nullptr)
    , endSurface(nullptr)
    , leadParameters(&sParameters)
    , leadVolume(nullptr)
    , leadLayer(nullptr)
    , leadLayerSurface(nullptr)
    , lastBoundaryParameters(0)
    , lastBoundarySurface(nullptr)
    , lastLeadParameters(&sParameters)
    , propDirection(pDir)
    , radialDirection(1)
    , nextGeometrySignature(Trk::Unsigned)
    , navigationStep(0)
    , pathLength(0.)
    , pathLimit(-1)
    , materialX0(0.)
    , materialLimitX0(-1.)
    , materialL0(0.)
    , materialLimitL0(-1.)
    , materialProcess(0)
    , pHypothesis(Trk::pion)
    , mFieldMode(Trk::MagneticFieldProperties(Trk::FullField))
    , materialUpdateMode(Trk::addNoise)
    , navigationCurvilinear(true)
    , sensitiveCurvilinear(false)
    , destinationCurvilinear(false)
    , extrapolationConfiguration(econfig)
    , eLoss(nullptr)
    , zOaTrX(0.)
    , zX(0.)
  {}
 
  void addConfigurationMode(ExtrapolationMode::eMode em)
  {
    // set the bit corresponding to this mode
    extrapolationConfiguration.addMode(em);
  }
 
  bool checkConfigurationMode(ExtrapolationMode::eMode em) const
  {
    // check if the bit is set or not
    return extrapolationConfiguration.checkMode(em);
  }
 
  bool onLastBoundary() const
  {
    return (leadParameters == lastBoundaryParameters);
  }
 
  void stepParameters(const T* pars, ExtrapolationMode::eMode fillMode);
 
  void stepTransport(const Surface& sf,
                     double pathLength = 0.,
                     const TransportJacobian* tjac = nullptr);
 
  void addMaterial(double sfactor, const MaterialProperties* mprop = nullptr);
  void addMaterial(double step, const Material* mat = nullptr);
 
  void stepMaterial(const Surface& sf,
                    const Layer* lay,
                    const Amg::Vector3D& position,
                    double sfactor,
                    const MaterialProperties* mprop = nullptr);
 
  bool initialVolume() const { return (leadVolume == startVolume); }
 
  bool finalVolumeReached() const
  {
    return (leadVolume == endVolume && endVolume);
  }
 
  bool pathLimitReached(double tolerance = 0.001) const
  {
    return (checkConfigurationMode(Trk::ExtrapolationMode::StopWithPathLimit) &&
            reachedLimit(pathLength, pathLimit, tolerance));
  }
 
  bool materialLimitReached(double tolerance = 0.001) const
  {
    return ((checkConfigurationMode(
               Trk::ExtrapolationMode::StopWithMaterialLimitX0) &&
             reachedLimit(materialX0, materialLimitX0, tolerance)) ||
            (checkConfigurationMode(
               Trk::ExtrapolationMode::StopWithMaterialLimitL0) &&
             reachedLimit(materialL0, materialLimitL0, tolerance)));
  }
 
  void restartAtDestination();
 
  void finalize(const ExtrapolationCode& ec);
 
  void emptyGarbageBin(const ExtrapolationCode& ec);
 
  void setParticleHypothesis(const ParticleHypothesis& hypo)
  {
    pHypothesis = hypo;
  }
 
  void setRadialDirection()
  {
    // in FATRAS extrapolation mode force radial direction to be outwards (+1)
    if (checkConfigurationMode(ExtrapolationMode::FATRAS))
      radialDirection = 1;
    else {
      // if the endSurface is given, it is used to evaluate the radial direction
      // else the leadParamenters are used
      if (endSurface) {
        if (leadParameters->position().perp() >
            endSurface->globalReferencePoint().perp())
          radialDirection = -1;
      } else {
        if (leadParameters->position().perp() >
            (leadParameters->position() +
             propDirection * leadParameters->momentum().unit())
              .perp())
          radialDirection = -1;
      }
    }
  }
 
  bool checkRadialCompatibility() const
  {
    // this checks the radial compatibility - not needed for outwards moving
    if (radialDirection > 0)
      return true;
    // this was radially inwards moving and stays like this
    if (leadParameters->position().perp() >
        (leadParameters->position() +
         propDirection * leadParameters->momentum().unit())
          .perp())
      return true;
    // radial direction changed
    return false;
  }
 
private:
  std::vector<const T*> m_garbageCollection;
};
 
template<class T>
void
ExtrapolationCell<T>::restartAtDestination()
{
  startParameters = endParameters;
  startVolume = endVolume;
  startLayer = endLayer;
  endParameters = 0;
  endVolume = 0;
  endLayer = 0;
  endSurface = 0;
  leadParameters = startParameters;
  leadVolume = 0;
  leadLayer = 0;
  leadLayerSurface = 0;
  lastBoundaryParameters = 0;
  lastBoundarySurface = 0;
  lastLeadParameters = startParameters;
  navigationStep = 0;
  pathLength = 0.;
  materialX0 = 0.;
  materialL0 = 0.;
  if (eLoss)
    eLoss->set(0., 0., 0., 0., 0., 0.);
  // clear the vector
  extrapolationSteps.clear();
}
 
template<class T>
void
ExtrapolationCell<T>::finalize(const ExtrapolationCode& ec)
{
  // set the leadParameters to the endParameters if anything happened here and
  // the code wass succesful
  if (ec.isSuccessOrRecovered() && leadParameters != startParameters) {
    // end parameters are the last lead parameters !< @TODO check if we need a
    // clone here! should not be necessary
    endParameters = leadParameters->clone();
  }
  // now do the cleanup - will delete the step content if eCode is failure
  emptyGarbageBin(ec);
  delete eLoss;
}
 
template<class T>
void
ExtrapolationCell<T>::emptyGarbageBin(const ExtrapolationCode& ec)
{
  for (auto bC : m_garbageCollection)
    delete bC;
  m_garbageCollection.clear();
  // in case of failure of the extrapolation stream, clear all the caches
  if (ec.isFailure()) {
    for (auto& es : extrapolationSteps) {
      delete es.parameters;
      delete es.transportJacobian;
    }
    // now clear the vector
    extrapolationSteps.clear();
  }
}
 
template<class T>
void
ExtrapolationCell<T>::stepParameters(const T* parameters,
                                     Trk::ExtrapolationMode::eMode fillMode)
{
  // this is the garbage bin collection
  if (!checkConfigurationMode(fillMode)) {
    // only dump them in the garbage bin
    m_garbageCollection.push_back(parameters);
    return;
  }
  // find out if you want to attach or you need a new one
  // current step surface
  const Surface* cssf = &(parameters->associatedSurface());
  // get the last step surface - if it is identical with the current one ->
  // attach information
  const Surface* lssf =
    extrapolationSteps.size()
      ? extrapolationSteps[extrapolationSteps.size() - 1].surface
      : 0;
  // create a new step
  if (cssf != lssf)
    extrapolationSteps.push_back(ExtrapolationStep<T>());
  // fill the parameters, the surface and add the mode
  extrapolationSteps[extrapolationSteps.size() - 1].parameters = parameters;
  extrapolationSteps[extrapolationSteps.size() - 1].surface = cssf;
  extrapolationSteps[extrapolationSteps.size() - 1].stepConfiguration.addMode(
    fillMode);
}
 
template<class T>
void
ExtrapolationCell<T>::stepTransport(const Surface& sf,
                                    double pLength,
                                    const TransportJacobian* tjac)
{
  // find out if you want to attach or you need a new one
  // current step surface
  const Surface* cssf = &sf;
  // get the last step surface - if it is identical with the current one ->
  // attach information
  const Surface* lssf =
    extrapolationSteps.size()
      ? extrapolationSteps[extrapolationSteps.size() - 1].surface
      : 0;
  // only create a new step for a transport jacobian
  if (tjac) {
    // create a new step
    if (cssf != lssf)
      extrapolationSteps.push_back(ExtrapolationStep<T>());
    // set the surface
    extrapolationSteps[extrapolationSteps.size() - 1].surface = cssf;
    // set the the transport information
    extrapolationSteps[extrapolationSteps.size() - 1].transportJacobian = tjac;
    extrapolationSteps[extrapolationSteps.size() - 1].stepConfiguration.addMode(
      Trk::ExtrapolationMode::CollectJacobians);
    // fill the step path length
    if (pLength > 0.) {
      extrapolationSteps[extrapolationSteps.size() - 1].pathLength = pLength;
      extrapolationSteps[extrapolationSteps.size() - 1]
        .stepConfiguration.addMode(Trk::ExtrapolationMode::CollectPathSteps);
    }
  } else {
    // let's just fill the pathLength information
    pathLength += pLength;
  }
}
 
template<class T>
void
ExtrapolationCell<T>::addMaterial(double sfactor,
                                  const MaterialProperties* mprop)
{
 
  // fill the material if there
  if (mprop) {
    // the overal material
    materialX0 += sfactor * mprop->thicknessInX0();
    materialL0 += sfactor * mprop->thicknessInL0();
    zOaTrX += mprop->zOverAtimesRho() * sfactor * mprop->thicknessInX0();
    zX += mprop->averageZ() * sfactor * mprop->thicknessInX0();
  }
}
 
template<class T>
void
ExtrapolationCell<T>::addMaterial(double step, const Material* mat)
{
 
  // fill the material if there
  if (mat && step > 0.) {
    // the overal material
    materialX0 += step / mat->X0;
    materialL0 += step / mat->L0;
    zOaTrX += mat->zOverAtimesRho() * step / mat->X0;
    zX += mat->averageZ() * step / mat->X0;
  }
}
 
template<class T>
void
ExtrapolationCell<T>::stepMaterial(const Surface& sf,
                                   const Layer* lay,
                                   const Amg::Vector3D& mposition,
                                   double sfactor,
                                   const MaterialProperties* mprop)
{
 
  // add material to the global counter
  addMaterial(sfactor, mprop);
  // find out if you want to attach or you need a new one
  // current step surface
  const Surface* cssf = &sf;
  // get the last step surface - if it is identical with the current one ->
  // attach information
  const Surface* lssf =
    extrapolationSteps.size()
      ? extrapolationSteps[extrapolationSteps.size() - 1].surface
      : 0;
  // create a new step
  if (cssf != lssf)
    extrapolationSteps.push_back(ExtrapolationStep<T>());
  // set the surface
  extrapolationSteps[extrapolationSteps.size() - 1].surface = cssf;
  extrapolationSteps[extrapolationSteps.size() - 1].layer = lay;
  // fill the material if there
  if (mprop) {
    // record the step information
    extrapolationSteps[extrapolationSteps.size() - 1].material = mprop;
    extrapolationSteps[extrapolationSteps.size() - 1].stepConfiguration.addMode(
      Trk::ExtrapolationMode::CollectMaterial);
    extrapolationSteps[extrapolationSteps.size() - 1].materialPosition =
      mposition;
    extrapolationSteps[extrapolationSteps.size() - 1].materialScaling = sfactor;
  }
}
 
} // end of namespace
 
#endif // TRKEXUTILS_EXTRAPOLATIONCELL_H