SteppingHelixPropagator/interface/SteppingHelixPropagator.h

0001 #ifndef SteppingHelixPropagator_SteppingHelixPropagator_h
0002 #define SteppingHelixPropagator_SteppingHelixPropagator_h
0003
0004 /** \class SteppingHelixPropagator
0005  *  Propagator implementation using steps along a helix.
0006  *  Minimal geometry navigation.
0007  *  Material effects (multiple scattering and energy loss) are based on tuning
0008  *  to MC and (eventually) data.
0009  *
0010  *  \author Vyacheslav Krutelyov (slava77)
0011  */
0012
0013 //
0014 // Original Author:  Vyacheslav Krutelyov
0015 //         Created:  Fri Mar  3 16:01:24 CST 2006
0016 //
0017 //
0018
0019 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0020 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0021 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0022 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0023 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0024 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0025
0026 #include "CLHEP/Matrix/SymMatrix.h"
0027 #include "CLHEP/Matrix/Matrix.h"
0028 #include "CLHEP/Vector/ThreeVector.h"
0029
0030 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixStateInfo.h"
0031
0032 class MagneticField;
0033 class VolumeBasedMagneticField;
0034 class MagVolume;
0035
0036 class SteppingHelixPropagator final : public Propagator {
0037 public:
0038   typedef CLHEP::Hep3Vector Vector;
0039   typedef CLHEP::Hep3Vector Point;
0040
0041   typedef SteppingHelixStateInfo StateInfo;
0042   typedef SteppingHelixStateInfo::Result Result;
0043
0044   struct Basis {
0045     Vector lX;
0046     Vector lY;
0047     Vector lZ;
0048   };
0049
0050   enum Pars { RADIUS_P = 0, Z_P = 0, PATHL_P = 0 };
0051
0052   enum DestType {
0053     RADIUS_DT = 0,
0054     Z_DT,
0055     PLANE_DT,
0056     CONE_DT,
0057     CYLINDER_DT,
0058     PATHL_DT,
0059     POINT_PCA_DT,
0060     LINE_PCA_DT,
0061     UNDEFINED_DT
0062   };
0063
0064   enum Fancy {
0065     HEL_AS_F = 0,  //simple analytical helix, eloss at end of step
0066     HEL_ALL_F,     //analytical helix with linear eloss
0067     POL_1_F,       //1st order approximation, straight line
0068     POL_2_F,       //2nd order
0069     POL_M_F        //highest available
0070   };
0071
0072   //! Constructors
0073   SteppingHelixPropagator();
0074   SteppingHelixPropagator(const MagneticField* field, PropagationDirection dir = alongMomentum);
0075
0076   SteppingHelixPropagator* clone() const override { return new SteppingHelixPropagator(*this); }
0077
0078   //! Destructor
0079   ~SteppingHelixPropagator() override;
0080
0081   const MagneticField* magneticField() const override { return field_; }
0082
0083   using Propagator::propagate;
0084   using Propagator::propagateWithPath;
0085
0086   //! Propagate to Plane given a starting point: return final
0087   //! TrajectoryState and path length from start to this point
0088   std::pair<TrajectoryStateOnSurface, double> propagateWithPath(const FreeTrajectoryState& ftsStart,
0089                                                                 const Plane& pDest) const override;
0090   //! Propagate to Cylinder given a starting point: return final TrajectoryState
0091   //!and path length from start to this point
0092   std::pair<TrajectoryStateOnSurface, double> propagateWithPath(const FreeTrajectoryState& ftsStart,
0093                                                                 const Cylinder& cDest) const override;
0094   //! Propagate to PCA to point given a starting point
0095   std::pair<FreeTrajectoryState, double> propagateWithPath(const FreeTrajectoryState& ftsStart,
0096                                                            const GlobalPoint& pDest) const override;
0097   //! Propagate to PCA to a line (given by 2 points) given a starting point
0098   std::pair<FreeTrajectoryState, double> propagateWithPath(const FreeTrajectoryState& ftsStart,
0099                                                            const GlobalPoint& pDest1,
0100                                                            const GlobalPoint& pDest2) const override;
0101   //! Propagate to PCA to a line (given by beamSpot position and slope) given a starting point
0102   std::pair<FreeTrajectoryState, double> propagateWithPath(const FreeTrajectoryState& ftsStart,
0103                                                            const reco::BeamSpot& beamSpot) const override;
0104
0105   //----------------------------------------------------------------------------------------------------------
0106
0107   //! Propagate to Plane given a starting point
0108   void propagate(const SteppingHelixStateInfo& ftsStart, const Surface& sDest, SteppingHelixStateInfo& out) const;
0109   void propagate(const SteppingHelixStateInfo& ftsStart, const Plane& pDest, SteppingHelixStateInfo& out) const;
0110   //! Propagate to Cylinder given a starting point (a Cylinder is assumed to be positioned at 0,0,0)
0111   void propagate(const SteppingHelixStateInfo& ftsStart, const Cylinder& cDest, SteppingHelixStateInfo& out) const;
0112   //! Propagate to PCA to point given a starting point
0113   void propagate(const SteppingHelixStateInfo& ftsStart, const GlobalPoint& pDest, SteppingHelixStateInfo& out) const;
0114   //! Propagate to PCA to a line (given by 2 points) given a starting point
0115   void propagate(const SteppingHelixStateInfo& ftsStart,
0116                  const GlobalPoint& pDest1,
0117                  const GlobalPoint& pDest2,
0118                  SteppingHelixStateInfo& out) const;
0119
0120   //! Switch debug printouts (to cout) .. very verbose
0121   void setDebug(bool debug) { debug_ = debug; }
0122
0123   //! Switch for material effects mode: no material effects if true
0124   void setMaterialMode(bool noMaterial) { noMaterialMode_ = noMaterial; }
0125
0126   //! Force no error propagation
0127   void setNoErrorPropagation(bool noErrorPropagation) { noErrorPropagation_ = noErrorPropagation; }
0128
0129   //! Apply radLength correction (1+0.036*ln(radX0+1)) to covariance matrix
0130   //! +1 is added for IR-safety
0131   //! Should be done with care: it's easy to make the end-point result dependent
0132   //! on the intermediate stop points
0133   void applyRadX0Correction(bool applyRadX0Correction) { applyRadX0Correction_ = applyRadX0Correction; }
0134
0135   //!Switch to using MagneticField Volumes .. as in VolumeBasedMagneticField
0136   void setUseMagVolumes(bool val) { useMagVolumes_ = val; }
0137
0138   //!Switch to using Material Volumes .. internally defined for now
0139   void setUseMatVolumes(bool val) { useMatVolumes_ = val; }
0140
0141   //! flag to return tangent plane for non-plane input
0142   void setReturnTangentPlane(bool val) { returnTangentPlane_ = val; }
0143
0144   //! flag to send LogWarning on failures
0145   void setSendLogWarning(bool val) { sendLogWarning_ = val; }
0146
0147   //! (temporary?) flag to identify if the volume is yoke based on MagVolume internals
0148   //! works if useMatVolumes_ is also true
0149   void setUseIsYokeFlag(bool val) { useIsYokeFlag_ = val; }
0150
0151   //! will use bigger steps ~tuned for good-enough L2/Standalone reco
0152   //! use this in hope of a speed-up
0153   void setUseTuningForL2Speed(bool val) { useTuningForL2Speed_ = val; }
0154
0155   //! set VolumBasedField pointer
0156   //! allows to have geometry description in uniformField scenario
0157   //! only important/relevant in the barrel region
0158   void setVBFPointer(const VolumeBasedMagneticField* val) { vbField_ = val; }
0159
0160   //! force getting field value from MagneticField, not the geometric one
0161   void setUseInTeslaFromMagField(bool val) { useInTeslaFromMagField_ = val; }
0162
0163   //! set shifts in Z for endcap pieces (includes EE, HE, ME, YE)
0164   void setEndcapShiftsInZPosNeg(double valPos, double valNeg) {
0165     ecShiftPos_ = valPos;
0166     ecShiftNeg_ = valNeg;
0167   }
0168
0169 protected:
0170   typedef SteppingHelixStateInfo::VolumeBounds MatBounds;
0171   static const int MAX_POINTS = 7;
0172   typedef StateInfo StateArray[MAX_POINTS + 1];
0173
0174   //! (Internals) set defaults needed for states used inside propagate methods
0175   void initStateArraySHPSpecific(StateArray& svBuf, bool flagsOnly) const;
0176
0177   //! (Internals) Init starting point
0178   void setIState(const SteppingHelixStateInfo& sStart, StateArray& svBuff, int& nPoints) const;
0179
0180   //! propagate: chose stop point by type argument
0181   //! propagate to fixed radius [ r = sqrt(x**2+y**2) ] with precision epsilon
0182   //! propagate to plane by [x0,y0,z0, n_x, n_y, n_z] parameters
0183   Result propagate(StateArray& svBuff,
0184                    int& nPoints,
0185                    SteppingHelixPropagator::DestType type,
0186                    const double pars[6],
0187                    double epsilon = 1e-3) const;
0188
0189   //! (Internals) compute transient values for initial point (resets step counter).
0190   //!  Called by setIState
0191   void loadState(SteppingHelixPropagator::StateInfo& svCurrent,
0192                  const SteppingHelixPropagator::Vector& p3,
0193                  const SteppingHelixPropagator::Point& r3,
0194                  int charge,
0195                  PropagationDirection dir,
0196                  const AlgebraicSymMatrix55& covCurv) const;
0197
0198   //! (Internals) compute transients for current point (increments step counter).
0199   //!  Called by makeAtomStep
0200   void getNextState(const SteppingHelixPropagator::StateInfo& svPrevious,
0201                     SteppingHelixPropagator::StateInfo& svNext,
0202                     double dP,
0203                     const SteppingHelixPropagator::Vector& tau,
0204                     const SteppingHelixPropagator::Vector& drVec,
0205                     double dS,
0206                     double dX0,
0207                     const AlgebraicMatrix55& dCovCurv) const;
0208
0209   //! Set/compute basis vectors for local coordinates at current step (called by incrementState)
0210   void setRep(SteppingHelixPropagator::Basis& rep, const SteppingHelixPropagator::Vector& tau) const;
0211
0212   //! main stepping function: compute next state vector after a step of length dS
0213   bool makeAtomStep(SteppingHelixPropagator::StateInfo& svCurrent,
0214                     SteppingHelixPropagator::StateInfo& svNext,
0215                     double dS,
0216                     PropagationDirection dir,
0217                     SteppingHelixPropagator::Fancy fancy) const;
0218
0219   //! estimate average (in fact smth. close to MPV and median) energy loss per unit path length
0220   double getDeDx(const SteppingHelixPropagator::StateInfo& sv,
0221                  double& dEdXPrime,
0222                  double& radX0,
0223                  MatBounds& rzLims) const;
0224
0225   //! (Internals) circular index for array of transients
0226   int cIndex_(int ind) const;
0227
0228   //! (Internals) determine distance and direction from the current position to the plane
0229   Result refToDest(DestType dest,
0230                    const SteppingHelixPropagator::StateInfo& sv,
0231                    const double pars[6],
0232                    double& dist,
0233                    double& tanDist,
0234                    PropagationDirection& refDirection,
0235                    double fastSkipDist = 1e12) const;
0236
0237   //! (Internals) determine distance and direction from the current position to the
0238   //! boundary of current mag volume
0239   Result refToMagVolume(const SteppingHelixPropagator::StateInfo& sv,
0240                         PropagationDirection dir,
0241                         double& dist,
0242                         double& tanDist,
0243                         double fastSkipDist = 1e12,
0244                         bool expectNewMagVolume = false,
0245                         double maxStep = 1e12) const;
0246
0247   Result refToMatVolume(const SteppingHelixPropagator::StateInfo& sv,
0248                         PropagationDirection dir,
0249                         double& dist,
0250                         double& tanDist,
0251                         double fastSkipDist = 1e12) const;
0252
0253   //! check if it's a yoke/iron based on this MagVol internals
0254   bool isYokeVolume(const MagVolume* vol) const;
0255
0256 private:
0257   typedef std::pair<TrajectoryStateOnSurface, double> TsosPP;
0258   typedef std::pair<FreeTrajectoryState, double> FtsPP;
0259   static const int MAX_STEPS = 10000;
0260   //mutable int nPoints_;
0261   //mutable StateInfo svBuf_[MAX_POINTS+1];
0262
0263   StateInfo invalidState_;
0264
0265   const MagneticField* field_;
0266   const VolumeBasedMagneticField* vbField_;
0267   const AlgebraicSymMatrix55 unit55_;
0268   bool debug_;
0269   bool noMaterialMode_;
0270   bool noErrorPropagation_;
0271   bool applyRadX0Correction_;
0272   bool useMagVolumes_;
0273   bool useIsYokeFlag_;
0274   bool useMatVolumes_;
0275   bool useInTeslaFromMagField_;
0276   bool returnTangentPlane_;
0277   bool sendLogWarning_;
0278   bool useTuningForL2Speed_;
0279
0280   double defaultStep_;
0281
0282   double ecShiftPos_;
0283   double ecShiftNeg_;
0284 };
0285
0286 #endif