HF/D0_helix/StHelix.cxx

0001 /***************************************************************************
0002  *
0003  * $Id: StHelix.cc,v 1.31 2015/07/20 17:30:51 jeromel Exp $
0004  *
0005  * Author: Thomas Ullrich, Sep 1997
0006  ***************************************************************************
0007  *
0008  * Description: Parametrization of a helix
0009  *
0010  ***************************************************************************
0011  *
0012  * $Log: StHelix.cc,v $
0013  * Revision 1.31  2015/07/20 17:30:51  jeromel
0014  * Use std::isnan to satisfy C++11
0015  *
0016  * Revision 1.30  2015/04/22 18:02:01  ullrich
0017  * Added two default argument to dca of two helices for HFT.
0018  *
0019  * Revision 1.29  2009/09/22 16:21:05  fine
0020  * Silence the compilation warning
0021  *
0022  * Revision 1.28  2008/09/11 20:34:31  ullrich
0023  * Fixed sign problem in seed calculation for helix-helix DCA.
0024  *
0025  * Revision 1.27  2007/08/20 23:25:39  perev
0026  * BugFix #1016
0027  *
0028  * Revision 1.26  2005/10/13 23:15:13  genevb
0029  * Save a few calculations
0030  *
0031  * Revision 1.25  2005/10/13 22:25:35  genevb
0032  * pathLength to plane now finds nearest approach to intersection regardless of # of loops
0033  *
0034  * Revision 1.24  2005/07/06 18:49:56  fisyak
0035  * Replace StHelixD, StLorentzVectorD,StLorentzVectorF,StMatrixD,StMatrixF,StPhysicalHelixD,StThreeVectorD,StThreeVectorF by templated version
0036  *
0037  * Revision 1.23  2004/12/02 02:51:16  ullrich
0038  * Added option to pathLenghth() and distance() to search for
0039  * DCA only within one period. Default stays as it was.
0040  *
0041  * Revision 1.22  2004/05/03 23:35:31  perev
0042  * Possible non init WarnOff
0043  *
0044  * Revision 1.21  2004/01/27 02:49:48  perev
0045  * Big value appropriate for float
0046  *
0047  * Revision 1.20  2003/12/22 18:59:36  perev
0048  * remove test for only +ve pathLeng added before
0049  *
0050  * Revision 1.19  2003/12/18 17:27:02  perev
0051  * Small bug fix in number of iters. i++ ==> i--
0052  *
0053  * Revision 1.18  2003/10/30 20:06:46  perev
0054  * Check of quality added
0055  *
0056  * Revision 1.17  2003/10/19 20:17:00  perev
0057  * Protection agains overfloat added into pathLength(StThreeVector,StThreeVector)
0058  *
0059  * Revision 1.16  2003/10/06 23:39:21  perev
0060  * sqrt(-ve) == no solution. infinity returns
0061  *
0062  * Revision 1.15  2003/09/02 17:59:34  perev
0063  * gcc 3.2 updates + WarnOff
0064  *
0065  * Revision 1.14  2003/06/26 17:15:56  ullrich
0066  * Changed local variable name in pathLenght.
0067  *
0068  * Revision 1.13  2002/06/21 17:49:25  genevb
0069  * Some minor speed improvements
0070  *
0071  * Revision 1.12  2002/04/24 02:40:25  ullrich
0072  * Restored old format and lost CVS statements.
0073  *
0074  * Revision 1.11  2002/02/12 19:37:51  jeromel
0075  * fix for Linux 7.2 (float.h). Took oportunity to doxygenize.
0076  *
0077  * Revision 1.10  2000/07/17 21:44:19  ullrich
0078  * Fixed problem in pathLength(cyl_rad).
0079  *
0080  * Revision 1.9  2000/05/22 21:38:28  ullrich
0081  * Add parenthesis to make Linux compiler happy.
0082  *
0083  * Revision 1.8  2000/05/22 21:11:21  ullrich
0084  * In pathLength(StThreeVector&): Increased number of max iteration
0085  * in Newton method from 10 to 100. Improved initial guess in case
0086  * it is off by n period.
0087  *
0088  * Revision 1.7  2000/03/06 20:24:25  ullrich
0089  * Parameter h for case B=0 correctly handled now.
0090  *
0091  * Revision 1.6  1999/12/22 15:14:39  ullrich
0092  * Added analytical solution for dca between two helices
0093  * in the case for B=0.
0094  *
0095  * Revision 1.5  1999/12/21 15:14:08  ullrich
0096  * Modified to cope with new compiler version on Sun (CC5.0).
0097  *
0098  * Revision 1.4  1999/11/29 21:45:38  fisyak
0099  * fix abs for HP
0100  *
0101  * Revision 1.3  1999/03/07 14:55:41  wenaus
0102  * fix scope problem
0103  *
0104  * Revision 1.2  1999/03/02 19:47:35  ullrich
0105  * Added method to find dca between two helices
0106  *
0107  * Revision 1.1  1999/01/30 03:59:02  fisyak
0108  * Root Version of StarClassLibrary
0109  *
0110  * Revision 1.1  1999/01/23 00:29:15  ullrich
0111  * Initial Revision
0112  *
0113  **************************************************************************/
0114 #if !defined(ST_NO_NUMERIC_LIMITS)
0115 #    include <limits>
0116 #    if !defined(ST_NO_NAMESPACES)
0117          using std::numeric_limits;
0118 #    endif
0119 #endif
0120 #define FOR_HELIX
0121 #include <float.h>
0122 #include "StHelix.h"
0123 #include "PhysicalConstants.h"
0124 #ifdef __ROOT__
0125 ClassImpT(StHelix,double);
0126 #endif
0127 const double StHelix::NoSolution = 3.e+33;
0128
0129 StHelix::StHelix(){ /*noop*/ }
0130
0131 StHelix::StHelix(double c, double d, double phase,
0132          const TVector3& o, int h)
0133 {
0134     setParameters(c, d, phase, o, h);
0135 }
0136
0137 StHelix::~StHelix() { /* noop */ };
0138
0139 void StHelix::setParameters(double c, double dip, double phase,
0140                 const TVector3& o, int h)
0141 {
0142     //
0143     //  The order in which the parameters are set is important
0144     //  since setCurvature might have to adjust the others.
0145     //
0146     mH = (h>=0) ? 1 : -1;    // Default is: positive particle
0147                              //             positive field
0148     mOrigin   = o;
0149     setDipAngle(dip);
0150     setPhase(phase);
0151
0152     //
0153     // Check for singularity and correct for negative curvature.
0154     // May change mH and mPhase. Must therefore be set last.
0155     //
0156     setCurvature(c);
0157
0158     //
0159     // For the case B=0, h is ill defined. In the following we
0160     // always assume h = +1. Since phase = psi - h * pi/2
0161     // we have to correct the phase in case h = -1.
0162     // This assumes that the user uses the same h for phase
0163     // as the one he passed to the constructor.
0164     //
0165     if (mSingularity && mH == -1) {
0166     mH = +1;
0167     setPhase(mPhase-M_PI);
0168     }
0169 }
0170
0171 void StHelix::setCurvature(double val)
0172 {
0173     if (val < 0) {
0174     mCurvature = -val;
0175     mH = -mH;
0176     setPhase(mPhase+M_PI);
0177     }
0178     else
0179     mCurvature = val;
0180
0181 #ifndef ST_NO_NUMERIC_LIMITS
0182     if (fabs(mCurvature) <= numeric_limits<double>::epsilon())
0183 #else
0184     if (fabs(mCurvature) <= static_cast<double>(0))
0185 #endif
0186     mSingularity = true;            // straight line
0187     else
0188     mSingularity = false;               // curved
0189 }
0190
0191 void StHelix::setPhase(double val)
0192 {
0193     mPhase       = val;
0194     mCosPhase    = cos(mPhase);
0195     mSinPhase    = sin(mPhase);
0196     if (fabs(mPhase) > M_PI)
0197     mPhase = atan2(mSinPhase, mCosPhase);  // force range [-pi,pi]
0198 }
0199
0200 void StHelix::setDipAngle(double val)
0201 {
0202     mDipAngle    = val;
0203     mCosDipAngle = cos(mDipAngle);
0204     mSinDipAngle = sin(mDipAngle);
0205 }
0206
0207 double StHelix::xcenter() const
0208 {
0209     if (mSingularity)
0210     return 0;
0211     else
0212     return mOrigin.x()-mCosPhase/mCurvature;
0213 }
0214
0215 double StHelix::ycenter() const
0216 {
0217     if (mSingularity)
0218     return 0;
0219     else
0220     return mOrigin.y()-mSinPhase/mCurvature;
0221 }
0222
0223 double StHelix::fudgePathLength(const TVector3& p) const
0224 {
0225     double s;
0226     double dx = p.x()-mOrigin.x();
0227     double dy = p.y()-mOrigin.y();
0228
0229     if (mSingularity) {
0230     s = (dy*mCosPhase - dx*mSinPhase)/mCosDipAngle;
0231     }
0232     else {
0233     s = atan2(dy*mCosPhase - dx*mSinPhase,
0234           1/mCurvature + dx*mCosPhase+dy*mSinPhase)/
0235         (mH*mCurvature*mCosDipAngle);
0236     }
0237     return s;
0238 }
0239
0240 double StHelix::distance(const TVector3& p, bool scanPeriods) const
0241 {
0242     return (this->at(pathLength(p,scanPeriods))-p).Mag();
0243 }
0244
0245 double StHelix::pathLength(const TVector3& p, bool scanPeriods) const
0246 {
0247     //
0248     //  Returns the path length at the distance of closest
0249     //  approach between the helix and point p.
0250     //  For the case of B=0 (straight line) the path length
0251     //  can be calculated analytically. For B>0 there is
0252     //  unfortunately no easy solution to the problem.
0253     //  Here we use the Newton method to find the root of the
0254     //  referring equation. The 'fudgePathLength' serves
0255     //  as a starting value.
0256     //
0257     double s;
0258     double dx = p.x()-mOrigin.x();
0259     double dy = p.y()-mOrigin.y();
0260     double dz = p.z()-mOrigin.z();
0261
0262     if (mSingularity) {
0263     s = mCosDipAngle*(mCosPhase*dy-mSinPhase*dx) +
0264         mSinDipAngle*dz;
0265     }
0266     else { //
0267 #ifndef ST_NO_NAMESPACES
0268     {
0269         using namespace units;
0270 #endif
0271         const double MaxPrecisionNeeded = micrometer;
0272         const int    MaxIterations      = 100;
0273
0274         //
0275         // The math is taken from Maple with C(expr,optimized) and
0276         // some hand-editing. It is not very nice but efficient.
0277         //
0278         double t34 = mCurvature*mCosDipAngle*mCosDipAngle;
0279         double t41 = mSinDipAngle*mSinDipAngle;
0280         double t6, t7, t11, t12, t19;
0281
0282         //
0283         // Get a first guess by using the dca in 2D. Since
0284         // in some extreme cases we might be off by n periods
0285         // we add (subtract) periods in case we get any closer.
0286         //
0287         s = fudgePathLength(p);
0288
0289         if (scanPeriods) {
0290             double ds = period();
0291             int    j, jmin = 0;
0292             double d, dmin = (at(s) - p).Mag();
0293             for(j=1; j<MaxIterations; j++) {
0294           if ((d = (at(s+j*ds) - p).Mag()) < dmin) {
0295               dmin = d;
0296               jmin = j;
0297           }
0298           else
0299               break;
0300             }
0301             for(j=-1; -j<MaxIterations; j--) {
0302           if ((d = (at(s+j*ds) - p).Mag()) < dmin) {
0303               dmin = d;
0304               jmin = j;
0305           }
0306           else
0307               break;
0308             }
0309             if (jmin) s += jmin*ds;
0310         }
0311
0312         //
0313         // Newtons method:
0314         // Stops after MaxIterations iterations or if the required
0315         // precision is obtained. Whatever comes first.
0316         //
0317         double sOld = s;
0318         for (int i=0; i<MaxIterations; i++) {
0319         t6  = mPhase+s*mH*mCurvature*mCosDipAngle;
0320         t7  = cos(t6);
0321         t11 = dx-(1/mCurvature)*(t7-mCosPhase);
0322         t12 = sin(t6);
0323         t19 = dy-(1/mCurvature)*(t12-mSinPhase);
0324         s  -= (t11*t12*mH*mCosDipAngle-t19*t7*mH*mCosDipAngle -
0325                (dz-s*mSinDipAngle)*mSinDipAngle)/
0326             (t12*t12*mCosDipAngle*mCosDipAngle+t11*t7*t34 +
0327              t7*t7*mCosDipAngle*mCosDipAngle +
0328              t19*t12*t34+t41);
0329         if (fabs(sOld-s) < MaxPrecisionNeeded) break;
0330         sOld = s;
0331         }
0332 #ifndef ST_NO_NAMESPACES
0333     }
0334 #endif
0335     }
0336     return s;
0337 }
0338
0339 double StHelix::period() const
0340 {
0341     if (mSingularity)
0342 #ifndef ST_NO_NUMERIC_LIMITS
0343             return numeric_limits<double>::max();
0344 #else
0345             return DBL_MAX;
0346 #endif
0347     else
0348     return fabs(2*M_PI/(mH*mCurvature*mCosDipAngle));
0349 }
0350
0351 pair<double, double> StHelix::pathLength(double r) const
0352 {
0353     pair<double,double> value;
0354     pair<double,double> VALUE(999999999.,999999999.);
0355     //
0356     // The math is taken from Maple with C(expr,optimized) and
0357     // some hand-editing. It is not very nice but efficient.
0358     // 'first' is the smallest of the two solutions (may be negative)
0359     // 'second' is the other.
0360     //
0361     if (mSingularity) {
0362     double t1 = mCosDipAngle*(mOrigin.x()*mSinPhase-mOrigin.y()*mCosPhase);
0363     double t12 = mOrigin.y()*mOrigin.y();
0364     double t13 = mCosPhase*mCosPhase;
0365     double t15 = r*r;
0366     double t16 = mOrigin.x()*mOrigin.x();
0367     double t20 = -mCosDipAngle*mCosDipAngle*(2.0*mOrigin.x()*mSinPhase*mOrigin.y()*mCosPhase +
0368                  t12-t12*t13-t15+t13*t16);
0369     if (t20<0.) return VALUE;
0370     t20 = ::sqrt(t20);
0371     value.first  = (t1-t20)/(mCosDipAngle*mCosDipAngle);
0372     value.second = (t1+t20)/(mCosDipAngle*mCosDipAngle);
0373     }
0374     else {
0375     double t1 = mOrigin.y()*mCurvature;
0376     double t2 = mSinPhase;
0377     double t3 = mCurvature*mCurvature;
0378     double t4 = mOrigin.y()*t2;
0379     double t5 = mCosPhase;
0380     double t6 = mOrigin.x()*t5;
0381     double t8 = mOrigin.x()*mOrigin.x();
0382     double t11 = mOrigin.y()*mOrigin.y();
0383     double t14 = r*r;
0384     double t15 = t14*mCurvature;
0385     double t17 = t8*t8;
0386     double t19 = t11*t11;
0387     double t21 = t11*t3;
0388     double t23 = t5*t5;
0389     double t32 = t14*t14;
0390     double t35 = t14*t3;
0391     double t38 = 8.0*t4*t6 - 4.0*t1*t2*t8 - 4.0*t11*mCurvature*t6 +
0392                  4.0*t15*t6 + t17*t3 + t19*t3 + 2.0*t21*t8 + 4.0*t8*t23 -
0393                  4.0*t8*mOrigin.x()*mCurvature*t5 - 4.0*t11*t23 -
0394                  4.0*t11*mOrigin.y()*mCurvature*t2 + 4.0*t11 - 4.0*t14 +
0395                  t32*t3 + 4.0*t15*t4 - 2.0*t35*t11 - 2.0*t35*t8;
0396     double t40 = (-t3*t38);
0397     if (t40<0.) return VALUE;
0398     t40 = ::sqrt(t40);
0399
0400     double t43 = mOrigin.x()*mCurvature;
0401     double t45 = 2.0*t5 - t35 + t21 + 2.0 - 2.0*t1*t2 -2.0*t43 - 2.0*t43*t5 + t8*t3;
0402     double t46 = mH*mCosDipAngle*mCurvature;
0403
0404     value.first = (-mPhase + 2.0*atan((-2.0*t1 + 2.0*t2 + t40)/t45))/t46;
0405     value.second = -(mPhase + 2.0*atan((2.0*t1 - 2.0*t2 + t40)/t45))/t46;
0406
0407     //
0408     //   Solution can be off by +/- one period, select smallest
0409     //
0410     double p = period();
0411     if (! std::isnan(value.first)) {
0412         if (fabs(value.first-p) < fabs(value.first)) value.first = value.first-p;
0413         else if (fabs(value.first+p) < fabs(value.first)) value.first = value.first+p;
0414     }
0415     if (! std::isnan(value.second)) {
0416         if (fabs(value.second-p) < fabs(value.second)) value.second = value.second-p;
0417         else if (fabs(value.second+p) < fabs(value.second)) value.second = value.second+p;
0418     }
0419     }
0420     if (value.first > value.second)
0421     swap(value.first,value.second);
0422     return(value);
0423 }
0424
0425 pair<double, double> StHelix::pathLength(double r, double x, double y)
0426 {
0427     double x0 = mOrigin.x();
0428     double y0 = mOrigin.y();
0429     mOrigin.SetX(x0-x);
0430     mOrigin.SetY(y0-y);
0431     pair<double, double> result = this->pathLength(r);
0432     mOrigin.SetX(x0);
0433     mOrigin.SetY(y0);
0434     return result;
0435 }
0436
0437 double StHelix::pathLength(const TVector3& r,
0438                    const TVector3& n) const
0439 {
0440     //
0441     // Vector 'r' defines the position of the center and
0442     // vector 'n' the normal vector of the plane.
0443     // For a straight line there is a simple analytical
0444     // solution. For curvatures > 0 the root is determined
0445     // by Newton method. In case no valid s can be found
0446     // the max. largest value for s is returned.
0447     //
0448     double s;
0449
0450     if (mSingularity) {
0451     double t = n.z()*mSinDipAngle +
0452                n.y()*mCosDipAngle*mCosPhase -
0453                n.x()*mCosDipAngle*mSinPhase;
0454     if (t == 0)
0455         s = NoSolution;
0456     else
0457         s = ((r - mOrigin)*n)/t;
0458     }
0459     else {
0460         const double MaxPrecisionNeeded = micrometer;
0461         const int    MaxIterations      = 20;
0462
0463     double A = mCurvature*((mOrigin - r)*n) -
0464                n.x()*mCosPhase -
0465                n.y()*mSinPhase;
0466     double t = mH*mCurvature*mCosDipAngle;
0467     double u = n.z()*mCurvature*mSinDipAngle;
0468
0469     double a, f, fp;
0470     double sOld = s = 0;
0471     double shiftOld = 0;
0472     double shift;
0473 //      (cos(angMax)-1)/angMax = 0.1
0474         const double angMax = 0.21;
0475         double deltas = fabs(angMax/(mCurvature*mCosDipAngle));
0476 //              dampingFactor = exp(-0.5);
0477 //  double dampingFactor = 0.60653;
0478     int i;
0479
0480     for (i=0; i<MaxIterations; i++) {
0481         a  = t*s+mPhase;
0482             double sina = sin(a);
0483             double cosa = cos(a);
0484         f  = A +
0485          n.x()*cosa +
0486          n.y()*sina +
0487          u*s;
0488         fp = -n.x()*sina*t +
0489           n.y()*cosa*t +
0490           u;
0491             if ( fabs(fp)*deltas <= fabs(f) ) { //too big step
0492                int sgn = 1;
0493                if (fp<0.) sgn = -sgn;
0494                if (f <0.) sgn = -sgn;
0495            shift = sgn*deltas;
0496                if (shift<0) shift*=0.9;  // don't get stuck shifting +/-deltas
0497             } else {
0498                shift = f/fp;
0499             }
0500         s -= shift;
0501         shiftOld = shift;
0502         if (fabs(sOld-s) < MaxPrecisionNeeded) break;
0503         sOld = s;
0504     }
0505         if (i == MaxIterations) return NoSolution;
0506     }
0507     return s;
0508 }
0509
0510 pair<double, double>
0511 StHelix::pathLengths(const StHelix& h, double minStepSize, double minRange) const
0512 {
0513     //
0514     //  Cannot handle case where one is a helix
0515     //  and the other one is a straight line.
0516     //
0517     if (mSingularity != h.mSingularity)
0518         return pair<double, double>(NoSolution, NoSolution);
0519
0520     double s1, s2;
0521
0522     if (mSingularity) {
0523         //
0524         //  Analytic solution
0525         //
0526         TVector3 dv = h.mOrigin - mOrigin;
0527         TVector3 a(-mCosDipAngle*mSinPhase,
0528                                 mCosDipAngle*mCosPhase,
0529                                 mSinDipAngle);
0530         TVector3 b(-h.mCosDipAngle*h.mSinPhase,
0531                                 h.mCosDipAngle*h.mCosPhase,
0532                                 h.mSinDipAngle);
0533         double ab = a*b;
0534         double g  = dv*a;
0535         double k  = dv*b;
0536         s2 = (k-ab*g)/(ab*ab-1.);
0537         s1 = g+s2*ab;
0538         return pair<double, double>(s1, s2);
0539     }
0540     else {
0541         //
0542         //  First step: get dca in the xy-plane as start value
0543         //
0544         double dx = h.xcenter() - xcenter();
0545         double dy = h.ycenter() - ycenter();
0546         double dd = ::sqrt(dx*dx + dy*dy);
0547         double r1 = 1/curvature();
0548         double r2 = 1/h.curvature();
0549
0550         double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd);
0551
0552         double s;
0553         double x, y;
0554         if (fabs(cosAlpha) < 1) {           // two solutions
0555             double sinAlpha = sin(acos(cosAlpha));
0556             x = xcenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd;
0557             y = ycenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd;
0558             s = pathLength(x, y);
0559             x = xcenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd;
0560             y = ycenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd;
0561             double a = pathLength(x, y);
0562             if (h.distance(at(a)) < h.distance(at(s))) s = a;
0563         }
0564         else {                              // no intersection (or exactly one)
0565             int rsign = ((r2-r1) > dd ? -1 : 1); // set -1 when *this* helix is
0566             // completely contained in the other
0567             x = xcenter() + rsign*r1*dx/dd;
0568             y = ycenter() + rsign*r1*dy/dd;
0569             s = pathLength(x, y);
0570         }
0571
0572         //
0573         //   Second step: scan in decreasing intervals around seed 's'
0574         //   minRange and minStepSize are passed as arguments to the method.
0575         //   They have default values defined in the header file.
0576         //
0577         double dmin              = h.distance(at(s));
0578         double range             = max(2*dmin, minRange);
0579         double ds                = range/10;
0580         double slast=-999999, ss, d;
0581         s1 = s - range/2.;
0582         s2 = s + range/2.;
0583
0584         while (ds > minStepSize) {
0585             for (ss=s1; ss<s2+ds; ss+=ds) {
0586                 d = h.distance(at(ss));
0587                 if (d < dmin) {
0588                     dmin = d;
0589                     s = ss;
0590                 }
0591                 slast = ss;
0592             }
0593             //
0594             //  In the rare cases where the minimum is at the
0595             //  the border of the current range we shift the range
0596             //  and start all over, i.e we do not decrease 'ds'.
0597             //  Else we decrease the search intervall around the
0598             //  current minimum and redo the scan in smaller steps.
0599             //
0600             if (s == s1) {
0601                 d = 0.8*(s2-s1);
0602                 s1 -= d;
0603                 s2 -= d;
0604             }
0605             else if (s == slast) {
0606                 d = 0.8*(s2-s1);
0607                 s1 += d;
0608                 s2 += d;
0609             }
0610             else {
0611                 s1 = s-ds;
0612                 s2 = s+ds;
0613                 ds /= 10;
0614             }
0615         }
0616         return pair<double, double>(s, h.pathLength(at(s)));
0617     }
0618 }
0619
0620
0621 void StHelix::moveOrigin(double s)
0622 {
0623     if (mSingularity)
0624     mOrigin = at(s);
0625     else {
0626     TVector3 newOrigin = at(s);
0627     double newPhase = atan2(newOrigin.y() - ycenter(),
0628                 newOrigin.x() - xcenter());
0629     mOrigin = newOrigin;
0630     setPhase(newPhase);
0631     }
0632 }
0633
0634 int operator== (const StHelix& a, const StHelix& b)
0635 {
0636     //
0637     // Checks for numerical identity only !
0638     //
0639     return (a.origin()    == b.origin()    &&
0640         a.dipAngle()  == b.dipAngle()  &&
0641         a.curvature() == b.curvature() &&
0642         a.phase()     == b.phase()     &&
0643         a.h()         == b.h());
0644 }
0645
0646 int operator!= (const StHelix& a, const StHelix& b) {return !(a == b);}
0647
0648 ostream& operator<<(ostream& os, const StHelix& h)
0649 {
0650     return os << "("
0651           << "curvature = "  << h.curvature() << ", "
0652           << "dip angle = "  << h.dipAngle()  << ", "
0653           << "phase = "      << h.phase()     << ", "
0654           << "h = "          << h.h()         << ", "
0655           << "origin = "     << h.origin().x() << " " << h.origin().y() << " " << h.origin().z() << ")";
0656 }
0657
0658
0659
0660