STEP_Propagator.h

Go to the documentation of this file.

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

//  Header file for class  STEP_Propagator
// (c) ATLAS Detector software
// Class for track status propagation through magnetic field
// Based on RungeKuttaPropagator Version 1.0 21/05/2004 I.Gavrilenko
// Version 0.1 16/2/2005 Esben Lund
// Extended by S.Todorova for propagation to multiple surfaces
 
#ifndef STEP_Propagator_H
#define STEP_Propagator_H
 
#include "AthenaBaseComps/AthAlgTool.h"
#include "AthenaKernel/IAthRNGSvc.h"
#include "GaudiKernel/ToolHandle.h"
#include "TrkEventPrimitives/PropDirection.h"
#include "TrkEventPrimitives/SurfaceTypes.h"
#include "TrkExInterfaces/IPropagator.h"
//
#include "TrkExUtils/MaterialInteraction.h"
#include "TrkExUtils/TrackSurfaceIntersection.h"
#include "TrkGeometry/BinnedMaterial.h" //Identified material typedef
#include "TrkMaterialOnTrack/EnergyLoss.h"
#include "TrkParameters/TrackParameters.h" //TrackParameters typedef
//#include "TrkExUtils/ExtrapolationCache.h"
// Amg
#include "EventPrimitives/EventPrimitives.h"
#include "GeoPrimitives/GeoPrimitives.h"
// MagField cache
#include "MagFieldConditions/AtlasFieldCacheCondObj.h"
#include "MagFieldElements/AtlasFieldCache.h"
//
#include "CxxUtils/restrict.h"
#include <list>
#include <vector>
 
namespace CLHEP{
  class HepRandomEngine;
}
 
