]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliHelix.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Additional protection in comparison macros (M.Ivanov)
[u/mrichter/AliRoot.git] / STEER / AliHelix.cxx
CommitLineData
7f572c00 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"
51ad6848 26#include "AliExternalTrackParam.h"
7f572c00 27#include "TMath.h"
28ClassImp(AliHelix)
29
30
31//_______________________________________________________________________
32AliHelix::AliHelix()
33{
34 //
35 // Default constructor
36 //
37 for (Int_t i =0;i<9;i++) fHelix[i]=0;
38}
39
40//_______________________________________________________________________
176aff27 41AliHelix::AliHelix(const AliHelix &t):TObject(t){
7f572c00 42 //
43 //
44 for (Int_t i=0;i<9;i++)
45 fHelix[i]=t.fHelix[i];
46}
47
48AliHelix::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
51ad6848 83
84AliHelix::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
7f572c00 121AliHelix::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
161void 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
173void 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 //
51ad6848 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;
7f572c00 203
51ad6848 204
7f572c00 205
206
207}
208
209
210void 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
232Double_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
176aff27 263Int_t AliHelix::GetPhase(Double_t /*r0*/, Double_t * /*t[2]*/)
7f572c00 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
275Double_t AliHelix::GetPhaseZ(Double_t z0)
276{
277 //
278 //
279 return (z0-fHelix[1])/fHelix[3];
280}
281
282
283Int_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]);
51ad6848 297 if (d<0.000000000001) return 0;
7f572c00 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];
51ad6848 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;
7f572c00 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
361Int_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
440Int_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/*
484Int_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
582Int_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