CMS 3D CMS Logo
SteppingHelixPropagator.cc
Go to the documentation of this file.
1 
11 //
12 // Original Author: Vyacheslav Krutelyov
13 // Created: Fri Mar 3 16:01:24 CST 2006
14 //
15 //
16 
17 #include "MagneticField/Engine/interface/MagneticField.h"
18 #include "MagneticField/VolumeBasedEngine/interface/VolumeBasedMagneticField.h"
19 #include "MagneticField/VolumeGeometry/interface/MagVolume.h"
20 #include "MagneticField/Interpolation/interface/MFGrid.h"
21 
22 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
23 #include "DataFormats/GeometrySurface/interface/Plane.h"
24 #include "DataFormats/GeometrySurface/interface/Cone.h"
25 
26 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
27 
28 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToCurvilinear.h"
29 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToCartesian.h"
30 
31 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
32 
33 #include "FWCore/MessageLogger/interface/MessageLogger.h"
34 
35 #include <sstream>
36 #include <typeinfo>
37 
38 void SteppingHelixPropagator::initStateArraySHPSpecific(StateArray& svBuf, bool flagsOnly) const {
39  for (int i = 0; i <= MAX_POINTS; i++) {
40  svBuf[i].isComplete = true;
41  svBuf[i].isValid_ = true;
42  svBuf[i].hasErrorPropagated_ = !noErrorPropagation_;
43  if (!flagsOnly) {
44  svBuf[i].p3 = Vector(0, 0, 0);
45  svBuf[i].r3 = Point(0, 0, 0);
46  svBuf[i].bf = Vector(0, 0, 0);
47  svBuf[i].bfGradLoc = Vector(0, 0, 0);
48  svBuf[i].covCurv = AlgebraicSymMatrix55();
49  svBuf[i].matDCovCurv = AlgebraicSymMatrix55();
50  }
51  }
52 }
53 
54 SteppingHelixPropagator::~SteppingHelixPropagator() {}
55 
56 SteppingHelixPropagator::SteppingHelixPropagator() : Propagator(anyDirection) { field_ = nullptr; }
57 
58 SteppingHelixPropagator::SteppingHelixPropagator(const MagneticField* field, PropagationDirection dir)
59  : Propagator(dir), unit55_(AlgebraicMatrixID()) {
60  field_ = field;
61  vbField_ = dynamic_cast<const VolumeBasedMagneticField*>(field_);
62  debug_ = false;
63  noMaterialMode_ = false;
64  noErrorPropagation_ = false;
65  applyRadX0Correction_ = true;
66  useMagVolumes_ = true;
67  useIsYokeFlag_ = true;
68  useMatVolumes_ = true;
69  useInTeslaFromMagField_ = false; //overrides behavior only if true
70  returnTangentPlane_ = true;
71  sendLogWarning_ = false;
72  useTuningForL2Speed_ = false;
73  defaultStep_ = 5.;
74 
75  ecShiftPos_ = 0;
76  ecShiftNeg_ = 0;
77 }
78 
79 std::pair<TrajectoryStateOnSurface, double> SteppingHelixPropagator::propagateWithPath(
80  const FreeTrajectoryState& ftsStart, const Plane& pDest) const {
81  StateArray svBuf;
82  initStateArraySHPSpecific(svBuf, true);
83  int nPoints = 0;
84  setIState(SteppingHelixStateInfo(ftsStart), svBuf, nPoints);
85 
86  StateInfo svCurrent;
87  propagate(svBuf[0], pDest, svCurrent);
88 
89  return TsosPP(svCurrent.getStateOnSurface(pDest), svCurrent.path());
90 }
91 
92 std::pair<TrajectoryStateOnSurface, double> SteppingHelixPropagator::propagateWithPath(
93  const FreeTrajectoryState& ftsStart, const Cylinder& cDest) const {
94  StateArray svBuf;
95  initStateArraySHPSpecific(svBuf, true);
96  int nPoints = 0;
97  setIState(SteppingHelixStateInfo(ftsStart), svBuf, nPoints);
98 
99  StateInfo svCurrent;
100  propagate(svBuf[0], cDest, svCurrent);
101 
102  return TsosPP(svCurrent.getStateOnSurface(cDest, returnTangentPlane_), svCurrent.path());
103 }
104 
105 std::pair<FreeTrajectoryState, double> SteppingHelixPropagator::propagateWithPath(const FreeTrajectoryState& ftsStart,
106  const GlobalPoint& pDest) const {
107  StateArray svBuf;
108  initStateArraySHPSpecific(svBuf, true);
109  int nPoints = 0;
110  setIState(SteppingHelixStateInfo(ftsStart), svBuf, nPoints);
111 
112  StateInfo svCurrent;
113  propagate(svBuf[0], pDest, svCurrent);
114 
115  FreeTrajectoryState ftsDest;
116  svCurrent.getFreeState(ftsDest);
117 
118  return FtsPP(ftsDest, svCurrent.path());
119 }
120 
121 std::pair<FreeTrajectoryState, double> SteppingHelixPropagator::propagateWithPath(const FreeTrajectoryState& ftsStart,
122  const GlobalPoint& pDest1,
123  const GlobalPoint& pDest2) const {
124  if ((pDest1 - pDest2).mag() < 1e-10) {
125  if (sendLogWarning_) {
126  edm::LogWarning("SteppingHelixPropagator")
127  << "Can't propagate: the points should be at a bigger distance" << std::endl;
128  }
129  return FtsPP();
130  }
131  StateArray svBuf;
132  initStateArraySHPSpecific(svBuf, true);
133  int nPoints = 0;
134  setIState(SteppingHelixStateInfo(ftsStart), svBuf, nPoints);
135 
136  StateInfo svCurrent;
137  propagate(svBuf[0], pDest1, pDest2, svCurrent);
138 
139  FreeTrajectoryState ftsDest;
140  svCurrent.getFreeState(ftsDest);
141 
142  return FtsPP(ftsDest, svCurrent.path());
143 }
144 
145 std::pair<FreeTrajectoryState, double> SteppingHelixPropagator::propagateWithPath(
146  const FreeTrajectoryState& ftsStart, const reco::BeamSpot& beamSpot) const {
147  GlobalPoint pDest1(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
148  GlobalPoint pDest2(pDest1.x() + beamSpot.dxdz() * 1e3, pDest1.y() + beamSpot.dydz() * 1e3, pDest1.z() + 1e3);
149  return propagateWithPath(ftsStart, pDest1, pDest2);
150 }
151 
152 void SteppingHelixPropagator::propagate(const SteppingHelixStateInfo& sStart,
153  const Surface& sDest,
154  SteppingHelixStateInfo& out) const {
155  if (!sStart.isValid()) {
156  if (sendLogWarning_) {
157  edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
158  }
159  out = invalidState_;
160  return;
161  }
162 
163  const Plane* pDest = dynamic_cast<const Plane*>(&sDest);
164  if (pDest != nullptr) {
165  propagate(sStart, *pDest, out);
166  return;
167  }
168 
169  const Cylinder* cDest = dynamic_cast<const Cylinder*>(&sDest);
170  if (cDest != nullptr) {
171  propagate(sStart, *cDest, out);
172  return;
173  }
174 
175  throw PropagationException("The surface is neither Cylinder nor Plane");
176 }
177 
178 void SteppingHelixPropagator::propagate(const SteppingHelixStateInfo& sStart,
179  const Plane& pDest,
180  SteppingHelixStateInfo& out) const {
181  if (!sStart.isValid()) {
182  if (sendLogWarning_) {
183  edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
184  }
185  out = invalidState_;
186  return;
187  }
188  StateArray svBuf;
189  initStateArraySHPSpecific(svBuf, true);
190  int nPoints = 0;
191  setIState(sStart, svBuf, nPoints);
192 
193  Point rPlane(pDest.position().x(), pDest.position().y(), pDest.position().z());
194  Vector nPlane(pDest.rotation().zx(), pDest.rotation().zy(), pDest.rotation().zz());
195  nPlane /= nPlane.mag();
196 
197  double pars[6] = {rPlane.x(), rPlane.y(), rPlane.z(), nPlane.x(), nPlane.y(), nPlane.z()};
198 
199  propagate(svBuf, nPoints, PLANE_DT, pars);
200 
201  //(re)set it before leaving: dir =1 (-1) if path increased (decreased) and 0 if it didn't change
202  //need to implement this somewhere else as a separate function
203  double lDir = 0;
204  if (sStart.path() < svBuf[cIndex_(nPoints - 1)].path())
205  lDir = 1.;
206  if (sStart.path() > svBuf[cIndex_(nPoints - 1)].path())
207  lDir = -1.;
208  svBuf[cIndex_(nPoints - 1)].dir = lDir;
209 
210  out = svBuf[cIndex_(nPoints - 1)];
211  return;
212 }
213 
214 void SteppingHelixPropagator::propagate(const SteppingHelixStateInfo& sStart,
215  const Cylinder& cDest,
216  SteppingHelixStateInfo& out) const {
217  if (!sStart.isValid()) {
218  if (sendLogWarning_) {
219  edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
220  }
221  out = invalidState_;
222  return;
223  }
224  StateArray svBuf;
225  initStateArraySHPSpecific(svBuf, true);
226  int nPoints = 0;
227  setIState(sStart, svBuf, nPoints);
228 
229  double pars[6] = {0, 0, 0, 0, 0, 0};
230  pars[RADIUS_P] = cDest.radius();
231 
232  propagate(svBuf, nPoints, RADIUS_DT, pars);
233 
234  //(re)set it before leaving: dir =1 (-1) if path increased (decreased) and 0 if it didn't change
235  //need to implement this somewhere else as a separate function
236  double lDir = 0;
237  if (sStart.path() < svBuf[cIndex_(nPoints - 1)].path())
238  lDir = 1.;
239  if (sStart.path() > svBuf[cIndex_(nPoints - 1)].path())
240  lDir = -1.;
241  svBuf[cIndex_(nPoints - 1)].dir = lDir;
242  out = svBuf[cIndex_(nPoints - 1)];
243  return;
244 }
245 
246 void SteppingHelixPropagator::propagate(const SteppingHelixStateInfo& sStart,
247  const GlobalPoint& pDest,
248  SteppingHelixStateInfo& out) const {
249  if (!sStart.isValid()) {
250  if (sendLogWarning_) {
251  edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
252  }
253  out = invalidState_;
254  return;
255  }
256  StateArray svBuf;
257  initStateArraySHPSpecific(svBuf, true);
258  int nPoints = 0;
259  setIState(sStart, svBuf, nPoints);
260 
261  double pars[6] = {pDest.x(), pDest.y(), pDest.z(), 0, 0, 0};
262 
263  propagate(svBuf, nPoints, POINT_PCA_DT, pars);
264 
265  out = svBuf[cIndex_(nPoints - 1)];
266  return;
267 }
268 
269 void SteppingHelixPropagator::propagate(const SteppingHelixStateInfo& sStart,
270  const GlobalPoint& pDest1,
271  const GlobalPoint& pDest2,
272  SteppingHelixStateInfo& out) const {
273  if ((pDest1 - pDest2).mag() < 1e-10 || !sStart.isValid()) {
274  if (sendLogWarning_) {
275  if ((pDest1 - pDest2).mag() < 1e-10)
276  edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: points are too close" << std::endl;
277  if (!sStart.isValid())
278  edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
279  }
280  out = invalidState_;
281  return;
282  }
283  StateArray svBuf;
284  initStateArraySHPSpecific(svBuf, true);
285  int nPoints = 0;
286  setIState(sStart, svBuf, nPoints);
287 
288  double pars[6] = {pDest1.x(), pDest1.y(), pDest1.z(), pDest2.x(), pDest2.y(), pDest2.z()};
289 
290  propagate(svBuf, nPoints, LINE_PCA_DT, pars);
291 
292  out = svBuf[cIndex_(nPoints - 1)];
293  return;
294 }
295 
296 void SteppingHelixPropagator::setIState(const SteppingHelixStateInfo& sStart, StateArray& svBuf, int& nPoints) const {
297  nPoints = 0;
298  svBuf[cIndex_(nPoints)] = sStart; //do this anyways to have a fresh start
299  if (sStart.isComplete) {
300  svBuf[cIndex_(nPoints)] = sStart;
301  nPoints++;
302  } else {
303  loadState(svBuf[cIndex_(nPoints)], sStart.p3, sStart.r3, sStart.q, propagationDirection(), sStart.covCurv);
304  nPoints++;
305  }
306  svBuf[cIndex_(0)].hasErrorPropagated_ = sStart.hasErrorPropagated_ & !noErrorPropagation_;
307 }
308 
309 SteppingHelixPropagator::Result SteppingHelixPropagator::propagate(StateArray& svBuf,
310  int& nPoints,
311  SteppingHelixPropagator::DestType type,
312  const double pars[6],
313  double epsilon) const {
314  static const std::string metname = "SteppingHelixPropagator";
315  StateInfo* svCurrent = &svBuf[cIndex_(nPoints - 1)];
316 
317  //check if it's going to work at all
318  double tanDist = 0;
319  double dist = 0;
320  PropagationDirection refDirection = anyDirection;
321  Result result = refToDest(type, (*svCurrent), pars, dist, tanDist, refDirection);
322 
323  if (result != SteppingHelixStateInfo::OK || fabs(dist) > 1e6) {
324  svCurrent->status_ = result;
325  if (fabs(dist) > 1e6)
326  svCurrent->status_ = SteppingHelixStateInfo::INACC;
327  svCurrent->isValid_ = false;
328  svCurrent->field = field_;
329  if (sendLogWarning_) {
330  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
331  << " Failed after first refToDest check with status "
332  << SteppingHelixStateInfo::ResultName[result] << std::endl;
333  }
334  return result;
335  }
336 
337  result = SteppingHelixStateInfo::UNDEFINED;
338  bool makeNextStep = true;
339  double dStep = defaultStep_;
340  PropagationDirection dir, oldDir;
341  dir = propagationDirection();
342  oldDir = dir;
343  int nOsc = 0;
344 
345  double distMag = 1e12;
346  double tanDistMag = 1e12;
347 
348  double distMat = 1e12;
349  double tanDistMat = 1e12;
350 
351  double tanDistNextCheck = -0.1; //just need a negative start val
352  double tanDistMagNextCheck = -0.1;
353  double tanDistMatNextCheck = -0.1;
354  double oldDStep = 0;
355  PropagationDirection oldRefDirection = propagationDirection();
356 
357  Result resultToMat = SteppingHelixStateInfo::UNDEFINED;
358  Result resultToMag = SteppingHelixStateInfo::UNDEFINED;
359 
360  bool isFirstStep = true;
361  bool expectNewMagVolume = false;
362 
363  int loopCount = 0;
364  while (makeNextStep) {
365  dStep = defaultStep_;
366  svCurrent = &svBuf[cIndex_(nPoints - 1)];
367  double curZ = svCurrent->r3.z();
368  double curR = svCurrent->r3.perp();
369  if (fabs(curZ) < 440 && curR < 260)
370  dStep = defaultStep_ * 2;
371 
372  //more such ifs might be scattered around
373  //even though dStep is large, it will still make stops at each volume boundary
374  if (useTuningForL2Speed_) {
375  dStep = 100.;
376  if (!useIsYokeFlag_ && fabs(curZ) < 667 && curR > 380 && curR < 850) {
377  dStep = 5 * (1 + 0.2 * svCurrent->p3.mag());
378  }
379  }
380 
381  // refDirection = propagationDirection();
382 
383  tanDistNextCheck -= oldDStep;
384  tanDistMagNextCheck -= oldDStep;
385  tanDistMatNextCheck -= oldDStep;
386 
387  if (tanDistNextCheck < 0) {
388  //use pre-computed values if it's the first step
389  if (!isFirstStep)
390  refToDest(type, (*svCurrent), pars, dist, tanDist, refDirection);
391  // constrain allowed path for a tangential approach
392  if (fabs(tanDist) > 4. * (fabs(dist)))
393  tanDist *= tanDist == 0 ? 0 : fabs(dist / tanDist * 4.);
394 
395  tanDistNextCheck = fabs(tanDist) * 0.5 - 0.5; //need a better guess (to-do)
396  //reasonable limit
397  if (tanDistNextCheck > defaultStep_ * 20.)
398  tanDistNextCheck = defaultStep_ * 20.;
399  oldRefDirection = refDirection;
400  } else {
401  tanDist = tanDist > 0. ? tanDist - oldDStep : tanDist + oldDStep;
402  refDirection = oldRefDirection;
403  if (debug_)
404  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
405  << "Skipped refToDest: guess tanDist = " << tanDist << " next check at " << tanDistNextCheck
406  << std::endl;
407  }
409  double fastSkipDist = fabs(dist) > fabs(tanDist) ? tanDist : dist;
410 
411  if (propagationDirection() == anyDirection) {
412  dir = refDirection;
413  } else {
414  dir = propagationDirection();
415  if (fabs(tanDist) < 0.1 && refDirection != dir) {
416  //how did it get here? nOsc++;
417  dir = refDirection;
418  if (debug_)
419  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
420  << "NOTE: overstepped last time: switch direction (can do it if within 1 mm)" << std::endl;
421  }
422  }
423 
424  if (useMagVolumes_ && !(fabs(dist) < fabs(epsilon))) { //need to know the general direction
425  if (tanDistMagNextCheck < 0) {
426  resultToMag = refToMagVolume(
427  (*svCurrent), dir, distMag, tanDistMag, fabs(fastSkipDist), expectNewMagVolume, fabs(tanDist));
428  // constrain allowed path for a tangential approach
429  if (fabs(tanDistMag) > 4. * (fabs(distMag)))
430  tanDistMag *= tanDistMag == 0 ? 0 : fabs(distMag / tanDistMag * 4.);
431 
432  tanDistMagNextCheck = fabs(tanDistMag) * 0.5 - 0.5; //need a better guess (to-do)
433  //reasonable limit; "turn off" checking if bounds are further than the destination
434  if (tanDistMagNextCheck > defaultStep_ * 20. || fabs(dist) < fabs(distMag) ||
435  resultToMag == SteppingHelixStateInfo::INACC)
436  tanDistMagNextCheck = defaultStep_ * 20 > fabs(fastSkipDist) ? fabs(fastSkipDist) : defaultStep_ * 20;
437  if (resultToMag != SteppingHelixStateInfo::INACC && resultToMag != SteppingHelixStateInfo::OK)
438  tanDistMagNextCheck = -1;
439  } else {
440  // resultToMag = SteppingHelixStateInfo::OK;
441  tanDistMag = tanDistMag > 0. ? tanDistMag - oldDStep : tanDistMag + oldDStep;
442  if (debug_)
443  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
444  << "Skipped refToMag: guess tanDistMag = " << tanDistMag << " next check at "
445  << tanDistMagNextCheck;
446  }
447  }
448 
449  if (useMatVolumes_ && !(fabs(dist) < fabs(epsilon))) { //need to know the general direction
450  if (tanDistMatNextCheck < 0) {
451  resultToMat = refToMatVolume((*svCurrent), dir, distMat, tanDistMat, fabs(fastSkipDist));
452  // constrain allowed path for a tangential approach
453  if (fabs(tanDistMat) > 4. * (fabs(distMat)))
454  tanDistMat *= tanDistMat == 0 ? 0 : fabs(distMat / tanDistMat * 4.);
455 
456  tanDistMatNextCheck = fabs(tanDistMat) * 0.5 - 0.5; //need a better guess (to-do)
457  //reasonable limit; "turn off" checking if bounds are further than the destination
458  if (tanDistMatNextCheck > defaultStep_ * 20. || fabs(dist) < fabs(distMat) ||
459  resultToMat == SteppingHelixStateInfo::INACC)
460  tanDistMatNextCheck = defaultStep_ * 20 > fabs(fastSkipDist) ? fabs(fastSkipDist) : defaultStep_ * 20;
461  if (resultToMat != SteppingHelixStateInfo::INACC && resultToMat != SteppingHelixStateInfo::OK)
462  tanDistMatNextCheck = -1;
463  } else {
464  // resultToMat = SteppingHelixStateInfo::OK;
465  tanDistMat = tanDistMat > 0. ? tanDistMat - oldDStep : tanDistMat + oldDStep;
466  if (debug_)
467  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
468  << "Skipped refToMat: guess tanDistMat = " << tanDistMat << " next check at "
469  << tanDistMatNextCheck;
470  }
471  }
472 
473  double rDotP = svCurrent->r3.dot(svCurrent->p3);
474  if ((fabs(curZ) > 1.5e3 || curR > 800.) &&
475  ((dir == alongMomentum && rDotP > 0) || (dir == oppositeToMomentum && rDotP < 0))) {
476  dStep = fabs(tanDist) - 1e-12;
477  }
478  double tanDistMin = fabs(tanDist);
479  if (tanDistMin > fabs(tanDistMag) + 0.05 &&
480  (resultToMag == SteppingHelixStateInfo::OK || resultToMag == SteppingHelixStateInfo::WRONG_VOLUME)) {
481  tanDistMin = fabs(tanDistMag) + 0.05; //try to step into the next volume
482  expectNewMagVolume = true;
483  } else
484  expectNewMagVolume = false;
485 
486  if (tanDistMin > fabs(tanDistMat) + 0.05 && resultToMat == SteppingHelixStateInfo::OK) {
487  tanDistMin = fabs(tanDistMat) + 0.05; //try to step into the next volume
488  if (expectNewMagVolume)
489  expectNewMagVolume = false;
490  }
491 
492  if (tanDistMin * fabs(svCurrent->dEdx) > 0.5 * svCurrent->p3.mag()) {
493  tanDistMin = 0.5 * svCurrent->p3.mag() / fabs(svCurrent->dEdx);
494  if (expectNewMagVolume)
495  expectNewMagVolume = false;
496  }
497 
498  double tanDistMinLazy = fabs(tanDistMin);
499  if ((type == POINT_PCA_DT || type == LINE_PCA_DT) && fabs(tanDist) < 2. * fabs(tanDistMin)) {
500  //being lazy here; the best is to take into account the curvature
501  tanDistMinLazy = fabs(tanDistMin) * 0.5;
502  }
503 
504  if (fabs(tanDistMinLazy) < dStep) {
505  dStep = fabs(tanDistMinLazy);
506  }
507 
508  //keep this path length for the next step
509  oldDStep = dStep;
510 
511  if (dStep > 1e-10 && !(fabs(dist) < fabs(epsilon))) {
512  StateInfo* svNext = &svBuf[cIndex_(nPoints)];
513  makeAtomStep((*svCurrent), (*svNext), dStep, dir, HEL_AS_F);
514  // if (useMatVolumes_ && expectNewMagVolume
515  // && svCurrent->magVol == svNext->magVol){
516  // double tmpDist=0;
517  // double tmpDistMag = 0;
518  // if (refToMagVolume((*svNext), dir, tmpDist, tmpDistMag, fabs(dist)) != SteppingHelixStateInfo::OK){
519  // //the point appears to be outside, but findVolume claims the opposite
520  // dStep += 0.05*fabs(tanDistMag/distMag); oldDStep = dStep; //do it again with a bigger step
521  // if (debug_) LogTrace(metname)
522  // <<"Failed to get into new mag volume: will try with new bigger step "<<dStep<<std::endl;
523  // makeAtomStep((*svCurrent), (*svNext), dStep, dir, HEL_AS_F);
524  // }
525  // }
526  nPoints++;
527  svCurrent = &svBuf[cIndex_(nPoints - 1)];
528  if (oldDir != dir) {
529  nOsc++;
530  tanDistNextCheck = -1; //check dist after osc
531  tanDistMagNextCheck = -1;
532  tanDistMatNextCheck = -1;
533  }
534  oldDir = dir;
535  }
536 
537  if (nOsc > 1 && fabs(dStep) > epsilon) {
538  if (debug_)
539  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Ooops" << std::endl;
540  }
541 
542  if (fabs(dist) < fabs(epsilon))
543  result = SteppingHelixStateInfo::OK;
544 
545  if ((type == POINT_PCA_DT || type == LINE_PCA_DT) && fabs(dStep) < fabs(epsilon)) {
546  //now check if it's not a branch point (peek ahead at 1 cm)
547  double nextDist = 0;
548  double nextTanDist = 0;
549  PropagationDirection nextRefDirection = anyDirection;
550  StateInfo* svNext = &svBuf[cIndex_(nPoints)];
551  makeAtomStep((*svCurrent), (*svNext), 1., dir, HEL_AS_F);
552  nPoints++;
553  svCurrent = &svBuf[cIndex_(nPoints - 1)];
554  refToDest(type, (*svCurrent), pars, nextDist, nextTanDist, nextRefDirection);
555  if (fabs(nextDist) > fabs(dist)) {
556  nPoints--;
557  svCurrent = &svBuf[cIndex_(nPoints - 1)];
558  result = SteppingHelixStateInfo::OK;
559  if (debug_) {
560  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
561  << "Found real local minimum in PCA" << std::endl;
562  }
563  } else {
564  //keep this trial point and continue
565  dStep = defaultStep_;
566  if (debug_) {
567  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Found branch point in PCA"
568  << std::endl;
569  }
570  }
571  }
572 
573  if (nPoints > MAX_STEPS * 1. / defaultStep_ || loopCount > MAX_STEPS * 100 || nOsc > 6)
574  result = SteppingHelixStateInfo::FAULT;
575 
576  if (svCurrent->p3.mag() < 0.1)
577  result = SteppingHelixStateInfo::RANGEOUT;
578 
579  curZ = svCurrent->r3.z();
580  curR = svCurrent->r3.perp();
581  if (curR > 20000 || fabs(curZ) > 20000)
582  result = SteppingHelixStateInfo::INACC;
583 
584  makeNextStep = result == SteppingHelixStateInfo::UNDEFINED;
585  svCurrent->status_ = result;
586  svCurrent->isValid_ = (result == SteppingHelixStateInfo::OK || makeNextStep);
587  svCurrent->field = field_;
588 
589  isFirstStep = false;
590  loopCount++;
591  }
592 
593  if (sendLogWarning_ && result != SteppingHelixStateInfo::OK) {
594  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
595  << " Propagation failed with status " << SteppingHelixStateInfo::ResultName[result]
596  << std::endl;
597  if (result == SteppingHelixStateInfo::RANGEOUT)
598  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
599  << " Momentum at last point is too low (<0.1) p_last = " << svCurrent->p3.mag()
600  << std::endl;
601  if (result == SteppingHelixStateInfo::INACC)
602  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
603  << " Went too far: the last point is at " << svCurrent->r3 << std::endl;
604  if (result == SteppingHelixStateInfo::FAULT && nOsc > 6)
605  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
606  << " Infinite loop condidtion detected: going in cycles. Break after 6 cycles"
607  << std::endl;
608  if (result == SteppingHelixStateInfo::FAULT && nPoints > MAX_STEPS * 1. / defaultStep_)
609  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
610  << " Tired to go farther. Made too many steps: more than "
611  << MAX_STEPS * 1. / defaultStep_ << std::endl;
612  }
613 
614  if (debug_) {
615  switch (type) {
616  case RADIUS_DT:
617  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to radius "
618  << pars[RADIUS_P] << std::endl;
619  break;
620  case Z_DT:
621  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to z " << pars[Z_P]
622  << std::endl;
623  break;
624  case PATHL_DT:
625  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to pathL "
626  << pars[PATHL_P] << std::endl;
627  break;
628  case PLANE_DT: {
629  Point rPlane(pars[0], pars[1], pars[2]);
630  Vector nPlane(pars[3], pars[4], pars[5]);
631  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to plane r0:" << rPlane
632  << " n:" << nPlane << std::endl;
633  } break;
634  case POINT_PCA_DT: {
635  Point rDest(pars[0], pars[1], pars[2]);
636  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to PCA to point "
637  << rDest << std::endl;
638  } break;
639  case LINE_PCA_DT: {
640  Point rDest1(pars[0], pars[1], pars[2]);
641  Point rDest2(pars[3], pars[4], pars[5]);
642  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to PCA to line "
643  << rDest1 << " - " << rDest2 << std::endl;
644  } break;
645  default:
646  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to NOT IMPLEMENTED"
647  << std::endl;
648  break;
649  }
650  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Made " << nPoints - 1
651  << " steps and stopped at(cur step) " << svCurrent->r3 << " nOsc " << nOsc << std::endl;
652  }
653 
654  return result;
655 }
656 
657 void SteppingHelixPropagator::loadState(SteppingHelixPropagator::StateInfo& svCurrent,
658  const SteppingHelixPropagator::Vector& p3,
659  const SteppingHelixPropagator::Point& r3,
660  int charge,
661  PropagationDirection dir,
662  const AlgebraicSymMatrix55& covCurv) const {
663  static const std::string metname = "SteppingHelixPropagator";
664 
665  svCurrent.q = charge;
666  svCurrent.p3 = p3;
667  svCurrent.r3 = r3;
668  svCurrent.dir = dir == alongMomentum ? 1. : -1.;
669 
670  svCurrent.path_ = 0; // this could've held the initial path
671  svCurrent.radPath_ = 0;
672 
673  GlobalPoint gPointNorZ(svCurrent.r3.x(), svCurrent.r3.y(), svCurrent.r3.z());
674 
675  float gpmag = gPointNorZ.mag2();
676  float pmag2 = p3.mag2();
677  if (gpmag > 1e20f) {
678  LogTrace(metname) << "Initial point is too far";
679  svCurrent.isValid_ = false;
680  return;
681  }
682  if (pmag2 < 1e-18f) {
683  LogTrace(metname) << "Initial momentum is too low";
684  svCurrent.isValid_ = false;
685  return;
686  }
687  if (!(gpmag == gpmag)) {
688  LogTrace(metname) << "Initial point is a nan";
689  edm::LogWarning("SteppingHelixPropagatorNAN") << "Initial point is a nan";
690  svCurrent.isValid_ = false;
691  return;
692  }
693  if (!(pmag2 == pmag2)) {
694  LogTrace(metname) << "Initial momentum is a nan";
695  edm::LogWarning("SteppingHelixPropagatorNAN") << "Initial momentum is a nan";
696  svCurrent.isValid_ = false;
697  return;
698  }
699 
700  GlobalVector bf(0, 0, 0);
701  // = field_->inTesla(gPoint);
702  if (useMagVolumes_) {
703  if (vbField_) {
704  svCurrent.magVol = vbField_->findVolume(gPointNorZ);
705  if (useIsYokeFlag_) {
706  double curRad = svCurrent.r3.perp();
707  if (curRad > 380 && curRad < 850 && fabs(svCurrent.r3.z()) < 667) {
708  svCurrent.isYokeVol = isYokeVolume(svCurrent.magVol);
709  } else {
710  svCurrent.isYokeVol = false;
711  }
712  }
713  } else {
714  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
715  << "Failed to cast into VolumeBasedMagneticField: fall back to the default behavior"
716  << std::endl;
717  svCurrent.magVol = nullptr;
718  }
719  if (debug_) {
720  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Got volume at "
721  << svCurrent.magVol << std::endl;
722  }
723  }
724 
725  if (useMagVolumes_ && svCurrent.magVol != nullptr && !useInTeslaFromMagField_) {
726  bf = svCurrent.magVol->inTesla(gPointNorZ);
727  if (debug_) {
728  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Loaded bfield float: " << bf
729  << " at global float " << gPointNorZ << " double " << svCurrent.r3 << std::endl;
730  LocalPoint lPoint(svCurrent.magVol->toLocal(gPointNorZ));
731  LocalVector lbf = svCurrent.magVol->fieldInTesla(lPoint);
732  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "\t cf in local locF: " << lbf
733  << " at " << lPoint << std::endl;
734  }
735  svCurrent.bf.set(bf.x(), bf.y(), bf.z());
736  } else {
737  GlobalPoint gPoint(r3.x(), r3.y(), r3.z());
738  bf = field_->inTesla(gPoint);
739  svCurrent.bf.set(bf.x(), bf.y(), bf.z());
740  }
741  if (svCurrent.bf.mag2() < 1e-32)
742  svCurrent.bf.set(0., 0., 1e-16);
743  if (debug_) {
744  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
745  << "Loaded bfield double: " << svCurrent.bf << " from float: " << bf << " at float "
746  << gPointNorZ << " double " << svCurrent.r3 << std::endl;
747  }
748 
749  double dEdXPrime = 0;
750  double dEdx = 0;
751  double radX0 = 0;
752  MatBounds rzLims;
753  dEdx = getDeDx(svCurrent, dEdXPrime, radX0, rzLims);
754  svCurrent.dEdx = dEdx;
755  svCurrent.dEdXPrime = dEdXPrime;
756  svCurrent.radX0 = radX0;
757  svCurrent.rzLims = rzLims;
758 
759  svCurrent.covCurv = covCurv;
760 
761  svCurrent.isComplete = true;
762  svCurrent.isValid_ = true;
763 
764  if (debug_) {
765  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
766  << "Loaded at path: " << svCurrent.path() << " radPath: " << svCurrent.radPath() << " p3 "
767  << " pt: " << svCurrent.p3.perp() << " phi: " << svCurrent.p3.phi()
768  << " eta: " << svCurrent.p3.eta() << " " << svCurrent.p3 << " r3: " << svCurrent.r3
769  << " bField: " << svCurrent.bf.mag() << std::endl;
770  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
771  << "Input Covariance in Curvilinear RF " << covCurv << std::endl;
772  }
773 }
774 
775 void SteppingHelixPropagator::getNextState(const SteppingHelixPropagator::StateInfo& svPrevious,
776  SteppingHelixPropagator::StateInfo& svNext,
777  double dP,
778  const SteppingHelixPropagator::Vector& tau,
779  const SteppingHelixPropagator::Vector& drVec,
780  double dS,
781  double dX0,
782  const AlgebraicMatrix55& dCovCurvTransform) const {
783  static const std::string metname = "SteppingHelixPropagator";
784  svNext.q = svPrevious.q;
785  svNext.dir = dS > 0.0 ? 1. : -1.;
786  svNext.p3 = tau;
787  svNext.p3 *= (svPrevious.p3.mag() - svNext.dir * fabs(dP));
788 
789  svNext.r3 = svPrevious.r3;
790  svNext.r3 += drVec;
791 
792  svNext.path_ = svPrevious.path() + dS;
793  svNext.radPath_ = svPrevious.radPath() + dX0;
794 
795  GlobalPoint gPointNorZ(svNext.r3.x(), svNext.r3.y(), svNext.r3.z());
796 
797  GlobalVector bf(0, 0, 0);
798 
799  if (useMagVolumes_) {
800  if (vbField_ != nullptr) {
801  svNext.magVol = vbField_->findVolume(gPointNorZ);
802  if (useIsYokeFlag_) {
803  double curRad = svNext.r3.perp();
804  if (curRad > 380 && curRad < 850 && fabs(svNext.r3.z()) < 667) {
805  svNext.isYokeVol = isYokeVolume(svNext.magVol);
806  } else {
807  svNext.isYokeVol = false;
808  }
809  }
810  } else {
811  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
812  << "Failed to cast into VolumeBasedMagneticField" << std::endl;
813  svNext.magVol = nullptr;
814  }
815  if (debug_) {
816  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Got volume at "
817  << svNext.magVol << std::endl;
818  }
819  }
820 
821  if (useMagVolumes_ && svNext.magVol != nullptr && !useInTeslaFromMagField_) {
822  bf = svNext.magVol->inTesla(gPointNorZ);
823  svNext.bf.set(bf.x(), bf.y(), bf.z());
824  } else {
825  GlobalPoint gPoint(svNext.r3.x(), svNext.r3.y(), svNext.r3.z());
826  bf = field_->inTesla(gPoint);
827  svNext.bf.set(bf.x(), bf.y(), bf.z());
828  }
829  if (svNext.bf.mag2() < 1e-32)
830  svNext.bf.set(0., 0., 1e-16);
831 
832  double dEdXPrime = 0;
833  double dEdx = 0;
834  double radX0 = 0;
835  MatBounds rzLims;
836  dEdx = getDeDx(svNext, dEdXPrime, radX0, rzLims);
837  svNext.dEdx = dEdx;
838  svNext.dEdXPrime = dEdXPrime;
839  svNext.radX0 = radX0;
840  svNext.rzLims = rzLims;
841 
842  //update Emat only if it's valid
843  svNext.hasErrorPropagated_ = svPrevious.hasErrorPropagated_;
844  if (svPrevious.hasErrorPropagated_) {
845  {
846  AlgebraicMatrix55 tmp = dCovCurvTransform * svPrevious.covCurv;
847  ROOT::Math::AssignSym::Evaluate(svNext.covCurv, tmp * ROOT::Math::Transpose(dCovCurvTransform));
848 
849  svNext.covCurv += svPrevious.matDCovCurv;
850  }
851  } else {
852  //could skip dragging along the unprop. cov later .. now
853  // svNext.cov = svPrevious.cov;
854  }
855 
856  if (debug_) {
857  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Now at path: " << svNext.path()
858  << " radPath: " << svNext.radPath() << " p3 "
859  << " pt: " << svNext.p3.perp() << " phi: " << svNext.p3.phi() << " eta: " << svNext.p3.eta()
860  << " " << svNext.p3 << " r3: " << svNext.r3
861  << " dPhi: " << acos(svNext.p3.unit().dot(svPrevious.p3.unit())) << " bField: " << svNext.bf.mag()
862  << " dedx: " << svNext.dEdx << std::endl;
863  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "New Covariance "
864  << svNext.covCurv << std::endl;
865  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Transf by dCovTransform "
866  << dCovCurvTransform << std::endl;
867  }
868 }
869 
870 void SteppingHelixPropagator::setRep(SteppingHelixPropagator::Basis& rep,
871  const SteppingHelixPropagator::Vector& tau) const {
872  Vector zRep(0., 0., 1.);
873  rep.lX = tau;
874  rep.lY = zRep.cross(tau);
875  rep.lY *= 1. / tau.perp();
876  rep.lZ = rep.lX.cross(rep.lY);
877 }
878 
879 bool SteppingHelixPropagator::makeAtomStep(SteppingHelixPropagator::StateInfo& svCurrent,
880  SteppingHelixPropagator::StateInfo& svNext,
881  double dS,
882  PropagationDirection dir,
883  SteppingHelixPropagator::Fancy fancy) const {
884  static const std::string metname = "SteppingHelixPropagator";
885  if (debug_) {
886  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Make atom step "
887  << svCurrent.path() << " with step " << dS << " in direction " << dir << std::endl;
888  }
889 
890  AlgebraicMatrix55 dCCurvTransform(unit55_);
891 
892  double dP = 0;
893  double curP = svCurrent.p3.mag();
894  Vector tau = svCurrent.p3;
895  tau *= 1. / curP;
896  Vector tauNext(tau);
897  Vector drVec(0, 0, 0);
898 
899  dS = dir == alongMomentum ? fabs(dS) : -fabs(dS);
900 
901  double radX0 = 1e24;
902 
903  switch (fancy) {
904  case HEL_AS_F:
905  case HEL_ALL_F: {
906  double p0 = curP; // see above = svCurrent.p3.mag();
907  double p0Inv = 1. / p0;
908  double b0 = svCurrent.bf.mag();
909 
910  //get to the mid-point first
911  double phi = (2.99792458e-3 * svCurrent.q * b0 * p0Inv) * dS / 2.;
912  bool phiSmall = fabs(phi) < 1e-4;
913 
914  double cosPhi = 0;
915  double sinPhi = 0;
916 
917  double oneLessCosPhi = 0;
918  double oneLessCosPhiOPhi = 0;
919  double phiLessSinPhiOPhi = 0;
920 
921  if (phiSmall) {
922  double phi2 = phi * phi;
923  double phi3 = phi2 * phi;
924  double phi4 = phi3 * phi;
925  oneLessCosPhi = phi2 / 2. - phi4 / 24. + phi2 * phi4 / 720.; // 0.5*phi*phi;//*(1.- phi*phi/12.);
926  oneLessCosPhiOPhi = 0.5 * phi - phi3 / 24. + phi2 * phi3 / 720.; //*(1.- phi*phi/12.);
927  phiLessSinPhiOPhi = phi * phi / 6. - phi4 / 120. + phi4 * phi2 / 5040.; //*(1. - phi*phi/20.);
928  } else {
929  cosPhi = cos(phi);
930  sinPhi = sin(phi);
931  oneLessCosPhi = 1. - cosPhi;
932  oneLessCosPhiOPhi = oneLessCosPhi / phi;
933  double sinPhiOPhi = sinPhi / phi;
934  phiLessSinPhiOPhi = 1 - sinPhiOPhi;
935  }
936 
937  Vector bHat = svCurrent.bf;
938  bHat *= 1. / b0; //bHat.mag();
939  // bool bAlongZ = fabs(bHat.z()) > 0.9999;
940 
941  Vector btVec(bHat.cross(tau)); // for balong z btVec.z()==0
942  double tauB = tau.dot(bHat);
943  Vector bbtVec(bHat * tauB - tau); // (-tau.x(), -tau.y(), 0)
944 
945  //don't need it here tauNext = tau + bbtVec*oneLessCosPhi - btVec*sinPhi;
946  drVec = bbtVec * phiLessSinPhiOPhi;
947  drVec -= btVec * oneLessCosPhiOPhi;
948  drVec += tau;
949  drVec *= dS / 2.;
950 
951  double dEdx = svCurrent.dEdx;
952  double dEdXPrime = svCurrent.dEdXPrime;
953  radX0 = svCurrent.radX0;
954  dP = dEdx * dS;
955 
956  //improve with above values:
957  drVec += svCurrent.r3;
958  GlobalVector bfGV(0, 0, 0);
959  Vector bf(0, 0, 0);
960  if (useMagVolumes_ && svCurrent.magVol != nullptr && !useInTeslaFromMagField_) {
961  bfGV = svCurrent.magVol->inTesla(GlobalPoint(drVec.x(), drVec.y(), drVec.z()));
962  bf.set(bfGV.x(), bfGV.y(), bfGV.z());
963  } else {
964  bfGV = field_->inTesla(GlobalPoint(drVec.x(), drVec.y(), drVec.z()));
965  bf.set(bfGV.x(), bfGV.y(), bfGV.z());
966  }
967  double b0Init = b0;
968  b0 = bf.mag();
969  if (b0 < 1e-16) {
970  b0 = 1e-16;
971  bf.set(0., 0., 1e-16);
972  }
973  if (debug_) {
974  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Improved b " << b0
975  << " at r3 " << drVec << std::endl;
976  }
977 
978  if (fabs((b0 - b0Init) * dS) > 1) {
979  //missed the mag volume boundary?
980  if (debug_) {
981  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Large bf*dS change "
982  << fabs((b0 - svCurrent.bf.mag()) * dS) << " --> recalc dedx" << std::endl;
983  }
984  svNext.r3 = drVec;
985  svNext.bf = bf;
986  svNext.p3 = svCurrent.p3;
987  svNext.isYokeVol = svCurrent.isYokeVol;
988  svNext.magVol = svCurrent.magVol;
989  MatBounds rzTmp;
990  dEdx = getDeDx(svNext, dEdXPrime, radX0, rzTmp);
991  dP = dEdx * dS;
992  if (debug_) {
993  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "New dEdX= " << dEdx
994  << " dP= " << dP << " for p0 " << p0 << std::endl;
995  }
996  }
997  //p0 is mid-way and b0 from mid-point
998  p0 += dP / 2.;
999  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;
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 }
JacobianCartesianToCurvilinear.h
Vector3DBase
Definition: Vector3DBase.h:8
SteppingHelixPropagator::Point
CLHEP::Hep3Vector Point
Definition: SteppingHelixPropagator.h:39
MagVolume.h
Cone::openingAngle
Geom::Theta< float > openingAngle() const
Angle of the cone.
Definition: Cone.h:50
SteppingHelixPropagator::CONE_DT
Definition: SteppingHelixPropagator.h:56
SteppingHelixPropagator::useIsYokeFlag_
bool useIsYokeFlag_
Definition: SteppingHelixPropagator.h:273
MagVolume::inTesla
::GlobalVector inTesla(const ::GlobalPoint &gp) const override
Definition: MagVolume.h:42
MagneticField::inTesla
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
SteppingHelixStateInfo::getFreeState
void getFreeState(FreeTrajectoryState &fts) const
convert internal structure into the fts
Definition: SteppingHelixStateInfo.cc:88
SteppingHelixPropagator::MatBounds
SteppingHelixStateInfo::VolumeBounds MatBounds
Definition: SteppingHelixPropagator.h:170
anyDirection
Definition: PropagationDirection.h:4
mps_fire.i
i
Definition: mps_fire.py:428
TkRotation::zy
T zy() const
Definition: extTkRotation.h:258
pwdgSkimBPark_cfi.beamSpot
beamSpot
Definition: pwdgSkimBPark_cfi.py:5
MessageLogger.h
Cylinder.h
Cylinder::radius
Scalar radius() const
Radius of the cylinder.
Definition: Cylinder.h:64
SteppingHelixPropagator.h
SteppingHelixStateInfo::RANGEOUT
Definition: SteppingHelixStateInfo.h:34
SteppingHelixStateInfo::radX0
double radX0
Definition: SteppingHelixStateInfo.h:101
SteppingHelixStateInfo::dEdXPrime
double dEdXPrime
Definition: SteppingHelixStateInfo.h:100
SteppingHelixPropagator::ecShiftPos_
double ecShiftPos_
Definition: SteppingHelixPropagator.h:282
JacobianCurvilinearToCartesian.h
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
SteppingHelixPropagator::initStateArraySHPSpecific
void initStateArraySHPSpecific(StateArray &svBuf, bool flagsOnly) const
(Internals) set defaults needed for states used inside propagate methods
Definition: SteppingHelixPropagator.cc:38
metsig::tau
Definition: SignAlgoResolutions.h:49
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
SteppingHelixPropagator::FtsPP
std::pair< FreeTrajectoryState, double > FtsPP
Definition: SteppingHelixPropagator.h:258
PropagationException
Common base class.
Definition: PropagationExceptions.h:14
SteppingHelixStateInfo::bf
Vector bf
Definition: SteppingHelixStateInfo.h:92
SteppingHelixStateInfo::hasErrorPropagated_
bool hasErrorPropagated_
Definition: SteppingHelixStateInfo.h:105
SteppingHelixPropagator::DestType
DestType
Definition: SteppingHelixPropagator.h:52
SteppingHelixPropagator::POINT_PCA_DT
Definition: SteppingHelixPropagator.h:59
SteppingHelixStateInfo::radPath_
double radPath_
Definition: SteppingHelixStateInfo.h:90
SteppingHelixPropagator::setIState
void setIState(const SteppingHelixStateInfo &sStart, StateArray &svBuff, int &nPoints) const
(Internals) Init starting point
Definition: SteppingHelixPropagator.cc:296
SteppingHelixStateInfo::ResultName
static const std::string ResultName[MAX_RESULT]
Definition: SteppingHelixStateInfo.h:36
SteppingHelixPropagator::field_
const MagneticField * field_
Definition: SteppingHelixPropagator.h:265
SteppingHelixPropagator::Vector
CLHEP::Hep3Vector Vector
Definition: SteppingHelixPropagator.h:38
oppositeToMomentum
Definition: PropagationDirection.h:4
CustomPhysics_cfi.gamma
gamma
Definition: CustomPhysics_cfi.py:17
Validation_hcalonly_cfi.sign
sign
Definition: Validation_hcalonly_cfi.py:32
SteppingHelixPropagator::refToDest
Result refToDest(DestType dest, const SteppingHelixPropagator::StateInfo &sv, const double pars[6], double &dist, double &tanDist, PropagationDirection &refDirection, double fastSkipDist=1e12) const
(Internals) determine distance and direction from the current position to the plane
Definition: SteppingHelixPropagator.cc:1781
Surface
Definition: Surface.h:36
Cone::vertex
GlobalPoint vertex() const
Global position of the cone vertex.
Definition: Cone.h:47
SteppingHelixStateInfo::p3
Vector p3
Definition: SteppingHelixStateInfo.h:85
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
PV3DBase::mag2
T mag2() const
Definition: PV3DBase.h:63
SteppingHelixStateInfo::dir
double dir
Definition: SteppingHelixStateInfo.h:91
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
SteppingHelixStateInfo::magVol
const MagVolume * magVol
Definition: SteppingHelixStateInfo.h:94
plot_hgcal_utils.dEdx
dEdx
Definition: plot_hgcal_utils.py:213
SteppingHelixStateInfo::OK
Definition: SteppingHelixStateInfo.h:34
SteppingHelixStateInfo::matDCovCurv
AlgebraicSymMatrix55 matDCovCurv
Definition: SteppingHelixStateInfo.h:88
SteppingHelixPropagator::SteppingHelixPropagator
SteppingHelixPropagator()
Constructors.
Definition: SteppingHelixPropagator.cc:56
geometryDiff.epsilon
int epsilon
Definition: geometryDiff.py:26
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
SteppingHelixStateInfo::getStateOnSurface
TrajectoryStateOnSurface getStateOnSurface(const Surface &surf, bool returnTangentPlane=false) const
Definition: SteppingHelixStateInfo.cc:60
SteppingHelixPropagator::returnTangentPlane_
bool returnTangentPlane_
Definition: SteppingHelixPropagator.h:276
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
SteppingHelixPropagator::Z_P
Definition: SteppingHelixPropagator.h:50
AlgebraicMatrixID
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
Definition: AlgebraicROOTObjects.h:72
SteppingHelixStateInfo::UNDEFINED
Definition: SteppingHelixStateInfo.h:34
Propagator
Definition: Propagator.h:44
Plane.h
SteppingHelixPropagator::makeAtomStep
bool makeAtomStep(SteppingHelixPropagator::StateInfo &svCurrent, SteppingHelixPropagator::StateInfo &svNext, double dS, PropagationDirection dir, SteppingHelixPropagator::Fancy fancy) const
main stepping function: compute next state vector after a step of length dS
Definition: SteppingHelixPropagator.cc:879
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
VolumeBasedMagneticField::findVolume
const MagVolume * findVolume(const GlobalPoint &gp) const
Definition: VolumeBasedMagneticField.cc:59
SteppingHelixPropagator::PATHL_P
Definition: SteppingHelixPropagator.h:50
AlgebraicMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepStd< double, 5, 5 > > AlgebraicMatrix55
Definition: AlgebraicROOTObjects.h:55
SteppingHelixPropagator::getNextState
void getNextState(const SteppingHelixPropagator::StateInfo &svPrevious, SteppingHelixPropagator::StateInfo &svNext, double dP, const SteppingHelixPropagator::Vector &tau, const SteppingHelixPropagator::Vector &drVec, double dS, double dX0, const AlgebraicMatrix55 &dCovCurv) const
Definition: SteppingHelixPropagator.cc:775
SteppingHelixStateInfo::isComplete
bool isComplete
Definition: SteppingHelixStateInfo.h:103
SteppingHelixStateInfo::INACC
Definition: SteppingHelixStateInfo.h:34
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
SteppingHelixPropagator::applyRadX0Correction_
bool applyRadX0Correction_
Definition: SteppingHelixPropagator.h:271
Propagator::propagationDirection
virtual PropagationDirection propagationDirection() const final
Definition: Propagator.h:139
pfDeepBoostedJetPreprocessParams_cfi.sv
sv
Definition: pfDeepBoostedJetPreprocessParams_cfi.py:352
reco::BeamSpot
Definition: BeamSpot.h:21
VolumeBasedMagneticField.h
Geom::pi
constexpr double pi()
Definition: Pi.h:31
SteppingHelixPropagator::LINE_PCA_DT
Definition: SteppingHelixPropagator.h:60
SteppingHelixStateInfo::q
int q
Definition: SteppingHelixStateInfo.h:84
SteppingHelixPropagator::~SteppingHelixPropagator
~SteppingHelixPropagator() override
Destructor.
Definition: SteppingHelixPropagator.cc:54
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
SteppingHelixPropagator::MAX_STEPS
static const int MAX_STEPS
Definition: SteppingHelixPropagator.h:259
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
SteppingHelixStateInfo::FAULT
Definition: SteppingHelixStateInfo.h:34
Point3DBase< float, GlobalTag >
MagVolume::inside
virtual bool inside(const GlobalPoint &gp, double tolerance=0.) const =0
SteppingHelixPropagator::defaultStep_
double defaultStep_
Definition: SteppingHelixPropagator.h:280
SteppingHelixPropagator::useTuningForL2Speed_
bool useTuningForL2Speed_
Definition: SteppingHelixPropagator.h:278
Cone
Definition: Cone.h:17
SteppingHelixPropagator::TsosPP
std::pair< TrajectoryStateOnSurface, double > TsosPP
Definition: SteppingHelixPropagator.h:257
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
SteppingHelixPropagator::invalidState_
StateInfo invalidState_
Definition: SteppingHelixPropagator.h:263
SteppingHelixStateInfo
Definition: SteppingHelixStateInfo.h:27
SteppingHelixPropagator::noMaterialMode_
bool noMaterialMode_
Definition: SteppingHelixPropagator.h:269
Cone.h
ALCARECOTkAlJpsiMuMu_cff.charge
charge
Definition: ALCARECOTkAlJpsiMuMu_cff.py:47
SteppingHelixStateInfo::VolumeBounds
Definition: SteppingHelixStateInfo.h:71
Point
Structure Point Contains parameters of Gaussian fits to DMRs.
Definition: DMRtrends.cc:57
SteppingHelixPropagator::RADIUS_DT
Definition: SteppingHelixPropagator.h:53
SteppingHelixPropagator::isYokeVolume
bool isYokeVolume(const MagVolume *vol) const
check if it's a yoke/iron based on this MagVol internals
Definition: SteppingHelixPropagator.cc:2432
SteppingHelixStateInfo::radPath
double radPath() const
Definition: SteppingHelixStateInfo.h:63
SteppingHelixStateInfo::isValid_
bool isValid_
Definition: SteppingHelixStateInfo.h:104
SteppingHelixPropagator::HEL_ALL_F
Definition: SteppingHelixPropagator.h:66
type
type
Definition: SiPixelVCal_PayloadInspector.cc:37
SteppingHelixStateInfo::path
double path() const
Definition: SteppingHelixStateInfo.h:62
SteppingHelixPropagator::setRep
void setRep(SteppingHelixPropagator::Basis &rep, const SteppingHelixPropagator::Vector &tau) const
Set/compute basis vectors for local coordinates at current step (called by incrementState)
Definition: SteppingHelixPropagator.cc:870
SteppingHelixPropagator::useInTeslaFromMagField_
bool useInTeslaFromMagField_
Definition: SteppingHelixPropagator.h:275
SteppingHelixPropagator::getDeDx
double getDeDx(const SteppingHelixPropagator::StateInfo &sv, double &dEdXPrime, double &radX0, MatBounds &rzLims) const
estimate average (in fact smth. close to MPV and median) energy loss per unit path length
Definition: SteppingHelixPropagator.cc:1317
SteppingHelixPropagator::PATHL_DT
Definition: SteppingHelixPropagator.h:58
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
MagVolume
Definition: MagVolume.h:13
SteppingHelixStateInfo::r3
Point r3
Definition: SteppingHelixStateInfo.h:86
SteppingHelixStateInfo::path_
double path_
Definition: SteppingHelixStateInfo.h:89
cuy.rep
rep
Definition: cuy.py:1190
SteppingHelixPropagator::propagateWithPath
std::pair< TrajectoryStateOnSurface, double > propagateWithPath(const FreeTrajectoryState &ftsStart, const Plane &pDest) const override
Definition: SteppingHelixPropagator.cc:79
SteppingHelixPropagator::vbField_
const VolumeBasedMagneticField * vbField_
Definition: SteppingHelixPropagator.h:266
SteppingHelixStateInfo::dEdx
double dEdx
Definition: SteppingHelixStateInfo.h:99
SteppingHelixPropagator::UNDEFINED_DT
Definition: SteppingHelixPropagator.h:61
SteppingHelixPropagator::useMagVolumes_
bool useMagVolumes_
Definition: SteppingHelixPropagator.h:272
SteppingHelixPropagator::useMatVolumes_
bool useMatVolumes_
Definition: SteppingHelixPropagator.h:274
SteppingHelixPropagator::loadState
void loadState(SteppingHelixPropagator::StateInfo &svCurrent, const SteppingHelixPropagator::Vector &p3, const SteppingHelixPropagator::Point &r3, int charge, PropagationDirection dir, const AlgebraicSymMatrix55 &covCurv) const
Definition: SteppingHelixPropagator.cc:657
submitPVResolutionJobs.q
q
Definition: submitPVResolutionJobs.py:84
b0
static constexpr float b0
Definition: L1EGammaCrystalsEmulatorProducer.cc:82
MagneticField.h
TkRotation::zz
T zz() const
Definition: extTkRotation.h:259
SteppingHelixStateInfo::Result
Result
Definition: SteppingHelixStateInfo.h:34
SteppingHelixStateInfo::isValid
bool isValid() const
Definition: SteppingHelixStateInfo.h:65
SteppingHelixPropagator::RADIUS_P
Definition: SteppingHelixPropagator.h:50
SteppingHelixPropagator::Fancy
Fancy
Definition: SteppingHelixPropagator.h:64
SiPixelPhase1Clusters_cfi.e3
e3
Definition: SiPixelPhase1Clusters_cfi.py:9
DDAxes::phi
GloballyPositioned::position
const PositionType & position() const
Definition: GloballyPositioned.h:36
MagVolume::faces
virtual const std::vector< VolumeSide > & faces() const =0
Access to volume faces.
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:127
SteppingHelixPropagator::refToMagVolume
Result refToMagVolume(const SteppingHelixPropagator::StateInfo &sv, PropagationDirection dir, double &dist, double &tanDist, double fastSkipDist=1e12, bool expectNewMagVolume=false, double maxStep=1e12) const
Definition: SteppingHelixPropagator.cc:2049
SteppingHelixPropagator::refToMatVolume
Result refToMatVolume(const SteppingHelixPropagator::StateInfo &sv, PropagationDirection dir, double &dist, double &tanDist, double fastSkipDist=1e12) const
Definition: SteppingHelixPropagator.cc:2308
SteppingHelixPropagator::StateArray
StateInfo StateArray[MAX_POINTS+1]
Definition: SteppingHelixPropagator.h:172
SteppingHelixPropagator::MAX_POINTS
static const int MAX_POINTS
Definition: SteppingHelixPropagator.h:171
FreeTrajectoryState
Definition: FreeTrajectoryState.h:27
SteppingHelixPropagator::Z_DT
Definition: SteppingHelixPropagator.h:54
GloballyPositioned::toLocal
LocalPoint toLocal(const GlobalPoint &gp) const
Definition: GloballyPositioned.h:98
MFGrid.h
PropagationDirection
PropagationDirection
Definition: PropagationDirection.h:4
Plane
Definition: Plane.h:16
SteppingHelixStateInfo::isYokeVol
bool isYokeVol
Definition: SteppingHelixStateInfo.h:95
MagVolume::isIron
bool isIron() const
Temporary hack to pass information on material. Will eventually be replaced!
Definition: MagVolume.h:45
SteppingHelixPropagator::ecShiftNeg_
double ecShiftNeg_
Definition: SteppingHelixPropagator.h:283
SteppingHelixPropagator::debug_
bool debug_
Definition: SteppingHelixPropagator.h:268
SteppingHelixPropagator::propagate
void propagate(const SteppingHelixStateInfo &ftsStart, const Surface &sDest, SteppingHelixStateInfo &out) const
Propagate to Plane given a starting point.
Definition: SteppingHelixPropagator.cc:152
SteppingHelixPropagator::noErrorPropagation_
bool noErrorPropagation_
Definition: SteppingHelixPropagator.h:270
p3
double p3[4]
Definition: TauolaWrapper.h:91
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
SteppingHelixPropagator::PLANE_DT
Definition: SteppingHelixPropagator.h:55
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
SteppingHelixStateInfo::covCurv
AlgebraicSymMatrix55 covCurv
Definition: SteppingHelixStateInfo.h:87
Cylinder
Definition: Cylinder.h:19
SteppingHelixPropagator::HEL_AS_F
Definition: SteppingHelixPropagator.h:65
mps_fire.result
result
Definition: mps_fire.py:311
SteppingHelixPropagator::sendLogWarning_
bool sendLogWarning_
Definition: SteppingHelixPropagator.h:277
SteppingHelixPropagator::unit55_
const AlgebraicSymMatrix55 unit55_
Definition: SteppingHelixPropagator.h:267
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:224
SteppingHelixPropagator::Basis
Definition: SteppingHelixPropagator.h:44
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
createTree.pp
pp
Definition: createTree.py:17
MagVolume::fieldInTesla
LocalVector fieldInTesla(const LocalPoint &lp) const
Definition: MagVolume.cc:11
GloballyPositioned::rotation
const RotationType & rotation() const
Definition: GloballyPositioned.h:38
SteppingHelixStateInfo::rzLims
VolumeBounds rzLims
Definition: SteppingHelixStateInfo.h:98
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
SteppingHelixStateInfo::NOT_IMPLEMENTED
Definition: SteppingHelixStateInfo.h:34
PropagationExceptions.h
SteppingHelixStateInfo::WRONG_VOLUME
Definition: SteppingHelixStateInfo.h:34
MagneticField
Definition: MagneticField.h:19
AlgebraicSymMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
Definition: AlgebraicROOTObjects.h:23
edm::Log
Definition: MessageLogger.h:70
TkRotation::zx
T zx() const
Definition: extTkRotation.h:257
PVValHelper::dx
Definition: PVValidationHelpers.h:48
alongMomentum
Definition: PropagationDirection.h:4
geometryCSVtoXML.xx
xx
Definition: geometryCSVtoXML.py:19
vertexPlots.e4
e4
Definition: vertexPlots.py:64
SteppingHelixStateInfo::field
const MagneticField * field
Definition: SteppingHelixStateInfo.h:96
SteppingHelixStateInfo::status_
Result status_
Definition: SteppingHelixStateInfo.h:107
mps_fire.dest
dest
Definition: mps_fire.py:179
DeadROC_duringRun.dir
dir
Definition: DeadROC_duringRun.py:23
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
metname
const std::string metname
Definition: MuonSeedOrcaPatternRecognition.cc:40
SteppingHelixPropagator::cIndex_
int cIndex_(int ind) const
(Internals) circular index for array of transients
Definition: SteppingHelixPropagator.cc:1773
Generated for CMSSW Reference Manual by doxygen 1.8.16