STEER/AliHelix.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //-------------------------------------------------------------------------
  19 //                Implementation of the AliHelix class
  20 //        Origin: Marian Ivanov, CERN, marian.ivanov@cern.ch
  21 //-------------------------------------------------------------------------
  22
  23
  24 #include "AliHelix.h"
  25 #include "AliKalmanTrack.h"
  26 #include "AliExternalTrackParam.h"
  27 #include "TMath.h"
  28 ClassImp(AliHelix)
  29
  30
  31 //_______________________________________________________________________
  32 AliHelix::AliHelix()
  33 {
  34   //
  35   // Default constructor
  36   //
  37   for (Int_t i =0;i<9;i++) fHelix[i]=0;
  38 }
  39
  40 //_______________________________________________________________________
  41 AliHelix::AliHelix(const AliHelix &t):TObject(t){
  42   //
  43   //
  44   for (Int_t i=0;i<9;i++)
  45     fHelix[i]=t.fHelix[i];
  46 }
  47
  48 AliHelix::AliHelix(const AliKalmanTrack &t)
  49 {
  50   //
  51   //
  52   Double_t alpha,x,cs,sn;
  53   t.GetExternalParameters(x,fHelix);
  54   alpha=t.GetAlpha();
  55   //
  56   //circle parameters
  57   fHelix[4]=fHelix[4]/t.GetConvConst();    // C
  58   cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
  59
  60   Double_t xc, yc, rc;
  61   rc  =  1/fHelix[4];
  62   xc  =  x-fHelix[2]*rc;
  63   yc  =  fHelix[0]+TMath::Sqrt(1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4])/fHelix[4];
  64
  65   fHelix[6] = xc*cs - yc*sn;
  66   fHelix[7] = xc*sn + yc*cs;
  67   fHelix[8] =  TMath::Abs(rc);
  68   //
  69   //
  70   fHelix[5]=x*cs - fHelix[0]*sn;            // x0
  71   fHelix[0]=x*sn + fHelix[0]*cs;            // y0
  72   //fHelix[1]=                               // z0
  73   fHelix[2]=TMath::ASin(fHelix[2]) + alpha; // phi0
  74   //fHelix[3]=                               // tgl
  75   //
  76   //
  77   fHelix[5]   = fHelix[6];
  78   fHelix[0]   = fHelix[7];
  79   //fHelix[5]-=TMath::Sin(fHelix[2])/fHelix[4];
  80   //fHelix[0]+=TMath::Cos(fHelix[2])/fHelix[4];
  81 }
  82
  83
  84 AliHelix::AliHelix(const AliExternalTrackParam &t)
  85 {
  86   //
  87   //
  88   Double_t alpha,x,cs,sn;
  89   const Double_t *param =t.GetParameter();
  90   for (Int_t i=0;i<5;i++) fHelix[i]=param[i];
  91   x = t.X();
  92   alpha=t.Alpha();
  93   //
  94   //circle parameters
  95   fHelix[4]=fHelix[4]/AliKalmanTrack::GetConvConst();    // C
  96   cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
  97
  98   Double_t xc, yc, rc;
  99   rc  =  1/fHelix[4];
 100   xc  =  x-fHelix[2]*rc;
 101   yc  =  fHelix[0]+TMath::Sqrt(1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4])/fHelix[4];
 102
 103   fHelix[6] = xc*cs - yc*sn;
 104   fHelix[7] = xc*sn + yc*cs;
 105   fHelix[8] =  TMath::Abs(rc);
 106   //
 107   //
 108   fHelix[5]=x*cs - fHelix[0]*sn;            // x0
 109   fHelix[0]=x*sn + fHelix[0]*cs;            // y0
 110   //fHelix[1]=                               // z0
 111   fHelix[2]=TMath::ASin(fHelix[2]) + alpha; // phi0
 112   //fHelix[3]=                               // tgl
 113   //
 114   //
 115   fHelix[5]   = fHelix[6];
 116   fHelix[0]   = fHelix[7];
 117   //fHelix[5]-=TMath::Sin(fHelix[2])/fHelix[4];
 118   //fHelix[0]+=TMath::Cos(fHelix[2])/fHelix[4];
 119 }
 120
 121 AliHelix::AliHelix(Double_t x[3], Double_t p[3], Double_t charge, Double_t conversion)
 122 {
 123   //
 124   //
 125   //
 126   Double_t pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
 127   if (TMath::Abs(conversion)<0.00000001)
 128     conversion = AliKalmanTrack::GetConvConst();
 129   //
 130   //
 131   fHelix[4] = charge/(conversion*pt); // C
 132   fHelix[3] = p[2]/pt;    // tgl
 133   //
 134   Double_t xc, yc, rc;
 135   rc  =  1/fHelix[4];
 136   xc  =  x[0]  -rc*p[1]/pt;
 137   yc  =  x[1]  +rc*p[0]/pt;
 138   //
 139   fHelix[5] = x[0];   // x0
 140   fHelix[0] = x[1];   // y0
 141   fHelix[1] = x[2];   // z0
 142   //
 143   fHelix[6] = xc;
 144   fHelix[7] = yc;
 145   fHelix[8] = TMath::Abs(rc);
 146   //
 147   fHelix[5]=xc;
 148   fHelix[0]=yc;
 149   //
 150   if (TMath::Abs(p[1])<TMath::Abs(p[0])){
 151     fHelix[2]=TMath::ASin(p[1]/pt);
 152     if (charge*yc<charge*x[1])  fHelix[2] = TMath::Pi()-fHelix[2];
 153   }
 154   else{
 155     fHelix[2]=TMath::ACos(p[0]/pt);
 156     if (charge*xc>charge*x[0])  fHelix[2] = -fHelix[2];
 157   }
 158
 159 }
 160
 161 void  AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion)
 162 {
 163   // return  momentum at given phase
 164   Double_t x[3],g[3],gg[3];
 165   Evaluate(phase,x,g,gg);
 166   if (TMath::Abs(conversion)<0.0001) conversion = AliKalmanTrack::GetConvConst();
 167   Double_t mt = TMath::Sqrt(g[0]*g[0]+g[1]*g[1]);
 168   p[0] = fHelix[8]*g[0]/(mt*conversion);
 169   p[1] = fHelix[8]*g[1]/(mt*conversion);
 170   p[2] = fHelix[8]*g[2]/(mt*conversion);
 171 }
 172
 173 void   AliHelix::GetAngle(Double_t t1, AliHelix &h, Double_t t2, Double_t angle[3])
 174 {
 175   //
 176   //
 177   //
 178   Double_t x1[3],g1[3],gg1[3];
 179   Double_t x2[3],g2[3],gg2[3];
 180   Evaluate(t1,x1,g1,gg1);
 181   h.Evaluate(t2,x2,g2,gg2);
 182
 183   //
 184   Double_t norm1r = g1[0]*g1[0]+g1[1]*g1[1];
 185   Double_t norm1  = TMath::Sqrt(norm1r+g1[2]*g1[2]);
 186   norm1r         = TMath::Sqrt(norm1r);
 187   //
 188   Double_t norm2r = g2[0]*g2[0]+g2[1]*g2[1];
 189   Double_t norm2  = TMath::Sqrt(norm2r+g2[2]*g2[2]);
 190   norm2r         = TMath::Sqrt(norm2r);
 191   //
 192   angle[0]  = (g1[0]*g2[0]+g1[1]*g2[1])/(norm1r*norm2r);   // angle in phi projection
 193   if (TMath::Abs(angle[0])<1.) angle[0] = TMath::ACos(angle[0]);
 194   else angle[0]=0;
 195   //
 196   angle[1]  = ((norm1r*norm2r)+g1[2]*g2[2])/(norm1*norm2); // angle in rz  projection
 197   if (TMath::Abs(angle[1])<1.) angle[1] = TMath::ACos(angle[1]);
 198   else angle[1]=0;
 199
 200   angle[2]  = (g1[0]*g2[0]+g1[1]*g2[1]+g1[2]*g2[2])/(norm1*norm2); //3D angle
 201   if (TMath::Abs(angle[2])<1.) angle[2] = TMath::ACos(angle[2]);
 202   else angle[2]=0;
 203
 204
 205
 206
 207 }
 208
 209
 210 void AliHelix::Evaluate(Double_t t,
 211                      Double_t r[3],  //radius vector
 212                      Double_t g[3],  //first defivatives
 213                      Double_t gg[3]) //second derivatives
 214 {
 215   //--------------------------------------------------------------------
 216   // Calculate position of a point on a track and some derivatives at given phase
 217   //--------------------------------------------------------------------
 218   Double_t phase=fHelix[4]*t+fHelix[2];
 219   Double_t sn=TMath::Sin(phase), cs=TMath::Cos(phase);
 220
 221   //r[0] = fHelix[5] + (sn - fHelix[6])/fHelix[4];
 222   //r[1] = fHelix[0] - (cs - fHelix[7])/fHelix[4];
 223   r[0] = fHelix[5] + sn/fHelix[4];
 224   r[1] = fHelix[0] - cs/fHelix[4];
 225   r[2] = fHelix[1] + fHelix[3]*t;
 226
 227   g[0] = cs; g[1]=sn; g[2]=fHelix[3];
 228
 229   gg[0]=-fHelix[4]*sn; gg[1]=fHelix[4]*cs; gg[2]=0.;
 230 }
 231
 232 Double_t  AliHelix::GetPhase(Double_t x, Double_t y )
 233
 234 {
 235   //
 236   //calculate helix param at given x,y  point
 237   //
 238   Double_t phase  =  (x-fHelix[5])*fHelix[4];
 239   if (TMath::Abs(phase)>=1)
 240     phase = TMath::Sign(0.99999999999,phase);
 241   phase = TMath::ASin(phase);
 242
 243   if ( ( ( fHelix[0]-y)*fHelix[4]) < 0.) {
 244     if (phase>0)
 245       phase = TMath::Pi()-phase;
 246     else
 247       phase = -(TMath::Pi()+phase);
 248   }
 249   if ( (phase-fHelix[2])>TMath::Pi())  phase = phase-2.*TMath::Pi();
 250   if ( (phase-fHelix[2])<-TMath::Pi()) phase = phase+2.*TMath::Pi();
 251
 252   Double_t t     = (phase-fHelix[2])/fHelix[4];
 253
 254   //  Double_t r[3];
 255   //Evaluate(t,r);
 256   //if  ( (TMath::Abs(r[0]-x)>0.01) || (TMath::Abs(r[1]-y)>0.01)){
 257   //     Double_t phase  =  (x-fHelix[5])*fHelix[4];
 258   //  printf("problem\n");
 259   //}
 260   return t;
 261 }
 262
 263 Int_t AliHelix::GetPhase(Double_t /*r0*/, Double_t * /*t[2]*/)
 264 {
 265   //
 266   //calculate helix param at given r  point - return nearest point ()
 267   //
 268   // not implemented yet
 269
 270
 271   return 0;
 272 }
 273
 274
 275 Double_t  AliHelix::GetPhaseZ(Double_t z0)
 276 {
 277   //
 278   //
 279   return (z0-fHelix[1])/fHelix[3];
 280 }
 281
 282
 283 Int_t    AliHelix::GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Double_t ri[2], Double_t cut)
 284 {
 285   //--------------------------------------------------------------------
 286   // This function returns  phase vectors with intesection between helix (0, 1 or 2)
 287   // in x-y plane projection
 288   //--------------------------------------------------------------------
 289   //
 290   //  Double_t * c1 = &fHelix[6];
 291   //Double_t * c2 = &(h.fHelix[6]);
 292   //  Double_t  c1[3] = {fHelix[5],fHelix[0],fHelix[8]};
 293   Double_t  c1[3] = {0,0,fHelix[8]};
 294   Double_t  c2[3] = {h.fHelix[5]-fHelix[5],h.fHelix[0]-fHelix[0],h.fHelix[8]};
 295
 296   Double_t d  = TMath::Sqrt(c2[0]*c2[0]+c2[1]*c2[1]);
 297   if (d<0.000000000001) return 0;
 298   //
 299   Double_t x0[2];
 300   Double_t y0[2];
 301   //
 302   if ( d>=(c1[2]+c2[2])){
 303     if (d>=(c1[2]+c2[2]+cut)) return 0;
 304     x0[0] = (d+c1[2]-c2[2])*c2[0]/(2*d)+ fHelix[5];
 305     y0[0] = (d+c1[2]-c2[2])*c2[1]/(2*d)+ fHelix[0];
 306     //    return 0;
 307     phase[0][0] = GetPhase(x0[0],y0[0]);
 308     phase[0][1] = h.GetPhase(x0[0],y0[0]);
 309     ri[0] = x0[0]*x0[0]+y0[0]*y0[0];
 310     return 1;
 311   }
 312   if ( (d+c2[2])<c1[2]){
 313     if ( (d+c2[2])+cut<c1[2]) return 0;
 314     //
 315     Double_t xx = c2[0]+ c2[0]*c2[2]/d+ fHelix[5];
 316     Double_t yy = c2[1]+ c2[1]*c2[2]/d+ fHelix[0];
 317     phase[0][1] = h.GetPhase(xx,yy);
 318     //
 319     Double_t xx2 = c2[0]*c1[2]/d+ fHelix[5];
 320     Double_t yy2 = c2[1]*c1[2]/d+ fHelix[0];
 321     phase[0][0] = GetPhase(xx2,yy2);
 322     ri[0] = xx*xx+yy*yy;
 323     return 1;
 324   }
 325
 326   if ( (d+c1[2])<c2[2]){
 327     if ( (d+c1[2])+cut<c2[2]) return 0;
 328     //
 329     Double_t xx = -c2[0]*c1[2]/d+ fHelix[5];
 330     Double_t yy = -c2[1]*c1[2]/d+ fHelix[0];
 331     phase[0][1] = GetPhase(xx,yy);
 332     //
 333     Double_t xx2 = c2[0]- c2[0]*c2[2]/d+ fHelix[5];
 334     Double_t yy2 = c2[1]- c2[1]*c2[2]/d+ fHelix[0];
 335     phase[0][0] = h.GetPhase(xx2,yy2);
 336     ri[0] = xx*xx+yy*yy;
 337     return 1;
 338   }
 339
 340   Double_t d1 = (d*d+c1[2]*c1[2]-c2[2]*c2[2])/(2.*d);
 341   Double_t v1 = c1[2]*c1[2]-d1*d1;
 342   if (v1<0) return 0;
 343   v1 = TMath::Sqrt(v1);
 344   //
 345   x0[0] = (c2[0]*d1+c2[1]*v1)/d + fHelix[5];
 346   y0[0] = (c2[1]*d1-c2[0]*v1)/d + fHelix[0];
 347   //
 348   x0[1] = (c2[0]*d1-c2[1]*v1)/d + fHelix[5];
 349   y0[1] = (c2[1]*d1+c2[0]*v1)/d + fHelix[0];
 350   //
 351   for (Int_t i=0;i<2;i++){
 352     phase[i][0] = GetPhase(x0[i],y0[i]);
 353     phase[i][1] = h.GetPhase(x0[i],y0[i]);
 354     ri[i] = x0[i]*x0[i]+y0[i]*y0[i];
 355   }
 356   return 2;
 357 }
 358
 359 /*
 360
 361 Int_t    AliHelix::GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Double_t ri[2], Double_t cut)
 362 {
 363   //--------------------------------------------------------------------
 364   // This function returns  phase vectors with intesection between helix (0, 1 or 2)
 365   // in x-y plane projection
 366   //--------------------------------------------------------------------
 367   //
 368   Double_t * c1 = &fHelix[6];
 369   Double_t * c2 = &(h.fHelix[6]);
 370   Double_t d  = TMath::Sqrt((c1[0]-c2[0])*(c1[0]-c2[0])+(c1[1]-c2[1])*(c1[1]-c2[1]));
 371   //
 372   Double_t x0[2];
 373   Double_t y0[2];
 374   //
 375   if ( d>=(c1[2]+c2[2])){
 376     if (d>=(c1[2]+c2[2]+cut)) return 0;
 377     x0[0] = c1[0]+ (d+c1[2]-c2[2])*(c2[0]-c1[0])/(2*d);
 378     y0[0] = c1[1]+ (d+c1[2]-c2[2])*(c2[1]-c1[1])/(2*d);
 379     return 0;
 380     phase[0][0] = GetPhase(x0[0],y0[0]);
 381     phase[0][1] = h.GetPhase(x0[0],y0[0]);
 382     ri[0] = x0[0]*x0[0]+y0[0]*y0[0];
 383     return 1;
 384   }
 385   if ( (d+c2[2])<c1[2]){
 386     if ( (d+c2[2])+cut<c1[2]) return 0;
 387     //
 388     Double_t xx = c2[0]+ (c2[0]-c1[0])*c2[2]/d;
 389     Double_t yy = c2[1]+ (c2[1]-c1[1])*c2[2]/d;
 390     phase[0][1] = h.GetPhase(xx,yy);
 391     //
 392     Double_t xx2 = c1[0]- (c1[0]-c2[0])*c1[2]/d;
 393     Double_t yy2 = c1[1]- (c1[1]-c2[1])*c1[2]/d;
 394     phase[0][0] = GetPhase(xx2,yy2);
 395     //if ( (TMath::Abs(xx2-xx)>cut)||(TMath::Abs(yy2-yy)>cut)){
 396     //  printf("problem\n");
 397     //}
 398     ri[0] = xx*xx+yy*yy;
 399     return 1;
 400   }
 401
 402   if ( (d+c1[2])<c2[2]){
 403     if ( (d+c1[2])+cut<c2[2]) return 0;
 404     //
 405     Double_t xx = c1[0]+ (c1[0]-c2[0])*c1[2]/d;
 406     Double_t yy = c1[1]+ (c1[1]-c2[1])*c1[2]/d;
 407     phase[0][1] = GetPhase(xx,yy);
 408     //
 409     Double_t xx2 = c2[0]- (c2[0]-c1[0])*c2[2]/d;
 410     Double_t yy2 = c2[1]- (c2[1]-c1[1])*c2[2]/d;
 411     phase[0][0] = h.GetPhase(xx2,yy2);
 412     //if ( (TMath::Abs(xx2-xx)>cut)||(TMath::Abs(yy2-yy)>cut)){
 413     //  printf("problem\n");
 414     //}
 415     ri[0] = xx*xx+yy*yy;
 416     return 1;
 417   }
 418
 419   Double_t d1 = (d*d+c1[2]*c1[2]-c2[2]*c2[2])/(2.*d);
 420   Double_t v1 = c1[2]*c1[2]-d1*d1;
 421   if (v1<0) return 0;
 422   v1 = TMath::Sqrt(v1);
 423   //
 424   x0[0] = c1[0]+ ((c2[0]-c1[0])*d1-(c1[1]-c2[1])*v1)/d;
 425   y0[0] = c1[1]+ ((c2[1]-c1[1])*d1+(c1[0]-c2[0])*v1)/d;
 426   //
 427   x0[1] = c1[0]+ ((c2[0]-c1[0])*d1+(c1[1]-c2[1])*v1)/d;
 428   y0[1] = c1[1]+ ((c2[1]-c1[1])*d1-(c1[0]-c2[0])*v1)/d;
 429   //
 430   for (Int_t i=0;i<2;i++){
 431     phase[i][0] = GetPhase(x0[i],y0[i]);
 432     phase[i][1] = h.GetPhase(x0[i],y0[i]);
 433     ri[i] = x0[i]*x0[i]+y0[i]*y0[i];
 434   }
 435   return 2;
 436 }
 437 */
 438
 439
 440 Int_t   AliHelix::LinearDCA(AliHelix &h, Double_t &t1, Double_t &t2,
 441                       Double_t &R, Double_t &dist)
 442 {
 443   //
 444   //
 445   // find intersection using linear approximation
 446   Double_t r1[3],g1[3],gg1[3];
 447   Double_t r2[3],g2[3],gg2[3];
 448   //
 449   Evaluate(t1,r1,g1,gg1);
 450   h.Evaluate(t2,r2,g2,gg2);
 451   //
 452   Double_t g1_2 = g1[0]*g1[0] +g1[1]*g1[1] +g1[2]*g1[2];
 453   Double_t g2_2 = g2[0]*g2[0] +g2[1]*g2[1] +g2[2]*g2[2];
 454   Double_t g1x2 = g1[0]*g2[0] +g1[1]*g2[1] +g1[2]*g2[2];
 455   Double_t det  = g1_2*g2_2   - g1x2*g1x2;
 456   //
 457   if (TMath::Abs(det)>0){
 458     //
 459     Double_t r1g1 = r1[0]*g1[0] +r1[1]*g1[1] +r1[2]*g1[2];
 460     Double_t r2g1 = r2[0]*g1[0] +r2[1]*g1[1] +r2[2]*g1[2];
 461     Double_t r1g2 = r1[0]*g2[0] +r1[1]*g2[1] +r1[2]*g2[2];
 462     Double_t r2g2 = r2[0]*g2[0] +r2[1]*g2[1] +r2[2]*g2[2];
 463     //
 464     Double_t dt    = - ( g2_2*(r1g1-r2g1) - g1x2*(r1g2-r2g2)) / det;
 465     Double_t dp    = - ( g1_2*(r2g2-r1g2) - g1x2*(r2g1-r1g1)) / det;
 466     //
 467     t1+=dt;
 468     t2+=dp;
 469     Evaluate(t1,r1);
 470     h.Evaluate(t2,r2);
 471     //
 472     dist = (r1[0]-r2[0])*(r1[0]-r2[0])+
 473                                   (r1[1]-r2[1])*(r1[1]-r2[1])+
 474                                   (r1[2]-r2[2])*(r1[2]-r2[2]);
 475     R = ((r1[0]+r2[0])*(r1[0]+r2[0])+(r1[1]+r2[1])*(r1[1]+r2[1]))/4.;
 476   }
 477   return 0;
 478 }
 479
 480
 481
 482
 483 /*
 484 Int_t  AliHelix::ParabolicDCA(AliHelix&h,  //helixes
 485                                Double_t &t1, Double_t &t2,
 486                                Double_t &R, Double_t &dist, Int_t iter)
 487 {
 488   //
 489   //
 490   // find intersection using linear fit
 491   Double_t r1[3],g1[3],gg1[3];
 492   Double_t r2[3],g2[3],gg2[3];
 493   //
 494   Evaluate(t1,r1,g1,gg1);
 495   h.Evaluate(t2,r2,g2,gg2);
 496
 497   //
 498   Double_t dx2=1.;
 499   Double_t dy2=1.;
 500   Double_t dz2=1.;
 501   //
 502   Double_t dx=r2[0]-r1[0], dy=r2[1]-r1[1], dz=r2[2]-r1[2];
 503   Double_t dm=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
 504   //
 505
 506  iter++;
 507  while (iter--) {
 508
 509      Double_t gt1=-(dx*g1[0]/dx2 + dy*g1[1]/dy2 + dz*g1[2]/dz2);
 510      Double_t gt2=+(dx*g2[0]/dx2 + dy*g2[1]/dy2 + dz*g2[2]/dz2);
 511      Double_t h11=(g1[0]*g1[0] - dx*gg1[0])/dx2 +
 512                   (g1[1]*g1[1] - dy*gg1[1])/dy2 +
 513                   (g1[2]*g1[2] - dz*gg1[2])/dz2;
 514      Double_t h22=(g2[0]*g2[0] + dx*gg2[0])/dx2 +
 515                   (g2[1]*g2[1] + dy*gg2[1])/dy2 +
 516                   (g2[2]*g2[2] + dz*gg2[2])/dz2;
 517      Double_t h12=-(g1[0]*g2[0]/dx2 + g1[1]*g2[1]/dy2 + g1[2]*g2[2]/dz2);
 518
 519      Double_t det=h11*h22-h12*h12;
 520
 521      Double_t dt1,dt2;
 522      if (TMath::Abs(det)<1.e-33) {
 523         //(quasi)singular Hessian
 524         dt1=-gt1; dt2=-gt2;
 525      } else {
 526         dt1=-(gt1*h22 - gt2*h12)/det;
 527         dt2=-(h11*gt2 - h12*gt1)/det;
 528      }
 529
 530      if ((dt1*gt1+dt2*gt2)>0) {dt1=-dt1; dt2=-dt2;}
 531
 532      //check delta(phase1) ?
 533      //check delta(phase2) ?
 534
 535      if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4)
 536      if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) {
 537        //if ((gt1*gt1+gt2*gt2) > 1.e-4/dy2/dy2)
 538        //  Warning("GetDCA"," stopped at not a stationary point !\n");
 539         Double_t lmb=h11+h22; lmb=lmb-TMath::Sqrt(lmb*lmb-4*det);
 540         if (lmb < 0.)
 541           //Warning("GetDCA"," stopped at not a minimum !\n");
 542         break;
 543      }
 544
 545      Double_t dd=dm;
 546      for (Int_t div=1 ; ; div*=2) {
 547         Evaluate(t1+dt1,r1,g1,gg1);
 548         h.Evaluate(t2+dt2,r2,g2,gg2);
 549         dx=r2[0]-r1[0]; dy=r2[1]-r1[1]; dz=r2[2]-r1[2];
 550         dd=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
 551         if (dd<dm) break;
 552         dt1*=0.5; dt2*=0.5;
 553         if (div>512) {
 554           //Warning("GetDCA"," overshoot !\n");
 555           break;
 556         }
 557      }
 558      dm=dd;
 559
 560      t1+=dt1;
 561      t2+=dt2;
 562
 563  }
 564
 565  Evaluate(t1,r1,g1,gg1);
 566  h.Evaluate(t2,r2,g2,gg2);
 567  //
 568  dist = (r1[0]-r2[0])*(r1[0]-r2[0])+
 569    (r1[1]-r2[1])*(r1[1]-r2[1])+
 570    (r1[2]-r2[2])*(r1[2]-r2[2]);
 571
 572  R = ((r1[0]+r2[0])*(r1[0]+r2[0])+(r1[1]+r2[1])*(r1[1]+r2[1]))/4;
 573
 574 }
 575 */
 576
 577
 578
 579
 580
 581
 582 Int_t  AliHelix::ParabolicDCA(AliHelix&h,  //helixes
 583                                Double_t &t1, Double_t &t2,
 584                                Double_t &R, Double_t &dist, Int_t iter)
 585 {
 586   //
 587   //
 588   // find intersection using linear fit
 589   Double_t r1[3],g1[3],gg1[3];
 590   Double_t r2[3],g2[3],gg2[3];
 591   //
 592   Evaluate(t1,r1,g1,gg1);
 593   h.Evaluate(t2,r2,g2,gg2);
 594
 595   //
 596   Double_t dx2=1.;
 597   Double_t dy2=1.;
 598   Double_t dz2=1.;
 599   //
 600   Double_t dx=r2[0]-r1[0], dy=r2[1]-r1[1], dz=r2[2]-r1[2];
 601   Double_t dm=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
 602   //
 603
 604  iter++;
 605  while (iter--) {
 606     Double_t gt1=-(dx*g1[0]/dx2 + dy*g1[1]/dy2 + dz*g1[2]/dz2);
 607     Double_t gt2=+(dx*g2[0]/dx2 + dy*g2[1]/dy2 + dz*g2[2]/dz2);
 608
 609     Double_t h11=(g1[0]*g1[0] - dx*gg1[0])/dx2 +
 610       (g1[1]*g1[1] - dy*gg1[1])/dy2 +
 611       (g1[2]*g1[2] - dz*gg1[2])/dz2;
 612     Double_t h22=(g2[0]*g2[0] + dx*gg2[0])/dx2 +
 613       (g2[1]*g2[1] + dy*gg2[1])/dy2 +
 614       (g2[2]*g2[2] + dz*gg2[2])/dz2;
 615     Double_t h12=-(g1[0]*g2[0]/dx2 + g1[1]*g2[1]/dy2 + g1[2]*g2[2]/dz2);
 616
 617     Double_t det=h11*h22-h12*h12;
 618
 619     Double_t dt1,dt2;
 620     if (TMath::Abs(det)<1.e-33) {
 621       //(quasi)singular Hessian
 622       dt1=-gt1; dt2=-gt2;
 623     } else {
 624       dt1=-(gt1*h22 - gt2*h12)/det;
 625       dt2=-(h11*gt2 - h12*gt1)/det;
 626     }
 627
 628     if ((dt1*gt1+dt2*gt2)>0) {dt1=-dt1; dt2=-dt2;}
 629
 630     //if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4)
 631     //  if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) {
 632     //  break;
 633     //  }
 634
 635     Double_t dd=dm;
 636     for (Int_t div=1 ; div<512 ; div*=2) {
 637       Evaluate(t1+dt1,r1,g1,gg1);
 638       h.Evaluate(t2+dt2,r2,g2,gg2);
 639       dx=r2[0]-r1[0]; dy=r2[1]-r1[1]; dz=r2[2]-r1[2];
 640       dd=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
 641       if (dd<dm) break;
 642       dt1*=0.5; dt2*=0.5;
 643       if (div==0){
 644         div =1;
 645       }
 646       if (div>512) {
 647         break;
 648       }
 649     }
 650     dm=dd;
 651     t1+=dt1;
 652     t2+=dt2;
 653  }
 654  Evaluate(t1,r1,g1,gg1);
 655  h.Evaluate(t2,r2,g2,gg2);
 656  //
 657  dist = (r1[0]-r2[0])*(r1[0]-r2[0])+
 658    (r1[1]-r2[1])*(r1[1]-r2[1])+
 659    (r1[2]-r2[2])*(r1[2]-r2[2]);
 660
 661  R = ((r1[0]+r2[0])*(r1[0]+r2[0])+(r1[1]+r2[1])*(r1[1]+r2[1]))/4;
 662  return 0;
 663
 664 }
 665