Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:45 

0001 /** \class SteppingHelixPropagator
0002  *  Propagator implementation using steps along a helix.
0003  *  Minimal geometry navigation.
0004  *  Material effects (multiple scattering and energy loss) are based on tuning
0005  *  to MC and (eventually) data. 
0006  *  Implementation file contents follow.
0007  *
0008  *  \author Vyacheslav Krutelyov (slava77)
0009  */
0010 
0011 //
0012 // Original Author:  Vyacheslav Krutelyov
0013 //         Created:  Fri Mar  3 16:01:24 CST 2006
0014 //
0015 //
0016 
0017 #include "MagneticField/Engine/interface/MagneticField.h"
0018 #include "MagneticField/VolumeBasedEngine/interface/VolumeBasedMagneticField.h"
0019 #include "MagneticField/VolumeGeometry/interface/MagVolume.h"
0020 #include "MagneticField/Interpolation/interface/MFGrid.h"
0021 
0022 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0023 #include "DataFormats/GeometrySurface/interface/Plane.h"
0024 #include "DataFormats/GeometrySurface/interface/Cone.h"
0025 
0026 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
0027 
0028 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToCurvilinear.h"
0029 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToCartesian.h"
0030 
0031 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0032 
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034 
0035 #include <sstream>
0036 #include <typeinfo>
0037 
0038 void SteppingHelixPropagator::initStateArraySHPSpecific(StateArray& svBuf, bool flagsOnly) const {
0039   for (int i = 0; i <= MAX_POINTS; i++) {
0040     svBuf[i].isComplete = true;
0041     svBuf[i].isValid_ = true;
0042     svBuf[i].hasErrorPropagated_ = !noErrorPropagation_;
0043     if (!flagsOnly) {
0044       svBuf[i].p3 = Vector(0, 0, 0);
0045       svBuf[i].r3 = Point(0, 0, 0);
0046       svBuf[i].bf = Vector(0, 0, 0);
0047       svBuf[i].bfGradLoc = Vector(0, 0, 0);
0048       svBuf[i].covCurv = AlgebraicSymMatrix55();
0049       svBuf[i].matDCovCurv = AlgebraicSymMatrix55();
0050     }
0051   }
0052 }
0053 
0054 SteppingHelixPropagator::~SteppingHelixPropagator() {}
0055 
0056 SteppingHelixPropagator::SteppingHelixPropagator() : Propagator(anyDirection) { field_ = nullptr; }
0057 
0058 SteppingHelixPropagator::SteppingHelixPropagator(const MagneticField* field, PropagationDirection dir)
0059     : Propagator(dir), unit55_(AlgebraicMatrixID()) {
0060   field_ = field;
0061   vbField_ = dynamic_cast<const VolumeBasedMagneticField*>(field_);
0062   debug_ = false;
0063   noMaterialMode_ = false;
0064   noErrorPropagation_ = false;
0065   applyRadX0Correction_ = true;
0066   useMagVolumes_ = true;
0067   useIsYokeFlag_ = true;
0068   useMatVolumes_ = true;
0069   useInTeslaFromMagField_ = false;  //overrides behavior only if true
0070   returnTangentPlane_ = true;
0071   sendLogWarning_ = false;
0072   useTuningForL2Speed_ = false;
0073   defaultStep_ = 5.;
0074 
0075   ecShiftPos_ = 0;
0076   ecShiftNeg_ = 0;
0077 }
0078 
0079 std::pair<TrajectoryStateOnSurface, double> SteppingHelixPropagator::propagateWithPath(
0080     const FreeTrajectoryState& ftsStart, const Plane& pDest) const {
0081   StateArray svBuf;
0082   initStateArraySHPSpecific(svBuf, true);
0083   int nPoints = 0;
0084   setIState(SteppingHelixStateInfo(ftsStart), svBuf, nPoints);
0085 
0086   StateInfo svCurrent;
0087   propagate(svBuf[0], pDest, svCurrent);
0088 
0089   return TsosPP(svCurrent.getStateOnSurface(pDest), svCurrent.path());
0090 }
0091 
0092 std::pair<TrajectoryStateOnSurface, double> SteppingHelixPropagator::propagateWithPath(
0093     const FreeTrajectoryState& ftsStart, const Cylinder& cDest) const {
0094   StateArray svBuf;
0095   initStateArraySHPSpecific(svBuf, true);
0096   int nPoints = 0;
0097   setIState(SteppingHelixStateInfo(ftsStart), svBuf, nPoints);
0098 
0099   StateInfo svCurrent;
0100   propagate(svBuf[0], cDest, svCurrent);
0101 
0102   return TsosPP(svCurrent.getStateOnSurface(cDest, returnTangentPlane_), svCurrent.path());
0103 }
0104 
0105 std::pair<FreeTrajectoryState, double> SteppingHelixPropagator::propagateWithPath(const FreeTrajectoryState& ftsStart,
0106                                                                                   const GlobalPoint& pDest) const {
0107   StateArray svBuf;
0108   initStateArraySHPSpecific(svBuf, true);
0109   int nPoints = 0;
0110   setIState(SteppingHelixStateInfo(ftsStart), svBuf, nPoints);
0111 
0112   StateInfo svCurrent;
0113   propagate(svBuf[0], pDest, svCurrent);
0114 
0115   FreeTrajectoryState ftsDest;
0116   svCurrent.getFreeState(ftsDest);
0117 
0118   return FtsPP(ftsDest, svCurrent.path());
0119 }
0120 
0121 std::pair<FreeTrajectoryState, double> SteppingHelixPropagator::propagateWithPath(const FreeTrajectoryState& ftsStart,
0122                                                                                   const GlobalPoint& pDest1,
0123                                                                                   const GlobalPoint& pDest2) const {
0124   if ((pDest1 - pDest2).mag() < 1e-10) {
0125     if (sendLogWarning_) {
0126       edm::LogWarning("SteppingHelixPropagator")
0127           << "Can't propagate: the points should be at a bigger distance" << std::endl;
0128     }
0129     return FtsPP();
0130   }
0131   StateArray svBuf;
0132   initStateArraySHPSpecific(svBuf, true);
0133   int nPoints = 0;
0134   setIState(SteppingHelixStateInfo(ftsStart), svBuf, nPoints);
0135 
0136   StateInfo svCurrent;
0137   propagate(svBuf[0], pDest1, pDest2, svCurrent);
0138 
0139   FreeTrajectoryState ftsDest;
0140   svCurrent.getFreeState(ftsDest);
0141 
0142   return FtsPP(ftsDest, svCurrent.path());
0143 }
0144 
0145 std::pair<FreeTrajectoryState, double> SteppingHelixPropagator::propagateWithPath(
0146     const FreeTrajectoryState& ftsStart, const reco::BeamSpot& beamSpot) const {
0147   GlobalPoint pDest1(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
0148   GlobalPoint pDest2(pDest1.x() + beamSpot.dxdz() * 1e3, pDest1.y() + beamSpot.dydz() * 1e3, pDest1.z() + 1e3);
0149   return propagateWithPath(ftsStart, pDest1, pDest2);
0150 }
0151 
0152 void SteppingHelixPropagator::propagate(const SteppingHelixStateInfo& sStart,
0153                                         const Surface& sDest,
0154                                         SteppingHelixStateInfo& out) const {
0155   if (!sStart.isValid()) {
0156     if (sendLogWarning_) {
0157       edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
0158     }
0159     out = invalidState_;
0160     return;
0161   }
0162 
0163   const Plane* pDest = dynamic_cast<const Plane*>(&sDest);
0164   if (pDest != nullptr) {
0165     propagate(sStart, *pDest, out);
0166     return;
0167   }
0168 
0169   const Cylinder* cDest = dynamic_cast<const Cylinder*>(&sDest);
0170   if (cDest != nullptr) {
0171     propagate(sStart, *cDest, out);
0172     return;
0173   }
0174 
0175   throw PropagationException("The surface is neither Cylinder nor Plane");
0176 }
0177 
0178 void SteppingHelixPropagator::propagate(const SteppingHelixStateInfo& sStart,
0179                                         const Plane& pDest,
0180                                         SteppingHelixStateInfo& out) const {
0181   if (!sStart.isValid()) {
0182     if (sendLogWarning_) {
0183       edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
0184     }
0185     out = invalidState_;
0186     return;
0187   }
0188   StateArray svBuf;
0189   initStateArraySHPSpecific(svBuf, true);
0190   int nPoints = 0;
0191   setIState(sStart, svBuf, nPoints);
0192 
0193   Point rPlane(pDest.position().x(), pDest.position().y(), pDest.position().z());
0194   Vector nPlane(pDest.rotation().zx(), pDest.rotation().zy(), pDest.rotation().zz());
0195   nPlane /= nPlane.mag();
0196 
0197   double pars[6] = {rPlane.x(), rPlane.y(), rPlane.z(), nPlane.x(), nPlane.y(), nPlane.z()};
0198 
0199   propagate(svBuf, nPoints, PLANE_DT, pars);
0200 
0201   //(re)set it before leaving: dir =1 (-1) if path increased (decreased) and 0 if it didn't change
0202   //need to implement this somewhere else as a separate function
0203   double lDir = 0;
0204   if (sStart.path() < svBuf[cIndex_(nPoints - 1)].path())
0205     lDir = 1.;
0206   if (sStart.path() > svBuf[cIndex_(nPoints - 1)].path())
0207     lDir = -1.;
0208   svBuf[cIndex_(nPoints - 1)].dir = lDir;
0209 
0210   out = svBuf[cIndex_(nPoints - 1)];
0211   return;
0212 }
0213 
0214 void SteppingHelixPropagator::propagate(const SteppingHelixStateInfo& sStart,
0215                                         const Cylinder& cDest,
0216                                         SteppingHelixStateInfo& out) const {
0217   if (!sStart.isValid()) {
0218     if (sendLogWarning_) {
0219       edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
0220     }
0221     out = invalidState_;
0222     return;
0223   }
0224   StateArray svBuf;
0225   initStateArraySHPSpecific(svBuf, true);
0226   int nPoints = 0;
0227   setIState(sStart, svBuf, nPoints);
0228 
0229   double pars[6] = {0, 0, 0, 0, 0, 0};
0230   pars[RADIUS_P] = cDest.radius();
0231 
0232   propagate(svBuf, nPoints, RADIUS_DT, pars);
0233 
0234   //(re)set it before leaving: dir =1 (-1) if path increased (decreased) and 0 if it didn't change
0235   //need to implement this somewhere else as a separate function
0236   double lDir = 0;
0237   if (sStart.path() < svBuf[cIndex_(nPoints - 1)].path())
0238     lDir = 1.;
0239   if (sStart.path() > svBuf[cIndex_(nPoints - 1)].path())
0240     lDir = -1.;
0241   svBuf[cIndex_(nPoints - 1)].dir = lDir;
0242   out = svBuf[cIndex_(nPoints - 1)];
0243   return;
0244 }
0245 
0246 void SteppingHelixPropagator::propagate(const SteppingHelixStateInfo& sStart,
0247                                         const GlobalPoint& pDest,
0248                                         SteppingHelixStateInfo& out) const {
0249   if (!sStart.isValid()) {
0250     if (sendLogWarning_) {
0251       edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
0252     }
0253     out = invalidState_;
0254     return;
0255   }
0256   StateArray svBuf;
0257   initStateArraySHPSpecific(svBuf, true);
0258   int nPoints = 0;
0259   setIState(sStart, svBuf, nPoints);
0260 
0261   double pars[6] = {pDest.x(), pDest.y(), pDest.z(), 0, 0, 0};
0262 
0263   propagate(svBuf, nPoints, POINT_PCA_DT, pars);
0264 
0265   out = svBuf[cIndex_(nPoints - 1)];
0266   return;
0267 }
0268 
0269 void SteppingHelixPropagator::propagate(const SteppingHelixStateInfo& sStart,
0270                                         const GlobalPoint& pDest1,
0271                                         const GlobalPoint& pDest2,
0272                                         SteppingHelixStateInfo& out) const {
0273   if ((pDest1 - pDest2).mag() < 1e-10 || !sStart.isValid()) {
0274     if (sendLogWarning_) {
0275       if ((pDest1 - pDest2).mag() < 1e-10)
0276         edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: points are too close" << std::endl;
0277       if (!sStart.isValid())
0278         edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
0279     }
0280     out = invalidState_;
0281     return;
0282   }
0283   StateArray svBuf;
0284   initStateArraySHPSpecific(svBuf, true);
0285   int nPoints = 0;
0286   setIState(sStart, svBuf, nPoints);
0287 
0288   double pars[6] = {pDest1.x(), pDest1.y(), pDest1.z(), pDest2.x(), pDest2.y(), pDest2.z()};
0289 
0290   propagate(svBuf, nPoints, LINE_PCA_DT, pars);
0291 
0292   out = svBuf[cIndex_(nPoints - 1)];
0293   return;
0294 }
0295 
0296 void SteppingHelixPropagator::setIState(const SteppingHelixStateInfo& sStart, StateArray& svBuf, int& nPoints) const {
0297   nPoints = 0;
0298   svBuf[cIndex_(nPoints)] = sStart;  //do this anyways to have a fresh start
0299   if (sStart.isComplete) {
0300     svBuf[cIndex_(nPoints)] = sStart;
0301     nPoints++;
0302   } else {
0303     loadState(svBuf[cIndex_(nPoints)], sStart.p3, sStart.r3, sStart.q, propagationDirection(), sStart.covCurv);
0304     nPoints++;
0305   }
0306   svBuf[cIndex_(0)].hasErrorPropagated_ = sStart.hasErrorPropagated_ & !noErrorPropagation_;
0307 }
0308 
0309 SteppingHelixPropagator::Result SteppingHelixPropagator::propagate(StateArray& svBuf,
0310                                                                    int& nPoints,
0311                                                                    SteppingHelixPropagator::DestType type,
0312                                                                    const double pars[6],
0313                                                                    double epsilon) const {
0314   static const std::string metname = "SteppingHelixPropagator";
0315   StateInfo* svCurrent = &svBuf[cIndex_(nPoints - 1)];
0316 
0317   //check if it's going to work at all
0318   double tanDist = 0;
0319   double dist = 0;
0320   PropagationDirection refDirection = anyDirection;
0321   Result result = refToDest(type, (*svCurrent), pars, dist, tanDist, refDirection);
0322 
0323   if (result != SteppingHelixStateInfo::OK || fabs(dist) > 1e6) {
0324     svCurrent->status_ = result;
0325     if (fabs(dist) > 1e6)
0326       svCurrent->status_ = SteppingHelixStateInfo::INACC;
0327     svCurrent->isValid_ = false;
0328     svCurrent->field = field_;
0329     if (sendLogWarning_) {
0330       edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0331                                << " Failed after first refToDest check with status "
0332                                << SteppingHelixStateInfo::ResultName[result] << std::endl;
0333     }
0334     return result;
0335   }
0336 
0337   result = SteppingHelixStateInfo::UNDEFINED;
0338   bool makeNextStep = true;
0339   double dStep = defaultStep_;
0340   PropagationDirection dir, oldDir;
0341   dir = propagationDirection();
0342   oldDir = dir;
0343   int nOsc = 0;
0344 
0345   double distMag = 1e12;
0346   double tanDistMag = 1e12;
0347 
0348   double distMat = 1e12;
0349   double tanDistMat = 1e12;
0350 
0351   double tanDistNextCheck = -0.1;  //just need a negative start val
0352   double tanDistMagNextCheck = -0.1;
0353   double tanDistMatNextCheck = -0.1;
0354   double oldDStep = 0;
0355   PropagationDirection oldRefDirection = propagationDirection();
0356 
0357   Result resultToMat = SteppingHelixStateInfo::UNDEFINED;
0358   Result resultToMag = SteppingHelixStateInfo::UNDEFINED;
0359 
0360   bool isFirstStep = true;
0361   bool expectNewMagVolume = false;
0362 
0363   int loopCount = 0;
0364   while (makeNextStep) {
0365     dStep = defaultStep_;
0366     svCurrent = &svBuf[cIndex_(nPoints - 1)];
0367     double curZ = svCurrent->r3.z();
0368     double curR = svCurrent->r3.perp();
0369     if (fabs(curZ) < 440 && curR < 260)
0370       dStep = defaultStep_ * 2;
0371 
0372     //more such ifs might be scattered around
0373     //even though dStep is large, it will still make stops at each volume boundary
0374     if (useTuningForL2Speed_) {
0375       dStep = 100.;
0376       if (!useIsYokeFlag_ && fabs(curZ) < 667 && curR > 380 && curR < 850) {
0377         dStep = 5 * (1 + 0.2 * svCurrent->p3.mag());
0378       }
0379     }
0380 
0381     //    refDirection = propagationDirection();
0382 
0383     tanDistNextCheck -= oldDStep;
0384     tanDistMagNextCheck -= oldDStep;
0385     tanDistMatNextCheck -= oldDStep;
0386 
0387     if (tanDistNextCheck < 0) {
0388       //use pre-computed values if it's the first step
0389       if (!isFirstStep)
0390         refToDest(type, (*svCurrent), pars, dist, tanDist, refDirection);
0391       // constrain allowed path for a tangential approach
0392       if (fabs(tanDist) > 4. * (fabs(dist)))
0393         tanDist *= tanDist == 0 ? 0 : fabs(dist / tanDist * 4.);
0394 
0395       tanDistNextCheck = fabs(tanDist) * 0.5 - 0.5;  //need a better guess (to-do)
0396       //reasonable limit
0397       if (tanDistNextCheck > defaultStep_ * 20.)
0398         tanDistNextCheck = defaultStep_ * 20.;
0399       oldRefDirection = refDirection;
0400     } else {
0401       tanDist = tanDist > 0. ? tanDist - oldDStep : tanDist + oldDStep;
0402       refDirection = oldRefDirection;
0403       if (debug_)
0404         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0405                           << "Skipped refToDest: guess tanDist = " << tanDist << " next check at " << tanDistNextCheck
0406                           << std::endl;
0407     }
0408     //! define a fast-skip distance: should be the shortest of a possible step or distance
0409     double fastSkipDist = fabs(dist) > fabs(tanDist) ? tanDist : dist;
0410 
0411     if (propagationDirection() == anyDirection) {
0412       dir = refDirection;
0413     } else {
0414       dir = propagationDirection();
0415       if (fabs(tanDist) < 0.1 && refDirection != dir) {
0416         //how did it get here?  nOsc++;
0417         dir = refDirection;
0418         if (debug_)
0419           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0420                             << "NOTE: overstepped last time: switch direction (can do it if within 1 mm)" << std::endl;
0421       }
0422     }
0423 
0424     if (useMagVolumes_ && !(fabs(dist) < fabs(epsilon))) {  //need to know the general direction
0425       if (tanDistMagNextCheck < 0) {
0426         resultToMag = refToMagVolume(
0427             (*svCurrent), dir, distMag, tanDistMag, fabs(fastSkipDist), expectNewMagVolume, fabs(tanDist));
0428         // constrain allowed path for a tangential approach
0429         if (fabs(tanDistMag) > 4. * (fabs(distMag)))
0430           tanDistMag *= tanDistMag == 0 ? 0 : fabs(distMag / tanDistMag * 4.);
0431 
0432         tanDistMagNextCheck = fabs(tanDistMag) * 0.5 - 0.5;  //need a better guess (to-do)
0433         //reasonable limit; "turn off" checking if bounds are further than the destination
0434         if (tanDistMagNextCheck > defaultStep_ * 20. || fabs(dist) < fabs(distMag) ||
0435             resultToMag == SteppingHelixStateInfo::INACC)
0436           tanDistMagNextCheck = defaultStep_ * 20 > fabs(fastSkipDist) ? fabs(fastSkipDist) : defaultStep_ * 20;
0437         if (resultToMag != SteppingHelixStateInfo::INACC && resultToMag != SteppingHelixStateInfo::OK)
0438           tanDistMagNextCheck = -1;
0439       } else {
0440         //  resultToMag = SteppingHelixStateInfo::OK;
0441         tanDistMag = tanDistMag > 0. ? tanDistMag - oldDStep : tanDistMag + oldDStep;
0442         if (debug_)
0443           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0444                             << "Skipped refToMag: guess tanDistMag = " << tanDistMag << " next check at "
0445                             << tanDistMagNextCheck;
0446       }
0447     }
0448 
0449     if (useMatVolumes_ && !(fabs(dist) < fabs(epsilon))) {  //need to know the general direction
0450       if (tanDistMatNextCheck < 0) {
0451         resultToMat = refToMatVolume((*svCurrent), dir, distMat, tanDistMat, fabs(fastSkipDist));
0452         // constrain allowed path for a tangential approach
0453         if (fabs(tanDistMat) > 4. * (fabs(distMat)))
0454           tanDistMat *= tanDistMat == 0 ? 0 : fabs(distMat / tanDistMat * 4.);
0455 
0456         tanDistMatNextCheck = fabs(tanDistMat) * 0.5 - 0.5;  //need a better guess (to-do)
0457         //reasonable limit; "turn off" checking if bounds are further than the destination
0458         if (tanDistMatNextCheck > defaultStep_ * 20. || fabs(dist) < fabs(distMat) ||
0459             resultToMat == SteppingHelixStateInfo::INACC)
0460           tanDistMatNextCheck = defaultStep_ * 20 > fabs(fastSkipDist) ? fabs(fastSkipDist) : defaultStep_ * 20;
0461         if (resultToMat != SteppingHelixStateInfo::INACC && resultToMat != SteppingHelixStateInfo::OK)
0462           tanDistMatNextCheck = -1;
0463       } else {
0464         //  resultToMat = SteppingHelixStateInfo::OK;
0465         tanDistMat = tanDistMat > 0. ? tanDistMat - oldDStep : tanDistMat + oldDStep;
0466         if (debug_)
0467           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0468                             << "Skipped refToMat: guess tanDistMat = " << tanDistMat << " next check at "
0469                             << tanDistMatNextCheck;
0470       }
0471     }
0472 
0473     double rDotP = svCurrent->r3.dot(svCurrent->p3);
0474     if ((fabs(curZ) > 1.5e3 || curR > 800.) &&
0475         ((dir == alongMomentum && rDotP > 0) || (dir == oppositeToMomentum && rDotP < 0))) {
0476       dStep = fabs(tanDist) - 1e-12;
0477     }
0478     double tanDistMin = fabs(tanDist);
0479     if (tanDistMin > fabs(tanDistMag) + 0.05 &&
0480         (resultToMag == SteppingHelixStateInfo::OK || resultToMag == SteppingHelixStateInfo::WRONG_VOLUME)) {
0481       tanDistMin = fabs(tanDistMag) + 0.05;  //try to step into the next volume
0482       expectNewMagVolume = true;
0483     } else
0484       expectNewMagVolume = false;
0485 
0486     if (tanDistMin > fabs(tanDistMat) + 0.05 && resultToMat == SteppingHelixStateInfo::OK) {
0487       tanDistMin = fabs(tanDistMat) + 0.05;  //try to step into the next volume
0488       if (expectNewMagVolume)
0489         expectNewMagVolume = false;
0490     }
0491 
0492     if (tanDistMin * fabs(svCurrent->dEdx) > 0.5 * svCurrent->p3.mag()) {
0493       tanDistMin = 0.5 * svCurrent->p3.mag() / fabs(svCurrent->dEdx);
0494       if (expectNewMagVolume)
0495         expectNewMagVolume = false;
0496     }
0497 
0498     double tanDistMinLazy = fabs(tanDistMin);
0499     if ((type == POINT_PCA_DT || type == LINE_PCA_DT) && fabs(tanDist) < 2. * fabs(tanDistMin)) {
0500       //being lazy here; the best is to take into account the curvature
0501       tanDistMinLazy = fabs(tanDistMin) * 0.5;
0502     }
0503 
0504     if (fabs(tanDistMinLazy) < dStep) {
0505       dStep = fabs(tanDistMinLazy);
0506     }
0507 
0508     //keep this path length for the next step
0509     oldDStep = dStep;
0510 
0511     if (dStep > 1e-10 && !(fabs(dist) < fabs(epsilon))) {
0512       StateInfo* svNext = &svBuf[cIndex_(nPoints)];
0513       makeAtomStep((*svCurrent), (*svNext), dStep, dir, HEL_AS_F);
0514       //       if (useMatVolumes_ && expectNewMagVolume
0515       //      && svCurrent->magVol == svNext->magVol){
0516       //    double tmpDist=0;
0517       //    double tmpDistMag = 0;
0518       //    if (refToMagVolume((*svNext), dir, tmpDist, tmpDistMag, fabs(dist)) != SteppingHelixStateInfo::OK){
0519       //    //the point appears to be outside, but findVolume claims the opposite
0520       //      dStep += 0.05*fabs(tanDistMag/distMag); oldDStep = dStep; //do it again with a bigger step
0521       //      if (debug_) LogTrace(metname)
0522       //        <<"Failed to get into new mag volume: will try with new bigger step "<<dStep<<std::endl;
0523       //      makeAtomStep((*svCurrent), (*svNext), dStep, dir, HEL_AS_F);
0524       //    }
0525       //       }
0526       nPoints++;
0527       svCurrent = &svBuf[cIndex_(nPoints - 1)];
0528       if (oldDir != dir) {
0529         nOsc++;
0530         tanDistNextCheck = -1;  //check dist after osc
0531         tanDistMagNextCheck = -1;
0532         tanDistMatNextCheck = -1;
0533       }
0534       oldDir = dir;
0535     }
0536 
0537     if (nOsc > 1 && fabs(dStep) > epsilon) {
0538       if (debug_)
0539         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Ooops" << std::endl;
0540     }
0541 
0542     if (fabs(dist) < fabs(epsilon))
0543       result = SteppingHelixStateInfo::OK;
0544 
0545     if ((type == POINT_PCA_DT || type == LINE_PCA_DT) && fabs(dStep) < fabs(epsilon)) {
0546       //now check if it's not a branch point (peek ahead at 1 cm)
0547       double nextDist = 0;
0548       double nextTanDist = 0;
0549       PropagationDirection nextRefDirection = anyDirection;
0550       StateInfo* svNext = &svBuf[cIndex_(nPoints)];
0551       makeAtomStep((*svCurrent), (*svNext), 1., dir, HEL_AS_F);
0552       nPoints++;
0553       svCurrent = &svBuf[cIndex_(nPoints - 1)];
0554       refToDest(type, (*svCurrent), pars, nextDist, nextTanDist, nextRefDirection);
0555       if (fabs(nextDist) > fabs(dist)) {
0556         nPoints--;
0557         svCurrent = &svBuf[cIndex_(nPoints - 1)];
0558         result = SteppingHelixStateInfo::OK;
0559         if (debug_) {
0560           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0561                             << "Found real local minimum in PCA" << std::endl;
0562         }
0563       } else {
0564         //keep this trial point and continue
0565         dStep = defaultStep_;
0566         if (debug_) {
0567           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Found branch point in PCA"
0568                             << std::endl;
0569         }
0570       }
0571     }
0572 
0573     if (nPoints > MAX_STEPS * 1. / defaultStep_ || loopCount > MAX_STEPS * 100 || nOsc > 6)
0574       result = SteppingHelixStateInfo::FAULT;
0575 
0576     if (svCurrent->p3.mag() < 0.1)
0577       result = SteppingHelixStateInfo::RANGEOUT;
0578 
0579     curZ = svCurrent->r3.z();
0580     curR = svCurrent->r3.perp();
0581     if (curR > 20000 || fabs(curZ) > 20000)
0582       result = SteppingHelixStateInfo::INACC;
0583 
0584     makeNextStep = result == SteppingHelixStateInfo::UNDEFINED;
0585     svCurrent->status_ = result;
0586     svCurrent->isValid_ = (result == SteppingHelixStateInfo::OK || makeNextStep);
0587     svCurrent->field = field_;
0588 
0589     isFirstStep = false;
0590     loopCount++;
0591   }
0592 
0593   if (sendLogWarning_ && result != SteppingHelixStateInfo::OK) {
0594     edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0595                              << " Propagation failed with status " << SteppingHelixStateInfo::ResultName[result]
0596                              << std::endl;
0597     if (result == SteppingHelixStateInfo::RANGEOUT)
0598       edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0599                                << " Momentum at last point is too low (<0.1) p_last = " << svCurrent->p3.mag()
0600                                << std::endl;
0601     if (result == SteppingHelixStateInfo::INACC)
0602       edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0603                                << " Went too far: the last point is at " << svCurrent->r3 << std::endl;
0604     if (result == SteppingHelixStateInfo::FAULT && nOsc > 6)
0605       edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0606                                << " Infinite loop condidtion detected: going in cycles. Break after 6 cycles"
0607                                << std::endl;
0608     if (result == SteppingHelixStateInfo::FAULT && nPoints > MAX_STEPS * 1. / defaultStep_)
0609       edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0610                                << " Tired to go farther. Made too many steps: more than "
0611                                << MAX_STEPS * 1. / defaultStep_ << std::endl;
0612   }
0613 
0614   if (debug_) {
0615     switch (type) {
0616       case RADIUS_DT:
0617         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to radius "
0618                           << pars[RADIUS_P] << std::endl;
0619         break;
0620       case Z_DT:
0621         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to z " << pars[Z_P]
0622                           << std::endl;
0623         break;
0624       case PATHL_DT:
0625         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to pathL "
0626                           << pars[PATHL_P] << std::endl;
0627         break;
0628       case PLANE_DT: {
0629         Point rPlane(pars[0], pars[1], pars[2]);
0630         Vector nPlane(pars[3], pars[4], pars[5]);
0631         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to plane r0:" << rPlane
0632                           << " n:" << nPlane << std::endl;
0633       } break;
0634       case POINT_PCA_DT: {
0635         Point rDest(pars[0], pars[1], pars[2]);
0636         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to PCA to point "
0637                           << rDest << std::endl;
0638       } break;
0639       case LINE_PCA_DT: {
0640         Point rDest1(pars[0], pars[1], pars[2]);
0641         Point rDest2(pars[3], pars[4], pars[5]);
0642         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to PCA to line "
0643                           << rDest1 << " - " << rDest2 << std::endl;
0644       } break;
0645       default:
0646         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to NOT IMPLEMENTED"
0647                           << std::endl;
0648         break;
0649     }
0650     LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Made " << nPoints - 1
0651                       << " steps and stopped at(cur step) " << svCurrent->r3 << " nOsc " << nOsc << std::endl;
0652   }
0653 
0654   return result;
0655 }
0656 
0657 void SteppingHelixPropagator::loadState(SteppingHelixPropagator::StateInfo& svCurrent,
0658                                         const SteppingHelixPropagator::Vector& p3,
0659                                         const SteppingHelixPropagator::Point& r3,
0660                                         int charge,
0661                                         PropagationDirection dir,
0662                                         const AlgebraicSymMatrix55& covCurv) const {
0663   static const std::string metname = "SteppingHelixPropagator";
0664 
0665   svCurrent.q = charge;
0666   svCurrent.p3 = p3;
0667   svCurrent.r3 = r3;
0668   svCurrent.dir = dir == alongMomentum ? 1. : -1.;
0669 
0670   svCurrent.path_ = 0;  // this could've held the initial path
0671   svCurrent.radPath_ = 0;
0672 
0673   GlobalPoint gPointNorZ(svCurrent.r3.x(), svCurrent.r3.y(), svCurrent.r3.z());
0674 
0675   float gpmag = gPointNorZ.mag2();
0676   float pmag2 = p3.mag2();
0677   if (gpmag > 1e20f) {
0678     LogTrace(metname) << "Initial point is too far";
0679     svCurrent.isValid_ = false;
0680     return;
0681   }
0682   if (pmag2 < 1e-18f) {
0683     LogTrace(metname) << "Initial momentum is too low";
0684     svCurrent.isValid_ = false;
0685     return;
0686   }
0687   if (!(gpmag == gpmag)) {
0688     LogTrace(metname) << "Initial point is a nan";
0689     edm::LogWarning("SteppingHelixPropagatorNAN") << "Initial point is a nan";
0690     svCurrent.isValid_ = false;
0691     return;
0692   }
0693   if (!(pmag2 == pmag2)) {
0694     LogTrace(metname) << "Initial momentum is a nan";
0695     edm::LogWarning("SteppingHelixPropagatorNAN") << "Initial momentum is a nan";
0696     svCurrent.isValid_ = false;
0697     return;
0698   }
0699 
0700   GlobalVector bf(0, 0, 0);
0701   // = field_->inTesla(gPoint);
0702   if (useMagVolumes_) {
0703     if (vbField_) {
0704       svCurrent.magVol = vbField_->findVolume(gPointNorZ);
0705       if (useIsYokeFlag_) {
0706         double curRad = svCurrent.r3.perp();
0707         if (curRad > 380 && curRad < 850 && fabs(svCurrent.r3.z()) < 667) {
0708           svCurrent.isYokeVol = isYokeVolume(svCurrent.magVol);
0709         } else {
0710           svCurrent.isYokeVol = false;
0711         }
0712       }
0713     } else {
0714       edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0715                                << "Failed to cast into VolumeBasedMagneticField: fall back to the default behavior"
0716                                << std::endl;
0717       svCurrent.magVol = nullptr;
0718     }
0719     if (debug_) {
0720       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Got volume at "
0721                         << svCurrent.magVol << std::endl;
0722     }
0723   }
0724 
0725   if (useMagVolumes_ && svCurrent.magVol != nullptr && !useInTeslaFromMagField_) {
0726     bf = svCurrent.magVol->inTesla(gPointNorZ);
0727     if (debug_) {
0728       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Loaded bfield float: " << bf
0729                         << " at global float " << gPointNorZ << " double " << svCurrent.r3 << std::endl;
0730       LocalPoint lPoint(svCurrent.magVol->toLocal(gPointNorZ));
0731       LocalVector lbf = svCurrent.magVol->fieldInTesla(lPoint);
0732       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "\t cf in local locF: " << lbf
0733                         << " at " << lPoint << std::endl;
0734     }
0735     svCurrent.bf.set(bf.x(), bf.y(), bf.z());
0736   } else {
0737     GlobalPoint gPoint(r3.x(), r3.y(), r3.z());
0738     bf = field_->inTesla(gPoint);
0739     svCurrent.bf.set(bf.x(), bf.y(), bf.z());
0740   }
0741   if (svCurrent.bf.mag2() < 1e-32)
0742     svCurrent.bf.set(0., 0., 1e-16);
0743   if (debug_) {
0744     LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0745                       << "Loaded bfield double: " << svCurrent.bf << "  from float: " << bf << " at float "
0746                       << gPointNorZ << " double " << svCurrent.r3 << std::endl;
0747   }
0748 
0749   double dEdXPrime = 0;
0750   double dEdx = 0;
0751   double radX0 = 0;
0752   MatBounds rzLims;
0753   dEdx = getDeDx(svCurrent, dEdXPrime, radX0, rzLims);
0754   svCurrent.dEdx = dEdx;
0755   svCurrent.dEdXPrime = dEdXPrime;
0756   svCurrent.radX0 = radX0;
0757   svCurrent.rzLims = rzLims;
0758 
0759   svCurrent.covCurv = covCurv;
0760 
0761   svCurrent.isComplete = true;
0762   svCurrent.isValid_ = true;
0763 
0764   if (debug_) {
0765     LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0766                       << "Loaded at  path: " << svCurrent.path() << " radPath: " << svCurrent.radPath() << " p3 "
0767                       << " pt: " << svCurrent.p3.perp() << " phi: " << svCurrent.p3.phi()
0768                       << " eta: " << svCurrent.p3.eta() << " " << svCurrent.p3 << " r3: " << svCurrent.r3
0769                       << " bField: " << svCurrent.bf.mag() << std::endl;
0770     LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0771                       << "Input Covariance in Curvilinear RF " << covCurv << std::endl;
0772   }
0773 }
0774 
0775 void SteppingHelixPropagator::getNextState(const SteppingHelixPropagator::StateInfo& svPrevious,
0776                                            SteppingHelixPropagator::StateInfo& svNext,
0777                                            double dP,
0778                                            const SteppingHelixPropagator::Vector& tau,
0779                                            const SteppingHelixPropagator::Vector& drVec,
0780                                            double dS,
0781                                            double dX0,
0782                                            const AlgebraicMatrix55& dCovCurvTransform) const {
0783   static const std::string metname = "SteppingHelixPropagator";
0784   svNext.q = svPrevious.q;
0785   svNext.dir = dS > 0.0 ? 1. : -1.;
0786   svNext.p3 = tau;
0787   svNext.p3 *= (svPrevious.p3.mag() - svNext.dir * fabs(dP));
0788 
0789   svNext.r3 = svPrevious.r3;
0790   svNext.r3 += drVec;
0791 
0792   svNext.path_ = svPrevious.path() + dS;
0793   svNext.radPath_ = svPrevious.radPath() + dX0;
0794 
0795   GlobalPoint gPointNorZ(svNext.r3.x(), svNext.r3.y(), svNext.r3.z());
0796 
0797   GlobalVector bf(0, 0, 0);
0798 
0799   if (useMagVolumes_) {
0800     if (vbField_ != nullptr) {
0801       svNext.magVol = vbField_->findVolume(gPointNorZ);
0802       if (useIsYokeFlag_) {
0803         double curRad = svNext.r3.perp();
0804         if (curRad > 380 && curRad < 850 && fabs(svNext.r3.z()) < 667) {
0805           svNext.isYokeVol = isYokeVolume(svNext.magVol);
0806         } else {
0807           svNext.isYokeVol = false;
0808         }
0809       }
0810     } else {
0811       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
0812                         << "Failed to cast into VolumeBasedMagneticField" << std::endl;
0813       svNext.magVol = nullptr;
0814     }
0815     if (debug_) {
0816       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Got volume at "
0817                         << svNext.magVol << std::endl;
0818     }
0819   }
0820 
0821   if (useMagVolumes_ && svNext.magVol != nullptr && !useInTeslaFromMagField_) {
0822     bf = svNext.magVol->inTesla(gPointNorZ);
0823     svNext.bf.set(bf.x(), bf.y(), bf.z());
0824   } else {
0825     GlobalPoint gPoint(svNext.r3.x(), svNext.r3.y(), svNext.r3.z());
0826     bf = field_->inTesla(gPoint);
0827     svNext.bf.set(bf.x(), bf.y(), bf.z());
0828   }
0829   if (svNext.bf.mag2() < 1e-32)
0830     svNext.bf.set(0., 0., 1e-16);
0831 
0832   double dEdXPrime = 0;
0833   double dEdx = 0;
0834   double radX0 = 0;
0835   MatBounds rzLims;
0836   dEdx = getDeDx(svNext, dEdXPrime, radX0, rzLims);
0837   svNext.dEdx = dEdx;
0838   svNext.dEdXPrime = dEdXPrime;
0839   svNext.radX0 = radX0;
0840   svNext.rzLims = rzLims;
0841 
0842   //update Emat only if it's valid
0843   svNext.hasErrorPropagated_ = svPrevious.hasErrorPropagated_;
0844   if (svPrevious.hasErrorPropagated_) {
0845     {
0846       AlgebraicMatrix55 tmp = dCovCurvTransform * svPrevious.covCurv;
0847       ROOT::Math::AssignSym::Evaluate(svNext.covCurv, tmp * ROOT::Math::Transpose(dCovCurvTransform));
0848 
0849       svNext.covCurv += svPrevious.matDCovCurv;
0850     }
0851   } else {
0852     //could skip dragging along the unprop. cov later .. now
0853     // svNext.cov = svPrevious.cov;
0854   }
0855 
0856   if (debug_) {
0857     LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Now at  path: " << svNext.path()
0858                       << " radPath: " << svNext.radPath() << " p3 "
0859                       << " pt: " << svNext.p3.perp() << " phi: " << svNext.p3.phi() << " eta: " << svNext.p3.eta()
0860                       << " " << svNext.p3 << " r3: " << svNext.r3
0861                       << " dPhi: " << acos(svNext.p3.unit().dot(svPrevious.p3.unit())) << " bField: " << svNext.bf.mag()
0862                       << " dedx: " << svNext.dEdx << std::endl;
0863     LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "New Covariance "
0864                       << svNext.covCurv << std::endl;
0865     LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Transf by dCovTransform "
0866                       << dCovCurvTransform << std::endl;
0867   }
0868 }
0869 
0870 void SteppingHelixPropagator::setRep(SteppingHelixPropagator::Basis& rep,
0871                                      const SteppingHelixPropagator::Vector& tau) const {
0872   Vector zRep(0., 0., 1.);
0873   rep.lX = tau;
0874   rep.lY = zRep.cross(tau);
0875   rep.lY *= 1. / tau.perp();
0876   rep.lZ = rep.lX.cross(rep.lY);
0877 }
0878 
0879 bool SteppingHelixPropagator::makeAtomStep(SteppingHelixPropagator::StateInfo& svCurrent,
0880                                            SteppingHelixPropagator::StateInfo& svNext,
0881                                            double dS,
0882                                            PropagationDirection dir,
0883                                            SteppingHelixPropagator::Fancy fancy) const {
0884   static const std::string metname = "SteppingHelixPropagator";
0885   if (debug_) {
0886     LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Make atom step "
0887                       << svCurrent.path() << " with step " << dS << " in direction " << dir << std::endl;
0888   }
0889 
0890   AlgebraicMatrix55 dCCurvTransform(unit55_);
0891 
0892   double dP = 0;
0893   double curP = svCurrent.p3.mag();
0894   Vector tau = svCurrent.p3;
0895   tau *= 1. / curP;
0896   Vector tauNext(tau);
0897   Vector drVec(0, 0, 0);
0898 
0899   dS = dir == alongMomentum ? fabs(dS) : -fabs(dS);
0900 
0901   double radX0 = 1e24;
0902 
0903   switch (fancy) {
0904     case HEL_AS_F:
0905     case HEL_ALL_F: {
0906       double p0 = curP;  // see above = svCurrent.p3.mag();
0907       double p0Inv = 1. / p0;
0908       double b0 = svCurrent.bf.mag();
0909 
0910       //get to the mid-point first
0911       double phi = (2.99792458e-3 * svCurrent.q * b0 * p0Inv) * dS / 2.;
0912       bool phiSmall = fabs(phi) < 1e-4;
0913 
0914       double cosPhi = 0;
0915       double sinPhi = 0;
0916 
0917       double oneLessCosPhi = 0;
0918       double oneLessCosPhiOPhi = 0;
0919       double phiLessSinPhiOPhi = 0;
0920 
0921       if (phiSmall) {
0922         double phi2 = phi * phi;
0923         double phi3 = phi2 * phi;
0924         double phi4 = phi3 * phi;
0925         oneLessCosPhi = phi2 / 2. - phi4 / 24. + phi2 * phi4 / 720.;             // 0.5*phi*phi;//*(1.- phi*phi/12.);
0926         oneLessCosPhiOPhi = 0.5 * phi - phi3 / 24. + phi2 * phi3 / 720.;         //*(1.- phi*phi/12.);
0927         phiLessSinPhiOPhi = phi * phi / 6. - phi4 / 120. + phi4 * phi2 / 5040.;  //*(1. - phi*phi/20.);
0928       } else {
0929         cosPhi = cos(phi);
0930         sinPhi = sin(phi);
0931         oneLessCosPhi = 1. - cosPhi;
0932         oneLessCosPhiOPhi = oneLessCosPhi / phi;
0933         double sinPhiOPhi = sinPhi / phi;
0934         phiLessSinPhiOPhi = 1 - sinPhiOPhi;
0935       }
0936 
0937       Vector bHat = svCurrent.bf;
0938       bHat *= 1. / b0;  //bHat.mag();
0939       //    bool bAlongZ = fabs(bHat.z()) > 0.9999;
0940 
0941       Vector btVec(bHat.cross(tau));  // for balong z btVec.z()==0
0942       double tauB = tau.dot(bHat);
0943       Vector bbtVec(bHat * tauB - tau);  // (-tau.x(), -tau.y(), 0)
0944 
0945       //don't need it here    tauNext = tau + bbtVec*oneLessCosPhi - btVec*sinPhi;
0946       drVec = bbtVec * phiLessSinPhiOPhi;
0947       drVec -= btVec * oneLessCosPhiOPhi;
0948       drVec += tau;
0949       drVec *= dS / 2.;
0950 
0951       double dEdx = svCurrent.dEdx;
0952       double dEdXPrime = svCurrent.dEdXPrime;
0953       radX0 = svCurrent.radX0;
0954       dP = dEdx * dS;
0955 
0956       //improve with above values:
0957       drVec += svCurrent.r3;
0958       GlobalVector bfGV(0, 0, 0);
0959       Vector bf(0, 0, 0);
0960       if (useMagVolumes_ && svCurrent.magVol != nullptr && !useInTeslaFromMagField_) {
0961         bfGV = svCurrent.magVol->inTesla(GlobalPoint(drVec.x(), drVec.y(), drVec.z()));
0962         bf.set(bfGV.x(), bfGV.y(), bfGV.z());
0963       } else {
0964         bfGV = field_->inTesla(GlobalPoint(drVec.x(), drVec.y(), drVec.z()));
0965         bf.set(bfGV.x(), bfGV.y(), bfGV.z());
0966       }
0967       double b0Init = b0;
0968       b0 = bf.mag();
0969       if (b0 < 1e-16) {
0970         b0 = 1e-16;
0971         bf.set(0., 0., 1e-16);
0972       }
0973       if (debug_) {
0974         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Improved b " << b0
0975                           << " at r3 " << drVec << std::endl;
0976       }
0977 
0978       if (fabs((b0 - b0Init) * dS) > 1) {
0979         //missed the mag volume boundary?
0980         if (debug_) {
0981           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Large bf*dS change "
0982                             << fabs((b0 - svCurrent.bf.mag()) * dS) << " --> recalc dedx" << std::endl;
0983         }
0984         svNext.r3 = drVec;
0985         svNext.bf = bf;
0986         svNext.p3 = svCurrent.p3;
0987         svNext.isYokeVol = svCurrent.isYokeVol;
0988         svNext.magVol = svCurrent.magVol;
0989         MatBounds rzTmp;
0990         dEdx = getDeDx(svNext, dEdXPrime, radX0, rzTmp);
0991         dP = dEdx * dS;
0992         if (debug_) {
0993           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "New dEdX= " << dEdx
0994                             << " dP= " << dP << " for p0 " << p0 << std::endl;
0995         }
0996       }
0997       //p0 is mid-way and b0 from mid-point
0998       p0 += dP / 2.;
0999       if (p0 < 1e-2)
1000         p0 = 1e-2;
1001       p0Inv = 1. / p0;
1002 
1003       phi = (2.99792458e-3 * svCurrent.q * b0 * p0Inv) * dS;
1004       phiSmall = fabs(phi) < 1e-4;
1005 
1006       if (phiSmall) {
1007         double phi2 = phi * phi;
1008         double phi3 = phi2 * phi;
1009         double phi4 = phi3 * phi;
1010         sinPhi = phi - phi3 / 6. + phi4 * phi / 120.;
1011         cosPhi = 1. - phi2 / 2. + phi4 / 24.;
1012         oneLessCosPhi = phi2 / 2. - phi4 / 24. + phi2 * phi4 / 720.;             // 0.5*phi*phi;//*(1.- phi*phi/12.);
1013         oneLessCosPhiOPhi = 0.5 * phi - phi3 / 24. + phi2 * phi3 / 720.;         //*(1.- phi*phi/12.);
1014         phiLessSinPhiOPhi = phi * phi / 6. - phi4 / 120. + phi4 * phi2 / 5040.;  //*(1. - phi*phi/20.);
1015       } else {
1016         cosPhi = cos(phi);
1017         sinPhi = sin(phi);
1018         oneLessCosPhi = 1. - cosPhi;
1019         oneLessCosPhiOPhi = oneLessCosPhi / phi;
1020         double sinPhiOPhi = sinPhi / phi;
1021         phiLessSinPhiOPhi = 1. - sinPhiOPhi;
1022       }
1023 
1024       bHat = bf;
1025       bHat *= 1. / b0;  // as above =1./bHat.mag();
1026       //    bAlongZ = fabs(bHat.z()) > 0.9999;
1027       btVec = bHat.cross(tau);  // for b||z (-tau.y(), tau.x() ,0)
1028       tauB = tau.dot(bHat);
1029       bbtVec = bHat * tauB - tau;  //bHat.cross(btVec); for b||z: (-tau.x(), -tau.y(), 0)
1030 
1031       tauNext = bbtVec * oneLessCosPhi;
1032       tauNext -= btVec * sinPhi;
1033       tauNext += tau;  //for b||z tauNext.z() == tau.z()
1034       double tauNextPerpInv = 1. / tauNext.perp();
1035       drVec = bbtVec * phiLessSinPhiOPhi;
1036       drVec -= btVec * oneLessCosPhiOPhi;
1037       drVec += tau;
1038       drVec *= dS;
1039 
1040       if (svCurrent.hasErrorPropagated_) {
1041         double theta02 = 0;
1042         double dX0 = fabs(dS) / radX0;
1043 
1044         if (applyRadX0Correction_) {
1045           // this provides the integrand for theta^2
1046           // if summed up along the path, should result in
1047           // theta_total^2 = Int_0^x0{ f(x)dX} = (13.6/p0)^2*x0*(1+0.036*ln(x0+1))
1048           // x0+1 above is to make the result infrared safe.
1049           double x0 = fabs(svCurrent.radPath());
1050           double alphaX0 = 13.6e-3 * p0Inv;
1051           alphaX0 *= alphaX0;
1052           double betaX0 = 0.038;
1053           double logx0p1 = log(x0 + 1);
1054           theta02 = dX0 * alphaX0 * (1 + betaX0 * logx0p1) * (1 + betaX0 * logx0p1 + 2. * betaX0 * x0 / (x0 + 1));
1055         } else {
1056           theta02 = 196e-6 * p0Inv * p0Inv *
1057                     dX0;  //14.e-3/p0*sqrt(fabs(dS)/radX0); // .. drop log term (this is non-additive)
1058         }
1059 
1060         double epsilonP0 = 1. + dP / (p0 - 0.5 * dP);
1061         //      double omegaP0 = -dP/(p0-0.5*dP) + dS*dEdXPrime;
1062         //      double dsp = dS/(p0-0.5*dP); //use the initial p0 (not the mid-point) to keep the transport properly additive
1063 
1064         Vector tbtVec(bHat - tauB * tau);  // for b||z tau.z()*(-tau.x(), -tau.y(), 1.-tau.z())
1065 
1066         {
1067           //Slightly modified copy of the curvilinear jacobian (don't use the original just because it's in float precision
1068           // and seems to have some assumptions about the field values
1069           // notation changes: p1--> tau, p2-->tauNext
1070           // theta --> phi
1071           //    Vector p1 = tau;
1072           //    Vector p2 = tauNext;
1073           Point xStart = svCurrent.r3;
1074           const Vector& dx = drVec;
1075           //GlobalVector h  = MagneticField::inInverseGeV(xStart);
1076           // Martijn: field is now given as parameter.. GlobalVector h  = globalParameters.magneticFieldInInverseGeV(xStart);
1077 
1078           //double qbp = fts.signedInverseMomentum();
1079           double qbp = svCurrent.q * p0Inv;
1080           //    double absS = dS;
1081 
1082           // calculate transport matrix
1083           // Origin: TRPRFN
1084           double t11 = tau.x();
1085           double t12 = tau.y();
1086           double t13 = tau.z();
1087           double t21 = tauNext.x();
1088           double t22 = tauNext.y();
1089           double t23 = tauNext.z();
1090           double cosl0 = tau.perp();
1091           //    double cosl1 = 1./tauNext.perp(); //not quite a cos .. it's a cosec--> change to cosecl1 below
1092           double cosecl1 = tauNextPerpInv;
1093           //AlgebraicMatrix a(5,5,1);
1094           // define average magnetic field and gradient
1095           // at initial point - inlike TRPRFN
1096           const Vector& hn = bHat;
1097           //    double qp = -2.99792458e-3*b0;
1098           //   double q = -h.mag()*qbp;
1099 
1100           double q = -phi / dS;  //qp*qbp; // -phi/dS
1101           //    double theta = -phi;
1102           double sint = -sinPhi;
1103           double cost = cosPhi;
1104           double hn1 = hn.x();
1105           double hn2 = hn.y();
1106           double hn3 = hn.z();
1107           double dx1 = dx.x();
1108           double dx2 = dx.y();
1109           double dx3 = dx.z();
1110           //    double hnDt1 = hn1*t11 + hn2*t12 + hn3*t13;
1111 
1112           double gamma = hn1 * t21 + hn2 * t22 + hn3 * t23;
1113           double an1 = hn2 * t23 - hn3 * t22;
1114           double an2 = hn3 * t21 - hn1 * t23;
1115           double an3 = hn1 * t22 - hn2 * t21;
1116           //      double auInv = sqrt(1.- t13*t13); double au = auInv>0 ? 1./auInv : 1e24;
1117           double auInv = cosl0;
1118           double au = auInv > 0 ? 1. / auInv : 1e24;
1119           //      double auInv = sqrt(t11*t11 + t12*t12); double au = auInv>0 ? 1./auInv : 1e24;
1120           double u11 = -au * t12;
1121           double u12 = au * t11;
1122           double v11 = -t13 * u12;
1123           double v12 = t13 * u11;
1124           double v13 = auInv;  //t11*u12 - t12*u11;
1125           auInv = sqrt(1. - t23 * t23);
1126           au = auInv > 0 ? 1. / auInv : 1e24;
1127           //      auInv = sqrt(t21*t21 + t22*t22); au = auInv>0 ? 1./auInv : 1e24;
1128           double u21 = -au * t22;
1129           double u22 = au * t21;
1130           double v21 = -t23 * u22;
1131           double v22 = t23 * u21;
1132           double v23 = auInv;  //t21*u22 - t22*u21;
1133           // now prepare the transport matrix
1134           // pp only needed in high-p case (WA)
1135           //   double pp = 1./qbp;
1136           ////    double pp = fts.momentum().mag();
1137           // moved up (where -h.mag() is needed()
1138           //   double qp = q*pp;
1139           double anv = -(hn1 * u21 + hn2 * u22);
1140           double anu = (hn1 * v21 + hn2 * v22 + hn3 * v23);
1141           double omcost = oneLessCosPhi;
1142           double tmsint = -phi * phiLessSinPhiOPhi;
1143 
1144           double hu1 = -hn3 * u12;
1145           double hu2 = hn3 * u11;
1146           double hu3 = hn1 * u12 - hn2 * u11;
1147 
1148           double hv1 = hn2 * v13 - hn3 * v12;
1149           double hv2 = hn3 * v11 - hn1 * v13;
1150           double hv3 = hn1 * v12 - hn2 * v11;
1151 
1152           //   1/p - doesn't change since |tau| = |tauNext| ... not. It changes now
1153           dCCurvTransform(0, 0) = 1. / (epsilonP0 * epsilonP0) * (1. + dS * dEdXPrime);
1154 
1155           //   lambda
1156 
1157           dCCurvTransform(1, 0) = phi * p0 / svCurrent.q * cosecl1 * (sinPhi * bbtVec.z() - cosPhi * btVec.z());
1158           //was dCCurvTransform(1,0) = -qp*anv*(t21*dx1 + t22*dx2 + t23*dx3); //NOTE (SK) this was found to have an opposite sign
1159           //from independent re-calculation ... in fact the tauNext.dot.dR piece isnt reproduced
1160 
1161           dCCurvTransform(1, 1) =
1162               cost * (v11 * v21 + v12 * v22 + v13 * v23) + sint * (hv1 * v21 + hv2 * v22 + hv3 * v23) +
1163               omcost * (hn1 * v11 + hn2 * v12 + hn3 * v13) * (hn1 * v21 + hn2 * v22 + hn3 * v23) +
1164               anv * (-sint * (v11 * t21 + v12 * t22 + v13 * t23) + omcost * (v11 * an1 + v12 * an2 + v13 * an3) -
1165                      tmsint * gamma * (hn1 * v11 + hn2 * v12 + hn3 * v13));
1166 
1167           dCCurvTransform(1, 2) = cost * (u11 * v21 + u12 * v22) + sint * (hu1 * v21 + hu2 * v22 + hu3 * v23) +
1168                                   omcost * (hn1 * u11 + hn2 * u12) * (hn1 * v21 + hn2 * v22 + hn3 * v23) +
1169                                   anv * (-sint * (u11 * t21 + u12 * t22) + omcost * (u11 * an1 + u12 * an2) -
1170                                          tmsint * gamma * (hn1 * u11 + hn2 * u12));
1171           dCCurvTransform(1, 2) *= cosl0;
1172 
1173           // Commented out in part for reproducibility purposes: these terms are zero in cart->curv
1174           //    dCCurvTransform(1,3) = -q*anv*(u11*t21 + u12*t22          ); //don't show up in cartesian setup-->curv
1175           //why would lambdaNext depend explicitely on initial position ? any arbitrary init point can be chosen not
1176           // affecting the final state's momentum direction ... is this the field gradient in curvilinear coord?
1177           //    dCCurvTransform(1,4) = -q*anv*(v11*t21 + v12*t22 + v13*t23); //don't show up in cartesian setup-->curv
1178 
1179           //   phi
1180 
1181           dCCurvTransform(2, 0) = -phi * p0 / svCurrent.q * cosecl1 * cosecl1 *
1182                                   (oneLessCosPhi * bHat.z() * btVec.mag2() + sinPhi * btVec.z() + cosPhi * tbtVec.z());
1183           //was     dCCurvTransform(2,0) = -qp*anu*(t21*dx1 + t22*dx2 + t23*dx3)*cosecl1;
1184 
1185           dCCurvTransform(2, 1) =
1186               cost * (v11 * u21 + v12 * u22) + sint * (hv1 * u21 + hv2 * u22) +
1187               omcost * (hn1 * v11 + hn2 * v12 + hn3 * v13) * (hn1 * u21 + hn2 * u22) +
1188               anu * (-sint * (v11 * t21 + v12 * t22 + v13 * t23) + omcost * (v11 * an1 + v12 * an2 + v13 * an3) -
1189                      tmsint * gamma * (hn1 * v11 + hn2 * v12 + hn3 * v13));
1190           dCCurvTransform(2, 1) *= cosecl1;
1191 
1192           dCCurvTransform(2, 2) = cost * (u11 * u21 + u12 * u22) + sint * (hu1 * u21 + hu2 * u22) +
1193                                   omcost * (hn1 * u11 + hn2 * u12) * (hn1 * u21 + hn2 * u22) +
1194                                   anu * (-sint * (u11 * t21 + u12 * t22) + omcost * (u11 * an1 + u12 * an2) -
1195                                          tmsint * gamma * (hn1 * u11 + hn2 * u12));
1196           dCCurvTransform(2, 2) *= cosecl1 * cosl0;
1197 
1198           // Commented out in part for reproducibility purposes: these terms are zero in cart->curv
1199           // dCCurvTransform(2,3) = -q*anu*(u11*t21 + u12*t22          )*cosecl1;
1200           //why would lambdaNext depend explicitely on initial position ? any arbitrary init point can be chosen not
1201           // affecting the final state's momentum direction ... is this the field gradient in curvilinear coord?
1202           // dCCurvTransform(2,4) = -q*anu*(v11*t21 + v12*t22 + v13*t23)*cosecl1;
1203 
1204           //   yt
1205 
1206           double pp = 1. / qbp;
1207           // (SK) these terms seem to consistently have a sign opp from private derivation
1208           dCCurvTransform(3, 0) = -pp * (u21 * dx1 + u22 * dx2);  //NB: modified from the original: changed the sign
1209           dCCurvTransform(4, 0) = -pp * (v21 * dx1 + v22 * dx2 + v23 * dx3);
1210 
1211           dCCurvTransform(3, 1) = (sint * (v11 * u21 + v12 * u22) + omcost * (hv1 * u21 + hv2 * u22) +
1212                                    tmsint * (hn1 * u21 + hn2 * u22) * (hn1 * v11 + hn2 * v12 + hn3 * v13)) /
1213                                   q;
1214 
1215           dCCurvTransform(3, 2) = (sint * (u11 * u21 + u12 * u22) + omcost * (hu1 * u21 + hu2 * u22) +
1216                                    tmsint * (hn1 * u21 + hn2 * u22) * (hn1 * u11 + hn2 * u12)) *
1217                                   cosl0 / q;
1218 
1219           dCCurvTransform(3, 3) = (u11 * u21 + u12 * u22);
1220 
1221           dCCurvTransform(3, 4) = (v11 * u21 + v12 * u22);
1222 
1223           //   zt
1224 
1225           dCCurvTransform(4, 1) =
1226               (sint * (v11 * v21 + v12 * v22 + v13 * v23) + omcost * (hv1 * v21 + hv2 * v22 + hv3 * v23) +
1227                tmsint * (hn1 * v21 + hn2 * v22 + hn3 * v23) * (hn1 * v11 + hn2 * v12 + hn3 * v13)) /
1228               q;
1229 
1230           dCCurvTransform(4, 2) = (sint * (u11 * v21 + u12 * v22) + omcost * (hu1 * v21 + hu2 * v22 + hu3 * v23) +
1231                                    tmsint * (hn1 * v21 + hn2 * v22 + hn3 * v23) * (hn1 * u11 + hn2 * u12)) *
1232                                   cosl0 / q;
1233 
1234           dCCurvTransform(4, 3) = (u11 * v21 + u12 * v22);
1235 
1236           dCCurvTransform(4, 4) = (v11 * v21 + v12 * v22 + v13 * v23);
1237           // end of TRPRFN
1238         }
1239 
1240         if (debug_) {
1241           Basis rep;
1242           setRep(rep, tauNext);
1243           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "rep X: " << rep.lX << " "
1244                             << rep.lX.mag() << "\t Y: " << rep.lY << " " << rep.lY.mag() << "\t Z: " << rep.lZ << " "
1245                             << rep.lZ.mag();
1246         }
1247         //mind the sign of dS and dP (dS*dP < 0 allways)
1248         //covariance should grow no matter which direction you propagate
1249         //==> take abs values.
1250         //reset not needed: fill all below  svCurrent.matDCov *= 0.;
1251         double mulRR = theta02 * dS * dS / 3.;
1252         double mulRP = theta02 * fabs(dS) * p0 / 2.;
1253         double mulPP = theta02 * p0 * p0;
1254         double losPP = dP * dP * 1.6 / fabs(dS) * (1.0 + p0 * 1e-3);
1255         //another guess .. makes sense for 1 cm steps 2./dS == 2 [cm] / dS [cm] at low pt
1256         //double it by 1TeV
1257         //not gaussian anyways
1258         // derived from the fact that sigma_p/eLoss ~ 0.08 after ~ 200 steps
1259 
1260         //curvilinear
1261         double sinThetaInv = tauNextPerpInv;
1262         double p0Mat =
1263             p0 + 0.5 * dP;  // FIXME change this to p0 after it's clear that there's agreement in everything else
1264         double p0Mat2 = p0Mat * p0Mat;
1265         // with 6x6 formulation
1266         svCurrent.matDCovCurv *= 0;
1267 
1268         svCurrent.matDCovCurv(0, 0) = losPP / (p0Mat2 * p0Mat2);
1269         //      svCurrent.matDCovCurv(0,1) = 0;
1270         //      svCurrent.matDCovCurv(0,2) = 0;
1271         //      svCurrent.matDCovCurv(0,3) = 0;
1272         //      svCurrent.matDCovCurv(0,4) = 0;
1273 
1274         //      svCurrent.matDCovCurv(1,0) = 0;
1275         svCurrent.matDCovCurv(1, 1) = mulPP / p0Mat2;
1276         //      svCurrent.matDCovCurv(1,2) = 0;
1277         //      svCurrent.matDCovCurv(1,3) = 0;
1278         svCurrent.matDCovCurv(1, 4) = mulRP / p0Mat;
1279 
1280         //      svCurrent.matDCovCurv(2,0) = 0;
1281         //      svCurrent.matDCovCurv(2,1) = 0;
1282         svCurrent.matDCovCurv(2, 2) = mulPP / p0Mat2 * (sinThetaInv * sinThetaInv);
1283         svCurrent.matDCovCurv(2, 3) = mulRP / p0Mat * sinThetaInv;
1284         //      svCurrent.matDCovCurv(2,4) = 0;
1285 
1286         //      svCurrent.matDCovCurv(3,0) = 0;
1287         //      svCurrent.matDCovCurv(3,1) = 0;
1288         svCurrent.matDCovCurv(3, 2) = mulRP / p0Mat * sinThetaInv;
1289         //      svCurrent.matDCovCurv(3,0) = 0;
1290         svCurrent.matDCovCurv(3, 3) = mulRR;
1291         //      svCurrent.matDCovCurv(3,4) = 0;
1292 
1293         //      svCurrent.matDCovCurv(4,0) = 0;
1294         svCurrent.matDCovCurv(4, 1) = mulRP / p0Mat;
1295         //      svCurrent.matDCovCurv(4,2) = 0;
1296         //      svCurrent.matDCovCurv(4,3) = 0;
1297         svCurrent.matDCovCurv(4, 4) = mulRR;
1298       }
1299       break;
1300     }
1301     default:
1302       break;
1303   }
1304 
1305   double pMag = curP;  //svCurrent.p3.mag();
1306 
1307   if (dir == oppositeToMomentum)
1308     dP = fabs(dP);
1309   else if (dP < 0) {  //the case of negative dP
1310     dP = -dP > pMag ? -pMag + 1e-5 : dP;
1311   }
1312 
1313   getNextState(svCurrent, svNext, dP, tauNext, drVec, dS, dS / radX0, dCCurvTransform);
1314   return true;
1315 }
1316 
1317 double SteppingHelixPropagator::getDeDx(const SteppingHelixPropagator::StateInfo& sv,
1318                                         double& dEdXPrime,
1319                                         double& radX0,
1320                                         MatBounds& rzLims) const {
1321   radX0 = 1.e24;
1322   dEdXPrime = 0.;
1323   rzLims = MatBounds();
1324   if (noMaterialMode_)
1325     return 0;
1326 
1327   double dEdx = 0.;
1328 
1329   double lR = sv.r3.perp();
1330   double lZ = fabs(sv.r3.z());
1331 
1332   //assume "Iron" .. seems to be quite the same for brass/iron/PbW04
1333   //good for Fe within 3% for 0.2 GeV to 10PeV
1334 
1335   double dEdX_HCal = 0.95;  //extracted from sim
1336   double dEdX_ECal = 0.45;
1337   double dEdX_coil = 0.35;  //extracted from sim .. closer to 40% in fact
1338   double dEdX_Fe = 1;
1339   double dEdX_MCh = 0.053;  //chambers on average
1340   double dEdX_Trk = 0.0114;
1341   double dEdX_Air = 2E-4;
1342   double dEdX_Vac = 0.0;
1343 
1344   double radX0_HCal = 1.44 / 0.8;  //guessing
1345   double radX0_ECal = 0.89 / 0.7;
1346   double radX0_coil = 4.;  //
1347   double radX0_Fe = 1.76;
1348   double radX0_MCh = 1e3;  //
1349   double radX0_Trk = 320.;
1350   double radX0_Air = 3.e4;
1351   double radX0_Vac = 3.e9;  //"big" number for vacuum
1352 
1353   //not all the boundaries are set below: this will be a default
1354   if (!(lR < 380 && lZ < 785)) {
1355     if (lZ > 785)
1356       rzLims = MatBounds(0, 1e4, 785, 1e4);
1357     if (lZ < 785)
1358       rzLims = MatBounds(380, 1e4, 0, 785);
1359   }
1360 
1361   //this def makes sense assuming we don't jump into endcap volume from the other z-side in one step
1362   //also, it is a positive shift only (at least for now): can't move ec further into the detector
1363   double ecShift = sv.r3.z() > 0 ? fabs(ecShiftPos_) : fabs(ecShiftNeg_);
1364 
1365   //this should roughly figure out where things are
1366   //(numbers taken from Fig1.1.2 TDR and from geom xmls)
1367   if (lR < 2.9) {  //inside beampipe
1368     dEdx = dEdX_Vac;
1369     radX0 = radX0_Vac;
1370     rzLims = MatBounds(0, 2.9, 0, 1E4);
1371   } else if (lR < 129) {
1372     if (lZ < 294) {
1373       rzLims = MatBounds(2.9, 129, 0, 294);  //tracker boundaries
1374       dEdx = dEdX_Trk;
1375       radX0 = radX0_Trk;
1376       //somewhat empirical formula that ~ matches the average if going from 0,0,0
1377       //assuming "uniform" tracker material
1378       //doesn't really track material layer to layer
1379       double lEtaDet = fabs(sv.r3.eta());
1380       double scaleRadX = lEtaDet > 1.5 ? 0.7724 : 1. / cosh(0.5 * lEtaDet);  //sin(2.*atan(exp(-0.5*lEtaDet)));
1381       scaleRadX *= scaleRadX;
1382       if (lEtaDet > 2 && lZ > 20)
1383         scaleRadX *= (lEtaDet - 1.);
1384       if (lEtaDet > 2.5 && lZ > 20)
1385         scaleRadX *= (lEtaDet - 1.);
1386       radX0 *= scaleRadX;
1387     }
1388     //endcap part begins here
1389     else if (lZ < 294 + ecShift) {
1390       //gap in front of EE here, piece inside 2.9<R<129
1391       rzLims = MatBounds(2.9, 129, 294, 294 + ecShift);
1392       dEdx = dEdX_Air;
1393       radX0 = radX0_Air;
1394     } else if (lZ < 372 + ecShift) {
1395       rzLims = MatBounds(2.9, 129, 294 + ecShift, 372 + ecShift);  //EE here, piece inside 2.9<R<129
1396       dEdx = dEdX_ECal;
1397       radX0 = radX0_ECal;
1398     }  //EE averaged out over a larger space
1399     else if (lZ < 398 + ecShift) {
1400       rzLims =
1401           MatBounds(2.9, 129, 372 + ecShift, 398 + ecShift);  //whatever goes behind EE 2.9<R<129 is air up to Z=398
1402       dEdx = dEdX_HCal * 0.05;
1403       radX0 = radX0_Air;
1404     }  //betw EE and HE
1405     else if (lZ < 555 + ecShift) {
1406       rzLims = MatBounds(2.9, 129, 398 + ecShift, 555 + ecShift);  //HE piece 2.9<R<129;
1407       dEdx = dEdX_HCal * 0.96;
1408       radX0 = radX0_HCal / 0.96;
1409     }  //HE calor abit less dense
1410     else {
1411       //      rzLims = MatBounds(2.9, 129, 555, 785);
1412       // set the boundaries first: they serve as stop-points here
1413       // the material is set below
1414       if (lZ < 568 + ecShift)
1415         rzLims = MatBounds(2.9, 129, 555 + ecShift, 568 + ecShift);  //a piece of HE support R<129, 555<Z<568
1416       else if (lZ < 625 + ecShift) {
1417         if (lR < 85 + ecShift)
1418           rzLims = MatBounds(2.9, 85, 568 + ecShift, 625 + ecShift);
1419         else
1420           rzLims = MatBounds(85, 129, 568 + ecShift, 625 + ecShift);
1421       } else if (lZ < 785 + ecShift)
1422         rzLims = MatBounds(2.9, 129, 625 + ecShift, 785 + ecShift);
1423       else if (lZ < 1500 + ecShift)
1424         rzLims = MatBounds(2.9, 129, 785 + ecShift, 1500 + ecShift);
1425       else
1426         rzLims = MatBounds(2.9, 129, 1500 + ecShift, 1E4);
1427 
1428       //iron .. don't care about no material in front of HF (too forward)
1429       if (!(lZ > 568 + ecShift && lZ < 625 + ecShift && lR > 85)  // HE support
1430           && !(lZ > 785 + ecShift && lZ < 850 + ecShift && lR > 118)) {
1431         dEdx = dEdX_Fe;
1432         radX0 = radX0_Fe;
1433       } else {
1434         dEdx = dEdX_MCh;
1435         radX0 = radX0_MCh;
1436       }  //ME at eta > 2.2
1437     }
1438   } else if (lR < 287) {
1439     if (lZ < 372 + ecShift && lR < 177) {  // 129<<R<177
1440       if (lZ < 304)
1441         rzLims = MatBounds(129, 177, 0, 304);  //EB  129<<R<177 0<Z<304
1442       else if (lZ < 343) {                     // 129<<R<177 304<Z<343
1443         if (lR < 135)
1444           rzLims = MatBounds(129, 135, 304, 343);  // tk piece 129<<R<135 304<Z<343
1445         else if (lR < 172) {                       //
1446           if (lZ < 311)
1447             rzLims = MatBounds(135, 172, 304, 311);
1448           else
1449             rzLims = MatBounds(135, 172, 311, 343);
1450         } else {
1451           if (lZ < 328)
1452             rzLims = MatBounds(172, 177, 304, 328);
1453           else
1454             rzLims = MatBounds(172, 177, 328, 343);
1455         }
1456       } else if (lZ < 343 + ecShift) {
1457         rzLims = MatBounds(129, 177, 343, 343 + ecShift);  //gap
1458       } else {
1459         if (lR < 156)
1460           rzLims = MatBounds(129, 156, 343 + ecShift, 372 + ecShift);
1461         else if ((lZ - 343 - ecShift) > (lR - 156) * 1.38)
1462           rzLims = MatBounds(156, 177, 127.72 + ecShift, 372 + ecShift, atan(1.38), 0);
1463         else
1464           rzLims = MatBounds(156, 177, 343 + ecShift, 127.72 + ecShift, 0, atan(1.38));
1465       }
1466 
1467       if (!(lR > 135 && lZ < 343 + ecShift && lZ > 304) &&
1468           !(lR > 156 && lZ < 372 + ecShift && lZ > 343 + ecShift && ((lZ - 343. - ecShift) < (lR - 156.) * 1.38))) {
1469         //the crystals are the same length, but they are not 100% of material
1470         double cosThetaEquiv = 0.8 / sqrt(1. + lZ * lZ / lR / lR) + 0.2;
1471         if (lZ > 343)
1472           cosThetaEquiv = 1.;
1473         dEdx = dEdX_ECal * cosThetaEquiv;
1474         radX0 = radX0_ECal / cosThetaEquiv;
1475       }  //EB
1476       else {
1477         if ((lZ > 304 && lZ < 328 && lR < 177 && lR > 135) && !(lZ > 311 && lR < 172)) {
1478           dEdx = dEdX_Fe;
1479           radX0 = radX0_Fe;
1480         }  //Tk_Support
1481         else if (lZ > 343 && lZ < 343 + ecShift) {
1482           dEdx = dEdX_Air;
1483           radX0 = radX0_Air;
1484         } else {
1485           dEdx = dEdX_ECal * 0.2;
1486           radX0 = radX0_Air;
1487         }  //cables go here <-- will be abit too dense for ecShift > 0
1488       }
1489     } else if (lZ < 554 + ecShift) {  // 129<R<177 372<Z<554  AND 177<R<287 0<Z<554
1490       if (lR < 177) {                 //  129<R<177 372<Z<554
1491         if (lZ > 372 + ecShift && lZ < 398 + ecShift)
1492           rzLims = MatBounds(129, 177, 372 + ecShift, 398 + ecShift);  // EE 129<R<177 372<Z<398
1493         else if (lZ < 548 + ecShift)
1494           rzLims = MatBounds(129, 177, 398 + ecShift, 548 + ecShift);  // HE 129<R<177 398<Z<548
1495         else
1496           rzLims = MatBounds(129, 177, 548 + ecShift, 554 + ecShift);  // HE gap 129<R<177 548<Z<554
1497       } else if (lR < 193) {                                           // 177<R<193 0<Z<554
1498         if ((lZ - 307) < (lR - 177.) * 1.739)
1499           rzLims = MatBounds(177, 193, 0, -0.803, 0, atan(1.739));
1500         else if (lZ < 389)
1501           rzLims = MatBounds(177, 193, -0.803, 389, atan(1.739), 0.);
1502         else if (lZ < 389 + ecShift)
1503           rzLims = MatBounds(177, 193, 389, 389 + ecShift);  // air insert
1504         else if (lZ < 548 + ecShift)
1505           rzLims = MatBounds(177, 193, 389 + ecShift, 548 + ecShift);  // HE 177<R<193 389<Z<548
1506         else
1507           rzLims = MatBounds(177, 193, 548 + ecShift, 554 + ecShift);  // HE gap 177<R<193 548<Z<554
1508       } else if (lR < 264) {                                           // 193<R<264 0<Z<554
1509         double anApex = 375.7278 - 193. / 1.327;                       // 230.28695599096
1510         if ((lZ - 375.7278) < (lR - 193.) / 1.327)
1511           rzLims = MatBounds(193, 264, 0, anApex, 0, atan(1. / 1.327));  //HB
1512         else if ((lZ - 392.7278) < (lR - 193.) / 1.327)
1513           rzLims = MatBounds(193, 264, anApex, anApex + 17., atan(1. / 1.327), atan(1. / 1.327));  // HB-HE gap
1514         else if ((lZ - 392.7278 - ecShift) < (lR - 193.) / 1.327)
1515           rzLims = MatBounds(193,
1516                              264,
1517                              anApex + 17.,
1518                              anApex + 17. + ecShift,
1519                              atan(1. / 1.327),
1520                              atan(1. / 1.327));  // HB-HE gap air insert
1521         // HE (372,193)-(517,193)-(517,264)-(417.8,264)
1522         else if (lZ < 517 + ecShift)
1523           rzLims = MatBounds(193, 264, anApex + 17. + ecShift, 517 + ecShift, atan(1. / 1.327), 0);
1524         else if (lZ < 548 + ecShift) {  // HE+gap 193<R<264 517<Z<548
1525           if (lR < 246)
1526             rzLims = MatBounds(193, 246, 517 + ecShift, 548 + ecShift);  // HE 193<R<246 517<Z<548
1527           else
1528             rzLims = MatBounds(246, 264, 517 + ecShift, 548 + ecShift);  // HE gap 246<R<264 517<Z<548
1529         } else
1530           rzLims = MatBounds(193, 264, 548 + ecShift, 554 + ecShift);  // HE gap  193<R<246 548<Z<554
1531       } else if (lR < 275) {                                           // 264<R<275 0<Z<554
1532         if (lZ < 433)
1533           rzLims = MatBounds(264, 275, 0, 433);  //HB
1534         else if (lZ < 554)
1535           rzLims = MatBounds(264, 275, 433, 554);  // HE gap 264<R<275 433<Z<554
1536         else
1537           rzLims = MatBounds(264, 275, 554, 554 + ecShift);  // HE gap 264<R<275 554<Z<554 air insert
1538       } else {                                               // 275<R<287 0<Z<554
1539         if (lZ < 402)
1540           rzLims = MatBounds(275, 287, 0, 402);  // HB 275<R<287 0<Z<402
1541         else if (lZ < 554)
1542           rzLims = MatBounds(275, 287, 402, 554);  //  //HE gap 275<R<287 402<Z<554
1543         else
1544           rzLims = MatBounds(275, 287, 554, 554 + ecShift);  //HE gap 275<R<287 554<Z<554 air insert
1545       }
1546 
1547       if ((lZ < 433 || lR < 264) && (lZ < 402 || lR < 275) &&
1548           (lZ < 517 + ecShift || lR < 246)  //notches
1549           //I should've made HE and HF different .. now need to shorten HE to match
1550           && lZ < 548 + ecShift && !(lZ < 389 + ecShift && lZ > 335 && lR < 193)        //not a gap EB-EE 129<R<193
1551           && !(lZ > 307 && lZ < 335 && lR < 193 && ((lZ - 307) > (lR - 177.) * 1.739))  //not a gap
1552           && !(lR < 177 && lZ < 398 + ecShift)                                          //under the HE nose
1553           && !(lR < 264 && lR > 193 &&
1554                fabs(441.5 + 0.5 * ecShift - lZ + (lR - 269.) / 1.327) < (8.5 + ecShift * 0.5)))  //not a gap
1555       {
1556         dEdx = dEdX_HCal;
1557         radX0 = radX0_HCal;  //hcal
1558       } else if ((lR < 193 && lZ > 389 && lZ < 389 + ecShift) ||
1559                  (lR > 264 && lR < 287 && lZ > 554 && lZ < 554 + ecShift) ||
1560                  (lR < 264 && lR > 193 &&
1561                   fabs(441.5 + 8.5 + 0.5 * ecShift - lZ + (lR - 269.) / 1.327) < ecShift * 0.5)) {
1562         dEdx = dEdX_Air;
1563         radX0 = radX0_Air;
1564       } else {
1565         dEdx = dEdX_HCal * 0.05;
1566         radX0 = radX0_Air;
1567       }  //endcap gap
1568     }
1569     //the rest is a tube of endcap volume  129<R<251 554<Z<largeValue
1570     else if (lZ < 564 + ecShift) {  // 129<R<287  554<Z<564
1571       if (lR < 251) {
1572         rzLims = MatBounds(129, 251, 554 + ecShift, 564 + ecShift);  // HE support 129<R<251  554<Z<564
1573         dEdx = dEdX_Fe;
1574         radX0 = radX0_Fe;
1575       }  //HE support
1576       else {
1577         rzLims = MatBounds(251, 287, 554 + ecShift, 564 + ecShift);  //HE/ME gap 251<R<287  554<Z<564
1578         dEdx = dEdX_MCh;
1579         radX0 = radX0_MCh;
1580       }
1581     } else if (lZ < 625 + ecShift) {  // ME/1/1 129<R<287  564<Z<625
1582       rzLims = MatBounds(129, 287, 564 + ecShift, 625 + ecShift);
1583       dEdx = dEdX_MCh;
1584       radX0 = radX0_MCh;
1585     } else if (lZ < 785 + ecShift) {  //129<R<287  625<Z<785
1586       if (lR < 275)
1587         rzLims = MatBounds(129, 275, 625 + ecShift, 785 + ecShift);  //YE/1 129<R<275 625<Z<785
1588       else {                                                         // 275<R<287  625<Z<785
1589         if (lZ < 720 + ecShift)
1590           rzLims = MatBounds(275, 287, 625 + ecShift, 720 + ecShift);  // ME/1/2 275<R<287  625<Z<720
1591         else
1592           rzLims = MatBounds(275, 287, 720 + ecShift, 785 + ecShift);  // YE/1 275<R<287  720<Z<785
1593       }
1594       if (!(lR > 275 && lZ < 720 + ecShift)) {
1595         dEdx = dEdX_Fe;
1596         radX0 = radX0_Fe;
1597       }  //iron
1598       else {
1599         dEdx = dEdX_MCh;
1600         radX0 = radX0_MCh;
1601       }
1602     } else if (lZ < 850 + ecShift) {
1603       rzLims = MatBounds(129, 287, 785 + ecShift, 850 + ecShift);
1604       dEdx = dEdX_MCh;
1605       radX0 = radX0_MCh;
1606     } else if (lZ < 910 + ecShift) {
1607       rzLims = MatBounds(129, 287, 850 + ecShift, 910 + ecShift);
1608       dEdx = dEdX_Fe;
1609       radX0 = radX0_Fe;
1610     }  //iron
1611     else if (lZ < 975 + ecShift) {
1612       rzLims = MatBounds(129, 287, 910 + ecShift, 975 + ecShift);
1613       dEdx = dEdX_MCh;
1614       radX0 = radX0_MCh;
1615     } else if (lZ < 1000 + ecShift) {
1616       rzLims = MatBounds(129, 287, 975 + ecShift, 1000 + ecShift);
1617       dEdx = dEdX_Fe;
1618       radX0 = radX0_Fe;
1619     }  //iron
1620     else if (lZ < 1063 + ecShift) {
1621       rzLims = MatBounds(129, 287, 1000 + ecShift, 1063 + ecShift);
1622       dEdx = dEdX_MCh;
1623       radX0 = radX0_MCh;
1624     } else if (lZ < 1073 + ecShift) {
1625       rzLims = MatBounds(129, 287, 1063 + ecShift, 1073 + ecShift);
1626       dEdx = dEdX_Fe;
1627       radX0 = radX0_Fe;
1628     } else if (lZ < 1E4) {
1629       rzLims = MatBounds(129, 287, 1073 + ecShift, 1E4);
1630       dEdx = dEdX_Air;
1631       radX0 = radX0_Air;
1632     } else {
1633       dEdx = dEdX_Air;
1634       radX0 = radX0_Air;
1635     }
1636   } else if (lR < 380 && lZ < 667) {
1637     if (lZ < 630) {
1638       if (lR < 315)
1639         rzLims = MatBounds(287, 315, 0, 630);
1640       else if (lR < 341)
1641         rzLims = MatBounds(315, 341, 0, 630);  //b-field ~linear rapid fall here
1642       else
1643         rzLims = MatBounds(341, 380, 0, 630);
1644     } else
1645       rzLims = MatBounds(287, 380, 630, 667);
1646 
1647     if (lZ < 630) {
1648       dEdx = dEdX_coil;
1649       radX0 = radX0_coil;
1650     }  //a guess for the solenoid average
1651     else {
1652       dEdx = dEdX_Air;
1653       radX0 = radX0_Air;
1654     }  //endcap gap
1655   } else {
1656     if (lZ < 667) {
1657       if (lR < 850) {
1658         bool isIron = false;
1659         //sanity check in addition to flags
1660         if (useIsYokeFlag_ && useMagVolumes_ && sv.magVol != nullptr) {
1661           isIron = sv.isYokeVol;
1662         } else {
1663           double bMag = sv.bf.mag();
1664           isIron = (bMag > 0.75 && !(lZ > 500 && lR < 500 && bMag < 1.15) && !(lZ < 450 && lR > 420 && bMag < 1.15));
1665         }
1666         //tell the material stepper where mat bounds are
1667         rzLims = MatBounds(380, 850, 0, 667);
1668         if (isIron) {
1669           dEdx = dEdX_Fe;
1670           radX0 = radX0_Fe;
1671         }  //iron
1672         else {
1673           dEdx = dEdX_MCh;
1674           radX0 = radX0_MCh;
1675         }
1676       } else {
1677         rzLims = MatBounds(850, 1E4, 0, 667);
1678         dEdx = dEdX_Air;
1679         radX0 = radX0_Air;
1680       }
1681     } else if (lR > 750) {
1682       rzLims = MatBounds(750, 1E4, 667, 1E4);
1683       dEdx = dEdX_Air;
1684       radX0 = radX0_Air;
1685     } else if (lZ < 667 + ecShift) {
1686       rzLims = MatBounds(287, 750, 667, 667 + ecShift);
1687       dEdx = dEdX_Air;
1688       radX0 = radX0_Air;
1689     }
1690     //the rest is endcap piece with 287<R<750 Z>667
1691     else if (lZ < 724 + ecShift) {
1692       if (lR < 380)
1693         rzLims = MatBounds(287, 380, 667 + ecShift, 724 + ecShift);
1694       else
1695         rzLims = MatBounds(380, 750, 667 + ecShift, 724 + ecShift);
1696       dEdx = dEdX_MCh;
1697       radX0 = radX0_MCh;
1698     } else if (lZ < 785 + ecShift) {
1699       if (lR < 380)
1700         rzLims = MatBounds(287, 380, 724 + ecShift, 785 + ecShift);
1701       else
1702         rzLims = MatBounds(380, 750, 724 + ecShift, 785 + ecShift);
1703       dEdx = dEdX_Fe;
1704       radX0 = radX0_Fe;
1705     }  //iron
1706     else if (lZ < 850 + ecShift) {
1707       rzLims = MatBounds(287, 750, 785 + ecShift, 850 + ecShift);
1708       dEdx = dEdX_MCh;
1709       radX0 = radX0_MCh;
1710     } else if (lZ < 910 + ecShift) {
1711       rzLims = MatBounds(287, 750, 850 + ecShift, 910 + ecShift);
1712       dEdx = dEdX_Fe;
1713       radX0 = radX0_Fe;
1714     }  //iron
1715     else if (lZ < 975 + ecShift) {
1716       rzLims = MatBounds(287, 750, 910 + ecShift, 975 + ecShift);
1717       dEdx = dEdX_MCh;
1718       radX0 = radX0_MCh;
1719     } else if (lZ < 1000 + ecShift) {
1720       rzLims = MatBounds(287, 750, 975 + ecShift, 1000 + ecShift);
1721       dEdx = dEdX_Fe;
1722       radX0 = radX0_Fe;
1723     }  //iron
1724     else if (lZ < 1063 + ecShift) {
1725       if (lR < 360) {
1726         rzLims = MatBounds(287, 360, 1000 + ecShift, 1063 + ecShift);
1727         dEdx = dEdX_MCh;
1728         radX0 = radX0_MCh;
1729       }
1730       //put empty air where me4/2 should be (future)
1731       else {
1732         rzLims = MatBounds(360, 750, 1000 + ecShift, 1063 + ecShift);
1733         dEdx = dEdX_Air;
1734         radX0 = radX0_Air;
1735       }
1736     } else if (lZ < 1073 + ecShift) {
1737       rzLims = MatBounds(287, 750, 1063 + ecShift, 1073 + ecShift);
1738       //this plate does not exist: air
1739       dEdx = dEdX_Air;
1740       radX0 = radX0_Air;
1741     } else if (lZ < 1E4) {
1742       rzLims = MatBounds(287, 750, 1073 + ecShift, 1E4);
1743       dEdx = dEdX_Air;
1744       radX0 = radX0_Air;
1745     } else {
1746       dEdx = dEdX_Air;
1747       radX0 = radX0_Air;
1748     }  //air
1749   }
1750 
1751   //dEdx so far is a relative number (relative to iron)
1752   //scale by what's expected for iron (the function was fit from pdg table)
1753   //0.065 (PDG) --> 0.044 to better match with MPV
1754   double p0 = sv.p3.mag();
1755   double logp0 = log(p0);
1756   double p0powN33 = 0;
1757   if (p0 > 3.) {
1758     // p0powN33 = exp(-0.33*logp0); //calculate for p>3GeV
1759     double xx = 1. / p0;
1760     xx = sqrt(sqrt(sqrt(sqrt(xx))));
1761     xx = xx * xx * xx * xx * xx;  // this is (p0)**(-5/16), close enough to -0.33
1762     p0powN33 = xx;
1763   }
1764   double dEdX_mat = -(11.4 + 0.96 * fabs(logp0 + log(2.8)) + 0.033 * p0 * (1.0 - p0powN33)) * 1e-3;
1765   //in GeV/cm .. 0.8 to get closer to the median or MPV
1766 
1767   dEdXPrime = dEdx == 0 ? 0 : -dEdx * (0.96 / p0 + 0.033 - 0.022 * p0powN33) * 1e-3;  //== d(dEdX)/dp
1768   dEdx *= dEdX_mat;
1769 
1770   return dEdx;
1771 }
1772 
1773 int SteppingHelixPropagator::cIndex_(int ind) const {
1774   int result = ind % MAX_POINTS;
1775   if (ind != 0 && result == 0) {
1776     result = MAX_POINTS;
1777   }
1778   return result;
1779 }
1780 
1781 SteppingHelixPropagator::Result SteppingHelixPropagator::refToDest(SteppingHelixPropagator::DestType dest,
1782                                                                    const SteppingHelixPropagator::StateInfo& sv,
1783                                                                    const double pars[6],
1784                                                                    double& dist,
1785                                                                    double& tanDist,
1786                                                                    PropagationDirection& refDirection,
1787                                                                    double fastSkipDist) const {
1788   static const std::string metname = "SteppingHelixPropagator";
1789   Result result = SteppingHelixStateInfo::NOT_IMPLEMENTED;
1790 
1791   switch (dest) {
1792     case RADIUS_DT: {
1793       double curR = sv.r3.perp();
1794       dist = pars[RADIUS_P] - curR;
1795       if (fabs(dist) > fastSkipDist) {
1796         result = SteppingHelixStateInfo::INACC;
1797         break;
1798       }
1799       double curP2 = sv.p3.mag2();
1800       double curPtPos2 = sv.p3.perp2();
1801       if (curPtPos2 < 1e-16)
1802         curPtPos2 = 1e-16;
1803 
1804       double cosDPhiPR =
1805           (sv.r3.x() * sv.p3.x() + sv.r3.y() * sv.p3.y());  //only the sign is needed cos((sv.r3.deltaPhi(sv.p3)));
1806       refDirection = dist * cosDPhiPR > 0 ? alongMomentum : oppositeToMomentum;
1807       tanDist = dist * sqrt(curP2 / curPtPos2);
1808       result = SteppingHelixStateInfo::OK;
1809     } break;
1810     case Z_DT: {
1811       double curZ = sv.r3.z();
1812       dist = pars[Z_P] - curZ;
1813       if (fabs(dist) > fastSkipDist) {
1814         result = SteppingHelixStateInfo::INACC;
1815         break;
1816       }
1817       double curP = sv.p3.mag();
1818       refDirection = sv.p3.z() * dist > 0. ? alongMomentum : oppositeToMomentum;
1819       tanDist = dist / sv.p3.z() * curP;
1820       result = SteppingHelixStateInfo::OK;
1821     } break;
1822     case PLANE_DT: {
1823       Point rPlane(pars[0], pars[1], pars[2]);
1824       Vector nPlane(pars[3], pars[4], pars[5]);
1825 
1826       // unfortunately this doesn't work: the numbers are too large
1827       //      bool pVertical = fabs(pars[5])>0.9999;
1828       //      double dRDotN = pVertical? (sv.r3.z() - rPlane.z())*nPlane.z() :(sv.r3 - rPlane).dot(nPlane);
1829       double dRDotN = (sv.r3.x() - rPlane.x()) * nPlane.x() + (sv.r3.y() - rPlane.y()) * nPlane.y() +
1830                       (sv.r3.z() - rPlane.z()) * nPlane.z();  //(sv.r3 - rPlane).dot(nPlane);
1831 
1832       dist = fabs(dRDotN);
1833       if (dist > fastSkipDist) {
1834         result = SteppingHelixStateInfo::INACC;
1835         break;
1836       }
1837       double curP = sv.p3.mag();
1838       double p0 = curP;
1839       double p0Inv = 1. / p0;
1840       Vector tau(sv.p3);
1841       tau *= p0Inv;
1842       double tN = tau.dot(nPlane);
1843       refDirection = tN * dRDotN < 0. ? alongMomentum : oppositeToMomentum;
1844       double b0 = sv.bf.mag();
1845       if (fabs(tN) > 1e-24) {
1846         tanDist = -dRDotN / tN;
1847       } else {
1848         tN = 1e-24;
1849         if (fabs(dRDotN) > 1e-24)
1850           tanDist = 1e6;
1851         else
1852           tanDist = 1;
1853       }
1854       if (fabs(tanDist) > 1e4)
1855         tanDist = 1e4;
1856       if (b0 > 1.5e-6) {
1857         double b0Inv = 1. / b0;
1858         double tNInv = 1. / tN;
1859         Vector bHat(sv.bf);
1860         bHat *= b0Inv;
1861         double bHatN = bHat.dot(nPlane);
1862         double cosPB = bHat.dot(tau);
1863         double kVal = 2.99792458e-3 * sv.q * p0Inv * b0;
1864         double aVal = tanDist * kVal;
1865         Vector lVec = bHat.cross(tau);
1866         double bVal = lVec.dot(nPlane) * tNInv;
1867         if (fabs(aVal * bVal) < 0.3) {
1868           double cVal =
1869               1. - bHatN * cosPB * tNInv;  // - sv.bf.cross(lVec).dot(nPlane)*b0Inv*tNInv; //1- bHat_n*bHat_tau/tau_n;
1870           double aacVal = cVal * aVal * aVal;
1871           if (fabs(aacVal) < 1) {
1872             double tanDCorr = bVal / 2. + (bVal * bVal / 2. + cVal / 6) * aVal;
1873             tanDCorr *= aVal;
1874             //+ (-bVal/24. + 0.625*bVal*bVal*bVal + 5./12.*bVal*cVal)*aVal*aVal*aVal
1875             if (debug_)
1876               LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << tanDist << " vs "
1877                                 << tanDist * (1. + tanDCorr) << " corr " << tanDist * tanDCorr << std::endl;
1878             tanDist *= (1. + tanDCorr);
1879           } else {
1880             if (debug_)
1881               LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "AACVal "
1882                                 << fabs(aacVal) << " = " << aVal << "**2 * " << cVal << " too large:: will not converge"
1883                                 << std::endl;
1884           }
1885         } else {
1886           if (debug_)
1887             LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "ABVal "
1888                               << fabs(aVal * bVal) << " = " << aVal << " * " << bVal << " too large:: will not converge"
1889                               << std::endl;
1890         }
1891       }
1892       result = SteppingHelixStateInfo::OK;
1893     } break;
1894     case CONE_DT: {
1895       Point cVertex(pars[0], pars[1], pars[2]);
1896       if (cVertex.perp2() < 1e-10) {
1897         //assumes the cone axis/vertex is along z
1898         Vector relV3 = sv.r3 - cVertex;
1899         double relV3mag = relV3.mag();
1900         //  double relV3Theta = relV3.theta();
1901         double theta(pars[3]);
1902         //  double dTheta = theta-relV3Theta;
1903         double sinTheta = sin(theta);
1904         double cosTheta = cos(theta);
1905         double cosV3Theta = relV3.z() / relV3mag;
1906         if (cosV3Theta > 1)
1907           cosV3Theta = 1;
1908         if (cosV3Theta < -1)
1909           cosV3Theta = -1;
1910         double sinV3Theta = sqrt(1. - cosV3Theta * cosV3Theta);
1911 
1912         double sinDTheta = sinTheta * cosV3Theta - cosTheta * sinV3Theta;  //sin(dTheta);
1913         double cosDTheta = cosTheta * cosV3Theta + sinTheta * sinV3Theta;  //cos(dTheta);
1914         bool isInside = sinTheta > sinV3Theta && cosTheta * cosV3Theta > 0;
1915         dist = isInside || cosDTheta > 0 ? relV3mag * sinDTheta : relV3mag;
1916         if (fabs(dist) > fastSkipDist) {
1917           result = SteppingHelixStateInfo::INACC;
1918           break;
1919         }
1920 
1921         double relV3phi = relV3.phi();
1922         double normPhi = isInside ? Geom::pi() + relV3phi : relV3phi;
1923         double normTheta = theta > Geom::pi() / 2. ? (isInside ? 1.5 * Geom::pi() - theta : theta - Geom::pi() / 2.)
1924                                                    : (isInside ? Geom::pi() / 2. - theta : theta + Geom::pi() / 2);
1925         //this is a normVector from the cone to the point
1926         Vector norm;
1927         norm.setRThetaPhi(fabs(dist), normTheta, normPhi);
1928         double curP = sv.p3.mag();
1929         double cosp3theta = sv.p3.z() / curP;
1930         if (cosp3theta > 1)
1931           cosp3theta = 1;
1932         if (cosp3theta < -1)
1933           cosp3theta = -1;
1934         double sineConeP = sinTheta * cosp3theta - cosTheta * sqrt(1. - cosp3theta * cosp3theta);
1935 
1936         double sinSolid = norm.dot(sv.p3) / (fabs(dist) * curP);
1937         tanDist = fabs(sinSolid) > fabs(sineConeP) ? dist / fabs(sinSolid) : dist / fabs(sineConeP);
1938 
1939         refDirection = sinSolid > 0 ? oppositeToMomentum : alongMomentum;
1940         if (debug_) {
1941           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Cone pars: theta " << theta
1942                             << " normTheta " << norm.theta() << " p3Theta " << sv.p3.theta() << " sinD: " << sineConeP
1943                             << " sinSolid " << sinSolid;
1944         }
1945         if (debug_) {
1946           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
1947                             << "refToDest:toCone the point is " << (isInside ? "in" : "out") << "side the cone"
1948                             << std::endl;
1949         }
1950       }
1951       result = SteppingHelixStateInfo::OK;
1952     } break;
1953       //   case CYLINDER_DT:
1954       //     break;
1955     case PATHL_DT: {
1956       double curS = fabs(sv.path());
1957       dist = pars[PATHL_P] - curS;
1958       if (fabs(dist) > fastSkipDist) {
1959         result = SteppingHelixStateInfo::INACC;
1960         break;
1961       }
1962       refDirection = pars[PATHL_P] > 0 ? alongMomentum : oppositeToMomentum;
1963       tanDist = dist;
1964       result = SteppingHelixStateInfo::OK;
1965     } break;
1966     case POINT_PCA_DT: {
1967       Point pDest(pars[0], pars[1], pars[2]);
1968       double curP = sv.p3.mag();
1969       dist = (sv.r3 - pDest).mag() + 1e-24;  //add a small number to avoid 1/0
1970       tanDist = (sv.r3.dot(sv.p3) - pDest.dot(sv.p3)) / curP;
1971       //account for bending in magnetic field (quite approximate)
1972       double b0 = sv.bf.mag();
1973       if (b0 > 1.5e-6) {
1974         double p0 = curP;
1975         double kVal = 2.99792458e-3 * sv.q / p0 * b0;
1976         double aVal = fabs(dist * kVal);
1977         tanDist *= 1. / (1. + aVal);
1978         if (debug_)
1979           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "corrected by aVal " << aVal
1980                             << " to " << tanDist;
1981       }
1982       refDirection = tanDist < 0 ? alongMomentum : oppositeToMomentum;
1983       result = SteppingHelixStateInfo::OK;
1984     } break;
1985     case LINE_PCA_DT: {
1986       Point rLine(pars[0], pars[1], pars[2]);
1987       Vector dLine(pars[3], pars[4], pars[5]);
1988       dLine = (dLine - rLine);
1989       dLine *= 1. / dLine.mag();
1990       if (debug_)
1991         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "dLine " << dLine;
1992 
1993       Vector dR = sv.r3 - rLine;
1994       double curP = sv.p3.mag();
1995       Vector dRPerp = dR - dLine * (dR.dot(dLine));
1996       if (debug_)
1997         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "dRperp " << dRPerp;
1998 
1999       dist = dRPerp.mag() + 1e-24;  //add a small number to avoid 1/0
2000       tanDist = dRPerp.dot(sv.p3) / curP;
2001       //angle wrt line
2002       double cosAlpha2 = dLine.dot(sv.p3) / curP;
2003       cosAlpha2 *= cosAlpha2;
2004       tanDist *= 1. / sqrt(fabs(1. - cosAlpha2) + 1e-96);
2005       //correct for dPhi in magnetic field: this isn't made quite right here
2006       //(the angle between the line and the trajectory plane is neglected .. conservative)
2007       double b0 = sv.bf.mag();
2008       if (b0 > 1.5e-6) {
2009         double p0 = curP;
2010         double kVal = 2.99792458e-3 * sv.q / p0 * b0;
2011         double aVal = fabs(dist * kVal);
2012         tanDist *= 1. / (1. + aVal);
2013         if (debug_)
2014           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "corrected by aVal " << aVal
2015                             << " to " << tanDist;
2016       }
2017       refDirection = tanDist < 0 ? alongMomentum : oppositeToMomentum;
2018       result = SteppingHelixStateInfo::OK;
2019     } break;
2020     default: {
2021       //some large number
2022       dist = 1e12;
2023       tanDist = 1e12;
2024       refDirection = anyDirection;
2025       result = SteppingHelixStateInfo::NOT_IMPLEMENTED;
2026     } break;
2027   }
2028 
2029   double tanDistConstrained = tanDist;
2030   if (fabs(tanDist) > 4. * fabs(dist))
2031     tanDistConstrained *= tanDist == 0 ? 0 : fabs(dist / tanDist * 4.);
2032 
2033   if (debug_) {
2034     LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "refToDest input: dest" << dest
2035                       << " pars[]: ";
2036     for (int i = 0; i < 6; i++) {
2037       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << ", " << i << " " << pars[i];
2038     }
2039     LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << std::endl;
2040     LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "refToDest output: "
2041                       << "\t dist" << dist << "\t tanDist " << tanDist << "\t tanDistConstr " << tanDistConstrained
2042                       << "\t refDirection " << refDirection << std::endl;
2043   }
2044   tanDist = tanDistConstrained;
2045 
2046   return result;
2047 }
2048 
2049 SteppingHelixPropagator::Result SteppingHelixPropagator::refToMagVolume(const SteppingHelixPropagator::StateInfo& sv,
2050                                                                         PropagationDirection dir,
2051                                                                         double& dist,
2052                                                                         double& tanDist,
2053                                                                         double fastSkipDist,
2054                                                                         bool expectNewMagVolume,
2055                                                                         double maxStep) const {
2056   static const std::string metname = "SteppingHelixPropagator";
2057   Result result = SteppingHelixStateInfo::NOT_IMPLEMENTED;
2058   const MagVolume* cVol = sv.magVol;
2059 
2060   if (cVol == nullptr)
2061     return result;
2062   const std::vector<VolumeSide>& cVolFaces(cVol->faces());
2063 
2064   double distToFace[6] = {0, 0, 0, 0, 0, 0};
2065   double tanDistToFace[6] = {0, 0, 0, 0, 0, 0};
2066   PropagationDirection refDirectionToFace[6] = {
2067       anyDirection, anyDirection, anyDirection, anyDirection, anyDirection, anyDirection};
2068   Result resultToFace[6] = {result, result, result, result, result, result};
2069   int iFDest = -1;
2070   int iDistMin = -1;
2071 
2072   unsigned int iFDestSorted[6] = {0, 0, 0, 0, 0, 0};
2073   int nDestSorted = 0;
2074   unsigned int nearParallels = 0;
2075 
2076   double curP = sv.p3.mag();
2077 
2078   if (debug_) {
2079     LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Trying volume "
2080                       << " with " << cVolFaces.size() << " faces" << std::endl;
2081   }
2082 
2083   unsigned int nFaces = cVolFaces.size();
2084   for (unsigned int iFace = 0; iFace < nFaces; ++iFace) {
2085     if (iFace > 5) {
2086       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Too many faces" << std::endl;
2087     }
2088     if (debug_) {
2089       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Start with face " << iFace
2090                         << std::endl;
2091     }
2092     //     const Plane* cPlane = dynamic_cast<const Plane*>(&cVolFaces[iFace].surface());
2093     //     const Cylinder* cCyl = dynamic_cast<const Cylinder*>(&cVolFaces[iFace].surface());
2094     //     const Cone* cCone = dynamic_cast<const Cone*>(&cVolFaces[iFace].surface());
2095     const Surface* cPlane = nullptr;  //only need to know the loc->glob transform
2096     const Cylinder* cCyl = nullptr;
2097     const Cone* cCone = nullptr;
2098     auto& iSurface = cVolFaces[iFace].surface();
2099     if (typeid(iSurface) == typeid(const Plane&)) {
2100       cPlane = &cVolFaces[iFace].surface();
2101     } else if (typeid(iSurface) == typeid(const Cylinder&)) {
2102       cCyl = dynamic_cast<const Cylinder*>(&cVolFaces[iFace].surface());
2103     } else if (typeid(iSurface) == typeid(const Cone&)) {
2104       cCone = dynamic_cast<const Cone*>(&cVolFaces[iFace].surface());
2105     } else {
2106       edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2107                                << "Could not cast a volume side surface to a known type" << std::endl;
2108     }
2109 
2110     if (debug_) {
2111       if (cPlane != nullptr)
2112         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "The face is a plane at "
2113                           << cPlane << std::endl;
2114       if (cCyl != nullptr)
2115         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "The face is a cylinder at "
2116                           << cCyl << std::endl;
2117     }
2118 
2119     double pars[6] = {0, 0, 0, 0, 0, 0};
2120     DestType dType = UNDEFINED_DT;
2121     if (cPlane != nullptr) {
2122       Point rPlane(cPlane->position().x(), cPlane->position().y(), cPlane->position().z());
2123       // = cPlane->toGlobal(LocalVector(0,0,1.)); nPlane = nPlane.unit();
2124       Vector nPlane(cPlane->rotation().zx(), cPlane->rotation().zy(), cPlane->rotation().zz());
2125       nPlane /= nPlane.mag();
2126 
2127       pars[0] = rPlane.x();
2128       pars[1] = rPlane.y();
2129       pars[2] = rPlane.z();
2130       pars[3] = nPlane.x();
2131       pars[4] = nPlane.y();
2132       pars[5] = nPlane.z();
2133       dType = PLANE_DT;
2134     } else if (cCyl != nullptr) {
2135       if (debug_) {
2136         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Cylinder at "
2137                           << cCyl->position() << " rorated by " << cCyl->rotation() << std::endl;
2138       }
2139       pars[RADIUS_P] = cCyl->radius();
2140       dType = RADIUS_DT;
2141     } else if (cCone != nullptr) {
2142       if (debug_) {
2143         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Cone at "
2144                           << cCone->position() << " rorated by " << cCone->rotation() << " vertex at "
2145                           << cCone->vertex() << " angle of " << cCone->openingAngle() << std::endl;
2146       }
2147       pars[0] = cCone->vertex().x();
2148       pars[1] = cCone->vertex().y();
2149       pars[2] = cCone->vertex().z();
2150       pars[3] = cCone->openingAngle();
2151       dType = CONE_DT;
2152     } else {
2153       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Unknown surface" << std::endl;
2154       resultToFace[iFace] = SteppingHelixStateInfo::UNDEFINED;
2155       continue;
2156     }
2157     resultToFace[iFace] =
2158         refToDest(dType, sv, pars, distToFace[iFace], tanDistToFace[iFace], refDirectionToFace[iFace], fastSkipDist);
2159 
2160     if (resultToFace[iFace] != SteppingHelixStateInfo::OK) {
2161       if (resultToFace[iFace] == SteppingHelixStateInfo::INACC)
2162         result = SteppingHelixStateInfo::INACC;
2163     }
2164 
2165     //keep those in right direction for later use
2166     if (resultToFace[iFace] == SteppingHelixStateInfo::OK) {
2167       double invDTFPosiF = 1. / (1e-32 + fabs(tanDistToFace[iFace]));
2168       double dSlope = fabs(distToFace[iFace]) * invDTFPosiF;
2169       double maxStepL = maxStep > 100 ? 100 : maxStep;
2170       if (maxStepL < 10)
2171         maxStepL = 10.;
2172       bool isNearParallel =
2173           fabs(tanDistToFace[iFace]) + 100. * curP * dSlope < maxStepL  //
2174           //a better choice is to use distance to next check of mag volume instead of 100cm; the last is for ~1.5arcLength(4T)+tandistance< maxStep
2175           && dSlope < 0.15;  //
2176       if (refDirectionToFace[iFace] == dir || isNearParallel) {
2177         if (isNearParallel)
2178           nearParallels++;
2179         iFDestSorted[nDestSorted] = iFace;
2180         nDestSorted++;
2181       }
2182     }
2183 
2184     //pick a shortest distance here (right dir only for now)
2185     if (refDirectionToFace[iFace] == dir) {
2186       if (iDistMin == -1)
2187         iDistMin = iFace;
2188       else if (fabs(distToFace[iFace]) < fabs(distToFace[iDistMin]))
2189         iDistMin = iFace;
2190     }
2191     if (debug_)
2192       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << cVol << " " << iFace << " "
2193                         << tanDistToFace[iFace] << " " << distToFace[iFace] << " " << refDirectionToFace[iFace] << " "
2194                         << dir << std::endl;
2195   }
2196 
2197   for (int i = 0; i < nDestSorted; ++i) {
2198     int iMax = nDestSorted - i - 1;
2199     for (int j = 0; j < nDestSorted - i; ++j) {
2200       if (fabs(tanDistToFace[iFDestSorted[j]]) > fabs(tanDistToFace[iFDestSorted[iMax]])) {
2201         iMax = j;
2202       }
2203     }
2204     int iTmp = iFDestSorted[nDestSorted - i - 1];
2205     iFDestSorted[nDestSorted - i - 1] = iFDestSorted[iMax];
2206     iFDestSorted[iMax] = iTmp;
2207   }
2208 
2209   if (debug_) {
2210     for (int i = 0; i < nDestSorted; ++i) {
2211       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << cVol << " " << i << " "
2212                         << iFDestSorted[i] << " " << tanDistToFace[iFDestSorted[i]] << std::endl;
2213     }
2214   }
2215 
2216   //now go from the shortest to the largest distance hoping to get a point in the volume.
2217   //other than in case of a near-parallel travel this should stop after the first try
2218 
2219   for (int i = 0; i < nDestSorted; ++i) {
2220     iFDest = iFDestSorted[i];
2221 
2222     double sign = dir == alongMomentum ? 1. : -1.;
2223     Point gPointEst(sv.r3);
2224     Vector lDelta(sv.p3);
2225     lDelta *= sign / curP * fabs(distToFace[iFDest]);
2226     gPointEst += lDelta;
2227     if (debug_) {
2228       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Linear est point " << gPointEst
2229                         << " for iFace " << iFDest << std::endl;
2230     }
2231     GlobalPoint gPointEstNorZ(gPointEst.x(), gPointEst.y(), gPointEst.z());
2232     if (cVol->inside(gPointEstNorZ)) {
2233       if (debug_) {
2234         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2235                           << "The point is inside the volume" << std::endl;
2236       }
2237 
2238       result = SteppingHelixStateInfo::OK;
2239       dist = distToFace[iFDest];
2240       tanDist = tanDistToFace[iFDest];
2241       if (debug_) {
2242         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2243                           << "Got a point near closest boundary -- face " << iFDest << std::endl;
2244       }
2245       break;
2246     }
2247   }
2248 
2249   if (result != SteppingHelixStateInfo::OK && expectNewMagVolume) {
2250     double sign = dir == alongMomentum ? 1. : -1.;
2251 
2252     //check if it's a wrong volume situation
2253     if (nDestSorted - nearParallels > 0)
2254       result = SteppingHelixStateInfo::WRONG_VOLUME;
2255     else {
2256       //get here if all faces in the corr direction were skipped
2257       Point gPointEst(sv.r3);
2258       double lDist = iDistMin == -1 ? fastSkipDist : fabs(distToFace[iDistMin]);
2259       if (lDist > fastSkipDist)
2260         lDist = fastSkipDist;
2261       Vector lDelta(sv.p3);
2262       lDelta *= sign / curP * lDist;
2263       gPointEst += lDelta;
2264       if (debug_) {
2265         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2266                           << "Linear est point to shortest dist " << gPointEst << " for iFace " << iDistMin
2267                           << " at distance " << lDist * sign << std::endl;
2268       }
2269       GlobalPoint gPointEstNorZ(gPointEst.x(), gPointEst.y(), gPointEst.z());
2270       if (cVol->inside(gPointEstNorZ)) {
2271         if (debug_) {
2272           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2273                             << "The point is inside the volume" << std::endl;
2274         }
2275 
2276       } else {
2277         result = SteppingHelixStateInfo::WRONG_VOLUME;
2278       }
2279     }
2280 
2281     if (result == SteppingHelixStateInfo::WRONG_VOLUME) {
2282       dist = sign * 0.05;
2283       tanDist = dist * 1.01;
2284       if (debug_) {
2285         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2286                           << "Wrong volume located: return small dist, tandist" << std::endl;
2287       }
2288     }
2289   }
2290 
2291   if (result == SteppingHelixStateInfo::INACC) {
2292     if (debug_)
2293       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "All faces are too far"
2294                         << std::endl;
2295   } else if (result == SteppingHelixStateInfo::WRONG_VOLUME) {
2296     if (debug_)
2297       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Appear to be in a wrong volume"
2298                         << std::endl;
2299   } else if (result != SteppingHelixStateInfo::OK) {
2300     if (debug_)
2301       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Something else went wrong"
2302                         << std::endl;
2303   }
2304 
2305   return result;
2306 }
2307 
2308 SteppingHelixPropagator::Result SteppingHelixPropagator::refToMatVolume(const SteppingHelixPropagator::StateInfo& sv,
2309                                                                         PropagationDirection dir,
2310                                                                         double& dist,
2311                                                                         double& tanDist,
2312                                                                         double fastSkipDist) const {
2313   static const std::string metname = "SteppingHelixPropagator";
2314   Result result = SteppingHelixStateInfo::NOT_IMPLEMENTED;
2315 
2316   double parLim[6] = {sv.rzLims.rMin, sv.rzLims.rMax, sv.rzLims.zMin, sv.rzLims.zMax, sv.rzLims.th1, sv.rzLims.th2};
2317 
2318   double distToFace[4] = {0, 0, 0, 0};
2319   double tanDistToFace[4] = {0, 0, 0, 0};
2320   PropagationDirection refDirectionToFace[4] = {anyDirection, anyDirection, anyDirection, anyDirection};
2321   Result resultToFace[4] = {result, result, result, result};
2322   int iFDest = -1;
2323 
2324   double curP = sv.p3.mag();
2325 
2326   for (unsigned int iFace = 0; iFace < 4; iFace++) {
2327     if (debug_) {
2328       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Start with mat face " << iFace
2329                         << std::endl;
2330     }
2331 
2332     double pars[6] = {0, 0, 0, 0, 0, 0};
2333     DestType dType = UNDEFINED_DT;
2334     if (iFace > 1) {
2335       if (fabs(parLim[iFace + 2]) < 1e-6) {  //plane
2336         if (sv.r3.z() < 0) {
2337           pars[0] = 0;
2338           pars[1] = 0;
2339           pars[2] = -parLim[iFace];
2340           pars[3] = 0;
2341           pars[4] = 0;
2342           pars[5] = 1;
2343         } else {
2344           pars[0] = 0;
2345           pars[1] = 0;
2346           pars[2] = parLim[iFace];
2347           pars[3] = 0;
2348           pars[4] = 0;
2349           pars[5] = 1;
2350         }
2351         dType = PLANE_DT;
2352       } else {
2353         if (sv.r3.z() > 0) {
2354           pars[0] = 0;
2355           pars[1] = 0;
2356           pars[2] = parLim[iFace];
2357           pars[3] = Geom::pi() / 2. - parLim[iFace + 2];
2358         } else {
2359           pars[0] = 0;
2360           pars[1] = 0;
2361           pars[2] = -parLim[iFace];
2362           pars[3] = Geom::pi() / 2. + parLim[iFace + 2];
2363         }
2364         dType = CONE_DT;
2365       }
2366     } else {
2367       pars[RADIUS_P] = parLim[iFace];
2368       dType = RADIUS_DT;
2369     }
2370 
2371     resultToFace[iFace] =
2372         refToDest(dType, sv, pars, distToFace[iFace], tanDistToFace[iFace], refDirectionToFace[iFace], fastSkipDist);
2373 
2374     if (resultToFace[iFace] != SteppingHelixStateInfo::OK) {
2375       if (resultToFace[iFace] == SteppingHelixStateInfo::INACC)
2376         result = SteppingHelixStateInfo::INACC;
2377       continue;
2378     }
2379     if (refDirectionToFace[iFace] == dir || fabs(distToFace[iFace]) < 2e-2 * fabs(tanDistToFace[iFace])) {
2380       double sign = dir == alongMomentum ? 1. : -1.;
2381       Point gPointEst(sv.r3);
2382       Vector lDelta(sv.p3);
2383       lDelta *= sign * fabs(distToFace[iFace]) / curP;
2384       gPointEst += lDelta;
2385       if (debug_) {
2386         LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Linear est point "
2387                           << gPointEst << std::endl;
2388       }
2389       double lZ = fabs(gPointEst.z());
2390       double lR = gPointEst.perp();
2391       double tan4 = parLim[4] == 0 ? 0 : tan(parLim[4]);
2392       double tan5 = parLim[5] == 0 ? 0 : tan(parLim[5]);
2393       if ((lZ - parLim[2]) > lR * tan4 && (lZ - parLim[3]) < lR * tan5 && lR > parLim[0] && lR < parLim[1]) {
2394         if (debug_) {
2395           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2396                             << "The point is inside the volume" << std::endl;
2397         }
2398         //OK, guessed a point still inside the volume
2399         if (iFDest == -1) {
2400           iFDest = iFace;
2401         } else {
2402           if (fabs(tanDistToFace[iFDest]) > fabs(tanDistToFace[iFace])) {
2403             iFDest = iFace;
2404           }
2405         }
2406       } else {
2407         if (debug_) {
2408           LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2409                             << "The point is NOT inside the volume" << std::endl;
2410         }
2411       }
2412     }
2413   }
2414   if (iFDest != -1) {
2415     result = SteppingHelixStateInfo::OK;
2416     dist = distToFace[iFDest];
2417     tanDist = tanDistToFace[iFDest];
2418     if (debug_) {
2419       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2420                         << "Got a point near closest boundary -- face " << iFDest << std::endl;
2421     }
2422   } else {
2423     if (debug_) {
2424       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2425                         << "Failed to find a dest point inside the volume" << std::endl;
2426     }
2427   }
2428 
2429   return result;
2430 }
2431 
2432 bool SteppingHelixPropagator::isYokeVolume(const MagVolume* vol) const {
2433   static const std::string metname = "SteppingHelixPropagator";
2434   if (vol == nullptr)
2435     return false;
2436   /*
2437   const MFGrid* mGrid = reinterpret_cast<const MFGrid*>(vol->provider());
2438   std::vector<int> dims(mGrid->dimensions());
2439   
2440   LocalVector lVCen(mGrid->nodeValue(dims[0]/2, dims[1]/2, dims[2]/2));
2441   LocalVector lVZLeft(mGrid->nodeValue(dims[0]/2, dims[1]/2, dims[2]/5));
2442   LocalVector lVZRight(mGrid->nodeValue(dims[0]/2, dims[1]/2, (dims[2]*4)/5));
2443 
2444   double mag2VCen = lVCen.mag2();
2445   double mag2VZLeft = lVZLeft.mag2();
2446   double mag2VZRight = lVZRight.mag2();
2447 
2448   bool result = false;
2449   if (mag2VCen > 0.6 && mag2VZLeft > 0.6 && mag2VZRight > 0.6){
2450     if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Volume is magnetic, located at "<<vol->position()<<std::endl;    
2451     result = true;
2452   } else {
2453     if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Volume is not magnetic, located at "<<vol->position()<<std::endl;
2454   }
2455 
2456   */
2457   bool result = vol->isIron();
2458   if (result) {
2459     if (debug_)
2460       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2461                         << "Volume is magnetic, located at " << vol->position() << std::endl;
2462   } else {
2463     if (debug_)
2464       LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2465                         << "Volume is not magnetic, located at " << vol->position() << std::endl;
2466   }
2467 
2468   return result;
2469 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!