AliRoot Core  edcc906 (edcc906)
AliExternalTrackParam.cxx
Go to the documentation of this file.
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 
19 // //
20 // Implementation of the external track parameterisation class. //
21 // //
22 // This parameterisation is used to exchange tracks between the detectors. //
23 // A set of functions returning the position and the momentum of tracks //
24 // in the global coordinate system as well as the track impact parameters //
25 // are implemented.
26 // Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch //
28 #include <cassert>
29 
30 #include <TVectorD.h>
31 #include <TMatrixDSym.h>
32 #include <TPolyMarker3D.h>
33 #include <TVector3.h>
34 #include <TMatrixD.h>
35 
36 #include "AliExternalTrackParam.h"
37 #include "AliVVertex.h"
38 #include "AliLog.h"
39 
40 ClassImp(AliExternalTrackParam)
41 
42 Double32_t AliExternalTrackParam::fgMostProbablePt=kMostProbablePt;
43 Bool_t AliExternalTrackParam::fgUseLogTermMS = kFALSE;;
44 //_____________________________________________________________________________
46  AliVTrack(),
47  fX(0),
48  fAlpha(0)
49 {
50  //
51  // default constructor
52  //
53  for (Int_t i = 0; i < 5; i++) fP[i] = 0;
54  for (Int_t i = 0; i < 15; i++) fC[i] = 0;
55 }
56 
57 //_____________________________________________________________________________
59  AliVTrack(track),
60  fX(track.fX),
61  fAlpha(track.fAlpha)
62 {
63  //
64  // copy constructor
65  //
66  for (Int_t i = 0; i < 5; i++) fP[i] = track.fP[i];
67  for (Int_t i = 0; i < 15; i++) fC[i] = track.fC[i];
69 }
70 
71 //_____________________________________________________________________________
73 {
74  //
75  // assignment operator
76  //
77 
78  if (this!=&trkPar) {
79  AliVTrack::operator=(trkPar);
80  fX = trkPar.fX;
81  fAlpha = trkPar.fAlpha;
82 
83  for (Int_t i = 0; i < 5; i++) fP[i] = trkPar.fP[i];
84  for (Int_t i = 0; i < 15; i++) fC[i] = trkPar.fC[i];
86  }
87 
88  return *this;
89 }
90 
91 //_____________________________________________________________________________
93  const Double_t param[5],
94  const Double_t covar[15]) :
95  AliVTrack(),
96  fX(x),
97  fAlpha(alpha)
98 {
99  //
100  // create external track parameters from given arguments
101  //
102  for (Int_t i = 0; i < 5; i++) fP[i] = param[i];
103  for (Int_t i = 0; i < 15; i++) fC[i] = covar[i];
104  CheckCovariance();
105 }
106 
107 //_____________________________________________________________________________
109 {
110  //
111  // Recreate TrackParams from VTrack
112  // This is not a copy contructor !
113  //
114  if (!vTrack) {
115  AliError("Source VTrack is NULL");
116  return;
117  }
118  if (this==vTrack) {
119  AliError("Copy of itself is requested");
120  return;
121  }
122  //
123  if (vTrack->InheritsFrom(AliExternalTrackParam::Class())) {
124  AliDebug(1,"Source itself is AliExternalTrackParam, using assignment operator");
125  *this = *(AliExternalTrackParam*)vTrack;
126  return;
127  }
128  //
129  AliVTrack::operator=( *vTrack );
130  //
131  Double_t xyz[3],pxpypz[3],cv[21];
132  vTrack->GetXYZ(xyz);
133  pxpypz[0]=vTrack->Px();
134  pxpypz[1]=vTrack->Py();
135  pxpypz[2]=vTrack->Pz();
136  vTrack->GetCovarianceXYZPxPyPz(cv);
137  Short_t sign = (Short_t)vTrack->Charge();
138  Set(xyz,pxpypz,cv,sign);
139 }
140 
141 //_____________________________________________________________________________
143  AliVTrack(),
144  fX(0.),
145  fAlpha(0.)
146 {
147  //
148  // Constructor from virtual track,
149  // This is not a copy contructor !
150  //
151 
152  if (vTrack->InheritsFrom("AliExternalTrackParam")) {
153  AliError("This is not a copy constructor. Use AliExternalTrackParam(const AliExternalTrackParam &) !");
154  AliWarning("Calling the default constructor...");
156  return;
157  }
158 
159  Double_t xyz[3],pxpypz[3],cv[21];
160  vTrack->GetXYZ(xyz);
161  pxpypz[0]=vTrack->Px();
162  pxpypz[1]=vTrack->Py();
163  pxpypz[2]=vTrack->Pz();
164  vTrack->GetCovarianceXYZPxPyPz(cv);
165  Short_t sign = (Short_t)vTrack->Charge();
166 
167  Set(xyz,pxpypz,cv,sign);
168 }
169 
170 //_____________________________________________________________________________
171 AliExternalTrackParam::AliExternalTrackParam(Double_t xyz[3],Double_t pxpypz[3],
172  Double_t cv[21],Short_t sign) :
173  AliVTrack(),
174  fX(0.),
175  fAlpha(0.)
176 {
177  //
178  // constructor from the global parameters
179  //
180 
181  Set(xyz,pxpypz,cv,sign);
182 }
183 
184 /*
185 //_____________________________________________________________________________
186 void AliExternalTrackParam::Set(Double_t xyz[3],Double_t pxpypz[3],
187  Double_t cv[21],Short_t sign)
188 {
189  //
190  // create external track parameters from the global parameters
191  // x,y,z,px,py,pz and their 6x6 covariance matrix
192  // A.Dainese 10.10.08
193 
194  // Calculate alpha: the rotation angle of the corresponding local system.
195  //
196  // For global radial position inside the beam pipe, alpha is the
197  // azimuthal angle of the momentum projected on (x,y).
198  //
199  // For global radial position outside the ITS, alpha is the
200  // azimuthal angle of the centre of the TPC sector in which the point
201  // xyz lies
202  //
203  const double kSafe = 1e-5;
204  Double_t radPos2 = xyz[0]*xyz[0]+xyz[1]*xyz[1];
205  Double_t radMax = 45.; // approximately ITS outer radius
206  if (radPos2 < radMax*radMax) { // inside the ITS
207  fAlpha = TMath::ATan2(pxpypz[1],pxpypz[0]);
208  } else { // outside the ITS
209  Float_t phiPos = TMath::Pi()+TMath::ATan2(-xyz[1], -xyz[0]);
210  fAlpha =
211  TMath::DegToRad()*(20*((((Int_t)(phiPos*TMath::RadToDeg()))/20))+10);
212  }
213  //
214  Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
215  // protection: avoid alpha being too close to 0 or +-pi/2
216  if (TMath::Abs(sn)<2*kSafe) {
217  if (fAlpha>0) fAlpha += fAlpha< TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
218  else fAlpha += fAlpha>-TMath::Pi()/2. ? -2*kSafe : 2*kSafe;
219  cs=TMath::Cos(fAlpha);
220  sn=TMath::Sin(fAlpha);
221  }
222  else if (TMath::Abs(cs)<2*kSafe) {
223  if (fAlpha>0) fAlpha += fAlpha> TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
224  else fAlpha += fAlpha>-TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
225  cs=TMath::Cos(fAlpha);
226  sn=TMath::Sin(fAlpha);
227  }
228  // Get the vertex of origin and the momentum
229  TVector3 ver(xyz[0],xyz[1],xyz[2]);
230  TVector3 mom(pxpypz[0],pxpypz[1],pxpypz[2]);
231  //
232  // avoid momenta along axis
233  if (TMath::Abs(mom[0])<kSafe) mom[0] = TMath::Sign(kSafe*TMath::Abs(mom[1]), mom[0]);
234  if (TMath::Abs(mom[1])<kSafe) mom[1] = TMath::Sign(kSafe*TMath::Abs(mom[0]), mom[1]);
235 
236  // Rotate to the local coordinate system
237  ver.RotateZ(-fAlpha);
238  mom.RotateZ(-fAlpha);
239 
240  //
241  // x of the reference plane
242  fX = ver.X();
243 
244  Double_t charge = (Double_t)sign;
245 
246  fP[0] = ver.Y();
247  fP[1] = ver.Z();
248  fP[2] = TMath::Sin(mom.Phi());
249  fP[3] = mom.Pz()/mom.Pt();
250  fP[4] = TMath::Sign(1/mom.Pt(),charge);
251  //
252  if (TMath::Abs( 1-fP[2]) < 3*kSafe) fP[2] = 1.- 3*kSafe; //Protection
253  else if (TMath::Abs(-1-fP[2]) < 3*kSafe) fP[2] =-1.+ 3*kSafe; //Protection
254  //
255  // Covariance matrix (formulas to be simplified)
256  Double_t pt=1./TMath::Abs(fP[4]);
257  // avoid alpha+phi being to close to +-pi/2 in the cov.matrix evaluation
258  double fp2 = fP[2];
259  Double_t r=TMath::Sqrt((1.-fp2)*(1.+fp2));
260  //
261  Double_t m00=-sn;// m10=cs;
262  Double_t m23=-pt*(sn + fp2*cs/r), m43=-pt*pt*(r*cs - fp2*sn);
263  Double_t m24= pt*(cs - fp2*sn/r), m44=-pt*pt*(r*sn + fp2*cs);
264  Double_t m35=pt, m45=-pt*pt*fP[3];
265 
266  m43*=GetSign();
267  m44*=GetSign();
268  m45*=GetSign();
269 
270  Double_t cv34 = TMath::Sqrt(cv[3 ]*cv[3 ]+cv[4 ]*cv[4 ]);
271  Double_t a1=cv[13]-cv[9]*(m23*m44+m43*m24)/m23/m43;
272  Double_t a2=m23*m24-m23*(m23*m44+m43*m24)/m43;
273  Double_t a3=m43*m44-m43*(m23*m44+m43*m24)/m23;
274  Double_t a4=cv[14]+2.*cv[9]; //cv[14]-2.*cv[9]*m24*m44/m23/m43;
275  Double_t a5=m24*m24-2.*m24*m44*m23/m43;
276  Double_t a6=m44*m44-2.*m24*m44*m43/m23;
277 
278  fC[0 ] = cv[0 ]+cv[2 ];
279  fC[1 ] = TMath::Sign(cv34,cv[3 ]/m00);
280  fC[2 ] = cv[5 ];
281  fC[3 ] = (cv[10]*m43-cv[6]*m44)/(m24*m43-m23*m44)/m00;
282  fC[10] = (cv[6]/m00-fC[3 ]*m23)/m43;
283  fC[6 ] = (cv[15]/m00-fC[10]*m45)/m35;
284  fC[4 ] = (cv[12]*m43-cv[8]*m44)/(m24*m43-m23*m44);
285  fC[11] = (cv[8]-fC[4]*m23)/m43;
286  fC[7 ] = cv[17]/m35-fC[11]*m45/m35;
287  fC[5 ] = TMath::Abs((a4*a3-a6*a1)/(a5*a3-a6*a2));
288  fC[14] = TMath::Abs((a1-a2*fC[5])/a3);
289  fC[12] = (cv[9]-fC[5]*m23*m23-fC[14]*m43*m43)/m23/m43;
290  Double_t b1=cv[18]-fC[12]*m23*m45-fC[14]*m43*m45;
291  Double_t b2=m23*m35;
292  Double_t b3=m43*m35;
293  Double_t b4=cv[19]-fC[12]*m24*m45-fC[14]*m44*m45;
294  Double_t b5=m24*m35;
295  Double_t b6=m44*m35;
296  fC[8 ] = (b4-b6*b1/b3)/(b5-b6*b2/b3);
297  fC[13] = b1/b3-b2*fC[8]/b3;
298  fC[9 ] = TMath::Abs((cv[20]-fC[14]*(m45*m45)-fC[13]*2.*m35*m45)/(m35*m35));
299 
300  CheckCovariance();
301 
302  return;
303 }
304 */
305 
306 //_____________________________________________________________________________
307 void AliExternalTrackParam::Set(Double_t xyz[3],Double_t pxpypz[3],
308  Double_t cv[21],Short_t sign)
309 {
310  //
311  // create external track parameters from the global parameters
312  // x,y,z,px,py,pz and their 6x6 covariance matrix
313  // A.Dainese 10.10.08
314 
315  // Calculate alpha: the rotation angle of the corresponding local system.
316  //
317  // For global radial position inside the beam pipe, alpha is the
318  // azimuthal angle of the momentum projected on (x,y).
319  //
320  // For global radial position outside the ITS, alpha is the
321  // azimuthal angle of the centre of the TPC sector in which the point
322  // xyz lies
323  //
324  const double kSafe = 1e-5;
325  Double_t radPos2 = xyz[0]*xyz[0]+xyz[1]*xyz[1];
326  Double_t radMax = 45.; // approximately ITS outer radius
327  if (radPos2 < radMax*radMax) { // inside the ITS
328  fAlpha = TMath::ATan2(pxpypz[1],pxpypz[0]);
329  } else { // outside the ITS
330  Float_t phiPos = TMath::Pi()+TMath::ATan2(-xyz[1], -xyz[0]);
331  fAlpha =
332  TMath::DegToRad()*(20*((((Int_t)(phiPos*TMath::RadToDeg()))/20))+10);
333  }
334  //
335  Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
336  // protection: avoid alpha being too close to 0 or +-pi/2
337  if (TMath::Abs(sn)<2*kSafe) {
338  if (fAlpha>0) fAlpha += fAlpha< TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
339  else fAlpha += fAlpha>-TMath::Pi()/2. ? -2*kSafe : 2*kSafe;
340  cs=TMath::Cos(fAlpha);
341  sn=TMath::Sin(fAlpha);
342  }
343  else if (TMath::Abs(cs)<2*kSafe) {
344  if (fAlpha>0) fAlpha += fAlpha> TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
345  else fAlpha += fAlpha>-TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
346  cs=TMath::Cos(fAlpha);
347  sn=TMath::Sin(fAlpha);
348  }
349  // Get the vertex of origin and the momentum
350  TVector3 ver(xyz[0],xyz[1],xyz[2]);
351  TVector3 mom(pxpypz[0],pxpypz[1],pxpypz[2]);
352  //
353  // Rotate to the local coordinate system
354  ver.RotateZ(-fAlpha);
355  mom.RotateZ(-fAlpha);
356 
357  //
358  // x of the reference plane
359  fX = ver.X();
360 
361  Double_t charge = (Double_t)sign;
362 
363  fP[0] = ver.Y();
364  fP[1] = ver.Z();
365  fP[2] = TMath::Sin(mom.Phi());
366  fP[3] = mom.Pz()/mom.Pt();
367  fP[4] = TMath::Sign(1/mom.Pt(),charge);
368  //
369  if (TMath::Abs( 1-fP[2]) < kSafe) fP[2] = 1.- kSafe; //Protection
370  else if (TMath::Abs(-1-fP[2]) < kSafe) fP[2] =-1.+ kSafe; //Protection
371  //
372  // Covariance matrix (formulas to be simplified)
373  Double_t pt=1./TMath::Abs(fP[4]);
374  Double_t r=TMath::Sqrt((1.-fP[2])*(1.+fP[2]));
375  //
376  Double_t cv34 = TMath::Sqrt(cv[3 ]*cv[3 ]+cv[4 ]*cv[4 ]);
377  //
378  Int_t special = 0;
379  double sgcheck = r*sn + fP[2]*cs;
380  if (TMath::Abs(sgcheck)>=1-kSafe) { // special case: lab phi is +-pi/2
381  special = 1;
382  sgcheck = TMath::Sign(1.0,sgcheck);
383  }
384  else if (TMath::Abs(sgcheck)<kSafe) {
385  sgcheck = TMath::Sign(1.0,cs);
386  special = 2; // special case: lab phi is 0
387  }
388  //
389  fC[0 ] = cv[0 ]+cv[2 ];
390  fC[1 ] = TMath::Sign(cv34,-cv[3 ]*sn);
391  fC[2 ] = cv[5 ];
392  //
393  if (special==1) {
394  double pti = 1./pt;
395  double pti2 = pti*pti;
396  int q = GetSign();
397  fC[3 ] = cv[6]*pti;
398  fC[4 ] = -sgcheck*cv[8]*r*pti;
399  fC[5 ] = TMath::Abs(cv[9]*r*r*pti2);
400  fC[6 ] = (cv[10]*fP[3]-sgcheck*cv[15])*pti/r;
401  fC[7 ] = (cv[17]-sgcheck*cv[12]*fP[3])*pti;
402  fC[8 ] = (-sgcheck*cv[18]+cv[13]*fP[3])*r*pti2;
403  fC[9 ] = TMath::Abs( cv[20]-2*sgcheck*cv[19]*fP[3]+cv[14]*fP[3]*fP[3])*pti2;
404  fC[10] = cv[10]*pti2/r*q;
405  fC[11] = -sgcheck*cv[12]*pti2*q;
406  fC[12] = cv[13]*r*pti*pti2*q;
407  fC[13] = (-sgcheck*cv[19]+cv[14]*fP[3])*r*pti2*pti;
408  fC[14] = TMath::Abs(cv[14]*pti2*pti2);
409  } else if (special==2) {
410  double pti = 1./pt;
411  double pti2 = pti*pti;
412  int q = GetSign();
413  fC[3 ] = -cv[10]*pti*cs/sn;
414  fC[4 ] = cv[12]*cs*pti;
415  fC[5 ] = TMath::Abs(cv[14]*cs*cs*pti2);
416  fC[6 ] = (sgcheck*cv[6]*fP[3]-cv[15])*pti/sn;
417  fC[7 ] = (cv[17]-sgcheck*cv[8]*fP[3])*pti;
418  fC[8 ] = (cv[19]-sgcheck*cv[13]*fP[3])*cs*pti2;
419  fC[9 ] = TMath::Abs( cv[20]-2*sgcheck*cv[18]*fP[3]+cv[9]*fP[3]*fP[3])*pti2;
420  fC[10] = sgcheck*cv[6]*pti2/sn*q;
421  fC[11] = -sgcheck*cv[8]*pti2*q;
422  fC[12] = -sgcheck*cv[13]*cs*pti*pti2*q;
423  fC[13] = (-sgcheck*cv[18]+cv[9]*fP[3])*pti2*pti*q;
424  fC[14] = TMath::Abs(cv[9]*pti2*pti2);
425  }
426  else {
427  Double_t m00=-sn;// m10=cs;
428  Double_t m23=-pt*(sn + fP[2]*cs/r), m43=-pt*pt*(r*cs - fP[2]*sn);
429  Double_t m24= pt*(cs - fP[2]*sn/r), m44=-pt*pt*(r*sn + fP[2]*cs);
430  Double_t m35=pt, m45=-pt*pt*fP[3];
431  //
432  m43*=GetSign();
433  m44*=GetSign();
434  m45*=GetSign();
435  //
436  Double_t a1=cv[13]-cv[9]*(m23*m44+m43*m24)/m23/m43;
437  Double_t a2=m23*m24-m23*(m23*m44+m43*m24)/m43;
438  Double_t a3=m43*m44-m43*(m23*m44+m43*m24)/m23;
439  Double_t a4=cv[14]+2.*cv[9]; //cv[14]-2.*cv[9]*m24*m44/m23/m43;
440  Double_t a5=m24*m24-2.*m24*m44*m23/m43;
441  Double_t a6=m44*m44-2.*m24*m44*m43/m23;
442  //
443  fC[3 ] = (cv[10]*m43-cv[6]*m44)/(m24*m43-m23*m44)/m00;
444  fC[10] = (cv[6]/m00-fC[3 ]*m23)/m43;
445  fC[6 ] = (cv[15]/m00-fC[10]*m45)/m35;
446  fC[4 ] = (cv[12]*m43-cv[8]*m44)/(m24*m43-m23*m44);
447  fC[11] = (cv[8]-fC[4]*m23)/m43;
448  fC[7 ] = cv[17]/m35-fC[11]*m45/m35;
449  fC[5 ] = TMath::Abs((a4*a3-a6*a1)/(a5*a3-a6*a2));
450  fC[14] = TMath::Abs((a1-a2*fC[5])/a3);
451  fC[12] = (cv[9]-fC[5]*m23*m23-fC[14]*m43*m43)/m23/m43;
452  Double_t b1=cv[18]-fC[12]*m23*m45-fC[14]*m43*m45;
453  Double_t b2=m23*m35;
454  Double_t b3=m43*m35;
455  Double_t b4=cv[19]-fC[12]*m24*m45-fC[14]*m44*m45;
456  Double_t b5=m24*m35;
457  Double_t b6=m44*m35;
458  fC[8 ] = (b4-b6*b1/b3)/(b5-b6*b2/b3);
459  fC[13] = b1/b3-b2*fC[8]/b3;
460  fC[9 ] = TMath::Abs((cv[20]-fC[14]*(m45*m45)-fC[13]*2.*m35*m45)/(m35*m35));
461  }
462  CheckCovariance();
463 
464  return;
465 }
466 
467 //_____________________________________________________________________________
469  //
470  // Resets all the parameters to 0
471  //
472  fX=fAlpha=0.;
473  for (Int_t i = 0; i < 5; i++) fP[i] = 0;
474  for (Int_t i = 0; i < 15; i++) fC[i] = 0;
475 }
476 
477 //_____________________________________________________________________________
478 void AliExternalTrackParam::AddCovariance(const Double_t c[15]) {
479  //
480  // Add "something" to the track covarince matrix.
481  // May be needed to account for unknown mis-calibration/mis-alignment
482  //
483  fC[0] +=c[0];
484  fC[1] +=c[1]; fC[2] +=c[2];
485  fC[3] +=c[3]; fC[4] +=c[4]; fC[5] +=c[5];
486  fC[6] +=c[6]; fC[7] +=c[7]; fC[8] +=c[8]; fC[9] +=c[9];
487  fC[10]+=c[10]; fC[11]+=c[11]; fC[12]+=c[12]; fC[13]+=c[13]; fC[14]+=c[14];
488  CheckCovariance();
489 }
490 
491 
492 Double_t AliExternalTrackParam::GetP() const {
493  //---------------------------------------------------------------------
494  // This function returns the track momentum
495  // Results for (nearly) straight tracks are meaningless !
496  //---------------------------------------------------------------------
497  if (TMath::Abs(fP[4])<=kAlmost0) return kVeryBig;
498  return TMath::Sqrt(1.+ fP[3]*fP[3])/TMath::Abs(fP[4]);
499 }
500 
502  //---------------------------------------------------------------------
503  // This function returns the 1/(track momentum)
504  //---------------------------------------------------------------------
505  return TMath::Abs(fP[4])/TMath::Sqrt(1.+ fP[3]*fP[3]);
506 }
507 
508 //_______________________________________________________________________
509 Double_t AliExternalTrackParam::GetD(Double_t x,Double_t y,Double_t b) const {
510  //------------------------------------------------------------------
511  // This function calculates the transverse impact parameter
512  // with respect to a point with global coordinates (x,y)
513  // in the magnetic field "b" (kG)
514  //------------------------------------------------------------------
515  if (TMath::Abs(b) < kAlmost0Field) return GetLinearD(x,y);
516  Double_t rp4=GetC(b);
517 
518  Double_t xt=fX, yt=fP[0];
519 
520  Double_t sn=TMath::Sin(fAlpha), cs=TMath::Cos(fAlpha);
521  Double_t a = x*cs + y*sn;
522  y = -x*sn + y*cs; x=a;
523  xt-=x; yt-=y;
524 
525  sn=rp4*xt - fP[2]; cs=rp4*yt + TMath::Sqrt((1.- fP[2])*(1.+fP[2]));
526  a=2*(xt*fP[2] - yt*TMath::Sqrt((1.-fP[2])*(1.+fP[2])))-rp4*(xt*xt + yt*yt);
527  return -a/(1 + TMath::Sqrt(sn*sn + cs*cs));
528 }
529 
530 //_______________________________________________________________________
532 GetDZ(Double_t x, Double_t y, Double_t z, Double_t b, Float_t dz[2]) const {
533  //------------------------------------------------------------------
534  // This function calculates the transverse and longitudinal impact parameters
535  // with respect to a point with global coordinates (x,y)
536  // in the magnetic field "b" (kG)
537  //------------------------------------------------------------------
538  Double_t f1 = fP[2], r1 = TMath::Sqrt((1.-f1)*(1.+f1));
539  Double_t xt=fX, yt=fP[0];
540  Double_t sn=TMath::Sin(fAlpha), cs=TMath::Cos(fAlpha);
541  Double_t a = x*cs + y*sn;
542  y = -x*sn + y*cs; x=a;
543  xt-=x; yt-=y;
544 
545  Double_t rp4=GetC(b);
546  if ((TMath::Abs(b) < kAlmost0Field) || (TMath::Abs(rp4) < kAlmost0)) {
547  dz[0] = -(xt*f1 - yt*r1);
548  dz[1] = fP[1] + (dz[0]*f1 - xt)/r1*fP[3] - z;
549  return;
550  }
551 
552  sn=rp4*xt - f1; cs=rp4*yt + r1;
553  a=2*(xt*f1 - yt*r1)-rp4*(xt*xt + yt*yt);
554  Double_t rr=TMath::Sqrt(sn*sn + cs*cs);
555  dz[0] = -a/(1 + rr);
556  Double_t f2 = -sn/rr, r2 = TMath::Sqrt((1.-f2)*(1.+f2));
557  dz[1] = fP[1] + fP[3]/rp4*TMath::ASin(f2*r1 - f1*r2) - z;
558 }
559 
560 //_______________________________________________________________________
561 Double_t AliExternalTrackParam::GetLinearD(Double_t xv,Double_t yv) const {
562  //------------------------------------------------------------------
563  // This function calculates the transverse impact parameter
564  // with respect to a point with global coordinates (xv,yv)
565  // neglecting the track curvature.
566  //------------------------------------------------------------------
567  Double_t sn=TMath::Sin(fAlpha), cs=TMath::Cos(fAlpha);
568  Double_t x= xv*cs + yv*sn;
569  Double_t y=-xv*sn + yv*cs;
570 
571  Double_t d = (fX-x)*fP[2] - (fP[0]-y)*TMath::Sqrt((1.-fP[2])*(1.+fP[2]));
572 
573  return -d;
574 }
575 
577 (Double_t xOverX0, Double_t xTimesRho, Double_t mass,
578  Double_t dEdx,
579  Bool_t anglecorr) {
580  //------------------------------------------------------------------
581  // This function corrects the track parameters for the crossed material.
582  // "xOverX0" - X/X0, the thickness in units of the radiation length.
583  // "xTimesRho" - is the product length*density (g/cm^2).
584  // It should be passed as negative when propagating tracks
585  // from the intreaction point to the outside of the central barrel.
586  // "mass" - the mass of this particle (GeV/c^2). Negative mass means charge=2 particle
587  // "dEdx" - mean enery loss (GeV/(g/cm^2)
588  // "anglecorr" - switch for the angular correction
589  //------------------------------------------------------------------
590  Double_t &fP2=fP[2];
591  Double_t &fP3=fP[3];
592  Double_t &fP4=fP[4];
593 
594  Double_t &fC22=fC[5];
595  Double_t &fC33=fC[9];
596  Double_t &fC43=fC[13];
597  Double_t &fC44=fC[14];
598 
599  //Apply angle correction, if requested
600  if(anglecorr) {
601  Double_t angle=TMath::Sqrt((1.+ fP3*fP3)/((1-fP2)*(1.+fP2)));
602  xOverX0 *=angle;
603  xTimesRho *=angle;
604  }
605 
606  Double_t p=GetP();
607  if (mass<0) p += p; // q=2 particle
608  Double_t p2=p*p;
609  Double_t beta2=p2/(p2 + mass*mass);
610 
611  //Calculating the multiple scattering corrections******************
612  Double_t cC22 = 0.;
613  Double_t cC33 = 0.;
614  Double_t cC43 = 0.;
615  Double_t cC44 = 0.;
616  if (xOverX0 != 0) {
617  //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*TMath::Abs(d)*9.36*2.33;
618  Double_t theta2=0.0136*0.0136/(beta2*p2)*TMath::Abs(xOverX0);
619  if (GetUseLogTermMS()) {
620  double lt = 1+0.038*TMath::Log(TMath::Abs(xOverX0));
621  if (lt>0) theta2 *= lt*lt;
622  }
623  if (mass<0) theta2 *= 4; // q=2 particle
624  if(theta2>TMath::Pi()*TMath::Pi()) return kFALSE;
625  cC22 = theta2*((1.-fP2)*(1.+fP2))*(1. + fP3*fP3);
626  cC33 = theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
627  cC43 = theta2*fP3*fP4*(1. + fP3*fP3);
628  cC44 = theta2*fP3*fP4*fP3*fP4;
629  }
630 
631  //Calculating the energy loss corrections************************
632  Double_t cP4=1.;
633  if ((xTimesRho != 0.) && (beta2 < 1.)) {
634  Double_t dE=dEdx*xTimesRho;
635  Double_t e=TMath::Sqrt(p2 + mass*mass);
636  if ( TMath::Abs(dE) > 0.3*e ) return kFALSE; //30% energy loss is too much!
637  if ( (1.+ dE/p2*(dE + 2*e)) < 0. ) return kFALSE;
638  cP4 = 1./TMath::Sqrt(1.+ dE/p2*(dE + 2*e)); //A precise formula by Ruben !
639  if (TMath::Abs(fP4*cP4)>100.) return kFALSE; //Do not track below 10 MeV/c
640 
641 
642  // Approximate energy loss fluctuation (M.Ivanov)
643  const Double_t knst=0.07; // To be tuned.
644  Double_t sigmadE=knst*TMath::Sqrt(TMath::Abs(dE));
645  cC44 += ((sigmadE*e/p2*fP4)*(sigmadE*e/p2*fP4));
646 
647  }
648 
649  //Applying the corrections*****************************
650  fC22 += cC22;
651  fC33 += cC33;
652  fC43 += cC43;
653  fC44 += cC44;
654  fP4 *= cP4;
655 
656  CheckCovariance();
657 
658  return kTRUE;
659 }
660 
662 (Double_t xOverX0, Double_t xTimesRho, Double_t mass,
663  Bool_t anglecorr,
664  Double_t (*Bethe)(Double_t)) {
665  //------------------------------------------------------------------
666  // This function corrects the track parameters for the crossed material.
667  // "xOverX0" - X/X0, the thickness in units of the radiation length.
668  // "xTimesRho" - is the product length*density (g/cm^2).
669  // It should be passed as negative when propagating tracks
670  // from the intreaction point to the outside of the central barrel.
671  // "mass" - the mass of this particle (GeV/c^2). mass<0 means charge=2
672  // "anglecorr" - switch for the angular correction
673  // "Bethe" - function calculating the energy loss (GeV/(g/cm^2))
674  //------------------------------------------------------------------
675 
676  Double_t bg=GetP()/mass;
677  if (mass<0) {
678  if (mass<-990) {
679  AliDebug(2,Form("Mass %f corresponds to unknown PID particle",mass));
680  return kFALSE;
681  }
682  bg = -2*bg;
683  }
684  Double_t dEdx=Bethe(bg);
685  if (mass<0) dEdx *= 4;
686 
687  return CorrectForMeanMaterialdEdx(xOverX0,xTimesRho,mass,dEdx,anglecorr);
688 }
689 
691 (Double_t xOverX0, Double_t xTimesRho, Double_t mass,
692  Double_t zOverA,
693  Double_t density,
694  Double_t exEnergy,
695  Double_t jp1,
696  Double_t jp2,
697  Bool_t anglecorr) {
698  //------------------------------------------------------------------
699  // This function corrects the track parameters for the crossed material
700  // using the full Geant-like Bethe-Bloch formula parameterization
701  // "xOverX0" - X/X0, the thickness in units of the radiation length.
702  // "xTimesRho" - is the product length*density (g/cm^2).
703  // It should be passed as negative when propagating tracks
704  // from the intreaction point to the outside of the central barrel.
705  // "mass" - the mass of this particle (GeV/c^2). mass<0 means charge=2 particle
706  // "density" - mean density (g/cm^3)
707  // "zOverA" - mean Z/A
708  // "exEnergy" - mean exitation energy (GeV)
709  // "jp1" - density effect first junction point
710  // "jp2" - density effect second junction point
711  // "anglecorr" - switch for the angular correction
712  //
713  // The default values of the parameters are for silicon
714  //
715  //------------------------------------------------------------------
716 
717  Double_t bg=GetP()/mass;
718  if (mass<0) {
719  if (mass<-990) {
720  AliDebug(2,Form("Mass %f corresponds to unknown PID particle",mass));
721  return kFALSE;
722  }
723  bg = -2*bg;
724  }
725  Double_t dEdx=BetheBlochGeant(bg,density,jp1,jp2,exEnergy,zOverA);
726 
727  if (mass<0) dEdx *= 4;
728  return CorrectForMeanMaterialdEdx(xOverX0,xTimesRho,mass,dEdx,anglecorr);
729 }
730 
731 
732 
734 (Double_t d, Double_t x0, Double_t mass, Double_t (*Bethe)(Double_t)) {
735  //------------------------------------------------------------------
736  // Deprecated function !
737  // Better use CorrectForMeanMaterial instead of it.
738  //
739  // This function corrects the track parameters for the crossed material
740  // "d" - the thickness (fraction of the radiation length)
741  // It should be passed as negative when propagating tracks
742  // from the intreaction point to the outside of the central barrel.
743  // "x0" - the radiation length (g/cm^2)
744  // "mass" - the mass of this particle (GeV/c^2)
745  //------------------------------------------------------------------
746 
747  return CorrectForMeanMaterial(d,x0*d,mass,kTRUE,Bethe);
748 
749 }
750 
752  Double_t kp1,
753  Double_t kp2,
754  Double_t kp3,
755  Double_t kp4,
756  Double_t kp5) {
757  //
758  // This is the empirical ALEPH parameterization of the Bethe-Bloch formula.
759  // It is normalized to 1 at the minimum.
760  //
761  // bg - beta*gamma
762  //
763  // The default values for the kp* parameters are for ALICE TPC.
764  // The returned value is in MIP units
765  //
766 
767  Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
768 
769  Double_t aa = TMath::Power(beta,kp4);
770  Double_t bb = TMath::Power(1./bg,kp5);
771 
772  bb=TMath::Log(kp3+bb);
773 
774  return (kp2-aa-bb)*kp1/aa;
775 }
776 
778  Double_t kp0,
779  Double_t kp1,
780  Double_t kp2,
781  Double_t kp3,
782  Double_t kp4) {
783  //
784  // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
785  //
786  // bg - beta*gamma
787  // kp0 - density [g/cm^3]
788  // kp1 - density effect first junction point
789  // kp2 - density effect second junction point
790  // kp3 - mean excitation energy [GeV]
791  // kp4 - mean Z/A
792  //
793  // The default values for the kp* parameters are for silicon.
794  // The returned value is in [GeV/(g/cm^2)].
795  //
796 
797  const Double_t mK = 0.307075e-3; // [GeV*cm^2/g]
798  const Double_t me = 0.511e-3; // [GeV/c^2]
799  const Double_t rho = kp0;
800  const Double_t x0 = kp1*2.303;
801  const Double_t x1 = kp2*2.303;
802  const Double_t mI = kp3;
803  const Double_t mZA = kp4;
804  const Double_t bg2 = bg*bg;
805  const Double_t maxT= 2*me*bg2; // neglecting the electron mass
806 
807  //*** Density effect
808  Double_t d2=0.;
809  const Double_t x=TMath::Log(bg);
810  const Double_t lhwI=TMath::Log(28.816*1e-9*TMath::Sqrt(rho*mZA)/mI);
811  if (x > x1) {
812  d2 = lhwI + x - 0.5;
813  } else if (x > x0) {
814  const Double_t r=(x1-x)/(x1-x0);
815  d2 = lhwI + x - 0.5 + (0.5 - lhwI - x0)*r*r*r;
816  }
817 
818  return mK*mZA*(1+bg2)/bg2*
819  (0.5*TMath::Log(2*me*bg2*maxT/(mI*mI)) - bg2/(1+bg2) - d2);
820 }
821 
823  //------------------------------------------------------------------
824  // This is an approximation of the Bethe-Bloch formula,
825  // reasonable for solid materials.
826  // All the parameters are, in fact, for Si.
827  // The returned value is in [GeV/(g/cm^2)]
828  //------------------------------------------------------------------
829 
830  return BetheBlochGeant(bg);
831 }
832 
834  //------------------------------------------------------------------
835  // This is an approximation of the Bethe-Bloch formula,
836  // reasonable for gas materials.
837  // All the parameters are, in fact, for Ne.
838  // The returned value is in [GeV/(g/cm^2)]
839  //------------------------------------------------------------------
840 
841  const Double_t rho = 0.9e-3;
842  const Double_t x0 = 2.;
843  const Double_t x1 = 4.;
844  const Double_t mI = 140.e-9;
845  const Double_t mZA = 0.49555;
846 
847  return BetheBlochGeant(bg,rho,x0,x1,mI,mZA);
848 }
849 
850 Bool_t AliExternalTrackParam::Rotate(Double_t alpha) {
851  //------------------------------------------------------------------
852  // Transform this track to the local coord. system rotated
853  // by angle "alpha" (rad) with respect to the global coord. system.
854  //------------------------------------------------------------------
855  if (TMath::Abs(fP[2]) >= kAlmost1) {
856  AliError(Form("Precondition is not satisfied: |sin(phi)|>1 ! %f",fP[2]));
857  return kFALSE;
858  }
859 
860  if (alpha < -TMath::Pi()) alpha += 2*TMath::Pi();
861  else if (alpha >= TMath::Pi()) alpha -= 2*TMath::Pi();
862 
863  Double_t &fP0=fP[0];
864  Double_t &fP2=fP[2];
865  Double_t &fC00=fC[0];
866  Double_t &fC10=fC[1];
867  Double_t &fC20=fC[3];
868  Double_t &fC21=fC[4];
869  Double_t &fC22=fC[5];
870  Double_t &fC30=fC[6];
871  Double_t &fC32=fC[8];
872  Double_t &fC40=fC[10];
873  Double_t &fC42=fC[12];
874 
875  Double_t x=fX;
876  Double_t ca=TMath::Cos(alpha-fAlpha), sa=TMath::Sin(alpha-fAlpha);
877  Double_t sf=fP2, cf=TMath::Sqrt((1.- fP2)*(1.+fP2)); // Improve precision
878  // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
879  // direction in local frame is along the X axis
880  if ((cf*ca+sf*sa)<0) {
881  AliDebug(1,Form("Rotation failed: local cos(phi) would become %.2f",cf*ca+sf*sa));
882  return kFALSE;
883  }
884  //
885  Double_t tmp=sf*ca - cf*sa;
886 
887  if (TMath::Abs(tmp) >= kAlmost1) {
888  if (TMath::Abs(tmp) > 1.+ Double_t(FLT_EPSILON))
889  AliWarning(Form("Rotation failed ! %.10e",tmp));
890  return kFALSE;
891  }
892  fAlpha = alpha;
893  fX = x*ca + fP0*sa;
894  fP0= -x*sa + fP0*ca;
895  fP2= tmp;
896 
897  if (TMath::Abs(cf)<kAlmost0) {
898  AliError(Form("Too small cosine value %f",cf));
899  cf = kAlmost0;
900  }
901 
902  Double_t rr=(ca+sf/cf*sa);
903 
904  fC00 *= (ca*ca);
905  fC10 *= ca;
906  fC20 *= ca*rr;
907  fC21 *= rr;
908  fC22 *= rr*rr;
909  fC30 *= ca;
910  fC32 *= rr;
911  fC40 *= ca;
912  fC42 *= rr;
913 
914  CheckCovariance();
915 
916  return kTRUE;
917 }
918 
919 //______________________________________________________
921 {
922  // rotate to new frame, ignore covariance
923  if (TMath::Abs(fP[2]) >= kAlmost1) {
924  AliError(Form("Precondition is not satisfied: |sin(phi)|>1 ! %f",fP[2]));
925  return kFALSE;
926  }
927  //
928  if (alpha < -TMath::Pi()) alpha += 2*TMath::Pi();
929  else if (alpha >= TMath::Pi()) alpha -= 2*TMath::Pi();
930  //
931  Double_t &fP0=fP[0];
932  Double_t &fP2=fP[2];
933  //
934  Double_t x=fX;
935  Double_t ca=TMath::Cos(alpha-fAlpha), sa=TMath::Sin(alpha-fAlpha);
936  Double_t sf=fP2, cf=TMath::Sqrt((1.- fP2)*(1.+fP2)); // Improve precision
937  // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
938  // direction in local frame is along the X axis
939  if ((cf*ca+sf*sa)<0) {
940  AliDebug(1,Form("Rotation failed: local cos(phi) would become %.2f",cf*ca+sf*sa));
941  return kFALSE;
942  }
943  //
944  Double_t tmp=sf*ca - cf*sa;
945 
946  if (TMath::Abs(tmp) >= kAlmost1) {
947  if (TMath::Abs(tmp) > 1.+ Double_t(FLT_EPSILON))
948  AliWarning(Form("Rotation failed ! %.10e",tmp));
949  return kFALSE;
950  }
951  fAlpha = alpha;
952  fX = x*ca + fP0*sa;
953  fP0= -x*sa + fP0*ca;
954  fP2= tmp;
955  return kTRUE;
956 }
957 
959  //------------------------------------------------------------------
960  // Transform this track to the local coord. system rotated by 180 deg.
961  //------------------------------------------------------------------
962  fX = -fX;
963  fAlpha += TMath::Pi();
964  while (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
965  while (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
966  //
967  fP[0] = -fP[0];
968  //fP[2] = -fP[2];
969  fP[3] = -fP[3];
970  fP[4] = -fP[4];
971  //
972  fC[1] = -fC[1]; // since the fP1 and fP2 are not inverted, their covariances with others change sign
973  fC[3] = -fC[3];
974  fC[7] = -fC[7];
975  fC[8] = -fC[8];
976  fC[11] = -fC[11];
977  fC[12] = -fC[12];
978  //
979  return kTRUE;
980 }
981 
982 Bool_t AliExternalTrackParam::PropagateTo(Double_t xk, Double_t b) {
983  //----------------------------------------------------------------
984  // Propagate this track to the plane X=xk (cm) in the field "b" (kG)
985  //----------------------------------------------------------------
986  Double_t dx=xk-fX;
987  if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
988 
989  Double_t crv=GetC(b);
990  if (TMath::Abs(b) < kAlmost0Field) crv=0.;
991 
992  Double_t x2r = crv*dx;
993  Double_t f1=fP[2], f2=f1 + x2r;
994  if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
995  if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
996  if (TMath::Abs(fP[4])< kAlmost0) return kFALSE;
997 
998  Double_t &fP0=fP[0], &fP1=fP[1], &fP2=fP[2], &fP3=fP[3], &fP4=fP[4];
999  Double_t
1000  &fC00=fC[0],
1001  &fC10=fC[1], &fC11=fC[2],
1002  &fC20=fC[3], &fC21=fC[4], &fC22=fC[5],
1003  &fC30=fC[6], &fC31=fC[7], &fC32=fC[8], &fC33=fC[9],
1004  &fC40=fC[10], &fC41=fC[11], &fC42=fC[12], &fC43=fC[13], &fC44=fC[14];
1005 
1006  Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
1007  if (TMath::Abs(r1)<kAlmost0) return kFALSE;
1008  if (TMath::Abs(r2)<kAlmost0) return kFALSE;
1009 
1010  fX=xk;
1011  double dy2dx = (f1+f2)/(r1+r2);
1012  fP0 += dx*dy2dx;
1013  fP2 += x2r;
1014  if (TMath::Abs(x2r)<0.05) fP1 += dx*(r2 + f2*dy2dx)*fP3; // Many thanks to P.Hristov !
1015  else {
1016  // for small dx/R the linear apporximation of the arc by the segment is OK,
1017  // but at large dx/R the error is very large and leads to incorrect Z propagation
1018  // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
1019  // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
1020  // double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
1021  // double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
1022  // fP1 += rot/crv*fP3;
1023  //
1024  double rot = TMath::ASin(r1*f2 - r2*f1); // more economic version from Yura.
1025  if (f1*f1+f2*f2>1 && f1*f2<0) { // special cases of large rotations or large abs angles
1026  if (f2>0) rot = TMath::Pi() - rot; //
1027  else rot = -TMath::Pi() - rot;
1028  }
1029  fP1 += fP3/crv*rot;
1030  }
1031 
1032  //f = F - 1
1033  /*
1034  Double_t f02= dx/(r1*r1*r1); Double_t cc=crv/fP4;
1035  Double_t f04=0.5*dx*dx/(r1*r1*r1); f04*=cc;
1036  Double_t f12= dx*fP3*f1/(r1*r1*r1);
1037  Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1); f14*=cc;
1038  Double_t f13= dx/r1;
1039  Double_t f24= dx; f24*=cc;
1040  */
1041  Double_t rinv = 1./r1;
1042  Double_t r3inv = rinv*rinv*rinv;
1043  Double_t f24= x2r/fP4;
1044  Double_t f02= dx*r3inv;
1045  Double_t f04=0.5*f24*f02;
1046  Double_t f12= f02*fP3*f1;
1047  Double_t f14=0.5*f24*f02*fP3*f1;
1048  Double_t f13= dx*rinv;
1049 
1050  //b = C*ft
1051  Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
1052  Double_t b02=f24*fC40;
1053  Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
1054  Double_t b12=f24*fC41;
1055  Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
1056  Double_t b22=f24*fC42;
1057  Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
1058  Double_t b42=f24*fC44;
1059  Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
1060  Double_t b32=f24*fC43;
1061 
1062  //a = f*b = f*C*ft
1063  Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
1064  Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
1065  Double_t a22=f24*b42;
1066 
1067  //F*C*Ft = C + (b + bt + a)
1068  fC00 += b00 + b00 + a00;
1069  fC10 += b10 + b01 + a01;
1070  fC20 += b20 + b02 + a02;
1071  fC30 += b30;
1072  fC40 += b40;
1073  fC11 += b11 + b11 + a11;
1074  fC21 += b21 + b12 + a12;
1075  fC31 += b31;
1076  fC41 += b41;
1077  fC22 += b22 + b22 + a22;
1078  fC32 += b32;
1079  fC42 += b42;
1080 
1081  CheckCovariance();
1082 
1083  return kTRUE;
1084 }
1085 
1086 Bool_t AliExternalTrackParam::PropagateParamOnlyTo(Double_t xk, Double_t b) {
1087  //----------------------------------------------------------------
1088  // Propagate this track to the plane X=xk (cm) in the field "b" (kG)
1089  // Only parameters are propagated, not the matrix. To be used for small
1090  // distances only (<mm, i.e. misalignment)
1091  //----------------------------------------------------------------
1092  Double_t dx=xk-fX;
1093  if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
1094 
1095  Double_t crv=GetC(b);
1096  if (TMath::Abs(b) < kAlmost0Field) crv=0.;
1097 
1098  Double_t x2r = crv*dx;
1099  Double_t f1=fP[2], f2=f1 + x2r;
1100  if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
1101  if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
1102  if (TMath::Abs(fP[4])< kAlmost0) return kFALSE;
1103 
1104  Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
1105  if (TMath::Abs(r1)<kAlmost0) return kFALSE;
1106  if (TMath::Abs(r2)<kAlmost0) return kFALSE;
1107 
1108  fX=xk;
1109  double dy2dx = (f1+f2)/(r1+r2);
1110  fP[0] += dx*dy2dx;
1111  fP[2] += x2r;
1112  if (TMath::Abs(x2r)<0.05) fP[1] += dx*(r2 + f2*dy2dx)*fP[3]; // Many thanks to P.Hristov !
1113  else {
1114  // for small dx/R the linear apporximation of the arc by the segment is OK,
1115  // but at large dx/R the error is very large and leads to incorrect Z propagation
1116  // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
1117  // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
1118  // double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
1119  // double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
1120  // fP1 += rot/crv*fP3;
1121  //
1122  double rot = TMath::ASin(r1*f2 - r2*f1); // more economic version from Yura.
1123  if (f1*f1+f2*f2>1 && f1*f2<0) { // special cases of large rotations or large abs angles
1124  if (f2>0) rot = TMath::Pi() - rot; //
1125  else rot = -TMath::Pi() - rot;
1126  }
1127  fP[1] += fP[3]/crv*rot;
1128  }
1129  return kTRUE;
1130 }
1131 
1132 Bool_t
1133 AliExternalTrackParam::Propagate(Double_t alpha, Double_t x, Double_t b) {
1134  //------------------------------------------------------------------
1135  // Transform this track to the local coord. system rotated
1136  // by angle "alpha" (rad) with respect to the global coord. system,
1137  // and propagate this track to the plane X=xk (cm) in the field "b" (kG)
1138  //------------------------------------------------------------------
1139 
1140  //Save the parameters
1141  Double_t as=fAlpha;
1142  Double_t xs=fX;
1143  Double_t ps[5], cs[15];
1144  for (Int_t i=0; i<5; i++) ps[i]=fP[i];
1145  for (Int_t i=0; i<15; i++) cs[i]=fC[i];
1146 
1147  if (Rotate(alpha))
1148  if (PropagateTo(x,b)) return kTRUE;
1149 
1150  //Restore the parameters, if the operation failed
1151  fAlpha=as;
1152  fX=xs;
1153  for (Int_t i=0; i<5; i++) fP[i]=ps[i];
1154  for (Int_t i=0; i<15; i++) fC[i]=cs[i];
1155  return kFALSE;
1156 }
1157 
1159 (Double_t alpha, Double_t x, Double_t b[3]) {
1160  //------------------------------------------------------------------
1161  // Transform this track to the local coord. system rotated
1162  // by angle "alpha" (rad) with respect to the global coord. system,
1163  // and propagate this track to the plane X=xk (cm),
1164  // taking into account all three components of the B field, "b[3]" (kG)
1165  //------------------------------------------------------------------
1166 
1167  //Save the parameters
1168  Double_t as=fAlpha;
1169  Double_t xs=fX;
1170  Double_t ps[5], cs[15];
1171  for (Int_t i=0; i<5; i++) ps[i]=fP[i];
1172  for (Int_t i=0; i<15; i++) cs[i]=fC[i];
1173 
1174  if (Rotate(alpha))
1175  if (PropagateToBxByBz(x,b)) return kTRUE;
1176 
1177  //Restore the parameters, if the operation failed
1178  fAlpha=as;
1179  fX=xs;
1180  for (Int_t i=0; i<5; i++) fP[i]=ps[i];
1181  for (Int_t i=0; i<15; i++) fC[i]=cs[i];
1182  return kFALSE;
1183 }
1184 
1185 
1186 void AliExternalTrackParam::Propagate(Double_t len, Double_t x[3],
1187 Double_t p[3], Double_t bz) const {
1188  //+++++++++++++++++++++++++++++++++++++++++
1189  // Origin: K. Shileev (Kirill.Shileev@cern.ch)
1190  // Extrapolate track along simple helix in magnetic field
1191  // Arguments: len -distance alogn helix, [cm]
1192  // bz - mag field, [kGaus]
1193  // Returns: x and p contain extrapolated positon and momentum
1194  // The momentum returned for straight-line tracks is meaningless !
1195  //+++++++++++++++++++++++++++++++++++++++++
1196  GetXYZ(x);
1197 
1198  if (OneOverPt() < kAlmost0 || TMath::Abs(bz) < kAlmost0Field || GetC(bz) < kAlmost0){ //straight-line tracks
1199  Double_t unit[3]; GetDirection(unit);
1200  x[0]+=unit[0]*len;
1201  x[1]+=unit[1]*len;
1202  x[2]+=unit[2]*len;
1203 
1204  p[0]=unit[0]/kAlmost0;
1205  p[1]=unit[1]/kAlmost0;
1206  p[2]=unit[2]/kAlmost0;
1207  } else {
1208  GetPxPyPz(p);
1209  Double_t pp=GetP();
1210  Double_t a = -kB2C*bz*GetSign();
1211  Double_t rho = a/pp;
1212  x[0] += p[0]*TMath::Sin(rho*len)/a - p[1]*(1-TMath::Cos(rho*len))/a;
1213  x[1] += p[1]*TMath::Sin(rho*len)/a + p[0]*(1-TMath::Cos(rho*len))/a;
1214  x[2] += p[2]*len/pp;
1215 
1216  Double_t p0=p[0];
1217  p[0] = p0 *TMath::Cos(rho*len) - p[1]*TMath::Sin(rho*len);
1218  p[1] = p[1]*TMath::Cos(rho*len) + p0 *TMath::Sin(rho*len);
1219  }
1220 }
1221 
1222 Bool_t AliExternalTrackParam::Intersect(Double_t pnt[3], Double_t norm[3],
1223 Double_t bz) const {
1224  //+++++++++++++++++++++++++++++++++++++++++
1225  // Origin: K. Shileev (Kirill.Shileev@cern.ch)
1226  // Finds point of intersection (if exists) of the helix with the plane.
1227  // Stores result in fX and fP.
1228  // Arguments: planePoint,planeNorm - the plane defined by any plane's point
1229  // and vector, normal to the plane
1230  // Returns: kTrue if helix intersects the plane, kFALSE otherwise.
1231  //+++++++++++++++++++++++++++++++++++++++++
1232  Double_t x0[3]; GetXYZ(x0); //get track position in MARS
1233 
1234  //estimates initial helix length up to plane
1235  Double_t s=
1236  (pnt[0]-x0[0])*norm[0] + (pnt[1]-x0[1])*norm[1] + (pnt[2]-x0[2])*norm[2];
1237  Double_t dist=99999,distPrev=dist;
1238  Double_t x[3],p[3];
1239  while(TMath::Abs(dist)>0.00001){
1240  //calculates helix at the distance s from x0 ALONG the helix
1241  Propagate(s,x,p,bz);
1242 
1243  //distance between current helix position and plane
1244  dist=(x[0]-pnt[0])*norm[0]+(x[1]-pnt[1])*norm[1]+(x[2]-pnt[2])*norm[2];
1245 
1246  if(TMath::Abs(dist) >= TMath::Abs(distPrev)) {return kFALSE;}
1247  distPrev=dist;
1248  s-=dist;
1249  }
1250  //on exit pnt is intersection point,norm is track vector at that point,
1251  //all in MARS
1252  for (Int_t i=0; i<3; i++) {pnt[i]=x[i]; norm[i]=p[i];}
1253  return kTRUE;
1254 }
1255 
1256 Double_t
1257 AliExternalTrackParam::GetPredictedChi2(const Double_t p[2],const Double_t cov[3]) const {
1258  //----------------------------------------------------------------
1259  // Estimate the chi2 of the space point "p" with the cov. matrix "cov"
1260  //----------------------------------------------------------------
1261  Double_t sdd = fC[0] + cov[0];
1262  Double_t sdz = fC[1] + cov[1];
1263  Double_t szz = fC[2] + cov[2];
1264  Double_t det = sdd*szz - sdz*sdz;
1265 
1266  if (TMath::Abs(det) < kAlmost0) return kVeryBig;
1267 
1268  Double_t d = fP[0] - p[0];
1269  Double_t z = fP[1] - p[1];
1270 
1271  return (d*szz*d - 2*d*sdz*z + z*sdd*z)/det;
1272 }
1273 
1274 Double_t AliExternalTrackParam::
1275 GetPredictedChi2(const Double_t p[3],const Double_t covyz[3],const Double_t covxyz[3]) const {
1276  //----------------------------------------------------------------
1277  // Estimate the chi2 of the 3D space point "p" and
1278  // the full covariance matrix "covyz" and "covxyz"
1279  //
1280  // Cov(x,x) ... : covxyz[0]
1281  // Cov(y,x) ... : covxyz[1] covyz[0]
1282  // Cov(z,x) ... : covxyz[2] covyz[1] covyz[2]
1283  //----------------------------------------------------------------
1284 
1285  Double_t res[3] = {
1286  GetX() - p[0],
1287  GetY() - p[1],
1288  GetZ() - p[2]
1289  };
1290 
1291  Double_t f=GetSnp();
1292  if (TMath::Abs(f) >= kAlmost1) return kVeryBig;
1293  Double_t r=TMath::Sqrt((1.-f)*(1.+f));
1294  Double_t a=f/r, b=GetTgl()/r;
1295 
1296  Double_t s2=333.*333.; //something reasonably big (cm^2)
1297 
1298  TMatrixDSym v(3);
1299  v(0,0)= s2; v(0,1)= a*s2; v(0,2)= b*s2;;
1300  v(1,0)=a*s2; v(1,1)=a*a*s2 + GetSigmaY2(); v(1,2)=a*b*s2 + GetSigmaZY();
1301  v(2,0)=b*s2; v(2,1)=a*b*s2 + GetSigmaZY(); v(2,2)=b*b*s2 + GetSigmaZ2();
1302 
1303  v(0,0)+=covxyz[0]; v(0,1)+=covxyz[1]; v(0,2)+=covxyz[2];
1304  v(1,0)+=covxyz[1]; v(1,1)+=covyz[0]; v(1,2)+=covyz[1];
1305  v(2,0)+=covxyz[2]; v(2,1)+=covyz[1]; v(2,2)+=covyz[2];
1306 
1307  v.Invert();
1308  if (!v.IsValid()) return kVeryBig;
1309 
1310  Double_t chi2=0.;
1311  for (Int_t i = 0; i < 3; i++)
1312  for (Int_t j = 0; j < 3; j++) chi2 += res[i]*res[j]*v(i,j);
1313 
1314  return chi2;
1315 }
1316 
1317 Double_t AliExternalTrackParam::
1319  //----------------------------------------------------------------
1320  // Estimate the chi2 (5 dof) of this track with respect to the track
1321  // given by the argument.
1322  // The two tracks must be in the same reference system
1323  // and estimated at the same reference plane.
1324  //----------------------------------------------------------------
1325 
1326  if (TMath::Abs(t->GetAlpha()-GetAlpha()) > FLT_EPSILON) {
1327  AliError("The reference systems of the tracks differ !");
1328  return kVeryBig;
1329  }
1330  if (TMath::Abs(t->GetX()-GetX()) > FLT_EPSILON) {
1331  AliError("The reference of the tracks planes differ !");
1332  return kVeryBig;
1333  }
1334 
1335  TMatrixDSym c(5);
1336  c(0,0)=GetSigmaY2();
1337  c(1,0)=GetSigmaZY(); c(1,1)=GetSigmaZ2();
1338  c(2,0)=GetSigmaSnpY(); c(2,1)=GetSigmaSnpZ(); c(2,2)=GetSigmaSnp2();
1339  c(3,0)=GetSigmaTglY(); c(3,1)=GetSigmaTglZ(); c(3,2)=GetSigmaTglSnp(); c(3,3)=GetSigmaTgl2();
1340  c(4,0)=GetSigma1PtY(); c(4,1)=GetSigma1PtZ(); c(4,2)=GetSigma1PtSnp(); c(4,3)=GetSigma1PtTgl(); c(4,4)=GetSigma1Pt2();
1341 
1342  c(0,0)+=t->GetSigmaY2();
1343  c(1,0)+=t->GetSigmaZY(); c(1,1)+=t->GetSigmaZ2();
1344  c(2,0)+=t->GetSigmaSnpY();c(2,1)+=t->GetSigmaSnpZ();c(2,2)+=t->GetSigmaSnp2();
1345  c(3,0)+=t->GetSigmaTglY();c(3,1)+=t->GetSigmaTglZ();c(3,2)+=t->GetSigmaTglSnp();c(3,3)+=t->GetSigmaTgl2();
1346  c(4,0)+=t->GetSigma1PtY();c(4,1)+=t->GetSigma1PtZ();c(4,2)+=t->GetSigma1PtSnp();c(4,3)+=t->GetSigma1PtTgl();c(4,4)+=t->GetSigma1Pt2();
1347  c(0,1)=c(1,0);
1348  c(0,2)=c(2,0); c(1,2)=c(2,1);
1349  c(0,3)=c(3,0); c(1,3)=c(3,1); c(2,3)=c(3,2);
1350  c(0,4)=c(4,0); c(1,4)=c(4,1); c(2,4)=c(4,2); c(3,4)=c(4,3);
1351 
1352  c.Invert();
1353  if (!c.IsValid()) return kVeryBig;
1354 
1355 
1356  Double_t res[5] = {
1357  GetY() - t->GetY(),
1358  GetZ() - t->GetZ(),
1359  GetSnp() - t->GetSnp(),
1360  GetTgl() - t->GetTgl(),
1361  GetSigned1Pt() - t->GetSigned1Pt()
1362  };
1363 
1364  Double_t chi2=0.;
1365  for (Int_t i = 0; i < 5; i++)
1366  for (Int_t j = 0; j < 5; j++) chi2 += res[i]*res[j]*c(i,j);
1367 
1368  return chi2;
1369 }
1370 
1372 PropagateTo(Double_t p[3],Double_t covyz[3],Double_t covxyz[3],Double_t bz) {
1373  //----------------------------------------------------------------
1374  // Propagate this track to the plane
1375  // the 3D space point "p" (with the covariance matrix "covyz" and "covxyz")
1376  // belongs to.
1377  // The magnetic field is "bz" (kG)
1378  //
1379  // The track curvature and the change of the covariance matrix
1380  // of the track parameters are negleted !
1381  // (So the "step" should be small compared with 1/curvature)
1382  //----------------------------------------------------------------
1383 
1384  Double_t f=GetSnp();
1385  if (TMath::Abs(f) >= kAlmost1) return kFALSE;
1386  Double_t r=TMath::Sqrt((1.-f)*(1.+f));
1387  Double_t a=f/r, b=GetTgl()/r;
1388 
1389  Double_t s2=333.*333.; //something reasonably big (cm^2)
1390 
1391  TMatrixDSym tV(3);
1392  tV(0,0)= s2; tV(0,1)= a*s2; tV(0,2)= b*s2;
1393  tV(1,0)=a*s2; tV(1,1)=a*a*s2; tV(1,2)=a*b*s2;
1394  tV(2,0)=b*s2; tV(2,1)=a*b*s2; tV(2,2)=b*b*s2;
1395 
1396  TMatrixDSym pV(3);
1397  pV(0,0)=covxyz[0]; pV(0,1)=covxyz[1]; pV(0,2)=covxyz[2];
1398  pV(1,0)=covxyz[1]; pV(1,1)=covyz[0]; pV(1,2)=covyz[1];
1399  pV(2,0)=covxyz[2]; pV(2,1)=covyz[1]; pV(2,2)=covyz[2];
1400 
1401  TMatrixDSym tpV(tV);
1402  tpV+=pV;
1403  tpV.Invert();
1404  if (!tpV.IsValid()) return kFALSE;
1405 
1406  TMatrixDSym pW(3),tW(3);
1407  for (Int_t i=0; i<3; i++)
1408  for (Int_t j=0; j<3; j++) {
1409  pW(i,j)=tW(i,j)=0.;
1410  for (Int_t k=0; k<3; k++) {
1411  pW(i,j) += tV(i,k)*tpV(k,j);
1412  tW(i,j) += pV(i,k)*tpV(k,j);
1413  }
1414  }
1415 
1416  Double_t t[3] = {GetX(), GetY(), GetZ()};
1417 
1418  Double_t x=0.;
1419  for (Int_t i=0; i<3; i++) x += (tW(0,i)*t[i] + pW(0,i)*p[i]);
1420  Double_t crv=GetC(bz);
1421  if (TMath::Abs(b) < kAlmost0Field) crv=0.;
1422  f += crv*(x-fX);
1423  if (TMath::Abs(f) >= kAlmost1) return kFALSE;
1424  fX=x;
1425 
1426  fP[0]=0.;
1427  for (Int_t i=0; i<3; i++) fP[0] += (tW(1,i)*t[i] + pW(1,i)*p[i]);
1428  fP[1]=0.;
1429  for (Int_t i=0; i<3; i++) fP[1] += (tW(2,i)*t[i] + pW(2,i)*p[i]);
1430 
1431  return kTRUE;
1432 }
1433 
1435 Double_t *p,Double_t *cov,Bool_t updated) const {
1436  //------------------------------------------------------------------
1437  // Returns the track residuals with the space point "p" having
1438  // the covariance matrix "cov".
1439  // If "updated" is kTRUE, the track parameters expected to be updated,
1440  // otherwise they must be predicted.
1441  //------------------------------------------------------------------
1442  static Double_t res[2];
1443 
1444  Double_t r00=cov[0], r01=cov[1], r11=cov[2];
1445  if (updated) {
1446  r00-=fC[0]; r01-=fC[1]; r11-=fC[2];
1447  } else {
1448  r00+=fC[0]; r01+=fC[1]; r11+=fC[2];
1449  }
1450  Double_t det=r00*r11 - r01*r01;
1451 
1452  if (TMath::Abs(det) < kAlmost0) return 0;
1453 
1454  Double_t tmp=r00; r00=r11/det; r11=tmp/det;
1455 
1456  if (r00 < 0.) return 0;
1457  if (r11 < 0.) return 0;
1458 
1459  Double_t dy = fP[0] - p[0];
1460  Double_t dz = fP[1] - p[1];
1461 
1462  res[0]=dy*TMath::Sqrt(r00);
1463  res[1]=dz*TMath::Sqrt(r11);
1464 
1465  return res;
1466 }
1467 
1468 Bool_t AliExternalTrackParam::Update(const Double_t p[2], const Double_t cov[3]) {
1469  //------------------------------------------------------------------
1470  // Update the track parameters with the space point "p" having
1471  // the covariance matrix "cov"
1472  //------------------------------------------------------------------
1473  Double_t &fP0=fP[0], &fP1=fP[1], &fP2=fP[2], &fP3=fP[3], &fP4=fP[4];
1474  Double_t
1475  &fC00=fC[0],
1476  &fC10=fC[1], &fC11=fC[2],
1477  &fC20=fC[3], &fC21=fC[4], &fC22=fC[5],
1478  &fC30=fC[6], &fC31=fC[7], &fC32=fC[8], &fC33=fC[9],
1479  &fC40=fC[10], &fC41=fC[11], &fC42=fC[12], &fC43=fC[13], &fC44=fC[14];
1480 
1481  Double_t r00=cov[0], r01=cov[1], r11=cov[2];
1482  r00+=fC00; r01+=fC10; r11+=fC11;
1483  Double_t det=r00*r11 - r01*r01;
1484 
1485  if (TMath::Abs(det) < kAlmost0) return kFALSE;
1486 
1487 
1488  Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
1489 
1490  Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
1491  Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
1492  Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
1493  Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
1494  Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
1495 
1496  Double_t dy=p[0] - fP0, dz=p[1] - fP1;
1497  Double_t sf=fP2 + k20*dy + k21*dz;
1498  if (TMath::Abs(sf) > kAlmost1) return kFALSE;
1499 
1500  fP0 += k00*dy + k01*dz;
1501  fP1 += k10*dy + k11*dz;
1502  fP2 = sf;
1503  fP3 += k30*dy + k31*dz;
1504  fP4 += k40*dy + k41*dz;
1505 
1506  Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
1507  Double_t c12=fC21, c13=fC31, c14=fC41;
1508 
1509  fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
1510  fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
1511  fC40-=k00*c04+k01*c14;
1512 
1513  fC11-=k10*c01+k11*fC11;
1514  fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
1515  fC41-=k10*c04+k11*c14;
1516 
1517  fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
1518  fC42-=k20*c04+k21*c14;
1519 
1520  fC33-=k30*c03+k31*c13;
1521  fC43-=k30*c04+k31*c14;
1522 
1523  fC44-=k40*c04+k41*c14;
1524 
1525  CheckCovariance();
1526 
1527  return kTRUE;
1528 }
1529 
1530 void
1531 AliExternalTrackParam::GetHelixParameters(Double_t hlx[6], Double_t b) const {
1532  //--------------------------------------------------------------------
1533  // External track parameters -> helix parameters
1534  // "b" - magnetic field (kG)
1535  //--------------------------------------------------------------------
1536  Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
1537 
1538  hlx[0]=fP[0]; hlx[1]=fP[1]; hlx[2]=fP[2]; hlx[3]=fP[3];
1539 
1540  hlx[5]=fX*cs - hlx[0]*sn; // x0
1541  hlx[0]=fX*sn + hlx[0]*cs; // y0
1542 //hlx[1]= // z0
1543  hlx[2]=TMath::ASin(hlx[2]) + fAlpha; // phi0
1544 //hlx[3]= // tgl
1545  hlx[4]=GetC(b); // C
1546 }
1547 
1548 
1549 static void Evaluate(const Double_t *h, Double_t t,
1550  Double_t r[3], //radius vector
1551  Double_t g[3], //first defivatives
1552  Double_t gg[3]) //second derivatives
1553 {
1554  //--------------------------------------------------------------------
1555  // Calculate position of a point on a track and some derivatives
1556  //--------------------------------------------------------------------
1557  Double_t phase=h[4]*t+h[2];
1558  Double_t sn=TMath::Sin(phase), cs=TMath::Cos(phase);
1559 
1560  r[0] = h[5];
1561  r[1] = h[0];
1562  if (TMath::Abs(h[4])>kAlmost0) {
1563  r[0] += (sn - h[6])/h[4];
1564  r[1] -= (cs - h[7])/h[4];
1565  } else {
1566  r[0] += t*cs;
1567  r[1] -= -t*sn;
1568  }
1569  r[2] = h[1] + h[3]*t;
1570 
1571  g[0] = cs; g[1]=sn; g[2]=h[3];
1572 
1573  gg[0]=-h[4]*sn; gg[1]=h[4]*cs; gg[2]=0.;
1574 }
1575 
1577 Double_t b, Double_t &xthis, Double_t &xp) const {
1578  //------------------------------------------------------------
1579  // Returns the (weighed !) distance of closest approach between
1580  // this track and the track "p".
1581  // Other returned values:
1582  // xthis, xt - coordinates of tracks' reference planes at the DCA
1583  //-----------------------------------------------------------
1584  Double_t dy2=GetSigmaY2() + p->GetSigmaY2();
1585  Double_t dz2=GetSigmaZ2() + p->GetSigmaZ2();
1586  Double_t dx2=dy2;
1587 
1588  Double_t p1[8]; GetHelixParameters(p1,b);
1589  p1[6]=TMath::Sin(p1[2]); p1[7]=TMath::Cos(p1[2]);
1590  Double_t p2[8]; p->GetHelixParameters(p2,b);
1591  p2[6]=TMath::Sin(p2[2]); p2[7]=TMath::Cos(p2[2]);
1592 
1593 
1594  Double_t r1[3],g1[3],gg1[3]; Double_t t1=0.;
1595  Evaluate(p1,t1,r1,g1,gg1);
1596  Double_t r2[3],g2[3],gg2[3]; Double_t t2=0.;
1597  Evaluate(p2,t2,r2,g2,gg2);
1598 
1599  Double_t dx=r2[0]-r1[0], dy=r2[1]-r1[1], dz=r2[2]-r1[2];
1600  Double_t dm=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
1601 
1602  Int_t max=27;
1603  while (max--) {
1604  Double_t gt1=-(dx*g1[0]/dx2 + dy*g1[1]/dy2 + dz*g1[2]/dz2);
1605  Double_t gt2=+(dx*g2[0]/dx2 + dy*g2[1]/dy2 + dz*g2[2]/dz2);
1606  Double_t h11=(g1[0]*g1[0] - dx*gg1[0])/dx2 +
1607  (g1[1]*g1[1] - dy*gg1[1])/dy2 +
1608  (g1[2]*g1[2] - dz*gg1[2])/dz2;
1609  Double_t h22=(g2[0]*g2[0] + dx*gg2[0])/dx2 +
1610  (g2[1]*g2[1] + dy*gg2[1])/dy2 +
1611  (g2[2]*g2[2] + dz*gg2[2])/dz2;
1612  Double_t h12=-(g1[0]*g2[0]/dx2 + g1[1]*g2[1]/dy2 + g1[2]*g2[2]/dz2);
1613 
1614  Double_t det=h11*h22-h12*h12;
1615 
1616  Double_t dt1,dt2;
1617  if (TMath::Abs(det)<1.e-33) {
1618  //(quasi)singular Hessian
1619  dt1=-gt1; dt2=-gt2;
1620  } else {
1621  dt1=-(gt1*h22 - gt2*h12)/det;
1622  dt2=-(h11*gt2 - h12*gt1)/det;
1623  }
1624 
1625  if ((dt1*gt1+dt2*gt2)>0) {dt1=-dt1; dt2=-dt2;}
1626 
1627  //check delta(phase1) ?
1628  //check delta(phase2) ?
1629 
1630  if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4)
1631  if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) {
1632  if ((gt1*gt1+gt2*gt2) > 1.e-4/dy2/dy2)
1633  AliDebug(1," stopped at not a stationary point !");
1634  Double_t lmb=h11+h22; lmb=lmb-TMath::Sqrt(lmb*lmb-4*det);
1635  if (lmb < 0.)
1636  AliDebug(1," stopped at not a minimum !");
1637  break;
1638  }
1639 
1640  Double_t dd=dm;
1641  for (Int_t div=1 ; ; div*=2) {
1642  Evaluate(p1,t1+dt1,r1,g1,gg1);
1643  Evaluate(p2,t2+dt2,r2,g2,gg2);
1644  dx=r2[0]-r1[0]; dy=r2[1]-r1[1]; dz=r2[2]-r1[2];
1645  dd=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
1646  if (dd<dm) break;
1647  dt1*=0.5; dt2*=0.5;
1648  if (div>512) {
1649  AliDebug(1," overshoot !"); break;
1650  }
1651  }
1652  dm=dd;
1653 
1654  t1+=dt1;
1655  t2+=dt2;
1656 
1657  }
1658 
1659  if (max<=0) AliDebug(1," too many iterations !");
1660 
1661  Double_t cs=TMath::Cos(GetAlpha());
1662  Double_t sn=TMath::Sin(GetAlpha());
1663  xthis=r1[0]*cs + r1[1]*sn;
1664 
1665  cs=TMath::Cos(p->GetAlpha());
1666  sn=TMath::Sin(p->GetAlpha());
1667  xp=r2[0]*cs + r2[1]*sn;
1668 
1669  return TMath::Sqrt(dm*TMath::Sqrt(dy2*dz2));
1670 }
1671 
1672 Double_t AliExternalTrackParam::
1674  //--------------------------------------------------------------
1675  // Propagates this track and the argument track to the position of the
1676  // distance of closest approach.
1677  // Returns the (weighed !) distance of closest approach.
1678  //--------------------------------------------------------------
1679  Double_t xthis,xp;
1680  Double_t dca=GetDCA(p,b,xthis,xp);
1681 
1682  if (!PropagateTo(xthis,b)) {
1683  //AliWarning(" propagation failed !");
1684  return 1e+33;
1685  }
1686 
1687  if (!p->PropagateTo(xp,b)) {
1688  //AliWarning(" propagation failed !";
1689  return 1e+33;
1690  }
1691 
1692  return dca;
1693 }
1694 
1695 
1697 Double_t b, Double_t maxd, Double_t dz[2], Double_t covar[3]) {
1698  //
1699  // Propagate this track to the DCA to vertex "vtx",
1700  // if the (rough) transverse impact parameter is not bigger then "maxd".
1701  // Magnetic field is "b" (kG).
1702  //
1703  // a) The track gets extapolated to the DCA to the vertex.
1704  // b) The impact parameters and their covariance matrix are calculated.
1705  //
1706  // In the case of success, the returned value is kTRUE
1707  // (otherwise, it's kFALSE)
1708  //
1709  Double_t alpha=GetAlpha();
1710  Double_t sn=TMath::Sin(alpha), cs=TMath::Cos(alpha);
1711  Double_t x=GetX(), y=GetParameter()[0], snp=GetParameter()[2];
1712  Double_t xv= vtx->GetX()*cs + vtx->GetY()*sn;
1713  Double_t yv=-vtx->GetX()*sn + vtx->GetY()*cs, zv=vtx->GetZ();
1714  x-=xv; y-=yv;
1715 
1716  //Estimate the impact parameter neglecting the track curvature
1717  Double_t d=TMath::Abs(x*snp - y*TMath::Sqrt((1.-snp)*(1.+snp)));
1718  if (d > maxd) return kFALSE;
1719 
1720  //Propagate to the DCA
1721  Double_t crv=GetC(b);
1722  if (TMath::Abs(b) < kAlmost0Field) crv=0.;
1723 
1724  Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt((1.-snp)*(1.+snp)));
1725  sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv); cs=TMath::Sqrt((1.-sn)*(1.+sn));
1726  if (TMath::Abs(tgfv)>0.) cs = sn/tgfv;
1727  else cs=1.;
1728 
1729  x = xv*cs + yv*sn;
1730  yv=-xv*sn + yv*cs; xv=x;
1731 
1732  if (!Propagate(alpha+TMath::ASin(sn),xv,b)) return kFALSE;
1733 
1734  if (dz==0) return kTRUE;
1735  dz[0] = GetParameter()[0] - yv;
1736  dz[1] = GetParameter()[1] - zv;
1737 
1738  if (covar==0) return kTRUE;
1739  Double_t cov[6]; vtx->GetCovarianceMatrix(cov);
1740 
1741  //***** Improvements by A.Dainese
1742  alpha=GetAlpha(); sn=TMath::Sin(alpha); cs=TMath::Cos(alpha);
1743  Double_t s2ylocvtx = cov[0]*sn*sn + cov[2]*cs*cs - 2.*cov[1]*cs*sn;
1744  covar[0] = GetCovariance()[0] + s2ylocvtx; // neglecting correlations
1745  covar[1] = GetCovariance()[1]; // between (x,y) and z
1746  covar[2] = GetCovariance()[2] + cov[5]; // in vertex's covariance matrix
1747  //*****
1748 
1749  return kTRUE;
1750 }
1751 
1753 Double_t b[3], Double_t maxd, Double_t dz[2], Double_t covar[3]) {
1754  //
1755  // Propagate this track to the DCA to vertex "vtx",
1756  // if the (rough) transverse impact parameter is not bigger then "maxd".
1757  //
1758  // This function takes into account all three components of the magnetic
1759  // field given by the b[3] arument (kG)
1760  //
1761  // a) The track gets extapolated to the DCA to the vertex.
1762  // b) The impact parameters and their covariance matrix are calculated.
1763  //
1764  // In the case of success, the returned value is kTRUE
1765  // (otherwise, it's kFALSE)
1766  //
1767  Double_t alpha=GetAlpha();
1768  Double_t sn=TMath::Sin(alpha), cs=TMath::Cos(alpha);
1769  Double_t x=GetX(), y=GetParameter()[0], snp=GetParameter()[2];
1770  Double_t xv= vtx->GetX()*cs + vtx->GetY()*sn;
1771  Double_t yv=-vtx->GetX()*sn + vtx->GetY()*cs, zv=vtx->GetZ();
1772  x-=xv; y-=yv;
1773 
1774  //Estimate the impact parameter neglecting the track curvature
1775  Double_t d=TMath::Abs(x*snp - y*TMath::Sqrt((1.-snp)*(1.+snp)));
1776  if (d > maxd) return kFALSE;
1777 
1778  //Propagate to the DCA
1779  Double_t crv=GetC(b[2]);
1780  if (TMath::Abs(b[2]) < kAlmost0Field) crv=0.;
1781 
1782  Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt((1.-snp)*(1.+snp)));
1783  sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv); cs=TMath::Sqrt((1.-sn)*(1.+sn));
1784  if (TMath::Abs(tgfv)>0.) cs = sn/tgfv;
1785  else cs=1.;
1786 
1787  x = xv*cs + yv*sn;
1788  yv=-xv*sn + yv*cs; xv=x;
1789 
1790  if (!PropagateBxByBz(alpha+TMath::ASin(sn),xv,b)) return kFALSE;
1791 
1792  if (dz==0) return kTRUE;
1793  dz[0] = GetParameter()[0] - yv;
1794  dz[1] = GetParameter()[1] - zv;
1795 
1796  if (covar==0) return kTRUE;
1797  Double_t cov[6]; vtx->GetCovarianceMatrix(cov);
1798 
1799  //***** Improvements by A.Dainese
1800  alpha=GetAlpha(); sn=TMath::Sin(alpha); cs=TMath::Cos(alpha);
1801  Double_t s2ylocvtx = cov[0]*sn*sn + cov[2]*cs*cs - 2.*cov[1]*cs*sn;
1802  covar[0] = GetCovariance()[0] + s2ylocvtx; // neglecting correlations
1803  covar[1] = GetCovariance()[1]; // between (x,y) and z
1804  covar[2] = GetCovariance()[2] + cov[5]; // in vertex's covariance matrix
1805  //*****
1806 
1807  return kTRUE;
1808 }
1809 
1810 void AliExternalTrackParam::GetDirection(Double_t d[3]) const {
1811  //----------------------------------------------------------------
1812  // This function returns a unit vector along the track direction
1813  // in the global coordinate system.
1814  //----------------------------------------------------------------
1815  Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
1816  Double_t snp=fP[2];
1817  Double_t csp =TMath::Sqrt((1.-snp)*(1.+snp));
1818  Double_t norm=TMath::Sqrt(1.+ fP[3]*fP[3]);
1819  d[0]=(csp*cs - snp*sn)/norm;
1820  d[1]=(snp*cs + csp*sn)/norm;
1821  d[2]=fP[3]/norm;
1822 }
1823 
1824 Bool_t AliExternalTrackParam::GetPxPyPz(Double_t p[3]) const {
1825  //---------------------------------------------------------------------
1826  // This function returns the global track momentum components
1827  // Results for (nearly) straight tracks are meaningless !
1828  //---------------------------------------------------------------------
1829  p[0]=fP[4]; p[1]=fP[2]; p[2]=fP[3];
1830  return Local2GlobalMomentum(p,fAlpha);
1831 }
1832 
1833 Double_t AliExternalTrackParam::Px() const {
1834  //---------------------------------------------------------------------
1835  // Returns x-component of momentum
1836  // Result for (nearly) straight tracks is meaningless !
1837  //---------------------------------------------------------------------
1838 
1839  Double_t p[3]={kVeryBig,kVeryBig,kVeryBig};
1840  GetPxPyPz(p);
1841 
1842  return p[0];
1843 }
1844 
1845 Double_t AliExternalTrackParam::Py() const {
1846  //---------------------------------------------------------------------
1847  // Returns y-component of momentum
1848  // Result for (nearly) straight tracks is meaningless !
1849  //---------------------------------------------------------------------
1850 
1851  Double_t p[3]={kVeryBig,kVeryBig,kVeryBig};
1852  GetPxPyPz(p);
1853 
1854  return p[1];
1855 }
1856 
1857 Double_t AliExternalTrackParam::Xv() const {
1858  //---------------------------------------------------------------------
1859  // Returns x-component of first track point
1860  //---------------------------------------------------------------------
1861 
1862  Double_t r[3]={0.,0.,0.};
1863  GetXYZ(r);
1864 
1865  return r[0];
1866 }
1867 
1868 Double_t AliExternalTrackParam::Yv() const {
1869  //---------------------------------------------------------------------
1870  // Returns y-component of first track point
1871  //---------------------------------------------------------------------
1872 
1873  Double_t r[3]={0.,0.,0.};
1874  GetXYZ(r);
1875 
1876  return r[1];
1877 }
1878 
1880  // return theta angle of momentum
1881 
1882  return 0.5*TMath::Pi() - TMath::ATan(fP[3]);
1883 }
1884 
1885 Double_t AliExternalTrackParam::Phi() const {
1886  //---------------------------------------------------------------------
1887  // Returns the azimuthal angle of momentum
1888  // 0 <= phi < 2*pi
1889  //---------------------------------------------------------------------
1890 
1891  Double_t phi=TMath::ASin(fP[2]) + fAlpha;
1892  if (phi<0.) phi+=2.*TMath::Pi();
1893  else if (phi>=2.*TMath::Pi()) phi-=2.*TMath::Pi();
1894 
1895  return phi;
1896 }
1897 
1899  //---------------------------------------------------------------------
1900  // Returns the azimuthal angle of position
1901  // 0 <= phi < 2*pi
1902  //---------------------------------------------------------------------
1903  Double_t r[3]={0.,0.,0.};
1904  GetXYZ(r);
1905  Double_t phi=TMath::ATan2(r[1],r[0]);
1906  if (phi<0.) phi+=2.*TMath::Pi();
1907 
1908  return phi;
1909 }
1910 
1911 Double_t AliExternalTrackParam::M() const {
1912  // return particle mass
1913 
1914  // No mass information available so far.
1915  // Redifine in derived class!
1916 
1917  return -999.;
1918 }
1919 
1920 Double_t AliExternalTrackParam::E() const {
1921  // return particle energy
1922 
1923  // No PID information available so far.
1924  // Redifine in derived class!
1925 
1926  return -999.;
1927 }
1928 
1929 Double_t AliExternalTrackParam::Eta() const {
1930  // return pseudorapidity
1931 
1932  return -TMath::Log(TMath::Tan(0.5 * Theta()));
1933 }
1934 
1935 Double_t AliExternalTrackParam::Y() const {
1936  // return rapidity
1937 
1938  // No PID information available so far.
1939  // Redifine in derived class!
1940 
1941  return -999.;
1942 }
1943 
1944 Bool_t AliExternalTrackParam::GetXYZ(Double_t *r) const {
1945  //---------------------------------------------------------------------
1946  // This function returns the global track position
1947  //---------------------------------------------------------------------
1948  r[0]=fX; r[1]=fP[0]; r[2]=fP[1];
1949  return Local2GlobalPosition(r,fAlpha);
1950 }
1951 
1952 Bool_t AliExternalTrackParam::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
1953  //---------------------------------------------------------------------
1954  // This function returns the global covariance matrix of the track params
1955  //
1956  // Cov(x,x) ... : cv[0]
1957  // Cov(y,x) ... : cv[1] cv[2]
1958  // Cov(z,x) ... : cv[3] cv[4] cv[5]
1959  // Cov(px,x)... : cv[6] cv[7] cv[8] cv[9]
1960  // Cov(py,x)... : cv[10] cv[11] cv[12] cv[13] cv[14]
1961  // Cov(pz,x)... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
1962  //
1963  // Results for (nearly) straight tracks are meaningless !
1964  //---------------------------------------------------------------------
1965  if (TMath::Abs(fP[4])<=kAlmost0) {
1966  for (Int_t i=0; i<21; i++) cv[i]=0.;
1967  return kFALSE;
1968  }
1969  if (TMath::Abs(fP[2]) > kAlmost1) {
1970  for (Int_t i=0; i<21; i++) cv[i]=0.;
1971  return kFALSE;
1972  }
1973  Double_t pt=1./TMath::Abs(fP[4]);
1974  Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
1975  Double_t r=TMath::Sqrt((1.-fP[2])*(1.+fP[2]));
1976 
1977  Double_t m00=-sn, m10=cs;
1978  Double_t m23=-pt*(sn + fP[2]*cs/r), m43=-pt*pt*(r*cs - fP[2]*sn);
1979  Double_t m24= pt*(cs - fP[2]*sn/r), m44=-pt*pt*(r*sn + fP[2]*cs);
1980  Double_t m35=pt, m45=-pt*pt*fP[3];
1981 
1982  m43*=GetSign();
1983  m44*=GetSign();
1984  m45*=GetSign();
1985 
1986  cv[0 ] = fC[0]*m00*m00;
1987  cv[1 ] = fC[0]*m00*m10;
1988  cv[2 ] = fC[0]*m10*m10;
1989  cv[3 ] = fC[1]*m00;
1990  cv[4 ] = fC[1]*m10;
1991  cv[5 ] = fC[2];
1992  cv[6 ] = m00*(fC[3]*m23 + fC[10]*m43);
1993  cv[7 ] = m10*(fC[3]*m23 + fC[10]*m43);
1994  cv[8 ] = fC[4]*m23 + fC[11]*m43;
1995  cv[9 ] = m23*(fC[5]*m23 + fC[12]*m43) + m43*(fC[12]*m23 + fC[14]*m43);
1996  cv[10] = m00*(fC[3]*m24 + fC[10]*m44);
1997  cv[11] = m10*(fC[3]*m24 + fC[10]*m44);
1998  cv[12] = fC[4]*m24 + fC[11]*m44;
1999  cv[13] = m23*(fC[5]*m24 + fC[12]*m44) + m43*(fC[12]*m24 + fC[14]*m44);
2000  cv[14] = m24*(fC[5]*m24 + fC[12]*m44) + m44*(fC[12]*m24 + fC[14]*m44);
2001  cv[15] = m00*(fC[6]*m35 + fC[10]*m45);
2002  cv[16] = m10*(fC[6]*m35 + fC[10]*m45);
2003  cv[17] = fC[7]*m35 + fC[11]*m45;
2004  cv[18] = m23*(fC[8]*m35 + fC[12]*m45) + m43*(fC[13]*m35 + fC[14]*m45);
2005  cv[19] = m24*(fC[8]*m35 + fC[12]*m45) + m44*(fC[13]*m35 + fC[14]*m45);
2006  cv[20] = m35*(fC[9]*m35 + fC[13]*m45) + m45*(fC[13]*m35 + fC[14]*m45);
2007 
2008  return kTRUE;
2009 }
2010 
2011 
2012 Bool_t
2013 AliExternalTrackParam::GetPxPyPzAt(Double_t x, Double_t b, Double_t *p) const {
2014  //---------------------------------------------------------------------
2015  // This function returns the global track momentum extrapolated to
2016  // the radial position "x" (cm) in the magnetic field "b" (kG)
2017  //---------------------------------------------------------------------
2018  p[0]=fP[4];
2019  p[1]=fP[2]+(x-fX)*GetC(b);
2020  p[2]=fP[3];
2021  return Local2GlobalMomentum(p,fAlpha);
2022 }
2023 
2024 Bool_t AliExternalTrackParam::GetYZAt(Double_t x, Double_t b, Double_t *yz) const
2025 {
2026  //---------------------------------------------------------------------
2027  // This function returns the local Y,Z-coordinates of the intersection
2028  // point between this track and the reference plane "x" (cm).
2029  // Magnetic field "b" (kG)
2030  //---------------------------------------------------------------------
2031  double dx=x-fX;
2032  if(TMath::Abs(dx)<=kAlmost0) {
2033  yz[0] = fP[0];
2034  yz[1] = fP[1];
2035  return kTRUE;
2036  }
2037  double crv=GetC(b);
2038  double f1=fP[2], x2r = crv*dx, f2=f1 + x2r;
2039  if (TMath::Abs(f1) >= kAlmost1 ||
2040  TMath::Abs(f2) >= kAlmost1 ||
2041  TMath::Abs(fP[4])< kAlmost0) return kFALSE;
2042  double r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
2043  double dy2dx = (f1+f2)/(r1+r2);
2044  yz[0] = fP[0] + dx*dy2dx;
2045  if (TMath::Abs(x2r)<0.05) yz[1] = fP[1] + dx*(r2 + f2*dy2dx)*fP[3];
2046  else {
2047  double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
2048  double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
2049  yz[1] = fP[1] + rot/crv*fP[3];
2050  }
2051  return kTRUE;
2052 }
2053 
2054 
2055 Bool_t
2056 AliExternalTrackParam::GetYAt(Double_t x, Double_t b, Double_t &y) const {
2057  //---------------------------------------------------------------------
2058  // This function returns the local Y-coordinate of the intersection
2059  // point between this track and the reference plane "x" (cm).
2060  // Magnetic field "b" (kG)
2061  //---------------------------------------------------------------------
2062  Double_t dx=x-fX;
2063  if(TMath::Abs(dx)<=kAlmost0) {y=fP[0]; return kTRUE;}
2064 
2065  Double_t f1=fP[2], f2=f1 + dx*GetC(b);
2066 
2067  if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
2068  if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
2069 
2070  Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
2071  y = fP[0] + dx*(f1+f2)/(r1+r2);
2072  return kTRUE;
2073 }
2074 
2075 Bool_t
2076 AliExternalTrackParam::GetZAt(Double_t x, Double_t b, Double_t &z) const {
2077  //---------------------------------------------------------------------
2078  // This function returns the local Z-coordinate of the intersection
2079  // point between this track and the reference plane "x" (cm).
2080  // Magnetic field "b" (kG)
2081  //---------------------------------------------------------------------
2082  Double_t dx=x-fX;
2083  if(TMath::Abs(dx)<=kAlmost0) {z=fP[1]; return kTRUE;}
2084 
2085  Double_t crv=GetC(b);
2086  Double_t x2r = crv*dx;
2087  Double_t f1=fP[2], f2=f1 + x2r;
2088 
2089  if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
2090  if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
2091 
2092  Double_t r1=sqrt((1.-f1)*(1.+f1)), r2=sqrt((1.-f2)*(1.+f2));
2093  double dy2dx = (f1+f2)/(r1+r2);
2094  if (TMath::Abs(x2r)<0.05) {
2095  z = fP[1] + dx*(r2 + f2*dy2dx)*fP[3]; // Many thanks to P.Hristov !
2096  }
2097  else {
2098  // for small dx/R the linear apporximation of the arc by the segment is OK,
2099  // but at large dx/R the error is very large and leads to incorrect Z propagation
2100  // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
2101  // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
2102  // Similarly, the rotation angle in linear in dx only for dx<<R
2103  double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
2104  double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
2105  z = fP[1] + rot/crv*fP[3];
2106  }
2107  return kTRUE;
2108 }
2109 
2110 Bool_t
2111 AliExternalTrackParam::GetXYZAt(Double_t x, Double_t b, Double_t *r) const {
2112  //---------------------------------------------------------------------
2113  // This function returns the global track position extrapolated to
2114  // the radial position "x" (cm) in the magnetic field "b" (kG)
2115  //---------------------------------------------------------------------
2116  Double_t dx=x-fX;
2117  if(TMath::Abs(dx)<=kAlmost0) return GetXYZ(r);
2118 
2119  Double_t crv=GetC(b);
2120  Double_t x2r = crv*dx;
2121  Double_t f1=fP[2], f2=f1 + dx*crv;
2122 
2123  if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
2124  if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
2125 
2126  Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
2127  double dy2dx = (f1+f2)/(r1+r2);
2128  r[0] = x;
2129  r[1] = fP[0] + dx*dy2dx;
2130  if (TMath::Abs(x2r)<0.05) {
2131  r[2] = fP[1] + dx*(r2 + f2*dy2dx)*fP[3];//Thanks to Andrea & Peter
2132  }
2133  else {
2134  // for small dx/R the linear apporximation of the arc by the segment is OK,
2135  // but at large dx/R the error is very large and leads to incorrect Z propagation
2136  // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
2137  // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
2138  // Similarly, the rotation angle in linear in dx only for dx<<R
2139  double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
2140  double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
2141  r[2] = fP[1] + rot/crv*fP[3];
2142  }
2143 
2144  return Local2GlobalPosition(r,fAlpha);
2145 }
2146 
2147 //_____________________________________________________________________________
2148 void AliExternalTrackParam::Print(Option_t* /*option*/) const
2149 {
2150 // print the parameters and the covariance matrix
2151 
2152  printf("AliExternalTrackParam: x = %-12g alpha = %-12g\n", fX, fAlpha);
2153  printf(" parameters: %12g %12g %12g %12g %12g\n",
2154  fP[0], fP[1], fP[2], fP[3], fP[4]);
2155  printf(" covariance: %12g\n", fC[0]);
2156  printf(" %12g %12g\n", fC[1], fC[2]);
2157  printf(" %12g %12g %12g\n", fC[3], fC[4], fC[5]);
2158  printf(" %12g %12g %12g %12g\n",
2159  fC[6], fC[7], fC[8], fC[9]);
2160  printf(" %12g %12g %12g %12g %12g\n",
2161  fC[10], fC[11], fC[12], fC[13], fC[14]);
2162 }
2163 
2164 Double_t AliExternalTrackParam::GetSnpAt(Double_t x,Double_t b) const {
2165  //
2166  // Get sinus at given x
2167  //
2168  Double_t crv=GetC(b);
2169  if (TMath::Abs(b) < kAlmost0Field) crv=0.;
2170  Double_t dx = x-fX;
2171  Double_t res = fP[2]+dx*crv;
2172  return res;
2173 }
2174 
2175 Bool_t AliExternalTrackParam::GetDistance(AliExternalTrackParam *param2, Double_t x, Double_t dist[3], Double_t bz){
2176  //------------------------------------------------------------------------
2177  // Get the distance between two tracks at the local position x
2178  // working in the local frame of this track.
2179  // Origin : Marian.Ivanov@cern.ch
2180  //-----------------------------------------------------------------------
2181  Double_t xyz[3];
2182  Double_t xyz2[3];
2183  xyz[0]=x;
2184  if (!GetYAt(x,bz,xyz[1])) return kFALSE;
2185  if (!GetZAt(x,bz,xyz[2])) return kFALSE;
2186  //
2187  //
2188  if (TMath::Abs(GetAlpha()-param2->GetAlpha())<kAlmost0){
2189  xyz2[0]=x;
2190  if (!param2->GetYAt(x,bz,xyz2[1])) return kFALSE;
2191  if (!param2->GetZAt(x,bz,xyz2[2])) return kFALSE;
2192  }else{
2193  //
2194  Double_t xyz1[3];
2195  Double_t dfi = param2->GetAlpha()-GetAlpha();
2196  Double_t ca = TMath::Cos(dfi), sa = TMath::Sin(dfi);
2197  xyz2[0] = xyz[0]*ca+xyz[1]*sa;
2198  xyz2[1] = -xyz[0]*sa+xyz[1]*ca;
2199  //
2200  xyz1[0]=xyz2[0];
2201  if (!param2->GetYAt(xyz2[0],bz,xyz1[1])) return kFALSE;
2202  if (!param2->GetZAt(xyz2[0],bz,xyz1[2])) return kFALSE;
2203  //
2204  xyz2[0] = xyz1[0]*ca-xyz1[1]*sa;
2205  xyz2[1] = +xyz1[0]*sa+xyz1[1]*ca;
2206  xyz2[2] = xyz1[2];
2207  }
2208  dist[0] = xyz[0]-xyz2[0];
2209  dist[1] = xyz[1]-xyz2[1];
2210  dist[2] = xyz[2]-xyz2[2];
2211 
2212  return kTRUE;
2213 }
2214 
2215 
2216 //
2217 // Draw functionality.
2218 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
2219 //
2220 
2221 void AliExternalTrackParam::DrawTrack(Float_t magf, Float_t minR, Float_t maxR, Float_t stepR){
2222  //
2223  // Draw track line
2224  //
2225  if (minR>maxR) return ;
2226  if (stepR<=0) return ;
2227  Int_t npoints = TMath::Nint((maxR-minR)/stepR)+1;
2228  if (npoints<1) return;
2229  TPolyMarker3D *polymarker = new TPolyMarker3D(npoints);
2230  FillPolymarker(polymarker, magf,minR,maxR,stepR);
2231  polymarker->Draw();
2232 }
2233 
2234 //
2235 void AliExternalTrackParam::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
2236  //
2237  // Fill points in the polymarker
2238  //
2239  Int_t counter=0;
2240  for (Double_t r=minR; r<maxR; r+=stepR){
2241  Double_t point[3];
2242  GetXYZAt(r,magF,point);
2243  pol->SetPoint(counter,point[0],point[1], point[2]);
2244  // printf("xyz\t%f\t%f\t%f\n",point[0], point[1],point[2]);
2245  counter++;
2246  }
2247 }
2248 
2249 Int_t AliExternalTrackParam::GetIndex(Int_t i, Int_t j){
2250  //
2251  Int_t min = TMath::Min(i,j);
2252  Int_t max = TMath::Max(i,j);
2253 
2254  return min+(max+1)*max/2;
2255 }
2256 
2257 
2258 void AliExternalTrackParam::g3helx3(Double_t qfield,
2259  Double_t step,
2260  Double_t vect[7]) {
2261 /******************************************************************
2262  * *
2263  * GEANT3 tracking routine in a constant field oriented *
2264  * along axis 3 *
2265  * Tracking is performed with a conventional *
2266  * helix step method *
2267  * *
2268  * Authors R.Brun, M.Hansroul ********* *
2269  * Rewritten V.Perevoztchikov *
2270  * *
2271  * Rewritten in C++ by I.Belikov *
2272  * *
2273  * qfield (kG) - particle charge times magnetic field *
2274  * step (cm) - step length along the helix *
2275  * vect[7](cm,GeV/c) - input/output x, y, z, px/p, py/p ,pz/p, p *
2276  * *
2277  ******************************************************************/
2278  const Int_t ix=0, iy=1, iz=2, ipx=3, ipy=4, ipz=5, ipp=6;
2279  const Double_t kOvSqSix=TMath::Sqrt(1./6.);
2280 
2281  Double_t cosx=vect[ipx], cosy=vect[ipy], cosz=vect[ipz];
2282 
2283  Double_t rho = qfield*kB2C/vect[ipp];
2284  Double_t tet = rho*step;
2285 
2286  Double_t tsint, sintt, sint, cos1t;
2287  if (TMath::Abs(tet) > 0.03) {
2288  sint = TMath::Sin(tet);
2289  sintt = sint/tet;
2290  tsint = (tet - sint)/tet;
2291  Double_t t=TMath::Sin(0.5*tet);
2292  cos1t = 2*t*t/tet;
2293  } else {
2294  tsint = tet*tet/6.;
2295  sintt = (1.-tet*kOvSqSix)*(1.+tet*kOvSqSix); // 1.- tsint;
2296  sint = tet*sintt;
2297  cos1t = 0.5*tet;
2298  }
2299 
2300  Double_t f1 = step*sintt;
2301  Double_t f2 = step*cos1t;
2302  Double_t f3 = step*tsint*cosz;
2303  Double_t f4 = -tet*cos1t;
2304  Double_t f5 = sint;
2305 
2306  vect[ix] += f1*cosx - f2*cosy;
2307  vect[iy] += f1*cosy + f2*cosx;
2308  vect[iz] += f1*cosz + f3;
2309 
2310  vect[ipx] += f4*cosx - f5*cosy;
2311  vect[ipy] += f4*cosy + f5*cosx;
2312 
2313 }
2314 
2315 Bool_t AliExternalTrackParam::PropagateToBxByBz(Double_t xk, const Double_t b[3]) {
2316  //----------------------------------------------------------------
2317  // Extrapolate this track to the plane X=xk in the field b[].
2318  //
2319  // X [cm] is in the "tracking coordinate system" of this track.
2320  // b[]={Bx,By,Bz} [kG] is in the Global coordidate system.
2321  //----------------------------------------------------------------
2322 
2323  Double_t dx=xk-fX;
2324  if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
2325  if (TMath::Abs(fP[4])<=kAlmost0) return kFALSE;
2326  // Do not propagate tracks outside the ALICE detector
2327  if (TMath::Abs(dx)>1e5 ||
2328  TMath::Abs(GetY())>1e5 ||
2329  TMath::Abs(GetZ())>1e5) {
2330  AliWarning(Form("Anomalous track, target X:%f",xk));
2331  Print();
2332  return kFALSE;
2333  }
2334 
2335  Double_t crv=GetC(b[2]);
2336  if (TMath::Abs(b[2]) < kAlmost0Field) crv=0.;
2337 
2338  Double_t x2r = crv*dx;
2339  Double_t f1=fP[2], f2=f1 + x2r;
2340  if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
2341  if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
2342 
2343 
2344  // Estimate the covariance matrix
2345  Double_t &fP3=fP[3], &fP4=fP[4];
2346  Double_t
2347  &fC00=fC[0],
2348  &fC10=fC[1], &fC11=fC[2],
2349  &fC20=fC[3], &fC21=fC[4], &fC22=fC[5],
2350  &fC30=fC[6], &fC31=fC[7], &fC32=fC[8], &fC33=fC[9],
2351  &fC40=fC[10], &fC41=fC[11], &fC42=fC[12], &fC43=fC[13], &fC44=fC[14];
2352 
2353  Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
2354 
2355  //f = F - 1
2356  /*
2357  Double_t f02= dx/(r1*r1*r1); Double_t cc=crv/fP4;
2358  Double_t f04=0.5*dx*dx/(r1*r1*r1); f04*=cc;
2359  Double_t f12= dx*fP3*f1/(r1*r1*r1);
2360  Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1); f14*=cc;
2361  Double_t f13= dx/r1;
2362  Double_t f24= dx; f24*=cc;
2363  */
2364  Double_t rinv = 1./r1;
2365  Double_t r3inv = rinv*rinv*rinv;
2366  Double_t f24= x2r/fP4;
2367  Double_t f02= dx*r3inv;
2368  Double_t f04=0.5*f24*f02;
2369  Double_t f12= f02*fP3*f1;
2370  Double_t f14=0.5*f24*f02*fP3*f1;
2371  Double_t f13= dx*rinv;
2372 
2373  //b = C*ft
2374  Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
2375  Double_t b02=f24*fC40;
2376  Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
2377  Double_t b12=f24*fC41;
2378  Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
2379  Double_t b22=f24*fC42;
2380  Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
2381  Double_t b42=f24*fC44;
2382  Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
2383  Double_t b32=f24*fC43;
2384 
2385  //a = f*b = f*C*ft
2386  Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
2387  Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
2388  Double_t a22=f24*b42;
2389 
2390  //F*C*Ft = C + (b + bt + a)
2391  fC00 += b00 + b00 + a00;
2392  fC10 += b10 + b01 + a01;
2393  fC20 += b20 + b02 + a02;
2394  fC30 += b30;
2395  fC40 += b40;
2396  fC11 += b11 + b11 + a11;
2397  fC21 += b21 + b12 + a12;
2398  fC31 += b31;
2399  fC41 += b41;
2400  fC22 += b22 + b22 + a22;
2401  fC32 += b32;
2402  fC42 += b42;
2403 
2404  CheckCovariance();
2405 
2406  // Appoximate step length
2407  double dy2dx = (f1+f2)/(r1+r2);
2408  Double_t step = (TMath::Abs(x2r)<0.05) ? dx*TMath::Abs(r2 + f2*dy2dx) // chord
2409  : 2.*TMath::ASin(0.5*dx*TMath::Sqrt(1.+dy2dx*dy2dx)*crv)/crv; // arc
2410  step *= TMath::Sqrt(1.+ GetTgl()*GetTgl());
2411 
2412  // Get the track's (x,y,z) and (px,py,pz) in the Global System
2413  Double_t r[3]; GetXYZ(r);
2414  Double_t p[3]; GetPxPyPz(p);
2415  Double_t pp=GetP();
2416  p[0] /= pp;
2417  p[1] /= pp;
2418  p[2] /= pp;
2419 
2420 
2421  // Rotate to the system where Bx=By=0.
2422  Double_t bt=TMath::Sqrt(b[0]*b[0] + b[1]*b[1]);
2423  Double_t cosphi=1., sinphi=0.;
2424  if (bt > kAlmost0) {cosphi=b[0]/bt; sinphi=b[1]/bt;}
2425  Double_t bb=TMath::Sqrt(b[0]*b[0] + b[1]*b[1] + b[2]*b[2]);
2426  Double_t costet=1., sintet=0.;
2427  if (bb > kAlmost0) {costet=b[2]/bb; sintet=bt/bb;}
2428  Double_t vect[7];
2429 
2430  vect[0] = costet*cosphi*r[0] + costet*sinphi*r[1] - sintet*r[2];
2431  vect[1] = -sinphi*r[0] + cosphi*r[1];
2432  vect[2] = sintet*cosphi*r[0] + sintet*sinphi*r[1] + costet*r[2];
2433 
2434  vect[3] = costet*cosphi*p[0] + costet*sinphi*p[1] - sintet*p[2];
2435  vect[4] = -sinphi*p[0] + cosphi*p[1];
2436  vect[5] = sintet*cosphi*p[0] + sintet*sinphi*p[1] + costet*p[2];
2437 
2438  vect[6] = pp;
2439 
2440 
2441  // Do the helix step
2442  g3helx3(GetSign()*bb,step,vect);
2443 
2444 
2445  // Rotate back to the Global System
2446  r[0] = cosphi*costet*vect[0] - sinphi*vect[1] + cosphi*sintet*vect[2];
2447  r[1] = sinphi*costet*vect[0] + cosphi*vect[1] + sinphi*sintet*vect[2];
2448  r[2] = -sintet*vect[0] + costet*vect[2];
2449 
2450  p[0] = cosphi*costet*vect[3] - sinphi*vect[4] + cosphi*sintet*vect[5];
2451  p[1] = sinphi*costet*vect[3] + cosphi*vect[4] + sinphi*sintet*vect[5];
2452  p[2] = -sintet*vect[3] + costet*vect[5];
2453 
2454 
2455  // Rotate back to the Tracking System
2456  Double_t cosalp = TMath::Cos(fAlpha);
2457  Double_t sinalp =-TMath::Sin(fAlpha);
2458 
2459  Double_t
2460  t = cosalp*r[0] - sinalp*r[1];
2461  r[1] = sinalp*r[0] + cosalp*r[1];
2462  r[0] = t;
2463 
2464  t = cosalp*p[0] - sinalp*p[1];
2465  p[1] = sinalp*p[0] + cosalp*p[1];
2466  p[0] = t;
2467 
2468 
2469  // Do the final correcting step to the target plane (linear approximation)
2470  Double_t x=r[0], y=r[1], z=r[2];
2471  if (TMath::Abs(dx) > kAlmost0) {
2472  if (TMath::Abs(p[0]) < kAlmost0) return kFALSE;
2473  dx = xk - r[0];
2474  x += dx;
2475  y += p[1]/p[0]*dx;
2476  z += p[2]/p[0]*dx;
2477  }
2478 
2479 
2480  // Calculate the track parameters
2481  t=TMath::Sqrt(p[0]*p[0] + p[1]*p[1]);
2482  fX = x;
2483  fP[0] = y;
2484  fP[1] = z;
2485  fP[2] = p[1]/t;
2486  fP[3] = p[2]/t;
2487  fP[4] = GetSign()/(t*pp);
2488 
2489  return kTRUE;
2490 }
2491 
2492 Bool_t AliExternalTrackParam::PropagateParamOnlyBxByBzTo(Double_t xk, const Double_t b[3]) {
2493  //----------------------------------------------------------------
2494  // Extrapolate this track params (w/o cov matrix) to the plane X=xk in the field b[].
2495  //
2496  // X [cm] is in the "tracking coordinate system" of this track.
2497  // b[]={Bx,By,Bz} [kG] is in the Global coordidate system.
2498  //----------------------------------------------------------------
2499 
2500  Double_t dx=xk-fX;
2501  if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
2502  if (TMath::Abs(fP[4])<=kAlmost0) return kFALSE;
2503  // Do not propagate tracks outside the ALICE detector
2504  if (TMath::Abs(dx)>1e5 ||
2505  TMath::Abs(GetY())>1e5 ||
2506  TMath::Abs(GetZ())>1e5) {
2507  AliWarning(Form("Anomalous track, target X:%f",xk));
2508  Print();
2509  return kFALSE;
2510  }
2511 
2512  Double_t crv=GetC(b[2]);
2513  if (TMath::Abs(b[2]) < kAlmost0Field) crv=0.;
2514 
2515  Double_t x2r = crv*dx;
2516  Double_t f1=fP[2], f2=f1 + x2r;
2517  if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
2518  if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
2519  //
2520  Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
2521  //
2522  // Appoximate step length
2523  double dy2dx = (f1+f2)/(r1+r2);
2524  Double_t step = (TMath::Abs(x2r)<0.05) ? dx*TMath::Abs(r2 + f2*dy2dx) // chord
2525  : 2.*TMath::ASin(0.5*dx*TMath::Sqrt(1.+dy2dx*dy2dx)*crv)/crv; // arc
2526  step *= TMath::Sqrt(1.+ GetTgl()*GetTgl());
2527 
2528  // Get the track's (x,y,z) and (px,py,pz) in the Global System
2529  Double_t r[3]; GetXYZ(r);
2530  Double_t p[3]; GetPxPyPz(p);
2531  Double_t pp=GetP();
2532  p[0] /= pp;
2533  p[1] /= pp;
2534  p[2] /= pp;
2535 
2536  // Rotate to the system where Bx=By=0.
2537  Double_t bt=TMath::Sqrt(b[0]*b[0] + b[1]*b[1]);
2538  Double_t cosphi=1., sinphi=0.;
2539  if (bt > kAlmost0) {cosphi=b[0]/bt; sinphi=b[1]/bt;}
2540  Double_t bb=TMath::Sqrt(b[0]*b[0] + b[1]*b[1] + b[2]*b[2]);
2541  Double_t costet=1., sintet=0.;
2542  if (bb > kAlmost0) {costet=b[2]/bb; sintet=bt/bb;}
2543  Double_t vect[7];
2544 
2545  vect[0] = costet*cosphi*r[0] + costet*sinphi*r[1] - sintet*r[2];
2546  vect[1] = -sinphi*r[0] + cosphi*r[1];
2547  vect[2] = sintet*cosphi*r[0] + sintet*sinphi*r[1] + costet*r[2];
2548 
2549  vect[3] = costet*cosphi*p[0] + costet*sinphi*p[1] - sintet*p[2];
2550  vect[4] = -sinphi*p[0] + cosphi*p[1];
2551  vect[5] = sintet*cosphi*p[0] + sintet*sinphi*p[1] + costet*p[2];
2552 
2553  vect[6] = pp;
2554 
2555  // Do the helix step
2556  g3helx3(GetSign()*bb,step,vect);
2557 
2558  // Rotate back to the Global System
2559  r[0] = cosphi*costet*vect[0] - sinphi*vect[1] + cosphi*sintet*vect[2];
2560  r[1] = sinphi*costet*vect[0] + cosphi*vect[1] + sinphi*sintet*vect[2];
2561  r[2] = -sintet*vect[0] + costet*vect[2];
2562 
2563  p[0] = cosphi*costet*vect[3] - sinphi*vect[4] + cosphi*sintet*vect[5];
2564  p[1] = sinphi*costet*vect[3] + cosphi*vect[4] + sinphi*sintet*vect[5];
2565  p[2] = -sintet*vect[3] + costet*vect[5];
2566 
2567  // Rotate back to the Tracking System
2568  Double_t cosalp = TMath::Cos(fAlpha);
2569  Double_t sinalp =-TMath::Sin(fAlpha);
2570 
2571  Double_t
2572  t = cosalp*r[0] - sinalp*r[1];
2573  r[1] = sinalp*r[0] + cosalp*r[1];
2574  r[0] = t;
2575 
2576  t = cosalp*p[0] - sinalp*p[1];
2577  p[1] = sinalp*p[0] + cosalp*p[1];
2578  p[0] = t;
2579 
2580  // Do the final correcting step to the target plane (linear approximation)
2581  Double_t x=r[0], y=r[1], z=r[2];
2582  if (TMath::Abs(dx) > kAlmost0) {
2583  if (TMath::Abs(p[0]) < kAlmost0) return kFALSE;
2584  dx = xk - r[0];
2585  x += dx;
2586  y += p[1]/p[0]*dx;
2587  z += p[2]/p[0]*dx;
2588  }
2589 
2590 
2591  // Calculate the track parameters
2592  t=TMath::Sqrt(p[0]*p[0] + p[1]*p[1]);
2593  fX = x;
2594  fP[0] = y;
2595  fP[1] = z;
2596  fP[2] = p[1]/t;
2597  fP[3] = p[2]/t;
2598  fP[4] = GetSign()/(t*pp);
2599 
2600  return kTRUE;
2601 }
2602 
2603 
2604 Bool_t AliExternalTrackParam::Translate(Double_t *vTrasl,Double_t *covV){
2605  //
2606  //Translation: in the event mixing, the tracks can be shifted
2607  //of the difference among primary vertices (vTrasl) and
2608  //the covariance matrix is changed accordingly
2609  //(covV = covariance of the primary vertex).
2610  //Origin: "Romita, Rossella" <R.Romita@gsi.de>
2611  //
2612  TVector3 translation;
2613  // vTrasl coordinates in the local system
2614  translation.SetXYZ(vTrasl[0],vTrasl[1],vTrasl[2]);
2615  translation.RotateZ(-fAlpha);
2616  translation.GetXYZ(vTrasl);
2617 
2618  //compute the new x,y,z of the track
2619  Double_t newX=fX-vTrasl[0];
2620  Double_t newY=fP[0]-vTrasl[1];
2621  Double_t newZ=fP[1]-vTrasl[2];
2622 
2623  //define the new parameters
2624  Double_t newParam[5];
2625  newParam[0]=newY;
2626  newParam[1]=newZ;
2627  newParam[2]=fP[2];
2628  newParam[3]=fP[3];
2629  newParam[4]=fP[4];
2630 
2631  // recompute the covariance matrix:
2632  // 1. covV in the local system
2633  Double_t cosRot=TMath::Cos(fAlpha), sinRot=TMath::Sin(fAlpha);
2634  TMatrixD qQi(3,3);
2635  qQi(0,0) = cosRot;
2636  qQi(0,1) = sinRot;
2637  qQi(0,2) = 0.;
2638  qQi(1,0) = -sinRot;
2639  qQi(1,1) = cosRot;
2640  qQi(1,2) = 0.;
2641  qQi(2,0) = 0.;
2642  qQi(2,1) = 0.;
2643  qQi(2,2) = 1.;
2644  TMatrixD uUi(3,3);
2645  uUi(0,0) = covV[0];
2646  uUi(0,0) = covV[0];
2647  uUi(1,0) = covV[1];
2648  uUi(0,1) = covV[1];
2649  uUi(2,0) = covV[3];
2650  uUi(0,2) = covV[3];
2651  uUi(1,1) = covV[2];
2652  uUi(2,2) = covV[5];
2653  uUi(1,2) = covV[4];
2654  if(uUi.Determinant() <= 0.) {return kFALSE;}
2655  TMatrixD uUiQi(uUi,TMatrixD::kMult,qQi);
2656  TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiQi);
2657 
2658  //2. compute the new covariance matrix of the track
2659  Double_t sigmaXX=m(0,0);
2660  Double_t sigmaXZ=m(2,0);
2661  Double_t sigmaXY=m(1,0);
2662  Double_t sigmaYY=GetSigmaY2()+m(1,1);
2663  Double_t sigmaYZ=fC[1]+m(1,2);
2664  Double_t sigmaZZ=fC[2]+m(2,2);
2665  Double_t covarianceYY=sigmaYY + (-1.)*((sigmaXY*sigmaXY)/sigmaXX);
2666  Double_t covarianceYZ=sigmaYZ-(sigmaXZ*sigmaXY/sigmaXX);
2667  Double_t covarianceZZ=sigmaZZ-((sigmaXZ*sigmaXZ)/sigmaXX);
2668 
2669  Double_t newCov[15];
2670  newCov[0]=covarianceYY;
2671  newCov[1]=covarianceYZ;
2672  newCov[2]=covarianceZZ;
2673  for(Int_t i=3;i<15;i++){
2674  newCov[i]=fC[i];
2675  }
2676 
2677  // set the new parameters
2678 
2679  Set(newX,fAlpha,newParam,newCov);
2680 
2681  return kTRUE;
2682  }
2683 
2685 
2686  // This function forces the diagonal elements of the covariance matrix to be positive.
2687  // In case the diagonal element is bigger than the maximal allowed value, it is set to
2688  // the limit and the off-diagonal elements that correspond to it are set to zero.
2689 
2690  fC[0] = TMath::Abs(fC[0]);
2691  if (fC[0]>kC0max) {
2692  double scl = TMath::Sqrt(kC0max/fC[0]);
2693  fC[0] = kC0max;
2694  fC[1] *= scl;
2695  fC[3] *= scl;
2696  fC[6] *= scl;
2697  fC[10] *= scl;
2698  }
2699  fC[2] = TMath::Abs(fC[2]);
2700  if (fC[2]>kC2max) {
2701  double scl = TMath::Sqrt(kC2max/fC[2]);
2702  fC[2] = kC2max;
2703  fC[1] *= scl;
2704  fC[4] *= scl;
2705  fC[7] *= scl;
2706  fC[11] *= scl;
2707  }
2708  fC[5] = TMath::Abs(fC[5]);
2709  if (fC[5]>kC5max) {
2710  double scl = TMath::Sqrt(kC5max/fC[5]);
2711  fC[5] = kC5max;
2712  fC[3] *= scl;
2713  fC[4] *= scl;
2714  fC[8] *= scl;
2715  fC[12] *= scl;
2716  }
2717  fC[9] = TMath::Abs(fC[9]);
2718  if (fC[9]>kC9max) {
2719  double scl = TMath::Sqrt(kC9max/fC[9]);
2720  fC[9] = kC9max;
2721  fC[6] *= scl;
2722  fC[7] *= scl;
2723  fC[8] *= scl;
2724  fC[13] *= scl;
2725  }
2726  fC[14] = TMath::Abs(fC[14]);
2727  if (fC[14]>kC14max) {
2728  double scl = TMath::Sqrt(kC14max/fC[14]);
2729  fC[14] = kC14max;
2730  fC[10] *= scl;
2731  fC[11] *= scl;
2732  fC[12] *= scl;
2733  fC[13] *= scl;
2734  }
2735 
2736  // The part below is used for tests and normally is commented out
2737 // TMatrixDSym m(5);
2738 // TVectorD eig(5);
2739 
2740 // m(0,0)=fC[0];
2741 // m(1,0)=fC[1]; m(1,1)=fC[2];
2742 // m(2,0)=fC[3]; m(2,1)=fC[4]; m(2,2)=fC[5];
2743 // m(3,0)=fC[6]; m(3,1)=fC[7]; m(3,2)=fC[8]; m(3,3)=fC[9];
2744 // m(4,0)=fC[10]; m(4,1)=fC[11]; m(4,2)=fC[12]; m(4,3)=fC[13]; m(4,4)=fC[14];
2745 
2746 // m(0,1)=m(1,0);
2747 // m(0,2)=m(2,0); m(1,2)=m(2,1);
2748 // m(0,3)=m(3,0); m(1,3)=m(3,1); m(2,3)=m(3,2);
2749 // m(0,4)=m(4,0); m(1,4)=m(4,1); m(2,4)=m(4,2); m(3,4)=m(4,3);
2750 // m.EigenVectors(eig);
2751 
2752 // // assert(eig(0)>=0 && eig(1)>=0 && eig(2)>=0 && eig(3)>=0 && eig(4)>=0);
2753 // if (!(eig(0)>=0 && eig(1)>=0 && eig(2)>=0 && eig(3)>=0 && eig(4)>=0)) {
2754 // AliWarning("Negative eigenvalues of the covariance matrix!");
2755 // this->Print();
2756 // eig.Print();
2757 // }
2758 }
2759 
2760 Bool_t AliExternalTrackParam::ConstrainToVertex(const AliVVertex* vtx, Double_t b[3])
2761 {
2762  // Constrain TPC inner params constrained
2763  //
2764  if (!vtx)
2765  return kFALSE;
2766 
2767  Double_t dz[2], cov[3];
2768  if (!PropagateToDCABxByBz(vtx, b, 3, dz, cov))
2769  return kFALSE;
2770 
2771  Double_t covar[6];
2772  vtx->GetCovarianceMatrix(covar);
2773 
2774  Double_t p[2]= { fP[0] - dz[0], fP[1] - dz[1] };
2775  Double_t c[3]= { covar[2], 0., covar[5] };
2776 
2777  Double_t chi2C = GetPredictedChi2(p,c);
2778  if (chi2C>kVeryBig)
2779  return kFALSE;
2780 
2781  if (!Update(p,c))
2782  return kFALSE;
2783 
2784  return kTRUE;
2785 }
2786 
2787 //___________________________________________________________________________________________
2788 Bool_t AliExternalTrackParam::GetXatLabR(Double_t r,Double_t &x, Double_t bz, Int_t dir) const
2789 {
2790  // Get local X of the track position estimated at the radius lab radius r.
2791  // The track curvature is accounted exactly
2792  //
2793  // The flag "dir" can be used to remove the ambiguity of which intersection to take (out of 2 possible)
2794  // 0 - take the intersection closest to the current track position
2795  // >0 - go along the track (increasing fX)
2796  // <0 - go backward (decreasing fX)
2797  //
2798  const Double_t &fy=fP[0], &sn = fP[2];
2799  const double kEps = 1.e-6;
2800  //
2801  double crv = GetC(bz);
2802  if (TMath::Abs(crv)>kAlmost0) { // helix
2803  // get center of the track circle
2804  double tR = 1./crv; // track radius (for the moment signed)
2805  double cs = TMath::Sqrt((1-sn)*(1+sn));
2806  double x0 = fX - sn*tR;
2807  double y0 = fy + cs*tR;
2808  double r0 = TMath::Sqrt(x0*x0+y0*y0);
2809  // printf("Xc:%+e Yc:%+e tR:%e r0:%e\n",x0,y0,tR,r0);
2810  //
2811  if (r0<=kAlmost0) return kFALSE; // the track is concentric to circle
2812  tR = TMath::Abs(tR);
2813  double tR2r0=1.,g=0,tmp=0;
2814  if (TMath::Abs(tR-r0)>kEps) {
2815  tR2r0 = tR/r0;
2816  g = 0.5*(r*r/(r0*tR) - tR2r0 - 1./tR2r0);
2817  tmp = 1.+g*tR2r0;
2818  }
2819  else {
2820  tR2r0 = 1.0;
2821  g = 0.5*r*r/(r0*tR) - 1;
2822  tmp = 0.5*r*r/(r0*r0);
2823  }
2824  double det = (1.-g)*(1.+g);
2825  if (det<0) return kFALSE; // does not reach raduis r
2826  det = TMath::Sqrt(det);
2827  //
2828  // the intersection happens in 2 points: {x0+tR*C,y0+tR*S}
2829  // with C=f*c0+-|s0|*det and S=f*s0-+c0 sign(s0)*det
2830  // where s0 and c0 make direction for the circle center (=x0/r0 and y0/r0)
2831  //
2832  x = x0*tmp;
2833  double y = y0*tmp;
2834  if (TMath::Abs(y0)>kAlmost0) { // when y0==0 the x,y is unique
2835  double dfx = tR2r0*TMath::Abs(y0)*det;
2836  double dfy = tR2r0*x0*TMath::Sign(det,y0);
2837  if (dir==0) { // chose the one which corresponds to smallest step
2838  double delta = (x-fX)*dfx-(y-fy)*dfy; // the choice of + in C will lead to smaller step if delta<0
2839  if (delta<0) x += dfx;
2840  else x -= dfx;
2841  }
2842  else if (dir>0) { // along track direction: x must be > fX
2843  x -= dfx; // try the smallest step (dfx is positive)
2844  double dfeps = fX-x; // handle special case of very small step
2845  if (dfeps<-kEps) return kTRUE;
2846  if (TMath::Abs(dfeps)<kEps && // are we already in right r?
2847  TMath::Abs(fX*fX+fy*fy - r*r)<kEps) return fX;
2848  x += dfx+dfx;
2849  if (x-fX>0) return kTRUE;
2850  if (x-fX<-kEps) return kFALSE;
2851  x = fX; // don't move
2852  }
2853  else { // backward: x must be < fX
2854  x += dfx; // try the smallest step (dfx is positive)
2855  double dfeps = x-fX; // handle special case of very small step
2856  if (dfeps<-kEps) return kTRUE;
2857  if (TMath::Abs(dfeps)<kEps && // are we already in right r?
2858  TMath::Abs(fX*fX+fy*fy - r*r)<kEps) return fX;
2859  x-=dfx+dfx;
2860  if (x-fX<0) return kTRUE;
2861  if (x-fX>kEps) return kFALSE;
2862  x = fX; // don't move
2863  }
2864  }
2865  else { // special case: track touching the circle just in 1 point
2866  if ( (dir>0&&x<fX) || (dir<0&&x>fX) ) return kFALSE;
2867  }
2868  }
2869  else { // this is a straight track
2870  if (TMath::Abs(sn)>=kAlmost1) { // || to Y axis
2871  double det = (r-fX)*(r+fX);
2872  if (det<0) return kFALSE; // does not reach raduis r
2873  x = fX;
2874  if (dir==0) return kTRUE;
2875  det = TMath::Sqrt(det);
2876  if (dir>0) { // along the track direction
2877  if (sn>0) {if (fy>det) return kFALSE;} // track is along Y axis and above the circle
2878  else {if (fy<-det) return kFALSE;} // track is against Y axis amd belo the circle
2879  }
2880  else if(dir>0) { // agains track direction
2881  if (sn>0) {if (fy<-det) return kFALSE;} // track is along Y axis
2882  else if (fy>det) return kFALSE; // track is against Y axis
2883  }
2884  }
2885  else if (TMath::Abs(sn)<=kAlmost0) { // || to X axis
2886  double det = (r-fy)*(r+fy);
2887  if (det<0) return kFALSE; // does not reach raduis r
2888  det = TMath::Sqrt(det);
2889  if (!dir) {
2890  x = fX>0 ? det : -det; // choose the solution requiring the smalest step
2891  return kTRUE;
2892  }
2893  else if (dir>0) { // along the track direction
2894  if (fX > det) return kFALSE; // current point is in on the right from the circle
2895  else if (fX <-det) x = -det; // on the left
2896  else x = det; // within the circle
2897  }
2898  else { // against the track direction
2899  if (fX <-det) return kFALSE;
2900  else if (fX > det) x = det;
2901  else x = -det;
2902  }
2903  }
2904  else { // general case of straight line
2905  double cs = TMath::Sqrt((1-sn)*(1+sn));
2906  double xsyc = fX*sn-fy*cs;
2907  double det = (r-xsyc)*(r+xsyc);
2908  if (det<0) return kFALSE; // does not reach raduis r
2909  det = TMath::Sqrt(det);
2910  double xcys = fX*cs+fy*sn;
2911  double t = -xcys;
2912  if (dir==0) t += t>0 ? -det:det; // chose the solution requiring the smalest step
2913  else if (dir>0) { // go in increasing fX direction. ( t+-det > 0)
2914  if (t>=-det) t += -det; // take minimal step giving t>0
2915  else return kFALSE; // both solutions have negative t
2916  }
2917  else { // go in increasing fX direction. (t+-det < 0)
2918  if (t<det) t -= det; // take minimal step giving t<0
2919  else return kFALSE; // both solutions have positive t
2920  }
2921  x = fX + cs*t;
2922  }
2923  }
2924  //
2925  return kTRUE;
2926 }
2927 //_________________________________________________________
2928 Bool_t AliExternalTrackParam::GetXYZatR(Double_t xr,Double_t bz, Double_t *xyz, Double_t* alpSect) const
2929 {
2930  // This method has 3 modes of behaviour
2931  // 1) xyz[3] array is provided but alpSect pointer is 0: calculate the position of track intersection
2932  // with circle of radius xr and fill it in xyz array
2933  // 2) alpSect pointer is provided: find alpha of the sector where the track reaches local coordinate xr
2934  // Note that in this case xr is NOT the radius but the local coordinate.
2935  // If the xyz array is provided, it will be filled by track lab coordinates at local X in this sector
2936  // 3) Neither alpSect nor xyz pointers are provided: just check if the track reaches radius xr
2937  //
2938  //
2939  double crv = GetC(bz);
2940  if ( (TMath::Abs(bz))<kAlmost0Field ) crv=0.;
2941  const double &fy = fP[0];
2942  const double &fz = fP[1];
2943  const double &sn = fP[2];
2944  const double &tgl = fP[3];
2945  //
2946  // general circle parameterization:
2947  // x = (r0+tR)cos(phi0) - tR cos(t+phi0)
2948  // y = (r0+tR)sin(phi0) - tR sin(t+phi0)
2949  // where qb is the sign of the curvature, tR is the track's signed radius and r0
2950  // is the DCA of helix to origin
2951  //
2952  double tR = 1./crv; // track radius signed
2953  double cs = TMath::Sqrt((1-sn)*(1+sn));
2954  double x0 = fX - sn*tR; // helix center coordinates
2955  double y0 = fy + cs*tR;
2956  double phi0 = TMath::ATan2(y0,x0); // angle of PCA wrt to the origin
2957  if (tR<0) phi0 += TMath::Pi();
2958  if (phi0 > TMath::Pi()) phi0 -= 2.*TMath::Pi();
2959  else if (phi0 <-TMath::Pi()) phi0 += 2.*TMath::Pi();
2960  double cs0 = TMath::Cos(phi0);
2961  double sn0 = TMath::Sin(phi0);
2962  double r0 = x0*cs0 + y0*sn0 - tR; // DCA to origin
2963  double r2R = 1.+r0/tR;
2964  //
2965  //
2966  if (r2R<kAlmost0) return kFALSE; // helix is centered at the origin, no specific intersection with other concetric circle
2967  if (!xyz && !alpSect) return kTRUE;
2968  double xr2R = xr/tR;
2969  double r2Ri = 1./r2R;
2970  // the intersection cos(t) = [1 + (r0/tR+1)^2 - (r0/tR)^2]/[2(1+r0/tR)]
2971  double cosT = 0.5*(r2R + (1-xr2R*xr2R)*r2Ri);
2972  if ( TMath::Abs(cosT)>kAlmost1 ) {
2973  // printf("Does not reach : %f %f\n",r0,tR);
2974  return kFALSE; // track does not reach the radius xr
2975  }
2976  //
2977  double t = TMath::ACos(cosT);
2978  if (tR<0) t = -t;
2979  // intersection point
2980  double xyzi[3];
2981  xyzi[0] = x0 - tR*TMath::Cos(t+phi0);
2982  xyzi[1] = y0 - tR*TMath::Sin(t+phi0);
2983  if (xyz) { // if postition is requested, then z is needed:
2984  double t0 = TMath::ATan2(cs,-sn) - phi0;
2985  double z0 = fz - t0*tR*tgl;
2986  xyzi[2] = z0 + tR*t*tgl;
2987  }
2988  else xyzi[2] = 0;
2989  //
2991  //
2992  if (xyz) {
2993  xyz[0] = xyzi[0];
2994  xyz[1] = xyzi[1];
2995  xyz[2] = xyzi[2];
2996  }
2997  //
2998  if (alpSect) {
2999  double &alp = *alpSect;
3000  // determine the sector of crossing
3001  double phiPos = TMath::Pi()+TMath::ATan2(-xyzi[1],-xyzi[0]);
3002  int sect = ((Int_t)(phiPos*TMath::RadToDeg()))/20;
3003  alp = TMath::DegToRad()*(20*sect+10);
3004  double x2r,f1,f2,r1,r2,dx,dy2dx,yloc=0, ylocMax = xr*TMath::Tan(TMath::Pi()/18); // min max Y within sector at given X
3005  //
3006  while(1) {
3007  Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
3008  if ((cs*ca+sn*sa)<0) {
3009  AliDebug(1,Form("Rotation to target sector impossible: local cos(phi) would become %.2f",cs*ca+sn*sa));
3010  return kFALSE;
3011  }
3012  //
3013  f1 = sn*ca - cs*sa;
3014  if (TMath::Abs(f1) >= kAlmost1) {
3015  AliDebug(1,Form("Rotation to target sector impossible: local sin(phi) would become %.2f",f1));
3016  return kFALSE;
3017  }
3018  //
3019  double tmpX = fX*ca + fy*sa;
3020  double tmpY = -fX*sa + fy*ca;
3021  //
3022  // estimate Y at X=xr
3023  dx=xr-tmpX;
3024  x2r = crv*dx;
3025  f2=f1 + x2r;
3026  if (TMath::Abs(f2) >= kAlmost1) {
3027  AliDebug(1,Form("Propagation in target sector failed ! %.10e",f2));
3028  return kFALSE;
3029  }
3030  r1 = TMath::Sqrt((1.-f1)*(1.+f1));
3031  r2 = TMath::Sqrt((1.-f2)*(1.+f2));
3032  dy2dx = (f1+f2)/(r1+r2);
3033  yloc = tmpY + dx*dy2dx;
3034  if (yloc>ylocMax) {alp += 2*TMath::Pi()/18; sect++;}
3035  else if (yloc<-ylocMax) {alp -= 2*TMath::Pi()/18; sect--;}
3036  else break;
3037  if (alp >= TMath::Pi()) alp -= 2*TMath::Pi();
3038  else if (alp < -TMath::Pi()) alp += 2*TMath::Pi();
3039  // if (sect>=18) sect = 0;
3040  // if (sect<=0) sect = 17;
3041  }
3042  //
3043  // if alpha was requested, then recalculate the position at intersection in sector
3044  if (xyz) {
3045  xyz[0] = xr;
3046  xyz[1] = yloc;
3047  if (TMath::Abs(x2r)<0.05) xyz[2] = fz + dx*(r2 + f2*dy2dx)*tgl;
3048  else {
3049  // for small dx/R the linear apporximation of the arc by the segment is OK,
3050  // but at large dx/R the error is very large and leads to incorrect Z propagation
3051  // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
3052  // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
3053  // Similarly, the rotation angle in linear in dx only for dx<<R
3054  double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
3055  double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
3056  xyz[2] = fz + rot/crv*tgl;
3057  }
3058  Local2GlobalPosition(xyz,alp);
3059  }
3060  }
3061  return kTRUE;
3062  //
3063 }
3064 
3065 
3066 Double_t AliExternalTrackParam::GetParameterAtRadius(Double_t r, Double_t bz, Int_t parType) const
3067 {
3068  //
3069  // Get track parameters at the radius of interest.
3070  // Given function is aimed to be used to interactivelly (tree->Draw())
3071  // access track properties at different radii
3072  //
3073  // TO BE USED WITH SPECICAL CARE -
3074  // it is aimed to be used for rough calculation as constant field and
3075  // no correction for material is used
3076  //
3077  // r - radius of interest
3078  // bz - magentic field
3079  // retun values dependens on parType:
3080  // parType = 0 -gx
3081  // parType = 1 -gy
3082  // parType = 2 -gz
3083  //
3084  // parType = 3 -pgx
3085  // parType = 4 -pgy
3086  // parType = 5 -pgz
3087  //
3088  // parType = 6 - r
3089  // parType = 7 - global position phi
3090  // parType = 8 - global direction phi
3091  // parType = 9 - direction phi- positionphi
3092  if (parType<0) {
3093  parType=-1;
3094  return 0;
3095  }
3096  Double_t xyz[3];
3097  Double_t pxyz[3];
3098  Double_t localX=0;
3099  Bool_t res = GetXatLabR(r,localX,bz,1);
3100  if (!res) {
3101  parType=-1;
3102  return 0;
3103  }
3104  //
3105  // position parameters
3106  //
3107  GetXYZAt(localX,bz,xyz);
3108  if (parType<3) {
3109  return xyz[parType];
3110  }
3111 
3112  if (parType==6) return TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
3113  if (parType==7) return TMath::ATan2(xyz[1],xyz[0]);
3114  //
3115  // momenta parameters
3116  //
3117  GetPxPyPzAt(localX,bz,pxyz);
3118  if (parType==8) return TMath::ATan2(pxyz[1],pxyz[0]);
3119  if (parType==9) {
3120  Double_t diff = TMath::ATan2(pxyz[1],pxyz[0])-TMath::ATan2(xyz[1],xyz[0]);
3121  if (diff>TMath::Pi()) diff-=TMath::TwoPi();
3122  if (diff<-TMath::Pi()) diff+=TMath::TwoPi();
3123  return diff;
3124  }
3125  if (parType>=3&&parType<6) {
3126  return pxyz[parType%3];
3127  }
3128  return 0;
3129 }
3130 
3131 //_______________________________________________________________________
3132 Bool_t AliExternalTrackParam::RelateToVVertexBxByBzDCA(const AliVVertex *vtx, Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam, Double_t dz[2], Double_t dzcov[3]) {
3133  //
3134  // Try to relate the track parameters to the vertex "vtx",
3135  // if the (rough) transverse impact parameter is not bigger then "maxd".
3136  //
3137  // All three components of the magnetic field ,"b[3]" (kG),
3138  // are taken into account.
3139  //
3140  // a) The paramters are extapolated to the DCA to the vertex.
3141  // b) The impact parameters and their covariance matrix are calculated.
3142  // c) An attempt to constrain the params to the vertex is done.
3143  // The constrained params are returned via "cParam".
3144  //
3145  // In the case of success, the returned value is kTRUE
3146  // otherwise, it's kFALSE)
3147  //
3148 
3149  if (!vtx) return kFALSE;
3150 
3151  if (!PropagateToDCABxByBz(vtx, b, maxd, dz, dzcov)) return kFALSE;
3152 
3153  Double_t covar[6]; vtx->GetCovarianceMatrix(covar);
3154  Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
3155  Double_t c[3]={covar[2],0.,covar[5]};
3156 
3157  Double_t chi2=GetPredictedChi2(p,c);
3158  if (chi2>kVeryBig) return kFALSE;
3159 
3160  if (!cParam) return kTRUE;
3161 
3162  *cParam = *this;
3163  if (!cParam->Update(p,c)) return kFALSE;
3164 
3165  return kTRUE;
3166 }
3167 
Double_t GetSigmaZY() const
TBrowser b
Definition: RunAnaESD.C:12
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
Bool_t GetPxPyPzAt(Double_t x, Double_t b, Double_t p[3]) const
void GetDZ(Double_t x, Double_t y, Double_t z, Double_t b, Float_t dz[2]) const
Bool_t Propagate(Double_t alpha, Double_t x, Double_t b)
Bool_t Intersect(Double_t pnt[3], Double_t norm[3], Double_t bz) const
Double_t GetSigma1PtSnp() const
void Print(Option_t *option="") const
Double_t GetSigmaZ2() const
Bool_t Rotate(Double_t alpha)
const Double_t kC2max
virtual Short_t Charge() const =0
Double_t GetSigmaTglZ() const
Bool_t Update(const Double_t p[2], const Double_t cov[3])
void Set(T x, T alpha, const T param[5], const T covar[15])
Bool_t GetXYZatR(Double_t xr, Double_t bz, Double_t *xyz=0, Double_t *alpSect=0) const
Bool_t GetYAt(Double_t x, Double_t b, Double_t &y) const
Bool_t GetXYZ(Double_t *p) const
static Int_t GetIndex(Int_t i, Int_t j)
virtual void DrawTrack(Float_t magF, Float_t minR, Float_t maxR, Float_t stepR)
virtual void GetCovarianceMatrix(Double_t covmatrix[6]) const =0
static Double_t BetheBlochGeant(Double_t bg, Double_t kp0=2.33, Double_t kp1=0.20, Double_t kp2=3.00, Double_t kp3=173e-9, Double_t kp4=0.49848)
Double_t * GetResiduals(Double_t *p, Double_t *cov, Bool_t updated=kTRUE) const
Bool_t RelateToVVertexBxByBzDCA(const AliVVertex *vtx, Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam=NULL, Double_t dz[2]=NULL, Double_t dzcov[3]=NULL)
Float_t p[]
Definition: kNNTest.C:133
static Double_t BetheBlochSolid(Double_t bg)
virtual Double_t Y() const
Double_t GetAlpha() const
const Double_t kC14max
Double_t GetSigmaSnpY() const
Bool_t Local2GlobalMomentum(Double_t p[3], Double_t alpha) const
Bool_t PropagateToDCABxByBz(const AliVVertex *vtx, Double_t b[3], Double_t maxd, Double_t dz[2]=0, Double_t cov[3]=0)
virtual Double_t GetZ() const =0
AliTPCfastTrack * track
Bool_t GetCovarianceXYZPxPyPz(Double_t cv[21]) const
Bool_t RotateParamOnly(Double_t alpha)
virtual Double_t GetTgl() const
Bool_t PropagateParamOnlyTo(Double_t xk, Double_t b)
npoints
Definition: driftITSTPC.C:85
Double_t PropagateToDCA(AliExternalTrackParam *p, Double_t b)
const Double_t * GetParameter() const
Double_t GetLinearD(Double_t xv, Double_t yv) const
const Float_t kAlmost0
void GetHelixParameters(Double_t h[6], Double_t b) const
virtual Bool_t GetCovarianceXYZPxPyPz(Double_t cv[21]) const =0
static Double_t BetheBlochGas(Double_t bg)
virtual Double_t Px() const =0
#define AliWarning(message)
Definition: AliLog.h:541
Double_t GetSigmaTglSnp() const
static void g3helx3(Double_t qfield, Double_t step, Double_t vect[7])
void CopyFromVTrack(const AliVTrack *vTrack)
Bool_t CorrectForMeanMaterialZA(Double_t xOverX0, Double_t xTimesRho, Double_t mass, Double_t zOverA=0.49848, Double_t density=2.33, Double_t exEnergy=173e-9, Double_t jp1=0.20, Double_t jp2=3.00, Bool_t anglecorr=kFALSE)
Bool_t bg
Definition: RunAnaESD.C:6
const Double_t kB2C
Definition: AliVParticle.h:25
TCut * cP4[4]
virtual Double_t GetX() const =0
Double_t chi2
Definition: AnalyzeLaser.C:7
Bool_t GetXatLabR(Double_t r, Double_t &x, Double_t bz, Int_t dir=0) const
Double_t GetSigned1Pt() const
virtual Double_t E() const
virtual Double_t Pz() const =0
Double_t GetSigmaY2() const
Double_t GetSigma1PtY() const
Bool_t ConstrainToVertex(const AliVVertex *vtx, Double_t b[3])
Double_t GetSigma1Pt2() const
const Double_t kC9max
Double_t GetSigmaTglY() const
virtual Bool_t GetXYZ(Double_t *p) const =0
Definition: AliVTrack.cxx:40
virtual Bool_t Translate(Double_t *vTrasl, Double_t *covV)
const Double_t * GetCovariance() const
Bool_t PropagateBxByBz(Double_t alpha, Double_t x, Double_t b[3])
const Double_t kVeryBig
Bool_t GetPxPyPz(Double_t *p) const
const Double_t kAlmost0Field
Definition: AliVParticle.h:26
Double_t GetSigma1PtZ() const
Double_t GetSigma1PtTgl() const
Bool_t PropagateParamOnlyBxByBzTo(Double_t xk, const Double_t b[3])
virtual Double_t M() const
TF1 * f
Definition: interpolTest.C:21
Bool_t GetZAt(Double_t x, Double_t b, Double_t &z) const
const Double_t kC0max
Bool_t CorrectForMaterial(Double_t d, Double_t x0, Double_t mass, Double_t(*f)(Double_t)=AliExternalTrackParam::BetheBlochSolid)
#define AliDebug(logLevel, message)
Definition: AliLog.h:300
Bool_t GetDistance(AliExternalTrackParam *param2, Double_t x, Double_t dist[3], Double_t b)
Bool_t CorrectForMeanMaterial(Double_t xOverX0, Double_t xTimesRho, Double_t mass, Bool_t anglecorr=kFALSE, Double_t(*f)(Double_t)=AliExternalTrackParam::BetheBlochSolid)
AliExternalTrackParam & operator=(const AliExternalTrackParam &trkPar)
Double_t GetPredictedChi2(const Double_t p[2], const Double_t cov[3]) const
Double_t GetSigmaTgl2() const
Double_t GetParameterAtRadius(Double_t r, Double_t bz, Int_t parType) const
virtual Double_t GetC(Double_t b) const
const Double_t kMostProbablePt
void AddCovariance(const Double_t cov[15])
Bool_t PropagateToBxByBz(Double_t x, const Double_t b[3])
#define AliError(message)
Definition: AliLog.h:591
virtual void FillPolymarker(TPolyMarker3D *pol, Float_t magf, Float_t minR, Float_t maxR, Float_t stepR)
Double_t GetSigmaSnp2() const
Bool_t Local2GlobalPosition(Double_t r[3], Double_t alpha) const
Double_t GetDCA(const AliExternalTrackParam *p, Double_t b, Double_t &xthis, Double_t &xp) const
Double_t GetSnpAt(Double_t x, Double_t b) const
Bool_t PropagateTo(Double_t p[3], Double_t covyz[3], Double_t covxyz[3], Double_t b)
Bool_t CorrectForMeanMaterialdEdx(Double_t xOverX0, Double_t xTimesRho, Double_t mass, Double_t dEdx, Bool_t anglecorr=kFALSE)
void res(Char_t i)
Definition: Resolution.C:2
for(Int_t itree=0;itree< arrInputTreesDistortionCalib->GetEntriesFast();++itree)
virtual Double_t GetY() const =0
static void Evaluate(const Double_t *h, Double_t t, Double_t r[3], Double_t g[3], Double_t gg[3])
Double_t GetSigmaSnpZ() const
void GetDirection(Double_t d[3]) const
static Double_t BetheBlochAleph(Double_t bg, Double_t kp1=0.76176e-1, Double_t kp2=10.632, Double_t kp3=0.13279e-4, Double_t kp4=1.8631, Double_t kp5=1.9479)
virtual Double_t Py() const =0
AliVTrack & operator=(const AliVTrack &vTrack)
Definition: AliVTrack.cxx:32
Double_t GetD(Double_t xv, Double_t yv, Double_t b) const
class TMatrixT< Double_t > TMatrixD
Bool_t GetXYZAt(Double_t x, Double_t b, Double_t r[3]) const
Bool_t GetYZAt(Double_t x, Double_t b, Double_t *yz) const
const Double_t kAlmost1
Definition: AliVParticle.h:22
const Double_t kC5max