PWG2/FEMTOSCOPY/AliFemto/Infrastructure/AliFmHelix.cxx

   1 /***************************************************************************
   2 *
   3 * $Id$
   4 *
   5 * Author: Thomas Ullrich, Sep 1997
   6 ***************************************************************************
   7 *
   8 * Description: Parametrization of a helix
   9 *
  10 ***************************************************************************
  11 *
  12 * $Log$
  13 * Revision 1.3  2007/04/27 07:24:34  akisiel
  14 * Make revisions needed for compilation from the main AliRoot tree
  15 *
  16 * Revision 1.1.1.1  2007/04/25 15:38:41  panos
  17 * Importing the HBT code dir
  18 *
  19 * Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
  20 * First version on CVS
  21 *
  22 * Revision 1.26  2005/10/13 23:15:13  genevb
  23 * Save a few calculations
  24 *
  25 * Revision 1.25  2005/10/13 22:25:35  genevb
  26 * pathLength to plane now finds nearest approach to intersection regardless of # of loops
  27 *
  28 * Revision 1.24  2005/07/06 18:49:56  fisyak
  29 * Replace AliFmHelixD, AliFmLorentzVectorD,AliFmLorentzVectorF,AliFmMatrixD,AliFmMatrixF,AliFmPhysicalHelixD,AliFmThreeVectorD,AliFmThreeVectorF by templated version
  30 *
  31 * Revision 1.23  2004/12/02 02:51:16  ullrich
  32 * Added option to pathLenghth() and distance() to search for
  33 * DCA only within one period. Default stays as it was.
  34 *
  35 * Revision 1.22  2004/05/03 23:35:31  perev
  36 * Possible non init WarnOff
  37 *
  38 * Revision 1.21  2004/01/27 02:49:48  perev
  39 * Big value appropriate for float
  40 *
  41 * Revision 1.20  2003/12/22 18:59:36  perev
  42 * remove test for only +ve pathLeng added before
  43 *
  44 * Revision 1.19  2003/12/18 17:27:02  perev
  45 * Small bug fix in number of iters. i++ ==> i--
  46 *
  47 * Revision 1.18  2003/10/30 20:06:46  perev
  48 * Check of quality added
  49 *
  50 * Revision 1.17  2003/10/19 20:17:00  perev
  51 * Protection agains overfloat added into pathLength(AliFmThreeVector,AliFmThreeVector)
  52 *
  53 * Revision 1.16  2003/10/06 23:39:21  perev
  54 * sqrt(-ve) == no solution. infinity returns
  55 *
  56 * Revision 1.15  2003/09/02 17:59:34  perev
  57 * gcc 3.2 updates + WarnOff
  58 *
  59 * Revision 1.14  2003/06/26 17:15:56  ullrich
  60 * Changed local variable name in pathLenght.
  61 *
  62 * Revision 1.13  2002/06/21 17:49:25  genevb
  63 * Some minor speed improvements
  64 *
  65 * Revision 1.12  2002/04/24 02:40:25  ullrich
  66 * Restored old format and lost CVS statements.
  67 *
  68 * Revision 1.11  2002/02/12 19:37:51  jeromel
  69 * fix for Linux 7.2 (float.h). Took oportunity to doxygenize.
  70 *
  71 * Revision 1.10  2000/07/17 21:44:19  ullrich
  72 * Fixed problem in pathLength(cyl_rad).
  73 *
  74 * Revision 1.9  2000/05/22 21:38:28  ullrich
  75 * Add parenthesis to make Linux compiler happy.
  76 *
  77 * Revision 1.8  2000/05/22 21:11:21  ullrich
  78 * In pathLength(AliFmThreeVector&): Increased number of max iteration
  79 * in Newton method from 10 to 100. Improved initial guess in case
  80 * it is off by n period.
  81 *
  82 * Revision 1.7  2000/03/06 20:24:25  ullrich
  83 * Parameter h for case B=0 correctly handled now.
  84 *
  85 * Revision 1.6  1999/12/22 15:14:39  ullrich
  86 * Added analytical solution for dca between two helices
  87 * in the case for B=0.
  88 *
  89 * Revision 1.5  1999/12/21 15:14:08  ullrich
  90 * Modified to cope with new compiler version on Sun (CC5.0).
  91 *
  92 * Revision 1.4  1999/11/29 21:45:38  fisyak
  93 * fix abs for HP
  94 *
  95 * Revision 1.3  1999/03/07 14:55:41  wenaus
  96 * fix scope problem
  97 *
  98 * Revision 1.2  1999/03/02 19:47:35  ullrich
  99 * Added method to find dca between two helices
 100 *
 101 * Revision 1.1  1999/01/30 03:59:02  fisyak
 102 * Root Version of AliFmarClassLibrary
 103 *
 104 * Revision 1.1  1999/01/23 00:29:15  ullrich
 105 * Initial Revision
 106 *
 107 **************************************************************************/
 108 #if !defined(ST_NO_NUMERIC_LIMITS)
 109 #    include <limits>
 110 #    if !defined(ST_NO_NAMESPACES)
 111 using std::numeric_limits;
 112 #    endif
 113 #endif
 114 #define FOR_HELIX
 115 #include <float.h>
 116 #include <assert.h>
 117
 118 #include "AliFmHelix.h"
 119 #include "Base/PhysicalConstants.h"
 120 #include "Base/SystemOfUnits.h"
 121 #ifdef __ROOT__
 122 ClassImpT(AliFmHelix,double);
 123 #endif
 124
 125 #ifdef WIN32
 126 #include "gcc2vs.h"
 127 #endif
 128
 129 #include <iostream>
 130 #include <fstream>
 131 using namespace std;
 132
 133 const double AliFmHelix::NoSolution = 3.e+33;
 134
 135 AliFmHelix::AliFmHelix() :
 136   mSingularity(0),
 137   mOrigin(0,0,0),
 138   mDipAngle(0),
 139   mCurvature(0),
 140   mPhase(0),
 141   mH(0),
 142   mCosDipAngle(0),
 143   mSinDipAngle(0),
 144   mCosPhase(0),
 145   mSinPhase(0)
 146 { /*noop*/ }
 147
 148 AliFmHelix::AliFmHelix(double c, double d, double phase,
 149                        const AliFmThreeVector<double>& o, int h) :
 150   mSingularity(0),
 151   mOrigin(0,0,0),
 152   mDipAngle(0),
 153   mCurvature(0),
 154   mPhase(0),
 155   mH(0),
 156   mCosDipAngle(0),
 157   mSinDipAngle(0),
 158   mCosPhase(0),
 159   mSinPhase(0)
 160 {
 161         setParameters(c, d, phase, o, h);
 162 }
 163
 164 AliFmHelix::~AliFmHelix() { /* noop */ };
 165
 166 void AliFmHelix::setParameters(double c, double dip, double phase,
 167                                                         const AliFmThreeVector<double>& o, int h)
 168 {
 169         //
 170         //  The order in which the parameters are set is important
 171         //  since setCurvature might have to adjust the others.
 172         //
 173         mH = (h>=0) ? 1 : -1;    // Default is: positive particle
 174         //             positive field
 175         mOrigin   = o;
 176         setDipAngle(dip);
 177         setPhase(phase);
 178
 179         //
 180         // Check for singularity and correct for negative curvature.
 181         // May change mH and mPhase. Must therefore be set last.
 182         //
 183         setCurvature(c);
 184
 185         //
 186         // For the case B=0, h is ill defined. In the following we
 187         // always assume h = +1. Since phase = psi - h * pi/2
 188         // we have to correct the phase in case h = -1.
 189         // This assumes that the user uses the same h for phase
 190         // as the one he passed to the constructor.
 191         //
 192         if (mSingularity && mH == -1) {
 193                 mH = +1;
 194                 setPhase(mPhase-M_PI);
 195         }
 196 }
 197
 198 void AliFmHelix::setCurvature(double val)
 199 {
 200         if (val < 0) {
 201                 mCurvature = -val;
 202                 mH = -mH;
 203                 setPhase(mPhase+M_PI);
 204         }
 205         else
 206                 mCurvature = val;
 207
 208 #ifndef ST_NO_NUMERIC_LIMITS
 209         if (fabs(mCurvature) <= numeric_limits<double>::epsilon())
 210 #else
 211         if (fabs(mCurvature) <= static_cast<double>(0))
 212 #endif
 213                 mSingularity = true;                    // straight line
 214         else
 215                 mSingularity = false;                   // curved
 216 }
 217
 218 void AliFmHelix::setPhase(double val)
 219 {
 220         mPhase       = val;
 221         mCosPhase    = cos(mPhase);
 222         mSinPhase    = sin(mPhase);
 223         if (fabs(mPhase) > M_PI)
 224                 mPhase = atan2(mSinPhase, mCosPhase);  // force range [-pi,pi]
 225 }
 226
 227 void AliFmHelix::setDipAngle(double val)
 228 {
 229         mDipAngle    = val;
 230         mCosDipAngle = cos(mDipAngle);
 231         mSinDipAngle = sin(mDipAngle);
 232 }
 233
 234 double AliFmHelix::xcenter() const
 235 {
 236         if (mSingularity)
 237                 return 0;
 238         else
 239                 return mOrigin.x()-mCosPhase/mCurvature;
 240 }
 241
 242 double AliFmHelix::ycenter() const
 243 {
 244         if (mSingularity)
 245                 return 0;
 246         else
 247                 return mOrigin.y()-mSinPhase/mCurvature;
 248 }
 249
 250 double AliFmHelix::fudgePathLength(const AliFmThreeVector<double>& p) const
 251 {
 252         double s;
 253         double dx = p.x()-mOrigin.x();
 254         double dy = p.y()-mOrigin.y();
 255
 256         if (mSingularity) {
 257                 s = (dy*mCosPhase - dx*mSinPhase)/mCosDipAngle;
 258         }
 259         else {
 260                 s = atan2(dy*mCosPhase - dx*mSinPhase,
 261                         1/mCurvature + dx*mCosPhase+dy*mSinPhase)/
 262                         (mH*mCurvature*mCosDipAngle);
 263         }
 264         return s;
 265 }
 266
 267 double AliFmHelix::distance(const AliFmThreeVector<double>& p, bool scanPeriods) const
 268 {
 269         return abs(this->at(pathLength(p,scanPeriods))-p);
 270 }
 271
 272 double AliFmHelix::pathLength(const AliFmThreeVector<double>& p, bool scanPeriods) const
 273 {
 274         //
 275         //  Returns the path length at the distance of closest
 276         //  approach between the helix and point p.
 277         //  For the case of B=0 (straight line) the path length
 278         //  can be calculated analytically. For B>0 there is
 279         //  unfortunately no easy solution to the problem.
 280         //  Here we use the Newton method to find the root of the
 281         //  referring equation. The 'fudgePathLength' serves
 282         //  as a starting value.
 283         //
 284         double s;
 285         double dx = p.x()-mOrigin.x();
 286         double dy = p.y()-mOrigin.y();
 287         double dz = p.z()-mOrigin.z();
 288
 289         if (mSingularity) {
 290                 s = mCosDipAngle*(mCosPhase*dy-mSinPhase*dx) +
 291                         mSinDipAngle*dz;
 292         }
 293         else { //
 294 #ifndef ST_NO_NAMESPACES
 295                 {
 296                         using namespace units;
 297 #endif
 298                         const double MaxPrecisionNeeded = micrometer;
 299                         const int    MaxIterations      = 100;
 300
 301                         //
 302                         // The math is taken from Maple with C(expr,optimized) and
 303                         // some hand-editing. It is not very nice but efficient.
 304                         //
 305                         double t34 = mCurvature*mCosDipAngle*mCosDipAngle;
 306                         double t41 = mSinDipAngle*mSinDipAngle;
 307                         double t6, t7, t11, t12, t19;
 308
 309                         //
 310                         // Get a first guess by using the dca in 2D. Since
 311                         // in some extreme cases we might be off by n periods
 312                         // we add (subtract) periods in case we get any closer.
 313                         //
 314                         s = fudgePathLength(p);
 315
 316                         if (scanPeriods) {
 317                                 double ds = period();
 318                                 int    j, jmin = 0;
 319                                 double d, dmin = abs(at(s) - p);
 320                                 for(j=1; j<MaxIterations; j++) {
 321                                         if ((d = abs(at(s+j*ds) - p)) < dmin) {
 322                                                 dmin = d;
 323                                                 jmin = j;
 324                                         }
 325                                         else
 326                                                 break;
 327                                 }
 328                                 for(j=-1; -j<MaxIterations; j--) {
 329                                         if ((d = abs(at(s+j*ds) - p)) < dmin) {
 330                                                 dmin = d;
 331                                                 jmin = j;
 332                                         }
 333                                         else
 334                                                 break;
 335                                 }
 336                                 if (jmin) s += jmin*ds;
 337                         }
 338
 339                         //
 340                         // Newtons method:
 341                         // Stops after MaxIterations iterations or if the required
 342                         // precision is obtained. Whatever comes first.
 343                         //
 344                         double sOld = s;
 345                         for (int i=0; i<MaxIterations; i++) {
 346                                 t6  = mPhase+s*mH*mCurvature*mCosDipAngle;
 347                                 t7  = cos(t6);
 348                                 t11 = dx-(1/mCurvature)*(t7-mCosPhase);
 349                                 t12 = sin(t6);
 350                                 t19 = dy-(1/mCurvature)*(t12-mSinPhase);
 351                                 s  -= (t11*t12*mH*mCosDipAngle-t19*t7*mH*mCosDipAngle -
 352                                         (dz-s*mSinDipAngle)*mSinDipAngle)/
 353                                         (t12*t12*mCosDipAngle*mCosDipAngle+t11*t7*t34 +
 354                                         t7*t7*mCosDipAngle*mCosDipAngle +
 355                                         t19*t12*t34+t41);
 356                                 if (fabs(sOld-s) < MaxPrecisionNeeded) break;
 357                                 sOld = s;
 358                         }
 359 #ifndef ST_NO_NAMESPACES
 360                 }
 361 #endif
 362         }
 363         return s;
 364 }
 365
 366 double AliFmHelix::period() const
 367 {
 368         if (mSingularity)
 369 #ifndef ST_NO_NUMERIC_LIMITS
 370                 return numeric_limits<double>::max();
 371 #else
 372                 return DBL_MAX;
 373 #endif
 374         else
 375                 return fabs(2*M_PI/(mH*mCurvature*mCosDipAngle));
 376 }
 377
 378 pair<double, double> AliFmHelix::pathLength(double r) const
 379 {
 380         pair<double,double> value;
 381         pair<double,double> VALUE(999999999.,999999999.);
 382         //
 383         // The math is taken from Maple with C(expr,optimized) and
 384         // some hand-editing. It is not very nice but efficient.
 385         // 'first' is the smallest of the two solutions (may be negative)
 386         // 'second' is the other.
 387         //
 388         if (mSingularity) {
 389                 double t1 = mCosDipAngle*(mOrigin.x()*mSinPhase-mOrigin.y()*mCosPhase);
 390                 double t12 = mOrigin.y()*mOrigin.y();
 391                 double t13 = mCosPhase*mCosPhase;
 392                 double t15 = r*r;
 393                 double t16 = mOrigin.x()*mOrigin.x();
 394                 double t20 = -mCosDipAngle*mCosDipAngle*(2.0*mOrigin.x()*mSinPhase*mOrigin.y()*mCosPhase +
 395                         t12-t12*t13-t15+t13*t16);
 396                 if (t20<0.) return VALUE;
 397                 t20 = ::sqrt(t20);
 398                 value.first  = (t1-t20)/(mCosDipAngle*mCosDipAngle);
 399                 value.second = (t1+t20)/(mCosDipAngle*mCosDipAngle);
 400         }
 401         else {
 402                 double t1 = mOrigin.y()*mCurvature;
 403                 double t2 = mSinPhase;
 404                 double t3 = mCurvature*mCurvature;
 405                 double t4 = mOrigin.y()*t2;
 406                 double t5 = mCosPhase;
 407                 double t6 = mOrigin.x()*t5;
 408                 double t8 = mOrigin.x()*mOrigin.x();
 409                 double t11 = mOrigin.y()*mOrigin.y();
 410                 double t14 = r*r;
 411                 double t15 = t14*mCurvature;
 412                 double t17 = t8*t8;
 413                 double t19 = t11*t11;
 414                 double t21 = t11*t3;
 415                 double t23 = t5*t5;
 416                 double t32 = t14*t14;
 417                 double t35 = t14*t3;
 418                 double t38 = 8.0*t4*t6 - 4.0*t1*t2*t8 - 4.0*t11*mCurvature*t6 +
 419                         4.0*t15*t6 + t17*t3 + t19*t3 + 2.0*t21*t8 + 4.0*t8*t23 -
 420                         4.0*t8*mOrigin.x()*mCurvature*t5 - 4.0*t11*t23 -
 421                         4.0*t11*mOrigin.y()*mCurvature*t2 + 4.0*t11 - 4.0*t14 +
 422                         t32*t3 + 4.0*t15*t4 - 2.0*t35*t11 - 2.0*t35*t8;
 423                 double t40 = (-t3*t38);
 424                 if (t40<0.) return VALUE;
 425                 t40 = ::sqrt(t40);
 426
 427                 double t43 = mOrigin.x()*mCurvature;
 428                 double t45 = 2.0*t5 - t35 + t21 + 2.0 - 2.0*t1*t2 -2.0*t43 - 2.0*t43*t5 + t8*t3;
 429                 double t46 = mH*mCosDipAngle*mCurvature;
 430
 431                 value.first = (-mPhase + 2.0*atan((-2.0*t1 + 2.0*t2 + t40)/t45))/t46;
 432                 value.second = -(mPhase + 2.0*atan((2.0*t1 - 2.0*t2 + t40)/t45))/t46;
 433
 434                 //
 435                 //   Solution can be off by +/- one period, select smallest
 436                 //
 437                 double p = period();
 438                 //      malisa deletes "isnan" check 22apr2006  if (!isnan(value.first)) {
 439                 if (fabs(value.first-p) < fabs(value.first)) value.first = value.first-p;
 440                 else if (fabs(value.first+p) < fabs(value.first)) value.first = value.first+p;
 441                 //      malisa  }
 442                 //      malisa deletes "isnan" check 22apr2006          if (!isnan(value.second)) {
 443                 if (fabs(value.second-p) < fabs(value.second)) value.second = value.second-p;
 444                 else if (fabs(value.second+p) < fabs(value.second)) value.second = value.second+p;
 445                 //      malisa }
 446         }
 447         if (value.first > value.second)
 448                 swap(value.first,value.second);
 449         return(value);
 450 }
 451
 452 pair<double, double> AliFmHelix::pathLength(double r, double x, double y, bool scanPeriods)
 453 {
 454         double x0 = mOrigin.x();
 455         double y0 = mOrigin.y();
 456         mOrigin.setX(x0-x);
 457         mOrigin.setY(y0-y);
 458         pair<double, double> result = this->pathLength(r);
 459         mOrigin.setX(x0);
 460         mOrigin.setY(y0);
 461         return result;
 462 }
 463
 464 double AliFmHelix::pathLength(const AliFmThreeVector<double>& r,
 465                                                    const AliFmThreeVector<double>& n) const
 466 {
 467         //
 468         // Vector 'r' defines the position of the center and
 469         // vector 'n' the normal vector of the plane.
 470         // For a straight line there is a simple analytical
 471         // solution. For curvatures > 0 the root is determined
 472         // by Newton method. In case no valid s can be found
 473         // the max. largest value for s is returned.
 474         //
 475         double s;
 476
 477         if (mSingularity) {
 478                 double t = n.z()*mSinDipAngle +
 479                         n.y()*mCosDipAngle*mCosPhase -
 480                         n.x()*mCosDipAngle*mSinPhase;
 481                 if (t == 0)
 482                         s = NoSolution;
 483                 else
 484                         s = ((r - mOrigin)*n)/t;
 485         }
 486         else {
 487                 const double MaxPrecisionNeeded = micrometer;
 488                 const int    MaxIterations      = 20;
 489
 490                 double A = mCurvature*((mOrigin - r)*n) -
 491                         n.x()*mCosPhase -
 492                         n.y()*mSinPhase;
 493                 double t = mH*mCurvature*mCosDipAngle;
 494                 double u = n.z()*mCurvature*mSinDipAngle;
 495
 496                 double a, f, fp;
 497                 double sOld = s = 0;
 498                 double shiftOld = 0;
 499                 double shift;
 500                 //              (cos(angMax)-1)/angMax = 0.1
 501                 const double angMax = 0.21;
 502                 double deltas = fabs(angMax/(mCurvature*mCosDipAngle));
 503                 //              dampingFactor = exp(-0.5);
 504                 double dampingFactor = 0.60653;
 505                 int i;
 506
 507                 for (i=0; i<MaxIterations; i++) {
 508                         a  = t*s+mPhase;
 509                         double sina = sin(a);
 510                         double cosa = cos(a);
 511                         f  = A +
 512                                 n.x()*cosa +
 513                                 n.y()*sina +
 514                                 u*s;
 515                         fp = -n.x()*sina*t +
 516                                 n.y()*cosa*t +
 517                                 u;
 518                         if ( fabs(fp)*deltas <= fabs(f) ) { //too big step
 519                                 int sgn = 1;
 520                                 if (fp<0.) sgn = -sgn;
 521                                 if (f <0.) sgn = -sgn;
 522                                 shift = sgn*deltas;
 523                                 if (shift == -shiftOld) { // don't get stuck shifting +/-deltas
 524                                         deltas *= dampingFactor; // dampen magnitude of shift
 525                                         shift = sgn*deltas;
 526                                         // allow iterations to run out
 527                                 } else {
 528                                         i--; // don't count against iterations
 529                                 }
 530                         } else {
 531                                 shift = f/fp;
 532                         }
 533                         s -= shift;
 534                         shiftOld = shift;
 535                         if (fabs(sOld-s) < MaxPrecisionNeeded) break;
 536                         sOld = s;
 537                 }
 538                 if (i == MaxIterations) return NoSolution;
 539         }
 540         return s;
 541 }
 542
 543 pair<double, double>
 544 AliFmHelix::pathLengths(const AliFmHelix& h, bool scanPeriods) const
 545 {
 546
 547         //
 548         //      Cannot handle case where one is a helix
 549         //  and the other one is a straight line.
 550         //
 551         if (mSingularity != h.mSingularity)
 552                 return pair<double, double>(NoSolution, NoSolution);
 553
 554         double s1, s2;
 555
 556         if (mSingularity) {
 557                 //
 558                 //  Analytic solution
 559                 //
 560                 AliFmThreeVector<double> dv = h.mOrigin - mOrigin;
 561                 AliFmThreeVector<double> a(-mCosDipAngle*mSinPhase,
 562                         mCosDipAngle*mCosPhase,
 563                         mSinDipAngle);
 564                 AliFmThreeVector<double> b(-h.mCosDipAngle*h.mSinPhase,
 565                         h.mCosDipAngle*h.mCosPhase,
 566                         h.mSinDipAngle);
 567                 double ab = a*b;
 568                 double g  = dv*a;
 569                 double k  = dv*b;
 570                 s2 = (k-ab*g)/(ab*ab-1.);
 571                 s1 = g+s2*ab;
 572                 return pair<double, double>(s1, s2);
 573         }
 574         else {
 575                 //
 576                 //  First step: get dca in the xy-plane as start value
 577                 //
 578                 double dx = h.xcenter() - xcenter();
 579                 double dy = h.ycenter() - ycenter();
 580                 double dd = ::sqrt(dx*dx + dy*dy);
 581                 double r1 = 1/curvature();
 582                 double r2 = 1/h.curvature();
 583
 584                 /* malisa 22apr2006 is commenting out the "isnan" check
 585                  *              if ( !finite(r1) || isnan(r1) || !finite(r2) || isnan(r2) ) {
 586                  *               cerr << __FUNCTION__ << " *** error *** ";
 587                  *                      cerr << "  r1=" << r1;
 588                  *                      cerr << "  r2=" << r2 << endl;
 589                  *                      return pair<double, double>(NoSolution, NoSolution);
 590                  *              }
 591                  */
 592
 593                 double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd);
 594
 595                 double s;
 596                 double x, y;
 597                 if (fabs(cosAlpha) < 1) {           // two solutions
 598                         double sinAlpha = sin(acos(cosAlpha));
 599                         x = xcenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd;
 600                         y = ycenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd;
 601                         s = pathLength(x, y);
 602                         x = xcenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd;
 603                         y = ycenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd;
 604                         double a = pathLength(x, y);
 605                         if (h.distance(at(a),scanPeriods) < h.distance(at(s),scanPeriods)) s = a;
 606                 }
 607                 else {                              // no intersection (or exactly one)
 608                         x = xcenter() + r1*dx/dd;
 609                         y = ycenter() + r1*dy/dd;
 610                         s = pathLength(x, y);
 611                 }
 612
 613                 //
 614                 //   Second step: scan in decreasing intervals around seed 's'
 615                 //
 616                 const double MinStepSize = 10*micrometer;
 617                 const double MinRange    = 10*centimeter;
 618                 double dmin              = h.distance(at(s),scanPeriods);
 619                 double range             = max(2*dmin, MinRange);
 620
 621                 /* malisa comments out the "isnan" check 22apr2006
 622                  *              if ( !finite(range) || isnan(range)) {
 623                  *                      cerr << __FUNCTION__ << " *** error *** ";
 624                  *                      cerr << "  range=" << range << endl;
 625                  *                      return pair<double, double>(NoSolution, NoSolution);
 626                  *              }
 627                  */
 628
 629                 double ds = range/10.;
 630                 double slast=-999999, d;
 631                 s1 = s - range/2.;
 632                 s2 = s + range/2.;
 633
 634                 double tmp1;
 635                 double tmp2;
 636
 637                 while (ds > MinStepSize) {
 638                         if ( !(s1<s2) ) {
 639                                 cerr << __FUNCTION__ << " *** error *** s1 = ";
 640                                 cerr << "    s2 = " << s2;
 641                                 cerr << "    range = " << range;
 642                                 cerr << "     ds = " << ds << endl;
 643                                 return pair<double, double>(NoSolution, NoSolution);
 644                         }
 645                         tmp1 = (s2-s1);
 646                         tmp2 = tmp1/10.0;
 647                         ds = tmp2;
 648                         if ( ds<0) {
 649                             cerr << __FUNCTION__ << " *** error *** ds = " << ds << endl;
 650                             return pair<double, double>(NoSolution, NoSolution);
 651                         }
 652                         int iterations = 0;
 653                         for (double ss=s1; ss<(s2+ds); ss+=ds) {
 654                             iterations++;
 655                             if ( iterations > 100 ) {
 656                                 cerr << __FUNCTION__ << " *** error *** iterations = " << iterations << endl;
 657                                 return pair<double, double>(NoSolution, NoSolution);
 658                             }
 659                             d = h.distance(at(ss),scanPeriods);
 660                             if (d < dmin) {
 661                                 dmin = d;
 662                                 s = ss;
 663                             }
 664                             slast = ss;
 665                         }
 666                         //
 667                         //  In the rare cases where the minimum is at the
 668                         //  the border of the current range we shift the range
 669                         //  and start all over, i.e we do not decrease 'ds'.
 670                         //  Else we decrease the search intervall around the
 671                         //  current minimum and redo the scan in smaller steps.
 672                         //
 673                         if (s == s1) {
 674                                 d = 0.8*(s2-s1);
 675                                 s1 -= d;
 676                                 s2 -= d;
 677                         }
 678                         else if (s == slast) {
 679                                 d = 0.8*(s2-s1);
 680                                 s1 += d;
 681                                 s2 += d;
 682                         }
 683                         else {
 684                                 s1 = s-ds;
 685                                 s2 = s+ds;
 686                         //      ds /= 10;
 687                         }
 688                 }
 689                 return pair<double, double>(s, h.pathLength(at(s),scanPeriods));
 690         }
 691 }
 692
 693
 694 void AliFmHelix::moveOrigin(double s)
 695 {
 696         if (mSingularity)
 697                 mOrigin = at(s);
 698         else {
 699                 AliFmThreeVector<double> newOrigin = at(s);
 700                 double newPhase = atan2(newOrigin.y() - ycenter(),
 701                         newOrigin.x() - xcenter());
 702                 mOrigin = newOrigin;
 703                 setPhase(newPhase);
 704         }
 705 }
 706
 707 int operator== (const AliFmHelix& a, const AliFmHelix& b)
 708 {
 709         //
 710         // Checks for numerical identity only !
 711         //
 712         return (a.origin()    == b.origin()    &&
 713                 a.dipAngle()  == b.dipAngle()  &&
 714                 a.curvature() == b.curvature() &&
 715                 a.phase()     == b.phase()     &&
 716                 a.h()         == b.h());
 717 }
 718
 719 int operator!= (const AliFmHelix& a, const AliFmHelix& b) {return !(a == b);}
 720
 721 ostream& operator<<(ostream& os, const AliFmHelix& h)
 722 {
 723         return os << '('
 724                 << "curvature = "  << h.curvature() << ", "
 725                 << "dip angle = "  << h.dipAngle()  << ", "
 726                 << "phase = "      << h.phase()     << ", "
 727                 << "h = "          << h.h()         << ", "
 728                 << "origin = "     << h.origin()    << ')';
 729 }
 730
 731
 732
 733
 734