CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SteppingHelixPropagator.cc
Go to the documentation of this file.
1 
13 //
14 // Original Author: Vyacheslav Krutelyov
15 // Created: Fri Mar 3 16:01:24 CST 2006
16 // $Id: SteppingHelixPropagator.cc,v 1.74 2010/10/03 17:23:11 elmer Exp $
17 //
18 //
19 
20 
25 
29 
31 
34 
36 
38 
39 #include <sstream>
40 #include <typeinfo>
41 
44 {
45  field_ = 0;
46 }
47 
50  Propagator(dir),
51  unit55_(AlgebraicMatrixID())
52 {
53  field_ = field;
54  vbField_ = dynamic_cast<const VolumeBasedMagneticField*>(field_);
57  debug_ = false;
58  noMaterialMode_ = false;
59  noErrorPropagation_ = false;
60  applyRadX0Correction_ = true;
61  useMagVolumes_ = true;
62  useIsYokeFlag_ = true;
63  useMatVolumes_ = true;
64  useInTeslaFromMagField_ = false; //overrides behavior only if true
65  returnTangentPlane_ = true;
66  sendLogWarning_ = false;
67  useTuningForL2Speed_ = false;
68  for (int i = 0; i <= MAX_POINTS; i++){
69  svBuf_[i].p3 = Vector(0,0,0);
70  svBuf_[i].r3 = Point(0,0,0);
71  svBuf_[i].bf = Vector(0,0,0);
72  svBuf_[i].bfGradLoc = Vector(0,0,0);
75  svBuf_[i].isComplete = true;
76  svBuf_[i].isValid_ = true;
78  }
79  defaultStep_ = 5.;
80 
81  ecShiftPos_ = 0;
82  ecShiftNeg_ = 0;
83 
84 }
85 
87 SteppingHelixPropagator::propagate(const FreeTrajectoryState& ftsStart, const Plane& pDest) const {
88  return propagateWithPath(ftsStart, pDest).first;
89 }
90 
93 {
94  return propagateWithPath(ftsStart, cDest).first;
95 }
96 
99 {
100  return propagateWithPath(ftsStart, pDest).first;
101 }
102 
105  const GlobalPoint& pDest1, const GlobalPoint& pDest2) const
106 {
107  return propagateWithPath(ftsStart, pDest1, pDest2).first;
108 }
109 
112  const reco::BeamSpot& beamSpot) const
113 {
114  return propagateWithPath(ftsStart, beamSpot).first;
115 }
116 
117 
118 std::pair<TrajectoryStateOnSurface, double>
120  const Plane& pDest) const {
121 
123 
124  const StateInfo& svCurrent = propagate(svBuf_[0], pDest);
125 
126  return TsosPP(svCurrent.getStateOnSurface(pDest), svCurrent.path());
127 }
128 
129 std::pair<TrajectoryStateOnSurface, double>
131  const Cylinder& cDest) const {
132 
134 
135  const StateInfo& svCurrent = propagate(svBuf_[0], cDest);
136 
137  return TsosPP(svCurrent.getStateOnSurface(cDest, returnTangentPlane_), svCurrent.path());
138 }
139 
140 
141 std::pair<FreeTrajectoryState, double>
143  const GlobalPoint& pDest) const {
145 
146  const StateInfo& svCurrent = propagate(svBuf_[0], pDest);
147 
148  FreeTrajectoryState ftsDest;
149  svCurrent.getFreeState(ftsDest);
150 
151  return FtsPP(ftsDest, svCurrent.path());
152 }
153 
154 std::pair<FreeTrajectoryState, double>
156  const GlobalPoint& pDest1, const GlobalPoint& pDest2) const {
157 
158  if ((pDest1-pDest2).mag() < 1e-10){
159  if (sendLogWarning_){
160  edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: the points should be at a bigger distance"
161  <<std::endl;
162  }
163  return FtsPP();
164  }
166 
167  const StateInfo& svCurrent = propagate(svBuf_[0], pDest1, pDest2);
168 
169  FreeTrajectoryState ftsDest;
170  svCurrent.getFreeState(ftsDest);
171 
172  return FtsPP(ftsDest, svCurrent.path());
173 }
174 
175 
176 std::pair<FreeTrajectoryState, double>
178  const reco::BeamSpot& beamSpot) const {
179  GlobalPoint pDest1(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
180  GlobalPoint pDest2(pDest1.x() + beamSpot.dxdz()*1e3,
181  pDest1.y() + beamSpot.dydz()*1e3,
182  pDest1.z() + 1e3);
183  return propagateWithPath(ftsStart, pDest1, pDest2);
184 }
185 
188  const Surface& sDest) const {
189 
190  if (! sStart.isValid()){
191  if (sendLogWarning_){
192  edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: invalid input state"
193  <<std::endl;
194  }
195  return invalidState_;
196  }
197 
198  const Plane* pDest = dynamic_cast<const Plane*>(&sDest);
199  if (pDest != 0) return propagate(sStart, *pDest);
200 
201  const Cylinder* cDest = dynamic_cast<const Cylinder*>(&sDest);
202  if (cDest != 0) return propagate(sStart, *cDest);
203 
204  throw PropagationException("The surface is neither Cylinder nor Plane");
205 
206 }
207 
210  const Plane& pDest) const {
211 
212  if (! sStart.isValid()){
213  if (sendLogWarning_){
214  edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: invalid input state"
215  <<std::endl;
216  }
217  return invalidState_;
218  }
219  setIState(sStart);
220 
221  Point rPlane(pDest.position().x(), pDest.position().y(), pDest.position().z());
222  Vector nPlane(pDest.rotation().zx(), pDest.rotation().zy(), pDest.rotation().zz()); nPlane /= nPlane.mag();
223 
224  double pars[6] = { rPlane.x(), rPlane.y(), rPlane.z(),
225  nPlane.x(), nPlane.y(), nPlane.z() };
226 
227  propagate(PLANE_DT, pars);
228 
229  //(re)set it before leaving: dir =1 (-1) if path increased (decreased) and 0 if it didn't change
230  //need to implement this somewhere else as a separate function
231  double lDir = 0;
232  if (sStart.path() < svBuf_[cIndex_(nPoints_-1)].path()) lDir = 1.;
233  if (sStart.path() > svBuf_[cIndex_(nPoints_-1)].path()) lDir = -1.;
234  svBuf_[cIndex_(nPoints_-1)].dir = lDir;
235  return svBuf_[cIndex_(nPoints_-1)];
236 }
237 
240  const Cylinder& cDest) const {
241 
242  if (! sStart.isValid()){
243  if (sendLogWarning_){
244  edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: invalid input state"
245  <<std::endl;
246  }
247  return invalidState_;
248  }
249  setIState(sStart);
250 
251  double pars[6] = {0,0,0,0,0,0};
252  pars[RADIUS_P] = cDest.radius();
253 
254 
255  propagate(RADIUS_DT, pars);
256 
257  //(re)set it before leaving: dir =1 (-1) if path increased (decreased) and 0 if it didn't change
258  //need to implement this somewhere else as a separate function
259  double lDir = 0;
260  if (sStart.path() < svBuf_[cIndex_(nPoints_-1)].path()) lDir = 1.;
261  if (sStart.path() > svBuf_[cIndex_(nPoints_-1)].path()) lDir = -1.;
262  svBuf_[cIndex_(nPoints_-1)].dir = lDir;
263  return svBuf_[cIndex_(nPoints_-1)];
264 }
265 
268  const GlobalPoint& pDest) const {
269 
270  if (! sStart.isValid()){
271  if (sendLogWarning_){
272  edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: invalid input state"
273  <<std::endl;
274  }
275  return invalidState_;
276  }
277  setIState(sStart);
278 
279  double pars[6] = {pDest.x(), pDest.y(), pDest.z(), 0, 0, 0};
280 
281 
282  propagate(POINT_PCA_DT, pars);
283 
284  return svBuf_[cIndex_(nPoints_-1)];
285 }
286 
289  const GlobalPoint& pDest1, const GlobalPoint& pDest2) const {
290 
291  if ((pDest1-pDest2).mag() < 1e-10 || !sStart.isValid()){
292  if (sendLogWarning_){
293  if ((pDest1-pDest2).mag() < 1e-10)
294  edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: points are too close"
295  <<std::endl;
296  if (!sStart.isValid())
297  edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: invalid input state"
298  <<std::endl;
299  }
300  return invalidState_;
301  }
302  setIState(sStart);
303 
304  double pars[6] = {pDest1.x(), pDest1.y(), pDest1.z(),
305  pDest2.x(), pDest2.y(), pDest2.z()};
306 
307  propagate(LINE_PCA_DT, pars);
308 
309  return svBuf_[cIndex_(nPoints_-1)];
310 }
311 
313  nPoints_ = 0;
314  svBuf_[cIndex_(nPoints_)] = sStart; //do this anyways to have a fresh start
315  if (sStart.isComplete ) {
316  svBuf_[cIndex_(nPoints_)] = sStart;
317  nPoints_++;
318  } else {
319  loadState(svBuf_[cIndex_(nPoints_)], sStart.p3, sStart.r3, sStart.q,
320  propagationDirection(), sStart.covCurv);
321  nPoints_++;
322  }
324 }
325 
328  const double pars[6], double epsilon) const{
329 
330  static const std::string metname = "SteppingHelixPropagator";
331  StateInfo* svCurrent = &svBuf_[cIndex_(nPoints_-1)];
332 
333  //check if it's going to work at all
334  double tanDist = 0;
335  double dist = 0;
336  PropagationDirection refDirection = anyDirection;
337  Result result = refToDest(type, (*svCurrent), pars, dist, tanDist, refDirection);
338 
339  if (result != SteppingHelixStateInfo::OK || fabs(dist) > 1e6){
340  svCurrent->status_ = result;
341  if (fabs(dist) > 1e6) svCurrent->status_ = SteppingHelixStateInfo::INACC;
342  svCurrent->isValid_ = false;
343  svCurrent->field = field_;
344  if (sendLogWarning_){
345  edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<" Failed after first refToDest check with status "
347  <<std::endl;
348  }
349  return result;
350  }
351 
353  bool makeNextStep = true;
354  double dStep = defaultStep_;
355  PropagationDirection dir,oldDir;
356  dir = propagationDirection();
357  oldDir = dir;
358  int nOsc = 0;
359 
360  double distMag = 1e12;
361  double tanDistMag = 1e12;
362 
363  double distMat = 1e12;
364  double tanDistMat = 1e12;
365 
366  double tanDistNextCheck = -0.1;//just need a negative start val
367  double tanDistMagNextCheck = -0.1;
368  double tanDistMatNextCheck = -0.1;
369  double oldDStep = 0;
370  PropagationDirection oldRefDirection = propagationDirection();
371 
374 
375  bool isFirstStep = true;
376  bool expectNewMagVolume = false;
377 
378  int loopCount = 0;
379  while (makeNextStep){
380  dStep = defaultStep_;
381  svCurrent = &svBuf_[cIndex_(nPoints_-1)];
382  double curZ = svCurrent->r3.z();
383  double curR = svCurrent->r3.perp();
384  if ( fabs(curZ) < 440 && curR < 260) dStep = defaultStep_*2;
385 
386  //more such ifs might be scattered around
387  //even though dStep is large, it will still make stops at each volume boundary
389  dStep = 100.;
390  if (! useIsYokeFlag_ && fabs(curZ) < 667 && curR > 380 && curR < 850){
391  dStep = 5*(1+0.2*svCurrent->p3.mag());
392  }
393  }
394 
395  // refDirection = propagationDirection();
396 
397  tanDistNextCheck -= oldDStep;
398  tanDistMagNextCheck -= oldDStep;
399  tanDistMatNextCheck -= oldDStep;
400 
401  if (tanDistNextCheck < 0){
402  //use pre-computed values if it's the first step
403  if (! isFirstStep) refToDest(type, (*svCurrent), pars, dist, tanDist, refDirection);
404  // constrain allowed path for a tangential approach
405  if (fabs(tanDist) > 4.*(fabs(dist)) ) tanDist *= tanDist == 0 ? 0 :fabs(dist/tanDist*4.);
406 
407  tanDistNextCheck = fabs(tanDist)*0.5 - 0.5; //need a better guess (to-do)
408  //reasonable limit
409  if (tanDistNextCheck > defaultStep_*20. ) tanDistNextCheck = defaultStep_*20.;
410  oldRefDirection = refDirection;
411  } else {
412  tanDist = tanDist > 0. ? tanDist - oldDStep : tanDist + oldDStep;
413  refDirection = oldRefDirection;
414  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Skipped refToDest: guess tanDist = "<<tanDist
415  <<" next check at "<<tanDistNextCheck<<std::endl;
416  }
418  double fastSkipDist = fabs(dist) > fabs(tanDist) ? tanDist : dist;
419 
421  dir = refDirection;
422  } else {
423  dir = propagationDirection();
424  if (fabs(tanDist)<0.1 && refDirection != dir ){
425  //how did it get here? nOsc++;
426  dir = refDirection;
427  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"NOTE: overstepped last time: switch direction (can do it if within 1 mm)"<<std::endl;
428  }
429  }
430 
431  if (useMagVolumes_ && ! (fabs(dist) < fabs(epsilon))){//need to know the general direction
432  if (tanDistMagNextCheck < 0){
433  resultToMag = refToMagVolume((*svCurrent), dir, distMag, tanDistMag, fabs(fastSkipDist), expectNewMagVolume, fabs(tanDist));
434  // constrain allowed path for a tangential approach
435  if (fabs(tanDistMag) > 4.*(fabs(distMag)) ) tanDistMag *= tanDistMag == 0 ? 0 : fabs(distMag/tanDistMag*4.);
436 
437  tanDistMagNextCheck = fabs(tanDistMag)*0.5-0.5; //need a better guess (to-do)
438  //reasonable limit; "turn off" checking if bounds are further than the destination
439  if (tanDistMagNextCheck > defaultStep_*20.
440  || fabs(dist) < fabs(distMag)
441  || resultToMag ==SteppingHelixStateInfo::INACC)
442  tanDistMagNextCheck = defaultStep_*20 > fabs(fastSkipDist) ? fabs(fastSkipDist) : defaultStep_*20;
443  if (resultToMag != SteppingHelixStateInfo::INACC
444  && resultToMag != SteppingHelixStateInfo::OK) tanDistMagNextCheck = -1;
445  } else {
446  // resultToMag = SteppingHelixStateInfo::OK;
447  tanDistMag = tanDistMag > 0. ? tanDistMag - oldDStep : tanDistMag + oldDStep;
448  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Skipped refToMag: guess tanDistMag = "<<tanDistMag
449  <<" next check at "<<tanDistMagNextCheck;
450  }
451  }
452 
453  if (useMatVolumes_ && ! (fabs(dist) < fabs(epsilon))){//need to know the general direction
454  if (tanDistMatNextCheck < 0){
455  resultToMat = refToMatVolume((*svCurrent), dir, distMat, tanDistMat, fabs(fastSkipDist));
456  // constrain allowed path for a tangential approach
457  if (fabs(tanDistMat) > 4.*(fabs(distMat)) ) tanDistMat *= tanDistMat == 0 ? 0 : fabs(distMat/tanDistMat*4.);
458 
459  tanDistMatNextCheck = fabs(tanDistMat)*0.5-0.5; //need a better guess (to-do)
460  //reasonable limit; "turn off" checking if bounds are further than the destination
461  if (tanDistMatNextCheck > defaultStep_*20.
462  || fabs(dist) < fabs(distMat)
463  || resultToMat ==SteppingHelixStateInfo::INACC )
464  tanDistMatNextCheck = defaultStep_*20 > fabs(fastSkipDist) ? fabs(fastSkipDist) : defaultStep_*20;
465  if (resultToMat != SteppingHelixStateInfo::INACC
466  && resultToMat != SteppingHelixStateInfo::OK) tanDistMatNextCheck = -1;
467  } else {
468  // resultToMat = SteppingHelixStateInfo::OK;
469  tanDistMat = tanDistMat > 0. ? tanDistMat - oldDStep : tanDistMat + oldDStep;
470  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Skipped refToMat: guess tanDistMat = "<<tanDistMat
471  <<" next check at "<<tanDistMatNextCheck;
472  }
473  }
474 
475  double rDotP = svCurrent->r3.dot(svCurrent->p3);
476  if ((fabs(curZ) > 1.5e3 || curR >800.)
477  && ((dir == alongMomentum && rDotP > 0)
478  || (dir == oppositeToMomentum && rDotP < 0) )
479  ){
480  dStep = fabs(tanDist) -1e-12;
481  }
482  double tanDistMin = fabs(tanDist);
483  if (tanDistMin > fabs(tanDistMag)+0.05 &&
484  (resultToMag == SteppingHelixStateInfo::OK || resultToMag == SteppingHelixStateInfo::WRONG_VOLUME)){
485  tanDistMin = fabs(tanDistMag)+0.05; //try to step into the next volume
486  expectNewMagVolume = true;
487  } else expectNewMagVolume = false;
488 
489  if (tanDistMin > fabs(tanDistMat)+0.05 && resultToMat == SteppingHelixStateInfo::OK){
490  tanDistMin = fabs(tanDistMat)+0.05; //try to step into the next volume
491  if (expectNewMagVolume) expectNewMagVolume = false;
492  }
493 
494  if (tanDistMin*fabs(svCurrent->dEdx) > 0.5*svCurrent->p3.mag()){
495  tanDistMin = 0.5*svCurrent->p3.mag()/fabs(svCurrent->dEdx);
496  if (expectNewMagVolume) expectNewMagVolume = false;
497  }
498 
499 
500 
501  double tanDistMinLazy = fabs(tanDistMin);
502  if ((type == POINT_PCA_DT || type == LINE_PCA_DT)
503  && fabs(tanDist) < 2.*fabs(tanDistMin) ){
504  //being lazy here; the best is to take into account the curvature
505  tanDistMinLazy = fabs(tanDistMin)*0.5;
506  }
507 
508  if (fabs(tanDistMinLazy) < dStep){
509  dStep = fabs(tanDistMinLazy);
510  }
511 
512  //keep this path length for the next step
513  oldDStep = dStep;
514 
515  if (dStep > 1e-10 && ! (fabs(dist) < fabs(epsilon))){
516  StateInfo* svNext = &svBuf_[cIndex_(nPoints_)];
517  makeAtomStep((*svCurrent), (*svNext), dStep, dir, HEL_AS_F);
518 // if (useMatVolumes_ && expectNewMagVolume
519 // && svCurrent->magVol == svNext->magVol){
520 // double tmpDist=0;
521 // double tmpDistMag = 0;
522 // if (refToMagVolume((*svNext), dir, tmpDist, tmpDistMag, fabs(dist)) != SteppingHelixStateInfo::OK){
523 // //the point appears to be outside, but findVolume claims the opposite
524 // dStep += 0.05*fabs(tanDistMag/distMag); oldDStep = dStep; //do it again with a bigger step
525 // if (debug_) LogTrace(metname)
526 // <<"Failed to get into new mag volume: will try with new bigger step "<<dStep<<std::endl;
527 // makeAtomStep((*svCurrent), (*svNext), dStep, dir, HEL_AS_F);
528 // }
529 // }
530  nPoints_++; svCurrent = &svBuf_[cIndex_(nPoints_-1)];
531  if (oldDir != dir){
532  nOsc++;
533  tanDistNextCheck = -1;//check dist after osc
534  tanDistMagNextCheck = -1;
535  tanDistMatNextCheck = -1;
536  }
537  oldDir = dir;
538  }
539 
540  if (nOsc>1 && fabs(dStep)>epsilon){
541  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Ooops"<<std::endl;
542  }
543 
544  if (fabs(dist) < fabs(epsilon) ) result = SteppingHelixStateInfo::OK;
545 
546  if ((type == POINT_PCA_DT || type == LINE_PCA_DT )
547  && fabs(dStep) < fabs(epsilon) ){
548  //now check if it's not a branch point (peek ahead at 1 cm)
549  double nextDist = 0;
550  double nextTanDist = 0;
551  PropagationDirection nextRefDirection = anyDirection;
552  StateInfo* svNext = &svBuf_[cIndex_(nPoints_)];
553  makeAtomStep((*svCurrent), (*svNext), 1., dir, HEL_AS_F);
554  nPoints_++; svCurrent = &svBuf_[cIndex_(nPoints_-1)];
555  refToDest(type, (*svCurrent), pars, nextDist, nextTanDist, nextRefDirection);
556  if ( fabs(nextDist) > fabs(dist)){
557  nPoints_--; svCurrent = &svBuf_[cIndex_(nPoints_-1)];
559  if (debug_){
560  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Found real local minimum in PCA"<<std::endl;
561  }
562  } else {
563  //keep this trial point and continue
564  dStep = defaultStep_;
565  if (debug_){
566  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Found branch point in PCA"<<std::endl;
567  }
568  }
569  }
570 
571  if (nPoints_ > MAX_STEPS*1./defaultStep_ || loopCount > MAX_STEPS*100
572  || nOsc > 6 ) result = SteppingHelixStateInfo::FAULT;
573 
574  if (svCurrent->p3.mag() < 0.1 ) result = SteppingHelixStateInfo::RANGEOUT;
575 
576  curZ = svCurrent->r3.z();
577  curR = svCurrent->r3.perp();
578  if ( curR > 20000 || fabs(curZ) > 20000 ) result = SteppingHelixStateInfo::INACC;
579 
580  makeNextStep = result == SteppingHelixStateInfo::UNDEFINED;
581  svCurrent->status_ = result;
582  svCurrent->isValid_ = (result == SteppingHelixStateInfo::OK || makeNextStep );
583  svCurrent->field = field_;
584 
585  isFirstStep = false;
586  loopCount++;
587  }
588 
590  edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<" Propagation failed with status "
592  <<std::endl;
593  if (result == SteppingHelixStateInfo::RANGEOUT)
594  edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<" Momentum at last point is too low (<0.1) p_last = "
595  <<svCurrent->p3.mag()
596  <<std::endl;
597  if (result == SteppingHelixStateInfo::INACC)
598  edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<" Went too far: the last point is at "<<svCurrent->r3
599  <<std::endl;
600  if (result == SteppingHelixStateInfo::FAULT && nOsc > 6)
601  edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<" Infinite loop condidtion detected: going in cycles. Break after 6 cycles"
602  <<std::endl;
604  edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<" Tired to go farther. Made too many steps: more than "
606  <<std::endl;
607 
608  }
609 
610  if (debug_){
611  switch (type) {
612  case RADIUS_DT:
613  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to radius "<<pars[RADIUS_P]<<std::endl;
614  break;
615  case Z_DT:
616  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to z "<<pars[Z_P]<<std::endl;
617  break;
618  case PATHL_DT:
619  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to pathL "<<pars[PATHL_P]<<std::endl;
620  break;
621  case PLANE_DT:
622  {
623  Point rPlane(pars[0], pars[1], pars[2]);
624  Vector nPlane(pars[3], pars[4], pars[5]);
625  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to plane r0:"<<rPlane<<" n:"<<nPlane<<std::endl;
626  }
627  break;
628  case POINT_PCA_DT:
629  {
630  Point rDest(pars[0], pars[1], pars[2]);
631  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to PCA to point "<<rDest<<std::endl;
632  }
633  break;
634  case LINE_PCA_DT:
635  {
636  Point rDest1(pars[0], pars[1], pars[2]);
637  Point rDest2(pars[3], pars[4], pars[5]);
638  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to PCA to line "<<rDest1<<" - "<<rDest2<<std::endl;
639  }
640  break;
641  default:
642  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to NOT IMPLEMENTED"<<std::endl;
643  break;
644  }
645  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Made "<<nPoints_-1<<" steps and stopped at(cur step) "<<svCurrent->r3<<" nOsc "<<nOsc<<std::endl;
646  }
647 
648  return result;
649 }
650 
653  const SteppingHelixPropagator::Point& r3, int charge,
655  const AlgebraicSymMatrix55& covCurv) const{
656  static const std::string metname = "SteppingHelixPropagator";
657 
658  svCurrent.q = charge;
659  svCurrent.p3 = p3;
660  svCurrent.r3 = r3;
661  svCurrent.dir = dir == alongMomentum ? 1.: -1.;
662 
663  svCurrent.path_ = 0; // this could've held the initial path
664  svCurrent.radPath_ = 0;
665 
666  GlobalPoint gPointNegZ(svCurrent.r3.x(), svCurrent.r3.y(), -fabs(svCurrent.r3.z()));
667  GlobalPoint gPointNorZ(svCurrent.r3.x(), svCurrent.r3.y(), svCurrent.r3.z());
668 
669  float gpmag = gPointNegZ.mag2();
670  if (gpmag > 1e20f ) {
671  LogTrace(metname)<<"Initial point is too far";
672  svCurrent.isValid_ = false;
673  return;
674  }
675  if (! (gpmag == gpmag) ) {
676  LogTrace(metname)<<"Initial point is a nan";
677  svCurrent.isValid_ = false;
678  return;
679  }
680 
681  GlobalVector bf(0,0,0);
682  // = field_->inTesla(gPoint);
683  if (useMagVolumes_){
684  if (vbField_ ){
685  if (vbField_->isZSymmetric()){
686  svCurrent.magVol = vbField_->findVolume(gPointNegZ);
687  } else {
688  svCurrent.magVol = vbField_->findVolume(gPointNorZ);
689  }
690  if (useIsYokeFlag_){
691  double curRad = svCurrent.r3.perp();
692  if (curRad > 380 && curRad < 850 && fabs(svCurrent.r3.z()) < 667){
693  svCurrent.isYokeVol = isYokeVolume(svCurrent.magVol);
694  } else {
695  svCurrent.isYokeVol = false;
696  }
697  }
698  } else {
699  edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Failed to cast into VolumeBasedMagneticField: fall back to the default behavior"<<std::endl;
700  svCurrent.magVol = 0;
701  }
702  if (debug_){
703  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Got volume at "<<svCurrent.magVol<<std::endl;
704  }
705  }
706 
707  if (useMagVolumes_ && svCurrent.magVol != 0 && ! useInTeslaFromMagField_){
708  if (vbField_->isZSymmetric()){
709  bf = svCurrent.magVol->inTesla(gPointNegZ);
710  if (debug_){
711  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Loaded bfield float: "<<bf
712  <<" at global float "<< gPointNegZ<<" double "<< svCurrent.r3<<std::endl;
713  LocalPoint lPoint(svCurrent.magVol->toLocal(gPointNegZ));
714  LocalVector lbf = svCurrent.magVol->fieldInTesla(lPoint);
715  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"\t cf in local locF: "<<lbf<<" at "<<lPoint<<std::endl;
716  }
717  } else {
718  bf = svCurrent.magVol->inTesla(gPointNorZ);
719  if (debug_){
720  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Loaded bfield float: "<<bf
721  <<" at global float "<< gPointNorZ<<" double "<< svCurrent.r3<<std::endl;
722  LocalPoint lPoint(svCurrent.magVol->toLocal(gPointNorZ));
723  LocalVector lbf = svCurrent.magVol->fieldInTesla(lPoint);
724  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"\t cf in local locF: "<<lbf<<" at "<<lPoint<<std::endl;
725  }
726  }
727  if (r3.z() > 0 && vbField_->isZSymmetric() ){
728  svCurrent.bf.set(-bf.x(), -bf.y(), bf.z());
729  } else {
730  svCurrent.bf.set(bf.x(), bf.y(), bf.z());
731  }
732  } else {
733  GlobalPoint gPoint(r3.x(), r3.y(), r3.z());
734  bf = field_->inTesla(gPoint);
735  svCurrent.bf.set(bf.x(), bf.y(), bf.z());
736  }
737  if (svCurrent.bf.mag2() < 1e-32) svCurrent.bf.set(0., 0., 1e-16);
738  if (debug_){
739  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific
740  <<"Loaded bfield double: "<<svCurrent.bf<<" from float: "<<bf
741  <<" at float "<< gPointNorZ<<" "<<gPointNegZ<<" double "<< svCurrent.r3<<std::endl;
742  }
743 
744 
745 
746  double dEdXPrime = 0;
747  double dEdx = 0;
748  double radX0 = 0;
749  MatBounds rzLims;
750  dEdx = getDeDx(svCurrent, dEdXPrime, radX0, rzLims);
751  svCurrent.dEdx = dEdx; svCurrent.dEdXPrime = dEdXPrime;
752  svCurrent.radX0 = radX0;
753  svCurrent.rzLims = rzLims;
754 
755  svCurrent.covCurv =covCurv;
756 
757  svCurrent.isComplete = true;
758  svCurrent.isValid_ = true;
759 
760  if (debug_){
761  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Loaded at path: "<<svCurrent.path()<<" radPath: "<<svCurrent.radPath()
762  <<" p3 "<<" pt: "<<svCurrent.p3.perp()<<" phi: "<<svCurrent.p3.phi()
763  <<" eta: "<<svCurrent.p3.eta()
764  <<" "<<svCurrent.p3
765  <<" r3: "<<svCurrent.r3
766  <<" bField: "<<svCurrent.bf.mag()
767  <<std::endl;
768  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Input Covariance in Curvilinear RF "<<covCurv<<std::endl;
769  }
770 }
771 
774  double dP, const SteppingHelixPropagator::Vector& tau,
775  const SteppingHelixPropagator::Vector& drVec, double dS, double dX0,
776  const AlgebraicMatrix55& dCovCurvTransform) const{
777  static const std::string metname = "SteppingHelixPropagator";
778  svNext.q = svPrevious.q;
779  svNext.dir = dS > 0.0 ? 1.: -1.;
780  svNext.p3 = tau; svNext.p3*=(svPrevious.p3.mag() - svNext.dir*fabs(dP));
781 
782  svNext.r3 = svPrevious.r3; svNext.r3 += drVec;
783 
784  svNext.path_ = svPrevious.path() + dS;
785  svNext.radPath_ = svPrevious.radPath() + dX0;
786 
787  GlobalPoint gPointNegZ(svNext.r3.x(), svNext.r3.y(), -fabs(svNext.r3.z()));
788  GlobalPoint gPointNorZ(svNext.r3.x(), svNext.r3.y(), svNext.r3.z());
789 
790  GlobalVector bf(0,0,0);
791 
792  if (useMagVolumes_){
793  if (vbField_ != 0){
794  if (vbField_->isZSymmetric()){
795  svNext.magVol = vbField_->findVolume(gPointNegZ);
796  } else {
797  svNext.magVol = vbField_->findVolume(gPointNorZ);
798  }
799  if (useIsYokeFlag_){
800  double curRad = svNext.r3.perp();
801  if (curRad > 380 && curRad < 850 && fabs(svNext.r3.z()) < 667){
802  svNext.isYokeVol = isYokeVolume(svNext.magVol);
803  } else {
804  svNext.isYokeVol = false;
805  }
806  }
807  } else {
808  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Failed to cast into VolumeBasedMagneticField"<<std::endl;
809  svNext.magVol = 0;
810  }
811  if (debug_){
812  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Got volume at "<<svNext.magVol<<std::endl;
813  }
814  }
815 
816  if (useMagVolumes_ && svNext.magVol != 0 && ! useInTeslaFromMagField_){
817  if (vbField_->isZSymmetric()){
818  bf = svNext.magVol->inTesla(gPointNegZ);
819  } else {
820  bf = svNext.magVol->inTesla(gPointNorZ);
821  }
822  if (svNext.r3.z() > 0 && vbField_->isZSymmetric() ){
823  svNext.bf.set(-bf.x(), -bf.y(), bf.z());
824  } else {
825  svNext.bf.set(bf.x(), bf.y(), bf.z());
826  }
827  } else {
828  GlobalPoint gPoint(svNext.r3.x(), svNext.r3.y(), svNext.r3.z());
829  bf = field_->inTesla(gPoint);
830  svNext.bf.set(bf.x(), bf.y(), bf.z());
831  }
832  if (svNext.bf.mag2() < 1e-32) svNext.bf.set(0., 0., 1e-16);
833 
834 
835  double dEdXPrime = 0;
836  double dEdx = 0;
837  double radX0 = 0;
838  MatBounds rzLims;
839  dEdx = getDeDx(svNext, dEdXPrime, radX0, rzLims);
840  svNext.dEdx = dEdx; svNext.dEdXPrime = dEdXPrime;
841  svNext.radX0 = radX0;
842  svNext.rzLims = rzLims;
843 
844  //update Emat only if it's valid
845  svNext.hasErrorPropagated_ = svPrevious.hasErrorPropagated_;
846  if (svPrevious.hasErrorPropagated_){
847  {
848  AlgebraicMatrix55 tmp = dCovCurvTransform*svPrevious.covCurv;
849  ROOT::Math::AssignSym::Evaluate(svNext.covCurv, tmp*ROOT::Math::Transpose(dCovCurvTransform));
850 
851  svNext.covCurv += svPrevious.matDCovCurv;
852  }
853  } else {
854  //could skip dragging along the unprop. cov later .. now
855  // svNext.cov = svPrevious.cov;
856  }
857 
858  if (debug_){
859  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Now at path: "<<svNext.path()<<" radPath: "<<svNext.radPath()
860  <<" p3 "<<" pt: "<<svNext.p3.perp()<<" phi: "<<svNext.p3.phi()
861  <<" eta: "<<svNext.p3.eta()
862  <<" "<<svNext.p3
863  <<" r3: "<<svNext.r3
864  <<" dPhi: "<<acos(svNext.p3.unit().dot(svPrevious.p3.unit()))
865  <<" bField: "<<svNext.bf.mag()
866  <<" dedx: "<<svNext.dEdx
867  <<std::endl;
868  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"New Covariance "<<svNext.covCurv<<std::endl;
869  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Transf by dCovTransform "<<dCovCurvTransform<<std::endl;
870  }
871 }
872 
874  const SteppingHelixPropagator::Vector& tau) const{
875  Vector zRep(0., 0., 1.);
876  rep.lX = tau;
877  rep.lY = zRep.cross(tau); rep.lY *= 1./tau.perp();
878  rep.lZ = rep.lX.cross(rep.lY);
879 }
880 
883  double dS,
885  SteppingHelixPropagator::Fancy fancy) const{
886  static const std::string metname = "SteppingHelixPropagator";
887  if (debug_){
888  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Make atom step "<<svCurrent.path()<<" with step "<<dS<<" in direction "<<dir<<std::endl;
889  }
890 
891  double dP = 0;
892  double curP = svCurrent.p3.mag();
893  Vector tau = svCurrent.p3; tau *= 1./curP;
894  Vector tauNext(tau);
895  Vector drVec(0,0,0);
896 
897  dS = dir == alongMomentum ? fabs(dS) : -fabs(dS);
898 
899 
900  double radX0 = 1e24;
901 
902  switch (fancy){
903  case HEL_AS_F:
904  case HEL_ALL_F:{
905  double p0 = curP;// see above = svCurrent.p3.mag();
906  double p0Inv = 1./p0;
907  double b0 = svCurrent.bf.mag();
908 
909  //get to the mid-point first
910  double phi = (2.99792458e-3*svCurrent.q*b0*p0Inv)*dS/2.;
911  bool phiSmall = fabs(phi) < 1e-4;
912 
913  double cosPhi = 0;
914  double sinPhi = 0;
915 
916  double oneLessCosPhi=0;
917  double oneLessCosPhiOPhi=0;
918  double sinPhiOPhi=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  sinPhi = phi - phi3/6. + phi4*phi/120.;
926  cosPhi = 1. -phi2/2. + phi4/24.;
927  oneLessCosPhi = phi2/2. - phi4/24. + phi2*phi4/720.; // 0.5*phi*phi;//*(1.- phi*phi/12.);
928  oneLessCosPhiOPhi = 0.5*phi - phi3/24. + phi2*phi3/720.;//*(1.- phi*phi/12.);
929  sinPhiOPhi = 1. - phi*phi/6. + phi4/120.;
930  phiLessSinPhiOPhi = phi*phi/6. - phi4/120. + phi4*phi2/5040.;//*(1. - phi*phi/20.);
931  } else {
932  cosPhi = cos(phi);
933  sinPhi = sin(phi);
934  oneLessCosPhi = 1.-cosPhi;
935  oneLessCosPhiOPhi = oneLessCosPhi/phi;
936  sinPhiOPhi = sinPhi/phi;
937  phiLessSinPhiOPhi = 1 - sinPhiOPhi;
938  }
939 
940  Vector bHat = svCurrent.bf; bHat *= 1./b0; //bHat.mag();
941  // bool bAlongZ = fabs(bHat.z()) > 0.9999;
942 
943  Vector btVec(bHat.cross(tau)); // for balong z btVec.z()==0
944  double tauB = tau.dot(bHat);
945  Vector bbtVec(bHat*tauB - tau); // (-tau.x(), -tau.y(), 0)
946 
947  //don't need it here tauNext = tau + bbtVec*oneLessCosPhi - btVec*sinPhi;
948  drVec = bbtVec*phiLessSinPhiOPhi; drVec -= btVec*oneLessCosPhiOPhi; 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); //(bfGV.x(), bfGV.y(), bfGV.z());
960  // = svCurrent.magVol->inTesla(GlobalPoint(drVec.x(), drVec.y(), -fabs(drVec.z())));
961  if (useMagVolumes_ && svCurrent.magVol != 0 && ! useInTeslaFromMagField_){
962  // this negative-z business will break at some point
963  if (vbField_->isZSymmetric()){
964  bfGV = svCurrent.magVol->inTesla(GlobalPoint(drVec.x(), drVec.y(), -fabs(drVec.z())));
965  } else {
966  bfGV = svCurrent.magVol->inTesla(GlobalPoint(drVec.x(), drVec.y(), drVec.z()));
967  }
968  if (drVec.z() > 0 && vbField_->isZSymmetric()){
969  bf.set(-bfGV.x(), -bfGV.y(), bfGV.z());
970  } else {
971  bf.set(bfGV.x(), bfGV.y(), bfGV.z());
972  }
973  } else {
974  bfGV = field_->inTesla(GlobalPoint(drVec.x(), drVec.y(), drVec.z()));
975  bf.set(bfGV.x(), bfGV.y(), bfGV.z());
976  }
977  double b0Init = b0;
978  b0 = bf.mag();
979  if (b0 < 1e-16) {
980  b0 = 1e-16;
981  bf.set(0., 0., 1e-16);
982  }
983  if (debug_){
984  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Improved b "<<b0
985  <<" at r3 "<<drVec<<std::endl;
986  }
987 
988  if (fabs((b0-b0Init)*dS) > 1){
989  //missed the mag volume boundary?
990  if (debug_){
991  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Large bf*dS change "<<fabs((b0-svCurrent.bf.mag())*dS)
992  <<" --> recalc dedx"<<std::endl;
993  }
994  svNext.r3 = drVec;
995  svNext.bf = bf;
996  svNext.p3 = svCurrent.p3;
997  svNext.isYokeVol = svCurrent.isYokeVol;
998  MatBounds rzTmp;
999  dEdx = getDeDx(svNext, dEdXPrime, radX0, rzTmp);
1000  dP = dEdx*dS;
1001  if (debug_){
1002  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"New dEdX= "<<dEdx
1003  <<" dP= "<<dP
1004  <<" for p0 "<<p0<<std::endl;
1005  }
1006  }
1007  //p0 is mid-way and b0 from mid-point
1008  p0 += dP/2.; if (p0 < 1e-2) p0 = 1e-2;
1009  p0Inv = 1./p0;
1010 
1011  phi = (2.99792458e-3*svCurrent.q*b0*p0Inv)*dS;
1012  phiSmall = fabs(phi) < 1e-4;
1013 
1014  if (phiSmall){
1015  double phi2 = phi*phi;
1016  double phi3 = phi2*phi;
1017  double phi4 = phi3*phi;
1018  sinPhi = phi - phi3/6. + phi4*phi/120.;
1019  cosPhi = 1. -phi2/2. + phi4/24.;
1020  oneLessCosPhi = phi2/2. - phi4/24. + phi2*phi4/720.; // 0.5*phi*phi;//*(1.- phi*phi/12.);
1021  oneLessCosPhiOPhi = 0.5*phi - phi3/24. + phi2*phi3/720.;//*(1.- phi*phi/12.);
1022  sinPhiOPhi = 1. - phi*phi/6. + phi4/120.;
1023  phiLessSinPhiOPhi = phi*phi/6. - phi4/120. + phi4*phi2/5040.;//*(1. - phi*phi/20.);
1024  }else {
1025  cosPhi = cos(phi);
1026  sinPhi = sin(phi);
1027  oneLessCosPhi = 1.-cosPhi;
1028  oneLessCosPhiOPhi = oneLessCosPhi/phi;
1029  sinPhiOPhi = sinPhi/phi;
1030  phiLessSinPhiOPhi = 1. - sinPhiOPhi;
1031  }
1032 
1033  bHat = bf; bHat *= 1./b0;// as above =1./bHat.mag();
1034  // bAlongZ = fabs(bHat.z()) > 0.9999;
1035  btVec = bHat.cross(tau); // for b||z (-tau.y(), tau.x() ,0)
1036  tauB = tau.dot(bHat);
1037  bbtVec = bHat*tauB - tau; //bHat.cross(btVec); for b||z: (-tau.x(), -tau.y(), 0)
1038 
1039  tauNext = bbtVec*oneLessCosPhi; tauNext -= btVec*sinPhi; tauNext += tau; //for b||z tauNext.z() == tau.z()
1040  double tauNextPerpInv = 1./tauNext.perp();
1041  drVec = bbtVec*phiLessSinPhiOPhi; drVec -= btVec*oneLessCosPhiOPhi; drVec += tau;
1042  drVec *= dS;
1043 
1044 
1045  if (svCurrent.hasErrorPropagated_){
1046  double theta02 = 0;
1047  double dX0 = fabs(dS)/radX0;
1048 
1049  if (applyRadX0Correction_){
1050  // this provides the integrand for theta^2
1051  // if summed up along the path, should result in
1052  // theta_total^2 = Int_0^x0{ f(x)dX} = (13.6/p0)^2*x0*(1+0.036*ln(x0+1))
1053  // x0+1 above is to make the result infrared safe.
1054  double x0 = fabs(svCurrent.radPath());
1055  double alphaX0 = 13.6e-3*p0Inv; alphaX0 *= alphaX0;
1056  double betaX0 = 0.038;
1057  double logx0p1 = log(x0+1);
1058  theta02 = dX0*alphaX0*(1+betaX0*logx0p1)*(1 + betaX0*logx0p1 + 2.*betaX0*x0/(x0+1) );
1059  } else {
1060  theta02 = 196e-6* p0Inv * p0Inv * dX0; //14.e-3/p0*sqrt(fabs(dS)/radX0); // .. drop log term (this is non-additive)
1061  }
1062 
1063  double epsilonP0 = 1.+ dP/(p0-0.5*dP);
1064  // double omegaP0 = -dP/(p0-0.5*dP) + dS*dEdXPrime;
1065  // double dsp = dS/(p0-0.5*dP); //use the initial p0 (not the mid-point) to keep the transport properly additive
1066 
1067  Vector tbtVec(bHat - tauB*tau); // for b||z tau.z()*(-tau.x(), -tau.y(), 1.-tau.z())
1068 
1070  {
1071  //Slightly modified copy of the curvilinear jacobian (don't use the original just because it's in float precision
1072  // and seems to have some assumptions about the field values
1073  // notation changes: p1--> tau, p2-->tauNext
1074  // theta --> phi
1075  // Vector p1 = tau;
1076  // Vector p2 = tauNext;
1077  Point xStart = svCurrent.r3;
1078  Vector dx = drVec;
1079  //GlobalVector h = MagneticField::inInverseGeV(xStart);
1080  // Martijn: field is now given as parameter.. GlobalVector h = globalParameters.magneticFieldInInverseGeV(xStart);
1081 
1082  //double qbp = fts.signedInverseMomentum();
1083  double qbp = svCurrent.q*p0Inv;
1084  // double absS = dS;
1085 
1086  // calculate transport matrix
1087  // Origin: TRPRFN
1088  double t11 = tau.x(); double t12 = tau.y(); double t13 = tau.z();
1089  double t21 = tauNext.x(); double t22 = tauNext.y(); 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  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; double cost = cosPhi;
1103  double hn1 = hn.x(); double hn2 = hn.y(); double hn3 = hn.z();
1104  double dx1 = dx.x(); double dx2 = dx.y(); double dx3 = dx.z();
1105  // double hnDt1 = hn1*t11 + hn2*t12 + hn3*t13;
1106 
1107  double gamma = hn1*t21 + hn2*t22 + hn3*t23;
1108  double an1 = hn2*t23 - hn3*t22;
1109  double an2 = hn3*t21 - hn1*t23;
1110  double an3 = hn1*t22 - hn2*t21;
1111  // double auInv = sqrt(1.- t13*t13); double au = auInv>0 ? 1./auInv : 1e24;
1112  double auInv = cosl0; double au = auInv>0 ? 1./auInv : 1e24;
1113  // double auInv = sqrt(t11*t11 + t12*t12); double au = auInv>0 ? 1./auInv : 1e24;
1114  double u11 = -au*t12; double u12 = au*t11;
1115  double v11 = -t13*u12; double v12 = t13*u11; double v13 = auInv;//t11*u12 - t12*u11;
1116  auInv = sqrt(1. - t23*t23); au = auInv>0 ? 1./auInv : 1e24;
1117  // auInv = sqrt(t21*t21 + t22*t22); au = auInv>0 ? 1./auInv : 1e24;
1118  double u21 = -au*t22; double u22 = au*t21;
1119  double v21 = -t23*u22; double v22 = t23*u21; double v23 = auInv;//t21*u22 - t22*u21;
1120  // now prepare the transport matrix
1121  // pp only needed in high-p case (WA)
1122  // double pp = 1./qbp;
1124  // moved up (where -h.mag() is needed()
1125  // double qp = q*pp;
1126  double anv = -(hn1*u21 + hn2*u22 );
1127  double anu = (hn1*v21 + hn2*v22 + hn3*v23);
1128  double omcost = oneLessCosPhi; double tmsint = -phi*phiLessSinPhiOPhi;
1129 
1130  double hu1 = - hn3*u12;
1131  double hu2 = hn3*u11;
1132  double hu3 = hn1*u12 - hn2*u11;
1133 
1134  double hv1 = hn2*v13 - hn3*v12;
1135  double hv2 = hn3*v11 - hn1*v13;
1136  double hv3 = hn1*v12 - hn2*v11;
1137 
1138  // 1/p - doesn't change since |tau| = |tauNext| ... not. It changes now
1139  dCCurvTransform_(0,0) = 1./(epsilonP0*epsilonP0)*(1. + dS*dEdXPrime);
1140 
1141  // lambda
1142 
1143  dCCurvTransform_(1,0) = phi*p0/svCurrent.q*cosecl1*
1144  (sinPhi*bbtVec.z() - cosPhi*btVec.z());
1145  //was dCCurvTransform_(1,0) = -qp*anv*(t21*dx1 + t22*dx2 + t23*dx3); //NOTE (SK) this was found to have an opposite sign
1146  //from independent re-calculation ... in fact the tauNext.dot.dR piece isnt reproduced
1147 
1148  dCCurvTransform_(1,1) = cost*(v11*v21 + v12*v22 + v13*v23) +
1149  sint*(hv1*v21 + hv2*v22 + hv3*v23) +
1150  omcost*(hn1*v11 + hn2*v12 + hn3*v13) * (hn1*v21 + hn2*v22 + hn3*v23) +
1151  anv*(-sint*(v11*t21 + v12*t22 + v13*t23) +
1152  omcost*(v11*an1 + v12*an2 + v13*an3) -
1153  tmsint*gamma*(hn1*v11 + hn2*v12 + hn3*v13) );
1154 
1155  dCCurvTransform_(1,2) = cost*(u11*v21 + u12*v22 ) +
1156  sint*(hu1*v21 + hu2*v22 + hu3*v23) +
1157  omcost*(hn1*u11 + hn2*u12 ) * (hn1*v21 + hn2*v22 + hn3*v23) +
1158  anv*(-sint*(u11*t21 + u12*t22 ) +
1159  omcost*(u11*an1 + u12*an2 ) -
1160  tmsint*gamma*(hn1*u11 + hn2*u12 ) );
1161  dCCurvTransform_(1,2) *= cosl0;
1162 
1163  // Commented out in part for reproducibility purposes: these terms are zero in cart->curv
1164  // dCCurvTransform_(1,3) = -q*anv*(u11*t21 + u12*t22 ); //don't show up in cartesian setup-->curv
1165  //why would lambdaNext depend explicitely on initial position ? any arbitrary init point can be chosen not
1166  // affecting the final state's momentum direction ... is this the field gradient in curvilinear coord?
1167  // dCCurvTransform_(1,4) = -q*anv*(v11*t21 + v12*t22 + v13*t23); //don't show up in cartesian setup-->curv
1168 
1169  // phi
1170 
1171  dCCurvTransform_(2,0) = - phi*p0/svCurrent.q*cosecl1*cosecl1*
1172  (oneLessCosPhi*bHat.z()*btVec.mag2() + sinPhi*btVec.z() + cosPhi*tbtVec.z()) ;
1173  //was dCCurvTransform_(2,0) = -qp*anu*(t21*dx1 + t22*dx2 + t23*dx3)*cosecl1;
1174 
1175  dCCurvTransform_(2,1) = cost*(v11*u21 + v12*u22 ) +
1176  sint*(hv1*u21 + hv2*u22 ) +
1177  omcost*(hn1*v11 + hn2*v12 + hn3*v13) *
1178  (hn1*u21 + hn2*u22 ) +
1179  anu*(-sint*(v11*t21 + v12*t22 + v13*t23) +
1180  omcost*(v11*an1 + v12*an2 + v13*an3) -
1181  tmsint*gamma*(hn1*v11 + hn2*v12 + hn3*v13) );
1182  dCCurvTransform_(2,1) *= cosecl1;
1183 
1184  dCCurvTransform_(2,2) = cost*(u11*u21 + u12*u22 ) +
1185  sint*(hu1*u21 + hu2*u22 ) +
1186  omcost*(hn1*u11 + hn2*u12 ) *
1187  (hn1*u21 + hn2*u22 ) +
1188  anu*(-sint*(u11*t21 + u12*t22 ) +
1189  omcost*(u11*an1 + u12*an2 ) -
1190  tmsint*gamma*(hn1*u11 + hn2*u12 ) );
1191  dCCurvTransform_(2,2) *= cosecl1*cosl0;
1192 
1193  // Commented out in part for reproducibility purposes: these terms are zero in cart->curv
1194  // dCCurvTransform_(2,3) = -q*anu*(u11*t21 + u12*t22 )*cosecl1;
1195  //why would lambdaNext depend explicitely on initial position ? any arbitrary init point can be chosen not
1196  // affecting the final state's momentum direction ... is this the field gradient in curvilinear coord?
1197  // dCCurvTransform_(2,4) = -q*anu*(v11*t21 + v12*t22 + v13*t23)*cosecl1;
1198 
1199  // yt
1200 
1201  double pp = 1./qbp;
1202  // (SK) these terms seem to consistently have a sign opp from private derivation
1203  dCCurvTransform_(3,0) = - pp*(u21*dx1 + u22*dx2 ); //NB: modified from the original: changed the sign
1204  dCCurvTransform_(4,0) = - pp*(v21*dx1 + v22*dx2 + v23*dx3);
1205 
1206 
1207  dCCurvTransform_(3,1) = (sint*(v11*u21 + v12*u22 ) +
1208  omcost*(hv1*u21 + hv2*u22 ) +
1209  tmsint*(hn1*u21 + hn2*u22 ) *
1210  (hn1*v11 + hn2*v12 + hn3*v13))/q;
1211 
1212  dCCurvTransform_(3,2) = (sint*(u11*u21 + u12*u22 ) +
1213  omcost*(hu1*u21 + hu2*u22 ) +
1214  tmsint*(hn1*u21 + hn2*u22 ) *
1215  (hn1*u11 + hn2*u12 ))*cosl0/q;
1216 
1217  dCCurvTransform_(3,3) = (u11*u21 + u12*u22 );
1218 
1219  dCCurvTransform_(3,4) = (v11*u21 + v12*u22 );
1220 
1221  // zt
1222 
1223  dCCurvTransform_(4,1) = (sint*(v11*v21 + v12*v22 + v13*v23) +
1224  omcost*(hv1*v21 + hv2*v22 + hv3*v23) +
1225  tmsint*(hn1*v21 + hn2*v22 + hn3*v23) *
1226  (hn1*v11 + hn2*v12 + hn3*v13))/q;
1227 
1228  dCCurvTransform_(4,2) = (sint*(u11*v21 + u12*v22 ) +
1229  omcost*(hu1*v21 + hu2*v22 + hu3*v23) +
1230  tmsint*(hn1*v21 + hn2*v22 + hn3*v23) *
1231  (hn1*u11 + hn2*u12 ))*cosl0/q;
1232 
1233  dCCurvTransform_(4,3) = (u11*v21 + u12*v22 );
1234 
1235  dCCurvTransform_(4,4) = (v11*v21 + v12*v22 + v13*v23);
1236  // end of TRPRFN
1237  }
1238 
1239  if (debug_){
1240  Basis rep; setRep(rep, tauNext);
1241  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"rep X: "<<rep.lX<<" "<<rep.lX.mag()
1242  <<"\t Y: "<<rep.lY<<" "<<rep.lY.mag()
1243  <<"\t Z: "<<rep.lZ<<" "<<rep.lZ.mag();
1244  }
1245  //mind the sign of dS and dP (dS*dP < 0 allways)
1246  //covariance should grow no matter which direction you propagate
1247  //==> take abs values.
1248  //reset not needed: fill all below svCurrent.matDCov *= 0.;
1249  double mulRR = theta02*dS*dS/3.;
1250  double mulRP = theta02*fabs(dS)*p0/2.;
1251  double mulPP = theta02*p0*p0;
1252  double losPP = dP*dP*1.6/fabs(dS)*(1.0 + p0*1e-3);
1253  //another guess .. makes sense for 1 cm steps 2./dS == 2 [cm] / dS [cm] at low pt
1254  //double it by 1TeV
1255  //not gaussian anyways
1256  // derived from the fact that sigma_p/eLoss ~ 0.08 after ~ 200 steps
1257 
1258  //curvilinear
1259  double sinThetaInv = tauNextPerpInv;
1260  double p0Mat = p0+ 0.5*dP; // FIXME change this to p0 after it's clear that there's agreement in everything else
1261  double p0Mat2 = p0Mat*p0Mat;
1262  // with 6x6 formulation
1263  svCurrent.matDCovCurv*=0;
1264 
1265  svCurrent.matDCovCurv(0,0) = losPP/(p0Mat2*p0Mat2);
1266  // svCurrent.matDCovCurv(0,1) = 0;
1267  // svCurrent.matDCovCurv(0,2) = 0;
1268  // svCurrent.matDCovCurv(0,3) = 0;
1269  // svCurrent.matDCovCurv(0,4) = 0;
1270 
1271  // svCurrent.matDCovCurv(1,0) = 0;
1272  svCurrent.matDCovCurv(1,1) = mulPP/p0Mat2;
1273  // svCurrent.matDCovCurv(1,2) = 0;
1274  // svCurrent.matDCovCurv(1,3) = 0;
1275  svCurrent.matDCovCurv(1,4) = mulRP/p0Mat;
1276 
1277  // svCurrent.matDCovCurv(2,0) = 0;
1278  // svCurrent.matDCovCurv(2,1) = 0;
1279  svCurrent.matDCovCurv(2,2) = mulPP/p0Mat2*(sinThetaInv*sinThetaInv);
1280  svCurrent.matDCovCurv(2,3) = mulRP/p0Mat*sinThetaInv;
1281  // svCurrent.matDCovCurv(2,4) = 0;
1282 
1283  // svCurrent.matDCovCurv(3,0) = 0;
1284  // svCurrent.matDCovCurv(3,1) = 0;
1285  svCurrent.matDCovCurv(3,2) = mulRP/p0Mat*sinThetaInv;
1286  // svCurrent.matDCovCurv(3,0) = 0;
1287  svCurrent.matDCovCurv(3,3) = mulRR;
1288  // svCurrent.matDCovCurv(3,4) = 0;
1289 
1290  // svCurrent.matDCovCurv(4,0) = 0;
1291  svCurrent.matDCovCurv(4,1) = mulRP/p0Mat;
1292  // svCurrent.matDCovCurv(4,2) = 0;
1293  // svCurrent.matDCovCurv(4,3) = 0;
1294  svCurrent.matDCovCurv(4,4) = mulRR;
1295  }
1296  break;
1297  }
1298  default:
1299  break;
1300  }
1301 
1302  double pMag = curP;//svCurrent.p3.mag();
1303 
1304  if (dir == oppositeToMomentum) dP = fabs(dP);
1305  else if( dP < 0) { //the case of negative dP
1306  dP = -dP > pMag ? -pMag+1e-5 : dP;
1307  }
1308 
1309  getNextState(svCurrent, svNext, dP, tauNext, drVec, dS, dS/radX0,
1311  return true;
1312 }
1313 
1315  double& dEdXPrime, double& radX0,
1316  MatBounds& rzLims) const{
1317  radX0 = 1.e24;
1318  dEdXPrime = 0.;
1319  rzLims = MatBounds();
1320  if (noMaterialMode_) return 0;
1321 
1322  double dEdx = 0.;
1323 
1324  double lR = sv.r3.perp();
1325  double lZ = fabs(sv.r3.z());
1326 
1327  //assume "Iron" .. seems to be quite the same for brass/iron/PbW04
1328  //good for Fe within 3% for 0.2 GeV to 10PeV
1329 
1330  double dEdX_HCal = 0.95; //extracted from sim
1331  double dEdX_ECal = 0.45;
1332  double dEdX_coil = 0.35; //extracted from sim .. closer to 40% in fact
1333  double dEdX_Fe = 1;
1334  double dEdX_MCh = 0.053; //chambers on average
1335  double dEdX_Trk = 0.0114;
1336  double dEdX_Air = 2E-4;
1337  double dEdX_Vac = 0.0;
1338 
1339  double radX0_HCal = 1.44/0.8; //guessing
1340  double radX0_ECal = 0.89/0.7;
1341  double radX0_coil = 4.; //
1342  double radX0_Fe = 1.76;
1343  double radX0_MCh = 1e3; //
1344  double radX0_Trk = 320.;
1345  double radX0_Air = 3.e4;
1346  double radX0_Vac = 3.e9; //"big" number for vacuum
1347 
1348 
1349  //not all the boundaries are set below: this will be a default
1350  if (! (lR < 380 && lZ < 785)){
1351  if (lZ > 785 ) rzLims = MatBounds(0, 1e4, 785, 1e4);
1352  if (lZ < 785 ) rzLims = MatBounds(380, 1e4, 0, 785);
1353  }
1354 
1355  //this def makes sense assuming we don't jump into endcap volume from the other z-side in one step
1356  //also, it is a positive shift only (at least for now): can't move ec further into the detector
1357  double ecShift = sv.r3.z() > 0 ? fabs(ecShiftPos_) : fabs(ecShiftNeg_);
1358 
1359  //this should roughly figure out where things are
1360  //(numbers taken from Fig1.1.2 TDR and from geom xmls)
1361  if (lR < 2.9){ //inside beampipe
1362  dEdx = dEdX_Vac; radX0 = radX0_Vac;
1363  rzLims = MatBounds(0, 2.9, 0, 1E4);
1364  }
1365  else if (lR < 129){
1366  if (lZ < 294){
1367  rzLims = MatBounds(2.9, 129, 0, 294); //tracker boundaries
1368  dEdx = dEdX_Trk; radX0 = radX0_Trk;
1369  //somewhat empirical formula that ~ matches the average if going from 0,0,0
1370  //assuming "uniform" tracker material
1371  //doesn't really track material layer to layer
1372  double lEtaDet = fabs(sv.r3.eta());
1373  double scaleRadX = lEtaDet > 1.5 ? 0.7724 : 1./cosh(0.5*lEtaDet);//sin(2.*atan(exp(-0.5*lEtaDet)));
1374  scaleRadX *= scaleRadX;
1375  if (lEtaDet > 2 && lZ > 20) scaleRadX *= (lEtaDet-1.);
1376  if (lEtaDet > 2.5 && lZ > 20) scaleRadX *= (lEtaDet-1.);
1377  radX0 *= scaleRadX;
1378  }
1379  //endcap part begins here
1380  else if ( lZ < 294 + ecShift ){
1381  //gap in front of EE here, piece inside 2.9<R<129
1382  rzLims = MatBounds(2.9, 129, 294, 294 + ecShift);
1383  dEdx = dEdX_Air; radX0 = radX0_Air;
1384  }
1385  else if (lZ < 372 + ecShift){
1386  rzLims = MatBounds(2.9, 129, 294 + ecShift, 372 + ecShift); //EE here, piece inside 2.9<R<129
1387  dEdx = dEdX_ECal; radX0 = radX0_ECal;
1388  }//EE averaged out over a larger space
1389  else if (lZ < 398 + ecShift){
1390  rzLims = MatBounds(2.9, 129, 372 + ecShift, 398 + ecShift); //whatever goes behind EE 2.9<R<129 is air up to Z=398
1391  dEdx = dEdX_HCal*0.05; radX0 = radX0_Air;
1392  }//betw EE and HE
1393  else if (lZ < 555 + ecShift){
1394  rzLims = MatBounds(2.9, 129, 398 + ecShift, 555 + ecShift); //HE piece 2.9<R<129;
1395  dEdx = dEdX_HCal*0.96; radX0 = radX0_HCal/0.96;
1396  } //HE calor abit less dense
1397  else {
1398  // rzLims = MatBounds(2.9, 129, 555, 785);
1399  // set the boundaries first: they serve as stop-points here
1400  // the material is set below
1401  if (lZ < 568 + ecShift) rzLims = MatBounds(2.9, 129, 555 + ecShift, 568 + ecShift); //a piece of HE support R<129, 555<Z<568
1402  else if (lZ < 625 + ecShift){
1403  if (lR < 85 + ecShift) rzLims = MatBounds(2.9, 85, 568 + ecShift, 625 + ecShift);
1404  else rzLims = MatBounds(85, 129, 568 + ecShift, 625 + ecShift);
1405  } else if (lZ < 785 + ecShift) rzLims = MatBounds(2.9, 129, 625 + ecShift, 785 + ecShift);
1406  else if (lZ < 1500 + ecShift) rzLims = MatBounds(2.9, 129, 785 + ecShift, 1500 + ecShift);
1407  else rzLims = MatBounds(2.9, 129, 1500 + ecShift, 1E4);
1408 
1409  //iron .. don't care about no material in front of HF (too forward)
1410  if (! (lZ > 568 + ecShift && lZ < 625 + ecShift && lR > 85 ) // HE support
1411  && ! (lZ > 785 + ecShift && lZ < 850 + ecShift && lR > 118)) {dEdx = dEdX_Fe; radX0 = radX0_Fe; }
1412  else { dEdx = dEdX_MCh; radX0 = radX0_MCh; } //ME at eta > 2.2
1413  }
1414  }
1415  else if (lR < 287){
1416  if (lZ < 372 + ecShift && lR < 177){ // 129<<R<177
1417  if (lZ < 304) rzLims = MatBounds(129, 177, 0, 304); //EB 129<<R<177 0<Z<304
1418  else if (lZ < 343){ // 129<<R<177 304<Z<343
1419  if (lR < 135 ) rzLims = MatBounds(129, 135, 304, 343);// tk piece 129<<R<135 304<Z<343
1420  else if (lR < 172 ){ //
1421  if (lZ < 311 ) rzLims = MatBounds(135, 172, 304, 311);
1422  else rzLims = MatBounds(135, 172, 311, 343);
1423  } else {
1424  if (lZ < 328) rzLims = MatBounds(172, 177, 304, 328);
1425  else rzLims = MatBounds(172, 177, 328, 343);
1426  }
1427  }
1428  else if ( lZ < 343 + ecShift){
1429  rzLims = MatBounds(129, 177, 343, 343 + ecShift); //gap
1430  }
1431  else {
1432  if (lR < 156 ) rzLims = MatBounds(129, 156, 343 + ecShift, 372 + ecShift);
1433  else if ( (lZ - 343 - ecShift) > (lR - 156)*1.38 )
1434  rzLims = MatBounds(156, 177, 127.72 + ecShift, 372 + ecShift, atan(1.38), 0);
1435  else rzLims = MatBounds(156, 177, 343 + ecShift, 127.72 + ecShift, 0, atan(1.38));
1436  }
1437 
1438  if (!(lR > 135 && lZ <343 + ecShift && lZ > 304 )
1439  && ! (lR > 156 && lZ < 372 + ecShift && lZ > 343 + ecShift && ((lZ-343. - ecShift )< (lR-156.)*1.38)))
1440  {
1441  //the crystals are the same length, but they are not 100% of material
1442  double cosThetaEquiv = 0.8/sqrt(1.+lZ*lZ/lR/lR) + 0.2;
1443  if (lZ > 343) cosThetaEquiv = 1.;
1444  dEdx = dEdX_ECal*cosThetaEquiv; radX0 = radX0_ECal/cosThetaEquiv;
1445  } //EB
1446  else {
1447  if ( (lZ > 304 && lZ < 328 && lR < 177 && lR > 135)
1448  && ! (lZ > 311 && lR < 172) ) {dEdx = dEdX_Fe; radX0 = radX0_Fe; } //Tk_Support
1449  else if ( lZ > 343 && lZ < 343 + ecShift) { dEdx = dEdX_Air; radX0 = radX0_Air; }
1450  else {dEdx = dEdX_ECal*0.2; radX0 = radX0_Air;} //cables go here <-- will be abit too dense for ecShift > 0
1451  }
1452  }
1453  else if (lZ < 554 + ecShift){ // 129<R<177 372<Z<554 AND 177<R<287 0<Z<554
1454  if (lR < 177){ // 129<R<177 372<Z<554
1455  if ( lZ > 372 + ecShift && lZ < 398 + ecShift )rzLims = MatBounds(129, 177, 372 + ecShift, 398 + ecShift); // EE 129<R<177 372<Z<398
1456  else if (lZ < 548 + ecShift) rzLims = MatBounds(129, 177, 398 + ecShift, 548 + ecShift); // HE 129<R<177 398<Z<548
1457  else rzLims = MatBounds(129, 177, 548 + ecShift, 554 + ecShift); // HE gap 129<R<177 548<Z<554
1458  }
1459  else if (lR < 193){ // 177<R<193 0<Z<554
1460  if ((lZ - 307) < (lR - 177.)*1.739) rzLims = MatBounds(177, 193, 0, -0.803, 0, atan(1.739));
1461  else if ( lZ < 389) rzLims = MatBounds(177, 193, -0.803, 389, atan(1.739), 0.);
1462  else if ( lZ < 389 + ecShift) rzLims = MatBounds(177, 193, 389, 389 + ecShift); // air insert
1463  else if ( lZ < 548 + ecShift ) rzLims = MatBounds(177, 193, 389 + ecShift, 548 + ecShift);// HE 177<R<193 389<Z<548
1464  else rzLims = MatBounds(177, 193, 548 + ecShift, 554 + ecShift); // HE gap 177<R<193 548<Z<554
1465  }
1466  else if (lR < 264){ // 193<R<264 0<Z<554
1467  double anApex = 375.7278 - 193./1.327; // 230.28695599096
1468  if ( (lZ - 375.7278) < (lR - 193.)/1.327) rzLims = MatBounds(193, 264, 0, anApex, 0, atan(1./1.327)); //HB
1469  else if ( (lZ - 392.7278 ) < (lR - 193.)/1.327)
1470  rzLims = MatBounds(193, 264, anApex, anApex+17., atan(1./1.327), atan(1./1.327)); // HB-HE gap
1471  else if ( (lZ - 392.7278 - ecShift ) < (lR - 193.)/1.327)
1472  rzLims = MatBounds(193, 264, anApex+17., anApex+17. + ecShift, atan(1./1.327), atan(1./1.327)); // HB-HE gap air insert
1473  // HE (372,193)-(517,193)-(517,264)-(417.8,264)
1474  else if ( lZ < 517 + ecShift ) rzLims = MatBounds(193, 264, anApex+17. + ecShift, 517 + ecShift, atan(1./1.327), 0);
1475  else if (lZ < 548 + ecShift){ // HE+gap 193<R<264 517<Z<548
1476  if (lR < 246 ) rzLims = MatBounds(193, 246, 517 + ecShift, 548 + ecShift); // HE 193<R<246 517<Z<548
1477  else rzLims = MatBounds(246, 264, 517 + ecShift, 548 + ecShift); // HE gap 246<R<264 517<Z<548
1478  }
1479  else rzLims = MatBounds(193, 264, 548 + ecShift, 554 + ecShift); // HE gap 193<R<246 548<Z<554
1480  }
1481  else if ( lR < 275){ // 264<R<275 0<Z<554
1482  if (lZ < 433) rzLims = MatBounds(264, 275, 0, 433); //HB
1483  else if (lZ < 554 ) rzLims = MatBounds(264, 275, 433, 554); // HE gap 264<R<275 433<Z<554
1484  else rzLims = MatBounds(264, 275, 554, 554 + ecShift); // HE gap 264<R<275 554<Z<554 air insert
1485  }
1486  else { // 275<R<287 0<Z<554
1487  if (lZ < 402) rzLims = MatBounds(275, 287, 0, 402);// HB 275<R<287 0<Z<402
1488  else if (lZ < 554) rzLims = MatBounds(275, 287, 402, 554);// //HE gap 275<R<287 402<Z<554
1489  else rzLims = MatBounds(275, 287, 554, 554 + ecShift); //HE gap 275<R<287 554<Z<554 air insert
1490  }
1491 
1492  if ((lZ < 433 || lR < 264) && (lZ < 402 || lR < 275) && (lZ < 517 + ecShift || lR < 246) //notches
1493  //I should've made HE and HF different .. now need to shorten HE to match
1494  && lZ < 548 + ecShift
1495  && ! (lZ < 389 + ecShift && lZ > 335 && lR < 193 ) //not a gap EB-EE 129<R<193
1496  && ! (lZ > 307 && lZ < 335 && lR < 193 && ((lZ - 307) > (lR - 177.)*1.739)) //not a gap
1497  && ! (lR < 177 && lZ < 398 + ecShift) //under the HE nose
1498  && ! (lR < 264 && lR > 193 && fabs(441.5+0.5*ecShift - lZ + (lR - 269.)/1.327) < (8.5 + ecShift*0.5)) ) //not a gap
1499  { dEdx = dEdX_HCal; radX0 = radX0_HCal; //hcal
1500  }
1501  else if( (lR < 193 && lZ > 389 && lZ < 389 + ecShift ) || (lR > 264 && lR < 287 && lZ > 554 && lZ < 554 + ecShift)
1502  ||(lR < 264 && lR > 193 && fabs(441.5+8.5+0.5*ecShift - lZ + (lR - 269.)/1.327) < ecShift*0.5) ) {
1503  dEdx = dEdX_Air; radX0 = radX0_Air;
1504  }
1505  else {dEdx = dEdX_HCal*0.05; radX0 = radX0_Air; }//endcap gap
1506  }
1507  //the rest is a tube of endcap volume 129<R<251 554<Z<largeValue
1508  else if (lZ < 564 + ecShift){ // 129<R<287 554<Z<564
1509  if (lR < 251) {
1510  rzLims = MatBounds(129, 251, 554 + ecShift, 564 + ecShift); // HE support 129<R<251 554<Z<564
1511  dEdx = dEdX_Fe; radX0 = radX0_Fe;
1512  }//HE support
1513  else {
1514  rzLims = MatBounds(251, 287, 554 + ecShift, 564 + ecShift); //HE/ME gap 251<R<287 554<Z<564
1515  dEdx = dEdX_MCh; radX0 = radX0_MCh;
1516  }
1517  }
1518  else if (lZ < 625 + ecShift){ // ME/1/1 129<R<287 564<Z<625
1519  rzLims = MatBounds(129, 287, 564 + ecShift, 625 + ecShift);
1520  dEdx = dEdX_MCh; radX0 = radX0_MCh;
1521  }
1522  else if (lZ < 785 + ecShift){ //129<R<287 625<Z<785
1523  if (lR < 275) rzLims = MatBounds(129, 275, 625 + ecShift, 785 + ecShift); //YE/1 129<R<275 625<Z<785
1524  else { // 275<R<287 625<Z<785
1525  if (lZ < 720 + ecShift) rzLims = MatBounds(275, 287, 625 + ecShift, 720 + ecShift); // ME/1/2 275<R<287 625<Z<720
1526  else rzLims = MatBounds(275, 287, 720 + ecShift, 785 + ecShift); // YE/1 275<R<287 720<Z<785
1527  }
1528  if (! (lR > 275 && lZ < 720 + ecShift)) { dEdx = dEdX_Fe; radX0 = radX0_Fe; }//iron
1529  else { dEdx = dEdX_MCh; radX0 = radX0_MCh; }
1530  }
1531  else if (lZ < 850 + ecShift){
1532  rzLims = MatBounds(129, 287, 785 + ecShift, 850 + ecShift);
1533  dEdx = dEdX_MCh; radX0 = radX0_MCh;
1534  }
1535  else if (lZ < 910 + ecShift){
1536  rzLims = MatBounds(129, 287, 850 + ecShift, 910 + ecShift);
1537  dEdx = dEdX_Fe; radX0 = radX0_Fe;
1538  }//iron
1539  else if (lZ < 975 + ecShift){
1540  rzLims = MatBounds(129, 287, 910 + ecShift, 975 + ecShift);
1541  dEdx = dEdX_MCh; radX0 = radX0_MCh;
1542  }
1543  else if (lZ < 1000 + ecShift){
1544  rzLims = MatBounds(129, 287, 975 + ecShift, 1000 + ecShift);
1545  dEdx = dEdX_Fe; radX0 = radX0_Fe;
1546  }//iron
1547  else if (lZ < 1063 + ecShift){
1548  rzLims = MatBounds(129, 287, 1000 + ecShift, 1063 + ecShift);
1549  dEdx = dEdX_MCh; radX0 = radX0_MCh;
1550  }
1551  else if ( lZ < 1073 + ecShift){
1552  rzLims = MatBounds(129, 287, 1063 + ecShift, 1073 + ecShift);
1553  dEdx = dEdX_Fe; radX0 = radX0_Fe;
1554  }
1555  else if (lZ < 1E4 ) {
1556  rzLims = MatBounds(129, 287, 1073 + ecShift, 1E4);
1557  dEdx = dEdX_Air; radX0 = radX0_Air;
1558  }
1559  else {
1560 
1561  dEdx = dEdX_Air; radX0 = radX0_Air;
1562  }
1563  }
1564  else if (lR <380 && lZ < 667){
1565  if (lZ < 630){
1566  if (lR < 315) rzLims = MatBounds(287, 315, 0, 630);
1567  else if (lR < 341 ) rzLims = MatBounds(315, 341, 0, 630); //b-field ~linear rapid fall here
1568  else rzLims = MatBounds(341, 380, 0, 630);
1569  } else rzLims = MatBounds(287, 380, 630, 667);
1570 
1571  if (lZ < 630) { dEdx = dEdX_coil; radX0 = radX0_coil; }//a guess for the solenoid average
1572  else {dEdx = dEdX_Air; radX0 = radX0_Air; }//endcap gap
1573  }
1574  else {
1575  if (lZ < 667) {
1576  if (lR < 850){
1577  bool isIron = false;
1578  //sanity check in addition to flags
1579  if (useIsYokeFlag_ && useMagVolumes_ && sv.magVol != 0){
1580  isIron = sv.isYokeVol;
1581  } else {
1582  double bMag = sv.bf.mag();
1583  isIron = (bMag > 0.75 && ! (lZ > 500 && lR <500 && bMag < 1.15)
1584  && ! (lZ < 450 && lR > 420 && bMag < 1.15 ) );
1585  }
1586  //tell the material stepper where mat bounds are
1587  rzLims = MatBounds(380, 850, 0, 667);
1588  if (isIron) { dEdx = dEdX_Fe; radX0 = radX0_Fe; }//iron
1589  else { dEdx = dEdX_MCh; radX0 = radX0_MCh; }
1590  } else {
1591  rzLims = MatBounds(850, 1E4, 0, 667);
1592  dEdx = dEdX_Air; radX0 = radX0_Air;
1593  }
1594  }
1595  else if (lR > 750 ){
1596  rzLims = MatBounds(750, 1E4, 667, 1E4);
1597  dEdx = dEdX_Air; radX0 = radX0_Air;
1598  }
1599  else if (lZ < 667 + ecShift){
1600  rzLims = MatBounds(287, 750, 667, 667 + ecShift);
1601  dEdx = dEdX_Air; radX0 = radX0_Air;
1602  }
1603  //the rest is endcap piece with 287<R<750 Z>667
1604  else if (lZ < 724 + ecShift){
1605  if (lR < 380 ) rzLims = MatBounds(287, 380, 667 + ecShift, 724 + ecShift);
1606  else rzLims = MatBounds(380, 750, 667 + ecShift, 724 + ecShift);
1607  dEdx = dEdX_MCh; radX0 = radX0_MCh;
1608  }
1609  else if (lZ < 785 + ecShift){
1610  if (lR < 380 ) rzLims = MatBounds(287, 380, 724 + ecShift, 785 + ecShift);
1611  else rzLims = MatBounds(380, 750, 724 + ecShift, 785 + ecShift);
1612  dEdx = dEdX_Fe; radX0 = radX0_Fe;
1613  }//iron
1614  else if (lZ < 850 + ecShift){
1615  rzLims = MatBounds(287, 750, 785 + ecShift, 850 + ecShift);
1616  dEdx = dEdX_MCh; radX0 = radX0_MCh;
1617  }
1618  else if (lZ < 910 + ecShift){
1619  rzLims = MatBounds(287, 750, 850 + ecShift, 910 + ecShift);
1620  dEdx = dEdX_Fe; radX0 = radX0_Fe;
1621  }//iron
1622  else if (lZ < 975 + ecShift){
1623  rzLims = MatBounds(287, 750, 910 + ecShift, 975 + ecShift);
1624  dEdx = dEdX_MCh; radX0 = radX0_MCh;
1625  }
1626  else if (lZ < 1000 + ecShift){
1627  rzLims = MatBounds(287, 750, 975 + ecShift, 1000 + ecShift);
1628  dEdx = dEdX_Fe; radX0 = radX0_Fe;
1629  }//iron
1630  else if (lZ < 1063 + ecShift){
1631  if (lR < 360){
1632  rzLims = MatBounds(287, 360, 1000 + ecShift, 1063 + ecShift);
1633  dEdx = dEdX_MCh; radX0 = radX0_MCh;
1634  }
1635  //put empty air where me4/2 should be (future)
1636  else {
1637  rzLims = MatBounds(360, 750, 1000 + ecShift, 1063 + ecShift);
1638  dEdx = dEdX_Air; radX0 = radX0_Air;
1639  }
1640  }
1641  else if ( lZ < 1073 + ecShift){
1642  rzLims = MatBounds(287, 750, 1063 + ecShift, 1073 + ecShift);
1643  //this plate does not exist: air
1644  dEdx = dEdX_Air; radX0 = radX0_Air;
1645  }
1646  else if (lZ < 1E4 ) {
1647  rzLims = MatBounds(287, 750, 1073 + ecShift, 1E4);
1648  dEdx = dEdX_Air; radX0 = radX0_Air;
1649  }
1650  else {dEdx = dEdX_Air; radX0 = radX0_Air; }//air
1651  }
1652 
1653  //dEdx so far is a relative number (relative to iron)
1654  //scale by what's expected for iron (the function was fit from pdg table)
1655  //0.065 (PDG) --> 0.044 to better match with MPV
1656  double p0 = sv.p3.mag();
1657  double logp0 = log(p0);
1658  double p0powN33 = 0;
1659  if (p0>3.) {
1660  // p0powN33 = exp(-0.33*logp0); //calculate for p>3GeV
1661  double xx=1./p0; xx=sqrt(sqrt(sqrt(sqrt(xx)))); xx=xx*xx*xx*xx*xx; // this is (p0)**(-5/16), close enough to -0.33
1662  p0powN33 = xx;
1663  }
1664  double dEdX_mat = -(11.4 + 0.96*fabs(logp0+log(2.8)) + 0.033*p0*(1.0 - p0powN33) )*1e-3;
1665  //in GeV/cm .. 0.8 to get closer to the median or MPV
1666 
1667  dEdXPrime = dEdx == 0 ? 0 : -dEdx*(0.96/p0 + 0.033 - 0.022*p0powN33)*1e-3; //== d(dEdX)/dp
1668  dEdx *= dEdX_mat;
1669 
1670  return dEdx;
1671 }
1672 
1673 
1675  int result = ind%MAX_POINTS;
1676  if (ind != 0 && result == 0){
1677  result = MAX_POINTS;
1678  }
1679  return result;
1680 }
1681 
1685  const double pars[6],
1686  double& dist, double& tanDist,
1687  PropagationDirection& refDirection,
1688  double fastSkipDist) const{
1689  static const std::string metname = "SteppingHelixPropagator";
1691 
1692  switch (dest){
1693  case RADIUS_DT:
1694  {
1695  double curR = sv.r3.perp();
1696  dist = pars[RADIUS_P] - curR;
1697  if (fabs(dist) > fastSkipDist){
1699  break;
1700  }
1701  double curP2 = sv.p3.mag2();
1702  double curPtPos2 = sv.p3.perp2(); if(curPtPos2< 1e-16) curPtPos2=1e-16;
1703 
1704  double cosDPhiPR = (sv.r3.x()*sv.p3.x()+sv.r3.y()*sv.p3.y());//only the sign is needed cos((sv.r3.deltaPhi(sv.p3)));
1705  refDirection = dist*cosDPhiPR > 0 ?
1707  tanDist = dist*sqrt(curP2/curPtPos2);
1708  result = SteppingHelixStateInfo::OK;
1709  }
1710  break;
1711  case Z_DT:
1712  {
1713  double curZ = sv.r3.z();
1714  dist = pars[Z_P] - curZ;
1715  if (fabs(dist) > fastSkipDist){
1717  break;
1718  }
1719  double curP = sv.p3.mag();
1720  refDirection = sv.p3.z()*dist > 0. ?
1722  tanDist = dist/sv.p3.z()*curP;
1723  result = SteppingHelixStateInfo::OK;
1724  }
1725  break;
1726  case PLANE_DT:
1727  {
1728  Point rPlane(pars[0], pars[1], pars[2]);
1729  Vector nPlane(pars[3], pars[4], pars[5]);
1730 
1731 
1732  // unfortunately this doesn't work: the numbers are too large
1733  // bool pVertical = fabs(pars[5])>0.9999;
1734  // double dRDotN = pVertical? (sv.r3.z() - rPlane.z())*nPlane.z() :(sv.r3 - rPlane).dot(nPlane);
1735  double dRDotN = (sv.r3.x()-rPlane.x())*nPlane.x() + (sv.r3.y()-rPlane.y())*nPlane.y() + (sv.r3.z()-rPlane.z())*nPlane.z();//(sv.r3 - rPlane).dot(nPlane);
1736 
1737  dist = fabs(dRDotN);
1738  if (dist > fastSkipDist){
1740  break;
1741  }
1742  double curP = sv.p3.mag();
1743  double p0 = curP;
1744  double p0Inv = 1./p0;
1745  Vector tau(sv.p3); tau *=p0Inv;
1746  double tN = tau.dot(nPlane);
1747  refDirection = tN*dRDotN < 0. ?
1749  double b0 = sv.bf.mag();
1750  if (fabs(tN)>1e-24){
1751  tanDist = -dRDotN/tN;
1752  } else {
1753  tN = 1e-24;
1754  if (fabs(dRDotN)>1e-24) tanDist = 1e6;
1755  else tanDist = 1;
1756  }
1757  if (fabs(tanDist) > 1e4) tanDist = 1e4;
1758  if (b0>1.5e-6){
1759  double b0Inv = 1./b0;
1760  double tNInv = 1./tN;
1761  Vector bHat(sv.bf); bHat *=b0Inv;
1762  double bHatN = bHat.dot(nPlane);
1763  double cosPB = bHat.dot(tau);
1764  double kVal = 2.99792458e-3*sv.q*p0Inv*b0;
1765  double aVal = tanDist*kVal;
1766  Vector lVec = bHat.cross(tau);
1767  double bVal = lVec.dot(nPlane)*tNInv;
1768  if (fabs(aVal*bVal)< 0.3){
1769  double cVal = 1. - bHatN*cosPB*tNInv; // - sv.bf.cross(lVec).dot(nPlane)*b0Inv*tNInv; //1- bHat_n*bHat_tau/tau_n;
1770  double aacVal = cVal*aVal*aVal;
1771  if (fabs(aacVal)<1){
1772  double tanDCorr = bVal/2. + (bVal*bVal/2. + cVal/6)*aVal;
1773  tanDCorr *= aVal;
1774  //+ (-bVal/24. + 0.625*bVal*bVal*bVal + 5./12.*bVal*cVal)*aVal*aVal*aVal
1775  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<tanDist<<" vs "
1776  <<tanDist*(1.+tanDCorr)<<" corr "<<tanDist*tanDCorr<<std::endl;
1777  tanDist *= (1.+tanDCorr);
1778  } else {
1779  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"AACVal "<< fabs(aacVal)
1780  <<" = "<<aVal<<"**2 * "<<cVal<<" too large:: will not converge"<<std::endl;
1781  }
1782  } else {
1783  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"ABVal "<< fabs(aVal*bVal)
1784  <<" = "<<aVal<<" * "<<bVal<<" too large:: will not converge"<<std::endl;
1785  }
1786  }
1787  result = SteppingHelixStateInfo::OK;
1788  }
1789  break;
1790  case CONE_DT:
1791  {
1792  Point cVertex(pars[0], pars[1], pars[2]);
1793  if (cVertex.perp2() < 1e-10){
1794  //assumes the cone axis/vertex is along z
1795  Vector relV3 = sv.r3 - cVertex;
1796  double relV3mag = relV3.mag();
1797  // double relV3Theta = relV3.theta();
1798  double theta(pars[3]);
1799  // double dTheta = theta-relV3Theta;
1800  double sinTheta = sin(theta);
1801  double cosTheta = cos(theta);
1802  double cosV3Theta = relV3.z()/relV3mag;
1803  if (cosV3Theta>1) cosV3Theta=1;
1804  if (cosV3Theta<-1) cosV3Theta=-1;
1805  double sinV3Theta = sqrt(1.-cosV3Theta*cosV3Theta);
1806 
1807  double sinDTheta = sinTheta*cosV3Theta - cosTheta*sinV3Theta;//sin(dTheta);
1808  double cosDTheta = cosTheta*cosV3Theta + sinTheta*sinV3Theta;//cos(dTheta);
1809  bool isInside = sinTheta > sinV3Theta && cosTheta*cosV3Theta > 0;
1810  dist = isInside || cosDTheta > 0 ? relV3mag*sinDTheta : relV3mag;
1811  if (fabs(dist) > fastSkipDist){
1813  break;
1814  }
1815 
1816  double relV3phi=relV3.phi();
1817  double normPhi = isInside ?
1818  Geom::pi() + relV3phi : relV3phi;
1819  double normTheta = theta > Geom::pi()/2. ?
1820  ( isInside ? 1.5*Geom::pi() - theta : theta - Geom::pi()/2. )
1821  : ( isInside ? Geom::pi()/2. - theta : theta + Geom::pi()/2 );
1822  //this is a normVector from the cone to the point
1823  Vector norm; norm.setRThetaPhi(fabs(dist), normTheta, normPhi);
1824  double curP = sv.p3.mag(); double cosp3theta = sv.p3.z()/curP;
1825  if (cosp3theta>1) cosp3theta=1;
1826  if (cosp3theta<-1) cosp3theta=-1;
1827  double sineConeP = sinTheta*cosp3theta - cosTheta*sqrt(1.-cosp3theta*cosp3theta);
1828 
1829  double sinSolid = norm.dot(sv.p3)/(fabs(dist)*curP);
1830  tanDist = fabs(sinSolid) > fabs(sineConeP) ? dist/fabs(sinSolid) : dist/fabs(sineConeP);
1831 
1832 
1833  refDirection = sinSolid > 0 ?
1835  if (debug_){
1836  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Cone pars: theta "<<theta
1837  <<" normTheta "<<norm.theta()
1838  <<" p3Theta "<<sv.p3.theta()
1839  <<" sinD: "<< sineConeP
1840  <<" sinSolid "<<sinSolid;
1841  }
1842  if (debug_){
1843  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"refToDest:toCone the point is "
1844  <<(isInside? "in" : "out")<<"side the cone"
1845  <<std::endl;
1846  }
1847  }
1848  result = SteppingHelixStateInfo::OK;
1849  }
1850  break;
1851  // case CYLINDER_DT:
1852  // break;
1853  case PATHL_DT:
1854  {
1855  double curS = fabs(sv.path());
1856  dist = pars[PATHL_P] - curS;
1857  if (fabs(dist) > fastSkipDist){
1859  break;
1860  }
1861  refDirection = pars[PATHL_P] > 0 ?
1863  tanDist = dist;
1864  result = SteppingHelixStateInfo::OK;
1865  }
1866  break;
1867  case POINT_PCA_DT:
1868  {
1869  Point pDest(pars[0], pars[1], pars[2]);
1870  double curP = sv.p3.mag();
1871  dist = (sv.r3 - pDest).mag()+ 1e-24;//add a small number to avoid 1/0
1872  tanDist = (sv.r3.dot(sv.p3) - pDest.dot(sv.p3))/curP;
1873  //account for bending in magnetic field (quite approximate)
1874  double b0 = sv.bf.mag();
1875  if (b0>1.5e-6){
1876  double p0 = curP;
1877  double kVal = 2.99792458e-3*sv.q/p0*b0;
1878  double aVal = fabs(dist*kVal);
1879  tanDist *= 1./(1.+ aVal);
1880  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"corrected by aVal "<<aVal<<" to "<<tanDist;
1881  }
1882  refDirection = tanDist < 0 ?
1884  result = SteppingHelixStateInfo::OK;
1885  }
1886  break;
1887  case LINE_PCA_DT:
1888  {
1889  Point rLine(pars[0], pars[1], pars[2]);
1890  Vector dLine(pars[3], pars[4], pars[5]);
1891  dLine = (dLine - rLine);
1892  dLine *= 1./dLine.mag(); if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"dLine "<<dLine;
1893 
1894  Vector dR = sv.r3 - rLine;
1895  double curP = sv.p3.mag();
1896  Vector dRPerp = dR - dLine*(dR.dot(dLine));
1897  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"dRperp "<<dRPerp;
1898 
1899  dist = dRPerp.mag() + 1e-24;//add a small number to avoid 1/0
1900  tanDist = dRPerp.dot(sv.p3)/curP;
1901  //angle wrt line
1902  double cosAlpha2 = dLine.dot(sv.p3)/curP; cosAlpha2 *= cosAlpha2;
1903  tanDist *= 1./sqrt(fabs(1.-cosAlpha2)+1e-96);
1904  //correct for dPhi in magnetic field: this isn't made quite right here
1905  //(the angle between the line and the trajectory plane is neglected .. conservative)
1906  double b0 = sv.bf.mag();
1907  if (b0>1.5e-6){
1908  double p0 = curP;
1909  double kVal = 2.99792458e-3*sv.q/p0*b0;
1910  double aVal = fabs(dist*kVal);
1911  tanDist *= 1./(1.+ aVal);
1912  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"corrected by aVal "<<aVal<<" to "<<tanDist;
1913  }
1914  refDirection = tanDist < 0 ?
1916  result = SteppingHelixStateInfo::OK;
1917  }
1918  break;
1919  default:
1920  {
1921  //some large number
1922  dist = 1e12;
1923  tanDist = 1e12;
1924  refDirection = anyDirection;
1926  }
1927  break;
1928  }
1929 
1930  double tanDistConstrained = tanDist;
1931  if (fabs(tanDist) > 4.*fabs(dist) ) tanDistConstrained *= tanDist == 0 ? 0 : fabs(dist/tanDist*4.);
1932 
1933  if (debug_){
1934  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"refToDest input: dest"<<dest<<" pars[]: ";
1935  for (int i = 0; i < 6; i++){
1936  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<", "<<i<<" "<<pars[i];
1937  }
1938  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<std::endl;
1939  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"refToDest output: "
1940  <<"\t dist" << dist
1941  <<"\t tanDist "<< tanDist
1942  <<"\t tanDistConstr "<< tanDistConstrained
1943  <<"\t refDirection "<< refDirection
1944  <<std::endl;
1945  }
1946  tanDist = tanDistConstrained;
1947 
1948  return result;
1949 }
1950 
1954  double& dist, double& tanDist,
1955  double fastSkipDist, bool expectNewMagVolume, double maxStep) const{
1956 
1957  static const std::string metname = "SteppingHelixPropagator";
1959  const MagVolume* cVol = sv.magVol;
1960 
1961  if (cVol == 0) return result;
1962  const std::vector<VolumeSide>& cVolFaces(cVol->faces());
1963 
1964  double distToFace[6] = {0,0,0,0,0,0};
1965  double tanDistToFace[6] = {0,0,0,0,0,0};
1966  PropagationDirection refDirectionToFace[6] = {anyDirection, anyDirection, anyDirection, anyDirection, anyDirection, anyDirection};
1967  Result resultToFace[6] = {result, result, result, result, result, result};
1968  int iFDest = -1;
1969  int iDistMin = -1;
1970 
1971  unsigned int iFDestSorted[6] = {0,0,0,0,0,0};
1972  unsigned int nDestSorted =0;
1973  unsigned int nearParallels = 0;
1974 
1975  double curP = sv.p3.mag();
1976 
1977  if (debug_){
1978  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Trying volume "<<DDSolidShapesName::name(cVol->shapeType())
1979  <<" with "<<cVolFaces.size()<<" faces"<<std::endl;
1980  }
1981 
1982  unsigned int nFaces = cVolFaces.size();
1983  for (unsigned int iFace = 0; iFace < nFaces; ++iFace){
1984  if (iFace > 5){
1985  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Too many faces"<<std::endl;
1986  }
1987  if (debug_){
1988  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Start with face "<<iFace<<std::endl;
1989  }
1990 // const Plane* cPlane = dynamic_cast<const Plane*>(&cVolFaces[iFace].surface());
1991 // const Cylinder* cCyl = dynamic_cast<const Cylinder*>(&cVolFaces[iFace].surface());
1992 // const Cone* cCone = dynamic_cast<const Cone*>(&cVolFaces[iFace].surface());
1993  const Surface* cPlane = 0; //only need to know the loc->glob transform
1994  const Cylinder* cCyl = 0;
1995  const Cone* cCone = 0;
1996  if (typeid(cVolFaces[iFace].surface()) == typeid(const Plane&)){
1997  cPlane = &cVolFaces[iFace].surface();
1998  } else if (typeid(cVolFaces[iFace].surface()) == typeid(const Cylinder&)){
1999  cCyl = dynamic_cast<const Cylinder*>(&cVolFaces[iFace].surface());
2000  } else if (typeid(cVolFaces[iFace].surface()) == typeid(const Cone&)){
2001  cCone = dynamic_cast<const Cone*>(&cVolFaces[iFace].surface());
2002  } else {
2003  edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Could not cast a volume side surface to a known type"<<std::endl;
2004  }
2005 
2006  if (debug_){
2007  if (cPlane!=0) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"The face is a plane at "<<cPlane<<std::endl;
2008  if (cCyl!=0) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"The face is a cylinder at "<<cCyl<<std::endl;
2009  }
2010 
2011  double pars[6] = {0,0,0,0,0,0};
2012  DestType dType = UNDEFINED_DT;
2013  if (cPlane != 0){
2014  Point rPlane(cPlane->position().x(),cPlane->position().y(),cPlane->position().z());
2015  // = cPlane->toGlobal(LocalVector(0,0,1.)); nPlane = nPlane.unit();
2016  Vector nPlane(cPlane->rotation().zx(), cPlane->rotation().zy(), cPlane->rotation().zz()); nPlane /= nPlane.mag();
2017 
2018  if (sv.r3.z() < 0 || !vbField_->isZSymmetric() ){
2019  pars[0] = rPlane.x(); pars[1] = rPlane.y(); pars[2] = rPlane.z();
2020  pars[3] = nPlane.x(); pars[4] = nPlane.y(); pars[5] = nPlane.z();
2021  } else {
2022  pars[0] = rPlane.x(); pars[1] = rPlane.y(); pars[2] = -rPlane.z();
2023  pars[3] = nPlane.x(); pars[4] = nPlane.y(); pars[5] = -nPlane.z();
2024  }
2025  dType = PLANE_DT;
2026  } else if (cCyl != 0){
2027  if (debug_){
2028  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Cylinder at "<<cCyl->position()
2029  <<" rorated by "<<cCyl->rotation()
2030  <<std::endl;
2031  }
2032  pars[RADIUS_P] = cCyl->radius();
2033  dType = RADIUS_DT;
2034  } else if (cCone != 0){
2035  if (debug_){
2036  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Cone at "<<cCone->position()
2037  <<" rorated by "<<cCone->rotation()
2038  <<" vertex at "<<cCone->vertex()
2039  <<" angle of "<<cCone->openingAngle()
2040  <<std::endl;
2041  }
2042  if (sv.r3.z() < 0 || !vbField_->isZSymmetric()){
2043  pars[0] = cCone->vertex().x(); pars[1] = cCone->vertex().y();
2044  pars[2] = cCone->vertex().z();
2045  pars[3] = cCone->openingAngle();
2046  } else {
2047  pars[0] = cCone->vertex().x(); pars[1] = cCone->vertex().y();
2048  pars[2] = -cCone->vertex().z();
2049  pars[3] = Geom::pi() - cCone->openingAngle();
2050  }
2051  dType = CONE_DT;
2052  } else {
2053  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Unknown surface"<<std::endl;
2054  resultToFace[iFace] = SteppingHelixStateInfo::UNDEFINED;
2055  continue;
2056  }
2057  resultToFace[iFace] =
2058  refToDest(dType, sv, pars,
2059  distToFace[iFace], tanDistToFace[iFace], refDirectionToFace[iFace], fastSkipDist);
2060 
2061 
2062  if (resultToFace[iFace] != SteppingHelixStateInfo::OK){
2063  if (resultToFace[iFace] == SteppingHelixStateInfo::INACC) result = SteppingHelixStateInfo::INACC;
2064  }
2065 
2066 
2067 
2068  //keep those in right direction for later use
2069  if (resultToFace[iFace] == SteppingHelixStateInfo::OK){
2070  double invDTFPosiF = 1./(1e-32+fabs(tanDistToFace[iFace]));
2071  double dSlope = fabs(distToFace[iFace])*invDTFPosiF;
2072  double maxStepL = maxStep> 100 ? 100 : maxStep; if (maxStepL < 10) maxStepL = 10.;
2073  bool isNearParallel = fabs(tanDistToFace[iFace]) + 100.*curP*dSlope < maxStepL //
2074  //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
2075  && dSlope < 0.15 ; //
2076  if (refDirectionToFace[iFace] == dir || isNearParallel){
2077  if (isNearParallel) nearParallels++;
2078  iFDestSorted[nDestSorted] = iFace;
2079  nDestSorted++;
2080  }
2081  }
2082 
2083  //pick a shortest distance here (right dir only for now)
2084  if (refDirectionToFace[iFace] == dir){
2085  if (iDistMin == -1) iDistMin = iFace;
2086  else if (fabs(distToFace[iFace]) < fabs(distToFace[iDistMin])) iDistMin = iFace;
2087  }
2088  if (debug_)
2089  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<cVol<<" "<<iFace<<" "
2090  <<tanDistToFace[iFace]<<" "<<distToFace[iFace]<<" "<<refDirectionToFace[iFace]<<" "<<dir<<std::endl;
2091 
2092  }
2093 
2094  for (unsigned int i = 0;i<nDestSorted; ++i){
2095  unsigned int iMax = nDestSorted-i-1;
2096  for (unsigned int j=0;j<nDestSorted-i; ++j){
2097  if (fabs(tanDistToFace[iFDestSorted[j]]) > fabs(tanDistToFace[iFDestSorted[iMax]]) ){
2098  iMax = j;
2099  }
2100  }
2101  unsigned int iTmp = iFDestSorted[nDestSorted-i-1];
2102  iFDestSorted[nDestSorted-i-1] = iFDestSorted[iMax];
2103  iFDestSorted[iMax] = iTmp;
2104  }
2105 
2106  if (debug_){
2107  for (unsigned int i=0;i<nDestSorted;++i){
2108  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<cVol<<" "<<i<<" "<<iFDestSorted[i]<<" "<<tanDistToFace[iFDestSorted[i]]<<std::endl;
2109  }
2110  }
2111 
2112  //now go from the shortest to the largest distance hoping to get a point in the volume.
2113  //other than in case of a near-parallel travel this should stop after the first try
2114 
2115  for (unsigned int i=0; i<nDestSorted;++i){
2116  iFDest = iFDestSorted[i];
2117 
2118  double sign = dir == alongMomentum ? 1. : -1.;
2119  Point gPointEst(sv.r3);
2120  Vector lDelta(sv.p3); lDelta *= sign/curP*fabs(distToFace[iFDest]);
2121  gPointEst += lDelta;
2122  if (debug_){
2123  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Linear est point "<<gPointEst
2124  <<" for iFace "<<iFDest<<std::endl;
2125  }
2126  GlobalPoint gPointEstNegZ(gPointEst.x(), gPointEst.y(), -fabs(gPointEst.z()));
2127  GlobalPoint gPointEstNorZ(gPointEst.x(), gPointEst.y(), gPointEst.z() );
2128  if ( (vbField_->isZSymmetric() && cVol->inside(gPointEstNegZ))
2129  || ( !vbField_->isZSymmetric() && cVol->inside(gPointEstNorZ) ) ){
2130  if (debug_){
2131  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"The point is inside the volume"<<std::endl;
2132  }
2133 
2134  result = SteppingHelixStateInfo::OK;
2135  dist = distToFace[iFDest];
2136  tanDist = tanDistToFace[iFDest];
2137  if (debug_){
2138  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Got a point near closest boundary -- face "<<iFDest<<std::endl;
2139  }
2140  break;
2141  }
2142  }
2143 
2144  if (result != SteppingHelixStateInfo::OK && expectNewMagVolume){
2145  double sign = dir == alongMomentum ? 1. : -1.;
2146 
2147  //check if it's a wrong volume situation
2148  if (nDestSorted-nearParallels > 0 ) result = SteppingHelixStateInfo::WRONG_VOLUME;
2149  else {
2150  //get here if all faces in the corr direction were skipped
2151  Point gPointEst(sv.r3);
2152  double lDist = fabs(distToFace[iDistMin]);
2153  if (lDist > fastSkipDist) lDist = fastSkipDist;
2154  Vector lDelta(sv.p3); lDelta *= sign/curP*lDist;
2155  gPointEst += lDelta;
2156  if (debug_){
2157  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Linear est point to shortest dist "<<gPointEst
2158  <<" for iFace "<<iDistMin<<" at distance "<<lDist*sign<<std::endl;
2159  }
2160  GlobalPoint gPointEstNegZ(gPointEst.x(), gPointEst.y(), -fabs(gPointEst.z()));
2161  GlobalPoint gPointEstNorZ(gPointEst.x(), gPointEst.y(), gPointEst.z() );
2162  if ( (vbField_->isZSymmetric() && cVol->inside(gPointEstNegZ))
2163  || ( !vbField_->isZSymmetric() && cVol->inside(gPointEstNorZ) ) ){
2164  if (debug_){
2165  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"The point is inside the volume"<<std::endl;
2166  }
2167 
2168  }else {
2170  }
2171  }
2172 
2173  if (result == SteppingHelixStateInfo::WRONG_VOLUME){
2174  dist = sign*0.05;
2175  tanDist = dist*1.01;
2176  if( debug_){
2177  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Wrong volume located: return small dist, tandist"<<std::endl;
2178  }
2179  }
2180  }
2181 
2182  if (result == SteppingHelixStateInfo::INACC){
2183  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"All faces are too far"<<std::endl;
2184  } else if (result == SteppingHelixStateInfo::WRONG_VOLUME){
2185  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Appear to be in a wrong volume"<<std::endl;
2186  } else if (result != SteppingHelixStateInfo::OK){
2187  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Something else went wrong"<<std::endl;
2188  }
2189 
2190 
2191  return result;
2192 }
2193 
2194 
2198  double& dist, double& tanDist,
2199  double fastSkipDist) const{
2200 
2201  static const std::string metname = "SteppingHelixPropagator";
2203 
2204  double parLim[6] = {sv.rzLims.rMin, sv.rzLims.rMax,
2205  sv.rzLims.zMin, sv.rzLims.zMax,
2206  sv.rzLims.th1, sv.rzLims.th2 };
2207 
2208  double distToFace[4] = {0,0,0,0};
2209  double tanDistToFace[4] = {0,0,0,0};
2210  PropagationDirection refDirectionToFace[4] = {anyDirection, anyDirection, anyDirection, anyDirection};
2211  Result resultToFace[4] = {result, result, result, result};
2212  int iFDest = -1;
2213 
2214  double curP = sv.p3.mag();
2215 
2216  for (unsigned int iFace = 0; iFace < 4; iFace++){
2217  if (debug_){
2218  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Start with mat face "<<iFace<<std::endl;
2219  }
2220 
2221  double pars[6] = {0,0,0,0,0,0};
2222  DestType dType = UNDEFINED_DT;
2223  if (iFace > 1){
2224  if (fabs(parLim[iFace+2])< 1e-6){//plane
2225  if (sv.r3.z() < 0){
2226  pars[0] = 0; pars[1] = 0; pars[2] = -parLim[iFace];
2227  pars[3] = 0; pars[4] = 0; pars[5] = 1;
2228  } else {
2229  pars[0] = 0; pars[1] = 0; pars[2] = parLim[iFace];
2230  pars[3] = 0; pars[4] = 0; pars[5] = 1;
2231  }
2232  dType = PLANE_DT;
2233  } else {
2234  if (sv.r3.z() > 0){
2235  pars[0] = 0; pars[1] = 0;
2236  pars[2] = parLim[iFace];
2237  pars[3] = Geom::pi()/2. - parLim[iFace+2];
2238  } else {
2239  pars[0] = 0; pars[1] = 0;
2240  pars[2] = - parLim[iFace];
2241  pars[3] = Geom::pi()/2. + parLim[iFace+2];
2242  }
2243  dType = CONE_DT;
2244  }
2245  } else {
2246  pars[RADIUS_P] = parLim[iFace];
2247  dType = RADIUS_DT;
2248  }
2249 
2250  resultToFace[iFace] =
2251  refToDest(dType, sv, pars,
2252  distToFace[iFace], tanDistToFace[iFace], refDirectionToFace[iFace], fastSkipDist);
2253 
2254  if (resultToFace[iFace] != SteppingHelixStateInfo::OK){
2255  if (resultToFace[iFace] == SteppingHelixStateInfo::INACC) result = SteppingHelixStateInfo::INACC;
2256  continue;
2257  }
2258  if (refDirectionToFace[iFace] == dir || fabs(distToFace[iFace]) < 2e-2*fabs(tanDistToFace[iFace]) ){
2259  double sign = dir == alongMomentum ? 1. : -1.;
2260  Point gPointEst(sv.r3);
2261  Vector lDelta(sv.p3); lDelta *= sign*fabs(distToFace[iFace])/curP;
2262  gPointEst += lDelta;
2263  if (debug_){
2264  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Linear est point "<<gPointEst
2265  <<std::endl;
2266  }
2267  double lZ = fabs(gPointEst.z());
2268  double lR = gPointEst.perp();
2269  double tan4 = parLim[4] == 0 ? 0 : tan(parLim[4]);
2270  double tan5 = parLim[5] == 0 ? 0 : tan(parLim[5]);
2271  if ( (lZ - parLim[2]) > lR*tan4
2272  && (lZ - parLim[3]) < lR*tan5
2273  && lR > parLim[0] && lR < parLim[1]){
2274  if (debug_){
2275  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"The point is inside the volume"<<std::endl;
2276  }
2277  //OK, guessed a point still inside the volume
2278  if (iFDest == -1){
2279  iFDest = iFace;
2280  } else {
2281  if (fabs(tanDistToFace[iFDest]) > fabs(tanDistToFace[iFace])){
2282  iFDest = iFace;
2283  }
2284  }
2285  } else {
2286  if (debug_){
2287  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"The point is NOT inside the volume"<<std::endl;
2288  }
2289  }
2290  }
2291 
2292  }
2293  if (iFDest != -1){
2294  result = SteppingHelixStateInfo::OK;
2295  dist = distToFace[iFDest];
2296  tanDist = tanDistToFace[iFDest];
2297  if (debug_){
2298  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Got a point near closest boundary -- face "<<iFDest<<std::endl;
2299  }
2300  } else {
2301  if (debug_){
2302  LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Failed to find a dest point inside the volume"<<std::endl;
2303  }
2304  }
2305 
2306  return result;
2307 }
2308 
2309 
2311  static const std::string metname = "SteppingHelixPropagator";
2312  if (vol == 0) return false;
2313  /*
2314  const MFGrid* mGrid = reinterpret_cast<const MFGrid*>(vol->provider());
2315  std::vector<int> dims(mGrid->dimensions());
2316 
2317  LocalVector lVCen(mGrid->nodeValue(dims[0]/2, dims[1]/2, dims[2]/2));
2318  LocalVector lVZLeft(mGrid->nodeValue(dims[0]/2, dims[1]/2, dims[2]/5));
2319  LocalVector lVZRight(mGrid->nodeValue(dims[0]/2, dims[1]/2, (dims[2]*4)/5));
2320 
2321  double mag2VCen = lVCen.mag2();
2322  double mag2VZLeft = lVZLeft.mag2();
2323  double mag2VZRight = lVZRight.mag2();
2324 
2325  bool result = false;
2326  if (mag2VCen > 0.6 && mag2VZLeft > 0.6 && mag2VZRight > 0.6){
2327  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Volume is magnetic, located at "<<vol->position()<<std::endl;
2328  result = true;
2329  } else {
2330  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Volume is not magnetic, located at "<<vol->position()<<std::endl;
2331  }
2332 
2333  */
2334  bool result = vol->isIron();
2335  if (result){
2336  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Volume is magnetic, located at "<<vol->position()<<std::endl;
2337  } else {
2338  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Volume is not magnetic, located at "<<vol->position()<<std::endl;
2339  }
2340 
2341  return result;
2342 }
type
Definition: HCALResponse.h:22
double z0() const
z coordinate
Definition: BeamSpot.h:69
int i
Definition: DBlmapReader.cc:9
T mag2() const
Definition: PV3DBase.h:60
Definition: Cone.h:20
DDSolidShape shapeType() const
Definition: MagVolume.h:31
virtual PropagationDirection propagationDirection() const
Definition: Propagator.h:143
tuple pp
Definition: createTree.py:15
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 ...
const std::string metname
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
void getFreeState(FreeTrajectoryState &fts) const
convert internal structure into the fts
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
bool isYokeVolume(const MagVolume *vol) const
check if it&#39;s a yoke/iron based on this MagVol internals
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
SteppingHelixStateInfo::VolumeBounds MatBounds
T y() const
Definition: PV3DBase.h:57
AlgebraicSymMatrix55 covCurv
const MagVolume * findVolume(const GlobalPoint &gp) const
std::pair< TrajectoryStateOnSurface, double > TsosPP
Geom::Theta< T > theta() const
int cIndex_(int ind) const
(Internals) circular index for array of transients
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
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 ...
PropagationDirection
std::pair< FreeTrajectoryState, double > FtsPP
Result refToMatVolume(const SteppingHelixPropagator::StateInfo &sv, PropagationDirection dir, double &dist, double &tanDist, double fastSkipDist=1e12) const
double charge(const std::vector< uint8_t > &Ampls)
Definition: Plane.h:17
void loadState(SteppingHelixPropagator::StateInfo &svCurrent, const SteppingHelixPropagator::Vector &p3, const SteppingHelixPropagator::Point &r3, int charge, PropagationDirection dir, const AlgebraicSymMatrix55 &covCurv) const
virtual ::GlobalVector inTesla(const ::GlobalPoint &gp) const
Definition: MagVolume.h:46
T zx() const
double dydz() const
dydz slope
Definition: BeamSpot.h:75
GlobalPoint vertex() const
Global position of the cone vertex.
Definition: Cone.h:50
T zz() const
T mag() const
Definition: PV3DBase.h:61
virtual bool inside(const GlobalPoint &gp, double tolerance=0.) const =0
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Scalar radius() const
Radius of the cylinder.
Definition: Cylinder.h:55
static const char * name(DDSolidShape s)
Definition: DDSolidShapes.h:20
TrajectoryStateOnSurface getStateOnSurface(const Surface &surf, bool returnTangentPlane=false) const
T sqrt(T t)
Definition: SSEVec.h:28
LocalPoint toLocal(const GlobalPoint &gp) const
tuple norm
Definition: lumiNorm.py:78
const VolumeBasedMagneticField * vbField_
T z() const
Definition: PV3DBase.h:58
tuple result
Definition: query.py:137
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
T zy() const
int j
Definition: DBlmapReader.cc:9
AlgebraicSymMatrix55 matDCovCurv
virtual TrajectoryStateOnSurface propagate(const FreeTrajectoryState &ftsStart, const Plane &pDest) const
Propagate to Plane given a starting point.
Common base class.
#define LogTrace(id)
double dxdz() const
dxdz slope
Definition: BeamSpot.h:73
const MagneticField * field
bool isIron() const
Temporary hack to pass information on material. Will eventually be replaced!
Definition: MagVolume.h:51
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
Result refToMagVolume(const SteppingHelixPropagator::StateInfo &sv, PropagationDirection dir, double &dist, double &tanDist, double fastSkipDist=1e12, bool expectNewMagVolume=false, double maxStep=1e12) const
Log< T >::type log(const T &t)
Definition: Log.h:22
virtual std::pair< TrajectoryStateOnSurface, double > propagateWithPath(const FreeTrajectoryState &ftsStart, const Plane &pDest) const
virtual const std::vector< VolumeSide > & faces() const =0
Access to volume faces.
const AlgebraicSymMatrix55 unit55_
void setRep(SteppingHelixPropagator::Basis &rep, const SteppingHelixPropagator::Vector &tau) const
Set/compute basis vectors for local coordinates at current step (called by incrementState) ...
LocalVector fieldInTesla(const LocalPoint &lp) const
Definition: MagVolume.cc:11
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
StateInfo svBuf_[MAX_POINTS+1]
double pi()
Definition: Pi.h:31
double y0() const
y coordinate
Definition: BeamSpot.h:67
void setIState(const SteppingHelixStateInfo &sStart) const
(Internals) Init starting point
const RotationType & rotation() const
dbl *** dir
Definition: mlp_gen.cc:35
const double epsilon
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
T x() const
Definition: PV3DBase.h:56
const PositionType & position() const
const MagneticField * field_
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepStd< double, 5, 5 > > AlgebraicMatrix55
static const std::string ResultName[MAX_RESULT]
Geom::Theta< float > openingAngle() const
Angle of the cone.
Definition: Cone.h:53
double p3[4]
Definition: TauolaWrapper.h:91
double x0() const
x coordinate
Definition: BeamSpot.h:65
Definition: DDAxes.h:10