Initialization of data members (fMass,fLocalConvConst). Additional protection in...
[u/mrichter/AliRoot.git] / STEER / AliExternalTrackParam.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 // track parameters in "external" format                                     //
21 //                                                                           //
22 // The track parameters are:                                                //
23 // - local y coordinate                                                      //
24 // - local z coordinate                                                      //
25 // - sin of azimutal angle                                                   //
26 // - tan of dip angle                                                        //
27 // - charge/pt                                                               //
28 // The parametrisation is given at the local x coordinate fX and the         //
29 // azimuthal angle fAlpha.                                                   //
30 //                                                                           //
31 // The external parametrisation can be used to exchange track parameters     //
32 // between different detectors.                                              //
33 //                                                                           //
34 ///////////////////////////////////////////////////////////////////////////////
35
36 #include <TMath.h>
37 #include <TVector3.h>
38
39 #include "AliExternalTrackParam.h"
40 #include "AliKalmanTrack.h"
41 #include "AliTrackReference.h"
42 #include "AliLog.h"
43
44 ClassImp(AliExternalTrackParam)
45
46 const AliMagF *AliExternalTrackParam::fgkFieldMap=0;
47 Double_t AliExternalTrackParam::fgConvConst=0.;
48
49
50 //_____________________________________________________________________________
51 AliExternalTrackParam::AliExternalTrackParam() :
52   fMass(-1),
53   fX(0),
54   fAlpha(0),
55   fLocalConvConst(0)
56 {
57   //
58   // default constructor
59   //
60   for (Int_t i = 0; i < 5; i++) fParam[i] = 0;
61   for (Int_t i = 0; i < 15; i++) fCovar[i] = 0;
62 }
63
64 //_____________________________________________________________________________
65 AliExternalTrackParam::AliExternalTrackParam(Double_t x, Double_t alpha, 
66                                              const Double_t param[5], 
67                                              const Double_t covar[15]) :
68   fMass(-1),
69   fX(x),
70   fAlpha(alpha),
71   fLocalConvConst(0)
72 {
73   //
74   // create external track parameters from given arguments
75   //
76   for (Int_t i = 0; i < 5; i++) fParam[i] = param[i];
77   for (Int_t i = 0; i < 15; i++) fCovar[i] = covar[i];
78 }
79
80 //_____________________________________________________________________________
81 AliExternalTrackParam::AliExternalTrackParam(const AliKalmanTrack& track) :
82   fMass(track.GetMass()),
83   fX(0),
84   fAlpha(track.GetAlpha()),
85   fLocalConvConst(0)
86 {
87   //
88   //
89   track.GetExternalParameters(fX,fParam);
90   track.GetExternalCovariance(fCovar);
91 }
92
93
94 AliExternalTrackParam * AliExternalTrackParam::MakeTrack(const AliTrackReference *ref, Double_t mass)
95 {
96   //
97   // Make dummy track from the track reference 
98   // negative mass means opposite charge 
99   //
100   Double_t xx[5];
101   Double_t cc[15];
102   for (Int_t i=0;i<15;i++) cc[i]=0;
103   Double_t x = ref->X(), y = ref->Y(), z = ref->Z();
104   Double_t alpha = TMath::ATan2(y,x);
105   Double_t xr = TMath::Sqrt(x*x+y*y);
106   xx[0] = 0;
107   xx[1] = z;
108   xx[3] = ref->Pz()/ref->Pt();
109   Float_t b[3];
110   Float_t xyz[3]={x,y,z};
111   Float_t convConst = 0;
112   (AliKalmanTrack::GetFieldMap())->Field(xyz,b);
113   convConst=1000/0.299792458/(1e-13 - b[2]);
114   xx[4] = 1./(convConst*ref->Pt()); // curvature rpresentation
115   if (mass<0) xx[4]*=-1.;  // negative mass - negative direction
116   Double_t alphap = TMath::ATan2(ref->Py(),ref->Px())-alpha;
117   if (alphap> TMath::Pi()) alphap-=TMath::Pi();
118   if (alphap<-TMath::Pi()) alphap+=TMath::Pi();
119   xx[2] = TMath::Sin(alphap);
120   xx[4]*=convConst;   // 1/pt representation 
121   //  AliExternalTrackParam * track = new  AliExternalTrackParam(xx,cc,xr,alpha);
122   AliExternalTrackParam * track = new  AliExternalTrackParam(xr,alpha,xx,cc);
123   track->fMass = TMath::Abs(mass);
124   //track->StartTimeIntegral();  
125   track->SaveLocalConvConst(); 
126   return track;
127 }
128
129
130
131 //_____________________________________________________________________________
132 const Double_t* AliExternalTrackParam::GetParameter() const
133 {
134 // get a pointer to the array of track parameters
135
136   return fParam;
137 }
138
139 //_____________________________________________________________________________
140 const Double_t* AliExternalTrackParam::GetCovariance() const
141 {
142 // get a pointer to the array of the track parameter covariance matrix
143
144   return fCovar;
145 }
146
147 //_____________________________________________________________________________
148 AliExternalTrackParam* AliExternalTrackParam::CreateExternalParam() const
149 {
150 // copy this instance
151
152   return new AliExternalTrackParam(fX, fAlpha, fParam, fCovar);
153 }
154
155 //_____________________________________________________________________________
156 void AliExternalTrackParam::ResetCovariance(Double_t factor,
157                                             Bool_t clearOffDiagonal)
158 {
159 // reset the covariance matrix ("forget" track history)
160
161   Int_t k = 0;
162   for (Int_t i = 0; i < 5; i++) {
163     for (Int_t j = 0; j < i; j++) {  // off diagonal elements
164       if (clearOffDiagonal) {
165         fCovar[k++] = 0;
166       } else {
167         fCovar[k++] *= factor;
168       }
169     }
170     fCovar[k++] *= factor;     // diagonal elements
171   }
172 }
173
174
175 //_____________________________________________________________________________
176 Bool_t AliExternalTrackParam::PropagateTo(Double_t xk, Double_t x0, Double_t rho)
177 {
178   //
179   // Propagate the track parameters to the given x coordinate assuming vacuum.
180   // If length is not NULL, the change of track length is added to it.
181   //
182   
183   Double_t lcc=GetLocalConvConst();  
184   Double_t cur = fParam[4]/lcc;
185   Double_t x1=fX, x2=xk, dx=x2-x1;
186   Double_t f1=fParam[2], f2=f1 + cur*dx;
187   if (TMath::Abs(f2) >= 0.98) {
188     // MI change  - don't propagate highly inclined tracks
189     //              covariance matrix distorted
190     return kFALSE;
191   }
192
193   // old position [SR, GSI, 17.02.2003]
194   Double_t oldX = fX, oldY = fParam[0], oldZ = fParam[1];
195   Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);  
196   fParam[0] += dx*(f1+f2)/(r1+r2);
197   fParam[1] += dx*(f1+f2)/(f1*r2 + f2*r1)*fParam[3];
198   fParam[2] += dx*cur;
199   // transform error matrix to the curvature
200   fCovar[10]/=lcc;
201   fCovar[11]/=lcc;
202   fCovar[12]/=lcc;
203   fCovar[13]/=lcc;
204   fCovar[14]/=lcc*lcc;
205
206   //f = F - 1
207   
208   Double_t f02=    dx/(r1*r1*r1);
209   Double_t f04=0.5*dx*dx/(r1*r1*r1);
210   Double_t f12=    dx*fParam[3]*f1/(r1*r1*r1);
211   Double_t f14=0.5*dx*dx*fParam[3]*f1/(r1*r1*r1);
212   Double_t f13=    dx/r1;
213   Double_t f24=    dx; 
214   
215   //b = C*ft
216   Double_t b00=f02*fCovar[3] + f04*fCovar[10], b01=f12*fCovar[3] + f14*fCovar[10] + f13*fCovar[6];
217   Double_t b02=f24*fCovar[10];
218   Double_t b10=f02*fCovar[4] + f04*fCovar[11], b11=f12*fCovar[4] + f14*fCovar[11] + f13*fCovar[7];
219   Double_t b12=f24*fCovar[11];
220   Double_t b20=f02*fCovar[5] + f04*fCovar[12], b21=f12*fCovar[5] + f14*fCovar[12] + f13*fCovar[8];
221   Double_t b22=f24*fCovar[12];
222   Double_t b40=f02*fCovar[12] + f04*fCovar[14], b41=f12*fCovar[12] + f14*fCovar[14] + f13*fCovar[13];
223   Double_t b42=f24*fCovar[14];
224   Double_t b30=f02*fCovar[8] + f04*fCovar[13], b31=f12*fCovar[8] + f14*fCovar[13] + f13*fCovar[9];
225   Double_t b32=f24*fCovar[13];
226   
227   //a = f*b = f*C*ft
228   Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
229   Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
230   Double_t a22=f24*b42;
231
232   //F*C*Ft = C + (b + bt + a)
233   fCovar[0] += b00 + b00 + a00;
234   fCovar[1] += b10 + b01 + a01; 
235   fCovar[2] += b11 + b11 + a11;
236   fCovar[3] += b20 + b02 + a02;
237   fCovar[4] += b21 + b12 + a12;
238   fCovar[5] += b22 + b22 + a22;
239   fCovar[6] += b30;
240   fCovar[7] += b31; 
241   fCovar[8] += b32;
242   fCovar[10] += b40;
243   fCovar[11] += b41;
244   fCovar[12] += b42;
245
246   fX=x2;
247
248   //Change of the magnetic field *************
249   SaveLocalConvConst();
250   // transform back error matrix from curvature to the 1/pt
251   fCovar[10]*=lcc;
252   fCovar[11]*=lcc;
253   fCovar[12]*=lcc;
254   fCovar[13]*=lcc;
255   fCovar[14]*=lcc*lcc;
256
257   Double_t dist = TMath::Sqrt((fX-oldX)*(fX-oldX)+(fParam[0]-oldY)*(fParam[0]-oldY)+
258                               (fParam[1]-oldZ)*(fParam[1]-oldZ));
259   if (!CorrectForMaterial(dist,x0,rho)) return 0;
260
261   // Integrated Time [SR, GSI, 17.02.2003]
262  //  if (IsStartedTimeIntegral() && fX>oldX) {
263 //     Double_t l2 = (fX-oldX)*(fX-oldX)+(fParam[0]-oldY)*(fParam[0]-oldY)+
264 //                   (fParam[1]-oldZ)*(fParam[1]-oldZ);
265 //     AddTimeStep(TMath::Sqrt(l2));
266 //   }
267   //
268
269   return kTRUE;
270 }
271
272 Bool_t     AliExternalTrackParam::PropagateToDCA(Double_t xd, Double_t yd,  Double_t x0, Double_t rho){
273   //
274   // Propagate the track parameters to the nearest point of given xv yv coordinate
275   //
276   Double_t a=fAlpha;
277   Double_t cs=TMath::Cos(a),sn=TMath::Sin(a);
278
279   Double_t xv= xd*cs + yd*sn;
280   Double_t yv=-xd*sn + yd*cs;   // vertex position in local frame
281   //  
282   Double_t c=fParam[4]/GetLocalConvConst(), snp=fParam[2];
283   //
284   Double_t x=fX, y=fParam[1];
285   Double_t tgfv=-(c*(x-xv)-snp)/(c*(y-yv) + TMath::Sqrt(1.-snp*snp));
286   Double_t fv=TMath::ATan(tgfv);
287   cs=TMath::Cos(fv); sn=TMath::Sin(fv);
288   x = xv*cs + yv*sn;
289   yv=-xv*sn + yv*cs; xv=x;
290   RotateTo(fv+a);
291   return PropagateTo(xv,x0,rho);  
292 }
293
294
295 //_____________________________________________________________________________
296 Bool_t AliExternalTrackParam::RotateTo(Double_t alp)
297 {
298   // Rotate the reference axis for the parametrisation to the given angle.
299   //
300   Double_t  x=fX;
301   Double_t p0=fParam[0];
302   //
303   if      (alp < -TMath::Pi()) alp += 2*TMath::Pi();
304   else if (alp >= TMath::Pi()) alp -= 2*TMath::Pi();
305   Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
306   Double_t sf=fParam[2], cf=TMath::Sqrt(1.- fParam[2]*fParam[2]);
307   // **** rotation **********************
308   
309   fAlpha = alp;
310   fX =  x*ca + p0*sa;
311   fParam[0]= -x*sa + p0*ca;
312   fParam[2]=  sf*ca - cf*sa;
313   Double_t rr=(ca+sf/cf*sa);  
314   //
315   fCovar[0] *= (ca*ca);
316   fCovar[1] *= ca; 
317   fCovar[3] *= ca*rr;
318   fCovar[6] *= ca;
319   fCovar[10] *= ca;
320   fCovar[4] *= rr;
321   fCovar[5] *= rr*rr;
322   fCovar[7] *= rr;
323   fCovar[11] *= rr;
324   return kTRUE;
325 }
326
327 //_____________________________________________________________________________
328 Bool_t AliExternalTrackParam::CorrectForMaterial(Double_t d, Double_t x0, Double_t rho)
329 {
330   //
331   // Take into account material effects assuming:
332   // x0  - mean rad length
333   // rho - mean density
334
335   //
336   // multiple scattering
337   //
338   if (fMass<=0) {
339     AliError("Non-positive mass");
340     return kFALSE;
341   }
342   Double_t p2=(1.+ fParam[3]*fParam[3])/(fParam[4]*fParam[4]);
343   Double_t beta2=p2/(p2 + fMass*fMass);
344   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
345   //
346   fCovar[5] += theta2*(1.- fParam[2]*fParam[2])*(1. + fParam[3]*fParam[3]);
347   fCovar[9] += theta2*(1. + fParam[3]*fParam[3])*(1. + fParam[3]*fParam[3]);
348   fCovar[13] += theta2*fParam[3]*fParam[4]*(1. + fParam[3]*fParam[3]);
349   fCovar[14] += theta2*fParam[3]*fParam[4]*fParam[3]*fParam[4];
350   
351   Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2+1e-10)) - beta2)*d*rho;  
352   fParam[4] *=(1.- TMath::Sqrt(p2+fMass*fMass)/p2*dE);
353   //
354   Double_t sigmade = 0.02*TMath::Sqrt(TMath::Abs(dE));   // energy loss fluctuation 
355   Double_t sigmac2 = sigmade*sigmade*fParam[4]*fParam[4]*(p2+fMass*fMass)/(p2*p2);
356   fCovar[14] += sigmac2;
357   //
358   //
359   
360
361
362   return kTRUE;
363 }
364
365 //_____________________________________________________________________________
366 Bool_t AliExternalTrackParam::GetProlongationAt(Double_t xk, 
367                                                 Double_t& y, 
368                                                 Double_t& z) const
369 {
370   //
371   // Get the local y and z coordinates at the given x value
372   //
373   Double_t lcc=GetLocalConvConst();  
374   Double_t cur = fParam[4]/lcc;
375   Double_t x1=fX, x2=xk, dx=x2-x1; 
376   Double_t f1=fParam[2], f2=f1 + cur*dx;
377   //  
378   if (TMath::Abs(f2) >= 0.98) {
379     // MI change  - don't propagate highly inclined tracks
380     //              covariance matrix distorted
381     return kFALSE;
382   }
383   Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
384   y = fParam[0] + dx*(f1+f2)/(r1+r2);
385   z = fParam[1] + dx*(f1+f2)/(f1*r2 + f2*r1)*fParam[3];
386   return kTRUE;
387 }
388
389 //_____________________________________________________________________________
390 Double_t AliExternalTrackParam::GetXAtVertex(Double_t /*x*/, 
391                                              Double_t /*y*/) const
392 {
393 // Get the x coordinate at the given vertex (x,y)
394 //
395 // NOT IMPLEMENTED for this class
396
397   return 0;
398 }
399
400
401 // //_____________________________________________________________________________
402 // Double_t AliExternalTrackParam::GetPredictedChi2(const AliCluster* /*cluster*/)
403 // {
404 // // calculate the chi2 contribution of the given cluster
405 // //
406 // // NOT IMPLEMENTED for this class
407
408 //   return -1;
409 // }
410
411 // //_____________________________________________________________________________
412 // Bool_t AliExternalTrackParam::Update(const AliCluster* /*cluster*/)
413 // {
414 // // update the track parameters using the position and error 
415 // // of the given cluster
416 // //
417 // // NOT IMPLEMENTED for this class
418
419 //   return kFALSE;
420 // }
421
422
423 //_____________________________________________________________________________
424 Double_t AliExternalTrackParam::SigmaPhi() const
425 {
426 // get the error of the azimuthal angle
427
428   return TMath::Sqrt(TMath::Abs(fCovar[5] / (1. - fParam[2]*fParam[2])));
429 }
430
431 //_____________________________________________________________________________
432 Double_t AliExternalTrackParam::SigmaTheta() const
433 {
434 // get the error of the polar angle
435
436   return TMath::Sqrt(TMath::Abs(fCovar[9])) / (1. + fParam[3]*fParam[3]);
437 }
438
439 //_____________________________________________________________________________
440 Double_t AliExternalTrackParam::SigmaPt() const
441 {
442 // get the error of the transversal component of the momentum
443
444   return TMath::Sqrt(fCovar[14]) / TMath::Abs(fParam[4]);
445 }
446
447 //_____________________________________________________________________________
448 TVector3 AliExternalTrackParam::Momentum() const
449 {
450 // get the momentum vector
451
452   Double_t phi = TMath::ASin(fParam[2]) + fAlpha;
453   Double_t pt = 1. / TMath::Abs(fParam[4]);
454   return TVector3(pt * TMath::Cos(phi), 
455                   pt * TMath::Sin(phi), 
456                   pt * fParam[3]);
457 }
458
459 //_____________________________________________________________________________
460 TVector3 AliExternalTrackParam::Position() const
461 {
462 // get the current spatial position in global coordinates
463
464   return TVector3(fX * TMath::Cos(fAlpha) - fParam[0] * TMath::Sin(fAlpha),
465                   fX * TMath::Sin(fAlpha) + fParam[0] * TMath::Cos(fAlpha),
466                   fParam[1]);
467 }
468
469
470 //_____________________________________________________________________________
471 void AliExternalTrackParam::Print(Option_t* /*option*/) const
472 {
473 // print the parameters and the covariance matrix
474
475   printf("AliExternalTrackParam: x = %-12g  alpha = %-12g\n", fX, fAlpha);
476   printf("  parameters: %12g %12g %12g %12g %12g\n",
477          fParam[0], fParam[1], fParam[2], fParam[3], fParam[4]);
478   printf("  covariance: %12g\n", fCovar[0]);
479   printf("              %12g %12g\n", fCovar[1], fCovar[2]);
480   printf("              %12g %12g %12g\n", fCovar[3], fCovar[4], fCovar[5]);
481   printf("              %12g %12g %12g %12g\n", 
482          fCovar[6], fCovar[7], fCovar[8], fCovar[9]);
483   printf("              %12g %12g %12g %12g %12g\n", 
484          fCovar[10], fCovar[11], fCovar[12], fCovar[13], fCovar[14]);
485 }