1 ///////////////////////////////////////////////////////////////////////////
3 // AliFmHelix: a helper helix class //
4 // Includes all the operations and specifications of the helix. Can be //
5 // used to determine path lengths, distance of closest approach etc. //
7 ///////////////////////////////////////////////////////////////////////////
8 #if !defined(ST_NO_NUMERIC_LIMITS)
10 # if !defined(ST_NO_NAMESPACES)
11 using std::numeric_limits;
18 #include "AliFmHelix.h"
19 #include "PhysicalConstants.h"
20 #include "SystemOfUnits.h"
22 ClassImpT(AliFmHelix,double);
33 const double AliFmHelix::fgkNoSolution = 3.e+33;
35 AliFmHelix::AliFmHelix() :
51 AliFmHelix::AliFmHelix(double c, double d, double phase,
52 const AliFmThreeVector<double>& o, int h) :
64 // Constructor with helix parameters
65 SetParameters(c, d, phase, o, h);
68 AliFmHelix::~AliFmHelix() {
73 void AliFmHelix::SetParameters(double c, double dip, double phase,
74 const AliFmThreeVector<double>& o, int h)
77 // The order in which the parameters are set is important
78 // since setCurvature might have to adjust the others.
80 fH = (h>=0) ? 1 : -1; // Default is: positive particle
87 // Check for singularity and correct for negative curvature.
88 // May change mH and mPhase. Must therefore be set last.
93 // For the case B=0, h is ill defined. In the following we
94 // always assume h = +1. Since phase = psi - h * pi/2
95 // we have to correct the phase in case h = -1.
96 // This assumes that the user uses the same h for phase
97 // as the one he passed to the constructor.
99 if (fSingularity && fH == -1) {
101 SetPhase(fPhase-M_PI);
105 void AliFmHelix::SetCurvature(double val)
107 // Set helix curvature
111 SetPhase(fPhase+M_PI);
116 #ifndef ST_NO_NUMERIC_LIMITS
117 if (fabs(fCurvature) <= numeric_limits<double>::epsilon())
119 if (fabs(fCurvature) <= static_cast<double>(0))
121 fSingularity = true; // straight line
123 fSingularity = false; // curved
126 void AliFmHelix::SetPhase(double val)
130 fCosPhase = cos(fPhase);
131 fSinPhase = sin(fPhase);
132 if (fabs(fPhase) > M_PI)
133 fPhase = atan2(fSinPhase, fCosPhase); // force range [-pi,pi]
136 void AliFmHelix::SetDipAngle(double val)
138 // Set helix dip angle
140 fCosDipAngle = cos(fDipAngle);
141 fSinDipAngle = sin(fDipAngle);
144 double AliFmHelix::XCenter() const
146 // Set helix center in X
150 return fOrigin.x()-fCosPhase/fCurvature;
153 double AliFmHelix::YCenter() const
155 // Set helix center in Y
159 return fOrigin.y()-fSinPhase/fCurvature;
162 double AliFmHelix::FudgePathLength(const AliFmThreeVector<double>& p) const
166 double dx = p.x()-fOrigin.x();
167 double dy = p.y()-fOrigin.y();
170 s = (dy*fCosPhase - dx*fSinPhase)/fCosDipAngle;
173 s = atan2(dy*fCosPhase - dx*fSinPhase,
174 1/fCurvature + dx*fCosPhase+dy*fSinPhase)/
175 (fH*fCurvature*fCosDipAngle);
180 double AliFmHelix::Distance(const AliFmThreeVector<double>& p, bool scanPeriods) const
182 // calculate distance between origin an p along the helix
183 return abs(this->At(PathLength(p,scanPeriods))-p);
186 double AliFmHelix::PathLength(const AliFmThreeVector<double>& p, bool scanPeriods) const
189 // Returns the path length at the distance of closest
190 // approach between the helix and point p.
191 // For the case of B=0 (straight line) the path length
192 // can be calculated analytically. For B>0 there is
193 // unfortunately no easy solution to the problem.
194 // Here we use the Newton method to find the root of the
195 // referring equation. The 'fudgePathLength' serves
196 // as a starting value.
199 double dx = p.x()-fOrigin.x();
200 double dy = p.y()-fOrigin.y();
201 double dz = p.z()-fOrigin.z();
204 s = fCosDipAngle*(fCosPhase*dy-fSinPhase*dx) +
208 #ifndef ST_NO_NAMESPACES
210 using namespace units;
212 const double ktMaxPrecisionNeeded = micrometer;
213 const int ktMaxIterations = 100;
216 // The math is taken from Maple with C(expr,optimized) and
217 // some hand-editing. It is not very nice but efficient.
219 double t34 = fCurvature*fCosDipAngle*fCosDipAngle;
220 double t41 = fSinDipAngle*fSinDipAngle;
221 double t6, t7, t11, t12, t19;
224 // Get a first guess by using the dca in 2D. Since
225 // in some extreme cases we might be off by n periods
226 // we add (subtract) periods in case we get any closer.
228 s = FudgePathLength(p);
231 double ds = Period();
233 double d, dmin = abs(At(s) - p);
234 for(j=1; j<ktMaxIterations; j++) {
235 if ((d = abs(At(s+j*ds) - p)) < dmin) {
242 for(j=-1; -j<ktMaxIterations; j--) {
243 if ((d = abs(At(s+j*ds) - p)) < dmin) {
250 if (jmin) s += jmin*ds;
255 // Stops after ktMaxIterations iterations or if the required
256 // precision is obtained. Whatever comes first.
259 for (int i=0; i<ktMaxIterations; i++) {
260 t6 = fPhase+s*fH*fCurvature*fCosDipAngle;
262 t11 = dx-(1/fCurvature)*(t7-fCosPhase);
264 t19 = dy-(1/fCurvature)*(t12-fSinPhase);
265 s -= (t11*t12*fH*fCosDipAngle-t19*t7*fH*fCosDipAngle -
266 (dz-s*fSinDipAngle)*fSinDipAngle)/
267 (t12*t12*fCosDipAngle*fCosDipAngle+t11*t7*t34 +
268 t7*t7*fCosDipAngle*fCosDipAngle +
270 if (fabs(sOld-s) < ktMaxPrecisionNeeded) break;
273 #ifndef ST_NO_NAMESPACES
280 double AliFmHelix::Period() const
284 #ifndef ST_NO_NUMERIC_LIMITS
285 return numeric_limits<double>::max();
290 return fabs(2*M_PI/(fH*fCurvature*fCosDipAngle));
293 pair<double, double> AliFmHelix::PathLength(double r) const
296 // The math is taken from Maple with C(expr,optimized) and
297 // some hand-editing. It is not very nice but efficient.
298 // 'first' is the smallest of the two solutions (may be negative)
299 // 'second' is the other.
301 pair<double,double> tvalue;
302 pair<double,double> tVALUE(999999999.,999999999.);
304 double t1 = fCosDipAngle*(fOrigin.x()*fSinPhase-fOrigin.y()*fCosPhase);
305 double t12 = fOrigin.y()*fOrigin.y();
306 double t13 = fCosPhase*fCosPhase;
308 double t16 = fOrigin.x()*fOrigin.x();
309 double t20 = -fCosDipAngle*fCosDipAngle*(2.0*fOrigin.x()*fSinPhase*fOrigin.y()*fCosPhase +
310 t12-t12*t13-t15+t13*t16);
311 if (t20<0.) return tVALUE;
313 tvalue.first = (t1-t20)/(fCosDipAngle*fCosDipAngle);
314 tvalue.second = (t1+t20)/(fCosDipAngle*fCosDipAngle);
317 double t1 = fOrigin.y()*fCurvature;
318 double t2 = fSinPhase;
319 double t3 = fCurvature*fCurvature;
320 double t4 = fOrigin.y()*t2;
321 double t5 = fCosPhase;
322 double t6 = fOrigin.x()*t5;
323 double t8 = fOrigin.x()*fOrigin.x();
324 double t11 = fOrigin.y()*fOrigin.y();
326 double t15 = t14*fCurvature;
328 double t19 = t11*t11;
331 double t32 = t14*t14;
333 double t38 = 8.0*t4*t6 - 4.0*t1*t2*t8 - 4.0*t11*fCurvature*t6 +
334 4.0*t15*t6 + t17*t3 + t19*t3 + 2.0*t21*t8 + 4.0*t8*t23 -
335 4.0*t8*fOrigin.x()*fCurvature*t5 - 4.0*t11*t23 -
336 4.0*t11*fOrigin.y()*fCurvature*t2 + 4.0*t11 - 4.0*t14 +
337 t32*t3 + 4.0*t15*t4 - 2.0*t35*t11 - 2.0*t35*t8;
338 double t40 = (-t3*t38);
339 if (t40<0.) return tVALUE;
342 double t43 = fOrigin.x()*fCurvature;
343 double t45 = 2.0*t5 - t35 + t21 + 2.0 - 2.0*t1*t2 -2.0*t43 - 2.0*t43*t5 + t8*t3;
344 double t46 = fH*fCosDipAngle*fCurvature;
346 tvalue.first = (-fPhase + 2.0*atan((-2.0*t1 + 2.0*t2 + t40)/t45))/t46;
347 tvalue.second = -(fPhase + 2.0*atan((2.0*t1 - 2.0*t2 + t40)/t45))/t46;
350 // Solution can be off by +/- one Period, select smallest
353 // malisa deletes "isnan" check 22apr2006 if (!isnan(value.first)) {
354 if (fabs(tvalue.first-p) < fabs(tvalue.first)) tvalue.first = tvalue.first-p;
355 else if (fabs(tvalue.first+p) < fabs(tvalue.first)) tvalue.first = tvalue.first+p;
357 // malisa deletes "isnan" check 22apr2006 if (!isnan(tvalue.second)) {
358 if (fabs(tvalue.second-p) < fabs(tvalue.second)) tvalue.second = tvalue.second-p;
359 else if (fabs(tvalue.second+p) < fabs(tvalue.second)) tvalue.second = tvalue.second+p;
362 if (tvalue.first > tvalue.second)
363 swap(tvalue.first,tvalue.second);
367 pair<double, double> AliFmHelix::PathLength(double r, double x, double y, bool scanPeriods)
370 double x0 = fOrigin.x();
371 double y0 = fOrigin.y();
374 pair<double, double> result = this->PathLength(r);
380 double AliFmHelix::PathLength(const AliFmThreeVector<double>& r,
381 const AliFmThreeVector<double>& n) const
384 // Vector 'r' defines the position of the center and
385 // vector 'n' the normal vector of the plane.
386 // For a straight line there is a simple analytical
387 // solution. For curvatures > 0 the root is determined
388 // by Newton method. In case no valid s can be found
389 // the max. largest value for s is returned.
394 double t = n.z()*fSinDipAngle +
395 n.y()*fCosDipAngle*fCosPhase -
396 n.x()*fCosDipAngle*fSinPhase;
400 s = ((r - fOrigin)*n)/t;
403 const double ktMaxPrecisionNeeded = micrometer;
404 const int ktMaxIterations = 20;
406 double tA = fCurvature*((fOrigin - r)*n) -
409 double t = fH*fCurvature*fCosDipAngle;
410 double u = n.z()*fCurvature*fSinDipAngle;
416 // (cos(kangMax)-1)/kangMax = 0.1
417 const double kangMax = 0.21;
418 double deltas = fabs(kangMax/(fCurvature*fCosDipAngle));
419 // dampingFactor = exp(-0.5);
420 double dampingFactor = 0.60653;
423 for (i=0; i<ktMaxIterations; i++) {
425 double sina = sin(a);
426 double cosa = cos(a);
434 if ( fabs(fp)*deltas <= fabs(f) ) { //too big step
436 if (fp<0.) sgn = -sgn;
437 if (f <0.) sgn = -sgn;
439 if (shift == -shiftOld) { // don't get stuck shifting +/-deltas
440 deltas *= dampingFactor; // dampen magnitude of shift
442 // allow iterations to run out
444 i--; // don't count against iterations
451 if (fabs(sOld-s) < ktMaxPrecisionNeeded) break;
454 if (i == ktMaxIterations) return fgkNoSolution;
460 AliFmHelix::PathLengths(const AliFmHelix& h, bool scanPeriods) const
463 // Cannot handle case where one is a helix
464 // and the other one is a straight line.
466 if (fSingularity != h.fSingularity)
467 return pair<double, double>(fgkNoSolution, fgkNoSolution);
475 AliFmThreeVector<double> dv = h.fOrigin - fOrigin;
476 AliFmThreeVector<double> a(-fCosDipAngle*fSinPhase,
477 fCosDipAngle*fCosPhase,
479 AliFmThreeVector<double> b(-h.fCosDipAngle*h.fSinPhase,
480 h.fCosDipAngle*h.fCosPhase,
485 s2 = (k-ab*g)/(ab*ab-1.);
487 return pair<double, double>(s1, s2);
491 // First step: get dca in the xy-plane as start value
493 double dx = h.XCenter() - XCenter();
494 double dy = h.YCenter() - YCenter();
495 double dd = ::sqrt(dx*dx + dy*dy);
496 double r1 = 1/Curvature();
497 double r2 = 1/h.Curvature();
499 /* malisa 22apr2006 is commenting out the "isnan" check
500 * if ( !finite(r1) || isnan(r1) || !finite(r2) || isnan(r2) ) {
501 * cerr << __FUNCTION__ << " *** error *** ";
502 * cerr << " r1=" << r1;
503 * cerr << " r2=" << r2 << endl;
504 * return pair<double, double>(fgkNoSolution, fgkNoSolution);
508 double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd);
512 if (fabs(cosAlpha) < 1) { // two solutions
513 double sinAlpha = sin(acos(cosAlpha));
514 x = XCenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd;
515 y = YCenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd;
516 s = PathLength(x, y);
517 x = XCenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd;
518 y = YCenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd;
519 double a = PathLength(x, y);
520 if (h.Distance(At(a),scanPeriods) < h.Distance(At(s),scanPeriods)) s = a;
522 else { // no intersection (or exactly one)
523 x = XCenter() + r1*dx/dd;
524 y = YCenter() + r1*dy/dd;
525 s = PathLength(x, y);
529 // Second step: scan in decreasing intervals around seed 's'
531 const double ktMinStepSize = 10*micrometer;
532 const double ktMinRange = 10*centimeter;
533 double dmin = h.Distance(At(s),scanPeriods);
534 double range = max(2*dmin, ktMinRange);
536 /* malisa comments out the "isnan" check 22apr2006
537 * if ( !finite(range) || isnan(range)) {
538 * cerr << __FUNCTION__ << " *** error *** ";
539 * cerr << " range=" << range << endl;
540 * return pair<double, double>(fgkNoSolution, fgkNoSolution);
544 double ds = range/10.;
545 double slast=-999999, d;
552 while (ds > ktMinStepSize) {
554 cerr << __FUNCTION__ << " *** error *** s1 = ";
555 cerr << " s2 = " << s2;
556 cerr << " range = " << range;
557 cerr << " ds = " << ds << endl;
558 return pair<double, double>(fgkNoSolution, fgkNoSolution);
564 cerr << __FUNCTION__ << " *** error *** ds = " << ds << endl;
565 return pair<double, double>(fgkNoSolution, fgkNoSolution);
568 for (double ss=s1; ss<(s2+ds); ss+=ds) {
570 if ( iterations > 100 ) {
571 cerr << __FUNCTION__ << " *** error *** iterations = " << iterations << endl;
572 return pair<double, double>(fgkNoSolution, fgkNoSolution);
574 d = h.Distance(At(ss),scanPeriods);
582 // In the rare cases where the minimum is at the
583 // the border of the current range we shift the range
584 // and start all over, i.e we do not decrease 'ds'.
585 // Else we decrease the search intervall around the
586 // current minimum and redo the scan in smaller steps.
593 else if (s == slast) {
604 return pair<double, double>(s, h.PathLength(At(s),scanPeriods));
609 void AliFmHelix::MoveOrigin(double s)
611 // Move the helix origin
615 AliFmThreeVector<double> newOrigin = At(s);
616 double newPhase = atan2(newOrigin.y() - YCenter(),
617 newOrigin.x() - XCenter());
623 int operator== (const AliFmHelix& a, const AliFmHelix& b)
626 // Checks for numerical identity only !
628 return (a.Origin() == b.Origin() &&
629 a.DipAngle() == b.DipAngle() &&
630 a.Curvature() == b.Curvature() &&
631 a.Phase() == b.Phase() &&
635 int operator!= (const AliFmHelix& a, const AliFmHelix& b) {return !(a == b);}
637 ostream& operator<<(ostream& os, const AliFmHelix& h)
640 << "curvature = " << h.Curvature() << ", "
641 << "dip angle = " << h.DipAngle() << ", "
642 << "phase = " << h.Phase() << ", "
643 << "h = " << h.H() << ", "
644 << "origin = " << h.Origin() << ')';