#include "TrkExInterfaces/ITimedMatEffUpdator.h"
namespace Trk {
class Surface;
class TrackingVolume;
class ScatteringAngles;
class AlignableTrackingVolume;
class ExtrapolationCache;
 
class STEP_Propagator final
  : public AthAlgTool
  , virtual public IPropagator
{
  // Public methods:
public:
  using IPropagator::propagate;
 
  STEP_Propagator(const std::string&, const std::string&, const IInterface*);
 
  virtual ~STEP_Propagator();

  virtual StatusCode initialize() override final;

  virtual StatusCode finalize() override final;

  virtual std::unique_ptr<Trk::NeutralParameters> propagate(
    const Trk::NeutralParameters&,
    const Trk::Surface&,
    Trk::PropDirection,
    const Trk::BoundaryCheck&,
    bool rC = false) const override final;

  virtual std::unique_ptr<Trk::TrackParameters> propagate(
    const EventContext& ctx,
    const Trk::TrackParameters& trackParameters,
    const Trk::Surface& targetSurface,
    Trk::PropDirection propagationDirection,
    const Trk::BoundaryCheck& boundaryCheck,
    const MagneticFieldProperties& magneticFieldProperties,
    ParticleHypothesis particle,
    bool returnCurv = false,
    const Trk::TrackingVolume* tVol = nullptr) const override final;

  virtual std::unique_ptr<Trk::TrackParameters> propagate(
    const EventContext& ctx,
    const Trk::TrackParameters& trackParameters,
    std::vector<Trk::DestSurf>& targetSurfaces,
    Trk::PropDirection propagationDirection,
    const MagneticFieldProperties& magneticFieldProperties,
    ParticleHypothesis particle,
    std::vector<unsigned int>& solutions,
    double& path,
    bool usePathLimit = false,
    bool returnCurv = false,
    const Trk::TrackingVolume* tVol = nullptr) const override final;

  virtual std::unique_ptr<Trk::TrackParameters> propagateT(
    const EventContext& ctx,
    const Trk::TrackParameters& trackParameters,
    std::vector<Trk::DestSurf>& targetSurfaces,
    Trk::PropDirection propagationDirection,
    const MagneticFieldProperties& magneticFieldProperties,
    ParticleHypothesis particle,
    std::vector<unsigned int>& solutions,
    Trk::PathLimit& path,
    Trk::TimeLimit& time,
    bool returnCurv,
    const Trk::TrackingVolume* tVol,
    std::vector<Trk::HitInfo>*& hitVector) const override final;

  virtual std::unique_ptr<Trk::TrackParameters> propagateM(
    const EventContext& ctx,
    const Trk::TrackParameters& trackParameters,
    std::vector<Trk::DestSurf>& targetSurfaces,
    Trk::PropDirection propagationDirection,
    const MagneticFieldProperties& magneticFieldProperties,
    ParticleHypothesis particle,
    std::vector<unsigned int>& solutions,
    std::vector<const Trk::TrackStateOnSurface*>* matstates,
    std::vector<std::pair<std::unique_ptr<Trk::TrackParameters>, int>>*
      intersections,
    double& path,
    bool usePathLimit = false,
    bool returnCurv = false,
    const Trk::TrackingVolume* tVol = nullptr,
    Trk::ExtrapolationCache* = nullptr) const override final;

  virtual std::unique_ptr<Trk::TrackParameters> propagate(
    const EventContext& ctx,
    const Trk::TrackParameters& trackParameters,
    const Trk::Surface& targetSurface,
    Trk::PropDirection propagationDirection,
    const Trk::BoundaryCheck& boundaryCheck,
    const MagneticFieldProperties& magneticFieldProperties,
    std::optional<Trk::TransportJacobian>& jacobian,
    double& pathLimit,
    ParticleHypothesis particle,
    bool returnCurv = false,
    const Trk::TrackingVolume* tVol = nullptr) const override final;

  virtual std::unique_ptr<Trk::TrackParameters> propagateParameters(
    const EventContext& ctx,
    const Trk::TrackParameters& trackParameters,
    const Trk::Surface& targetSurface,
    Trk::PropDirection propagationDirection,
    const Trk::BoundaryCheck& boundaryCheck,
    const MagneticFieldProperties& magneticFieldProperties,
    ParticleHypothesis particle,
    bool returnCurv = false,
    const Trk::TrackingVolume* tVol = nullptr) const override final;

  virtual std::unique_ptr<Trk::TrackParameters> propagateParameters(
    const EventContext& ctx,
    const Trk::TrackParameters& trackParameters,
    const Trk::Surface& targetSurface,
    Trk::PropDirection propagationDirection,
    const Trk::BoundaryCheck& boundaryCheck,
    const MagneticFieldProperties& magneticFieldProperties,
    std::optional<Trk::TransportJacobian>& jacobian,
    ParticleHypothesis particle,
    bool returnCurv = false,
    const Trk::TrackingVolume* tVol = nullptr) const override final;

  virtual std::optional<TrackSurfaceIntersection> intersect(
    const EventContext& ctx,
    const Trk::TrackParameters& trackParameters,
    const Trk::Surface& targetSurface,
    const Trk::MagneticFieldProperties& magneticFieldProperties,
    ParticleHypothesis particle,
    const Trk::TrackingVolume* tVol = nullptr) const override final;

  virtual std::optional<TrackSurfaceIntersection> intersectSurface(
    const EventContext& ctx,
    const Surface& surface,
    const TrackSurfaceIntersection& trackIntersection,
    const double qOverP,
    const MagneticFieldProperties& mft,
    ParticleHypothesis particle) const override final;

  virtual void globalPositions(
    const EventContext& ctx,
    std::deque<Amg::Vector3D>& positionsList,
    const TrackParameters& trackParameters,
    const MagneticFieldProperties& magneticFieldProperties,
    const CylinderBounds& cylinderBounds,
    double maxStepSize,
    ParticleHypothesis particle,
    const Trk::TrackingVolume* tVol = 0) const override final;

  virtual Trk::MultiComponentState multiStatePropagate(
    const EventContext&,
    const MultiComponentState&,
    const Surface&,
    const MagneticFieldProperties&,
    const PropDirection,
    const BoundaryCheck&,
    const ParticleHypothesis) const override final
  {
    ATH_MSG_ERROR("Call to non-implemented multiStatePropagate");
    return {};
  }

  // Private methods:

  struct Cache
  {
    bool m_energyLoss{ true };
    bool m_detailedElossFlag{ true };
    bool m_straggling{ true };
    bool m_solenoid{ false };
    bool m_matPropOK{ true };
    bool m_brem{ false };
    bool m_includeBgradients{ true };
    bool m_includeGgradient{ false };
    bool m_MPV{ false };
    bool m_multipleScattering{ true };
    int m_propagateWithPathLimit{};
    size_t m_currentLayerBin{};
    double m_matupd_lastmom{};
    double m_matupd_lastpath{};
    double m_matdump_lastpath{};
    double m_delRad{ 0 };    // deRad/dl;
    double m_delIoni{ 0 };   // deIoni/dl;
    double m_sigmaIoni{ 0 }; // dsigma(ioni)/dl;
    double m_kazL{ 0 };      // kazL constant;
    double m_sigmaRad{ 0 };  // dsigma(rad)/dl;
    // cache for input variance
    double m_inputThetaVariance{};
    double m_stragglingVariance{};
    double m_pathLimit{};
    double m_timeOfFlight{};
    double m_timeStep{};
    double m_particleMass{ 0 }; 
    double m_charge{};
    double m_combinedThickness{};
    double m_tolerance{ 1e-05 };
    double m_momentumCutOff{ 50. };
    double m_scatteringScale{ 1. };
    double m_maxPath{ 100000. };
    int m_maxSteps{ 10000 };
    double m_layXmax{ 1. };
    // secondary interactions
    double m_timeIn{};
    double m_bremMom{ 0. };
    double m_bremEmitThreshold{ 0. };
    double m_bremSampleThreshold{ 0. };
    double m_P[45]{};
 
    const Trk::BinnedMaterial* m_binMat{ nullptr };
    std::vector<const Trk::TrackStateOnSurface*>* m_matstates{ nullptr };
    std::vector<std::pair<std::unique_ptr<Trk::TrackParameters>, int>>*
      m_identifiedParameters{ nullptr };

    std::vector<Trk::HitInfo>* m_hitVector{ nullptr };
 
    ParticleHypothesis m_particle{};
    const TrackingVolume* m_trackingVolume{ nullptr };
    const Material* m_material{ nullptr };
    Trk::ExtrapolationCache* m_extrapolationCache{ nullptr };
    // cache for combined covariance matrix contribution
    AmgSymMatrix(5) m_combinedCovariance;
    // cache for differential covariance matrix contribution ( partial material
    // dump )
    AmgSymMatrix(5) m_covariance;
    Trk::EnergyLoss m_combinedEloss;
    std::vector<std::pair<int, std::pair<double, double>>> m_currentDist;
    MagField::AtlasFieldCache m_fieldCache;
    CLHEP::HepRandomEngine* m_randomEngine { nullptr };
    const EventContext& m_ctx;
 
    Cache (const EventContext& ctx) : m_ctx (ctx) {}
  };
 
private:
  // the configured properties.
  void setCacheFromProperties(Cache& cache) const
  {
    cache.m_includeBgradients = m_includeBgradients;
    cache.m_includeGgradient = m_includeGgradient;
    cache.m_energyLoss = m_energyLoss;
    cache.m_detailedElossFlag = m_detailedEloss;
    cache.m_MPV = m_MPV;
    cache.m_multipleScattering = m_multipleScattering;
    cache.m_straggling = m_straggling;
    cache.m_tolerance = m_tolerance;
    cache.m_momentumCutOff = m_momentumCutOff;
    cache.m_scatteringScale = m_scatteringScale;
    cache.m_maxPath = m_maxPath;
    cache.m_maxSteps = m_maxSteps;
    cache.m_layXmax = m_layXmax;
  }

  // Main function for propagation
  // with search of closest surface (ST)
  std::unique_ptr<Trk::TrackParameters> propagateRungeKutta(
    Cache& cache,
    bool errorPropagation,
    const Trk::TrackParameters& trackParameters,
    std::vector<DestSurf>& targetSurfaces,
    Trk::PropDirection propagationDirection,
    const MagneticFieldProperties& magneticFieldProperties,
    ParticleHypothesis particle,
    std::vector<unsigned int>& solutions,
    double& path,
    bool returnCurv = false) const;

  // Method for propagation with search of closest surface (ST)
  bool propagateWithJacobian(Cache& cache,
                             bool errorPropagation,
                             std::vector<DestSurf>& sfs,
                             double* P,
                             Trk::PropDirection propDir,
                             std::vector<unsigned int>& solutions,
                             double& path,
                             double sumPath) const;

  // dump material effects
  void dumpMaterialEffects(Cache& cache,
                           const Trk::CurvilinearParameters* trackParameters,
                           double path) const;

  // Momentum smearing (simulation mode)
  void smear(Cache& cache,
             double& phi,
             double& theta,
             const Trk::TrackParameters* parms,
             double radDist) const;
  // Bremstrahlung (simulation mode)
  void sampleBrem(Cache& cache, double mom) const;
 
  DoubleProperty m_tolerance{this, "Tolerance", 1e-05,
    "Error tolerance. Low tolerance gives igh accuracy"};
  BooleanProperty m_materialEffects{this, "MaterialEffects", true,
    "Switch material effects on or off"};
  BooleanProperty m_includeBgradients{this, "IncludeBgradients", true,
    "Include B-field gradients in the error propagation"};
  BooleanProperty m_includeGgradient{this, "IncludeGgradient", false,
    "Include dg/dlambda into the error propagation? Only relevant when energy loss is true."};
  DoubleProperty m_momentumCutOff{this, "MomentumCutOff", 50.,
    "Stop propagation below this momentum in MeV"};
  BooleanProperty m_multipleScattering{this, "MultipleScattering", true,
    "Add multiple scattering to the covariance matrix?"};
  BooleanProperty m_energyLoss{this, "EnergyLoss", true, "Include energy loss?"};
  BooleanProperty m_detailedEloss{this, "DetailedEloss", true,
    "Provide the extended EnergyLoss object with MopIonization etc."};
  BooleanProperty m_straggling{this, "Straggling", true,
    "Add energy loss fluctuations (straggling) to the covariance matrix?"};
  BooleanProperty m_MPV{this, "MostProbableEnergyLoss", false,
    "Use the most probable value of the energy loss, else use the mean energy loss."};
  DoubleProperty m_stragglingScale{this, "StragglingScale", 1.,
    "Scale for adjusting the width of the energy loss fluctuations."};
  DoubleProperty m_scatteringScale{this, "MultipleScatteringScale", 1.,
    "Scale for adjusting the multiple scattering contribution to the covariance matrix."};
  DoubleProperty m_maxPath{this, "MaxPath", 100000.,
    "Maximum propagation length in mm."};
  IntegerProperty m_maxSteps{this, "MaxSteps", 10000,
    "Maximum number of allowed steps (to avoid infinite loops)."};
  DoubleProperty m_layXmax{this, "MSstepMax", 1.,
    "maximal layer thickness for multiple scattering calculations"};
 
  // The following are needed for simulation
  // enable simulation mode
  BooleanProperty m_simulation{this, "SimulationMode", false,
    "flag for simulation mode"};

  ToolHandle<ITimedMatEffUpdator> m_simMatUpdator{this, "SimMatEffUpdator", ""};
  ServiceHandle<IAthRNGSvc> m_rndGenSvc{this, "RandomNumberService", "AthRNGSvc",
    "Random number generator"};

  ATHRNG::RNGWrapper* m_rngWrapper = nullptr;
  StringProperty m_randomEngineName{this, "RandomStreamName", "FatrasRnd",
    "Name of the random number stream"};
 
  // Read handle for conditions object to get the field cache
  SG::ReadCondHandleKey<AtlasFieldCacheCondObj> m_fieldCacheCondObjInputKey{
    this,
    "AtlasFieldCacheCondObj",
    "fieldCondObj",
    "Name of the Magnetic Field conditions object key"
  };
  void getFieldCacheObject(Cache& cache, const EventContext& ctx) const;
 
  CLHEP::HepRandomEngine* getRandomEngine (const EventContext& ctx) const;
};
 
} // end of namespace Trk
 
#endif // STEP_Propagator_H