Bug corrected.
[u/mrichter/AliRoot.git] / STEER / AliTrackerBase.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: AliTrackerBase.cxx 38069 2009-12-24 16:56:18Z belikov $ */
17
18 //-------------------------------------------------------------------------
19 //               Implementation of the AliTrackerBase class
20 //                that is the base for the AliTracker class    
21 //                     Origin: Marian.Ivanov@cern.ch
22 //-------------------------------------------------------------------------
23 #include <TClass.h>
24 #include <TMath.h>
25 #include <TGeoManager.h>
26
27 #include "AliLog.h"
28 #include "AliTrackerBase.h"
29 #include "AliExternalTrackParam.h"
30 #include "AliTrackPointArray.h"
31 #include "TVectorD.h"
32
33 extern TGeoManager *gGeoManager;
34
35 ClassImp(AliTrackerBase)
36
37 AliTrackerBase::AliTrackerBase():
38   TObject(),
39   fX(0),
40   fY(0),
41   fZ(0),
42   fSigmaX(0.005),
43   fSigmaY(0.005),
44   fSigmaZ(0.010)
45 {
46   //--------------------------------------------------------------------
47   // The default constructor.
48   //--------------------------------------------------------------------
49   if (!TGeoGlobalMagField::Instance()->GetField())
50     AliWarning("Field map is not set.");
51 }
52
53 //__________________________________________________________________________
54 AliTrackerBase::AliTrackerBase(const AliTrackerBase &atr):
55   TObject(atr),
56   fX(atr.fX),
57   fY(atr.fY),
58   fZ(atr.fZ),
59   fSigmaX(atr.fSigmaX),
60   fSigmaY(atr.fSigmaY),
61   fSigmaZ(atr.fSigmaZ)
62 {
63   //--------------------------------------------------------------------
64   // The default constructor.
65   //--------------------------------------------------------------------
66   if (!TGeoGlobalMagField::Instance()->GetField())
67     AliWarning("Field map is not set.");
68 }
69
70 //__________________________________________________________________________
71 Double_t AliTrackerBase::GetBz()
72 {
73   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
74   if (!fld) return 0.5*kAlmost0Field;
75   Double_t bz = fld->SolenoidField();
76   return TMath::Sign(0.5*kAlmost0Field,bz) + bz;
77 }
78
79 //__________________________________________________________________________
80 Double_t AliTrackerBase::GetBz(const Double_t *r) {
81   //------------------------------------------------------------------
82   // Returns Bz (kG) at the point "r" .
83   //------------------------------------------------------------------
84   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
85   if (!fld) return  0.5*kAlmost0Field;
86   Double_t bz = fld->GetBz(r);
87   return  TMath::Sign(0.5*kAlmost0Field,bz) + bz;
88 }
89
90 //__________________________________________________________________________
91 void AliTrackerBase::GetBxByBz(const Double_t r[3], Double_t b[3]) {
92   //------------------------------------------------------------------
93   // Returns Bx, By and Bz (kG) at the point "r" .
94   //------------------------------------------------------------------
95   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
96   if (!fld) {
97      b[0] = b[1] = 0.;
98      b[2] = 0.5*kAlmost0Field;
99      return;
100   }
101
102   if (fld->IsUniform()) {
103      b[0] = b[1] = 0.;
104      b[2] = fld->SolenoidField();
105   }  else {
106      fld->Field(r,b);
107   }
108   b[2] = (TMath::Sign(0.5*kAlmost0Field,b[2]) + b[2]);
109   return;
110 }
111
112 Double_t AliTrackerBase::MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam)
113 {
114   // 
115   // Calculate mean material budget and material properties between 
116   //    the points "start" and "end".
117   //
118   // "mparam" - parameters used for the energy and multiple scattering
119   //  corrections: 
120   //
121   // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3]
122   // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional]
123   // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional]
124   // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional]
125   // mparam[4] - length: sum(x_i) [cm]
126   // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional]
127   // mparam[6] - number of boundary crosses
128   //
129   //  Origin:  Marian Ivanov, Marian.Ivanov@cern.ch
130   //
131   //  Corrections and improvements by
132   //        Andrea Dainese, Andrea.Dainese@lnl.infn.it,
133   //        Andrei Gheata,  Andrei.Gheata@cern.ch
134   //
135
136   mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0;
137   mparam[4]=0; mparam[5]=0; mparam[6]=0;
138   //
139   Double_t bparam[6]; // total parameters
140   Double_t lparam[6]; // local parameters
141
142   for (Int_t i=0;i<6;i++) bparam[i]=0;
143
144   if (!gGeoManager) {
145     AliErrorClass("No TGeo\n");
146     return 0.;
147   }
148   //
149   Double_t length;
150   Double_t dir[3];
151   length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+
152                        (end[1]-start[1])*(end[1]-start[1])+
153                        (end[2]-start[2])*(end[2]-start[2]));
154   mparam[4]=length;
155   if (length<TGeoShape::Tolerance()) return 0.0;
156   Double_t invlen = 1./length;
157   dir[0] = (end[0]-start[0])*invlen;
158   dir[1] = (end[1]-start[1])*invlen;
159   dir[2] = (end[2]-start[2])*invlen;
160
161   // Initialize start point and direction
162   TGeoNode *currentnode = 0;
163   TGeoNode *startnode = gGeoManager->InitTrack(start, dir);
164   if (!startnode) {
165     AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f",
166                   start[0],start[1],start[2]));
167     return 0.0;
168   }
169   TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial();
170   lparam[0]   = material->GetDensity();
171   lparam[1]   = material->GetRadLen();
172   lparam[2]   = material->GetA();
173   lparam[3]   = material->GetZ();
174   lparam[4]   = length;
175   lparam[5]   = lparam[3]/lparam[2];
176   if (material->IsMixture()) {
177     TGeoMixture * mixture = (TGeoMixture*)material;
178     lparam[5] =0;
179     Double_t sum =0;
180     for (Int_t iel=0;iel<mixture->GetNelements();iel++){
181       sum  += mixture->GetWmixt()[iel];
182       lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
183     }
184     lparam[5]/=sum;
185   }
186
187   // Locate next boundary within length without computing safety.
188   // Propagate either with length (if no boundary found) or just cross boundary
189   gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
190   Double_t step = 0.0; // Step made
191   Double_t snext = gGeoManager->GetStep();
192   // If no boundary within proposed length, return current density
193   if (!gGeoManager->IsOnBoundary()) {
194     mparam[0] = lparam[0];
195     mparam[1] = lparam[4]/lparam[1];
196     mparam[2] = lparam[2];
197     mparam[3] = lparam[3];
198     mparam[4] = lparam[4];
199     return lparam[0];
200   }
201   // Try to cross the boundary and see what is next
202   Int_t nzero = 0;
203   while (length>TGeoShape::Tolerance()) {
204     currentnode = gGeoManager->GetCurrentNode();
205     if (snext<2.*TGeoShape::Tolerance()) nzero++;
206     else nzero = 0;
207     if (nzero>3) {
208       // This means navigation has problems on one boundary
209       // Try to cross by making a small step
210       AliErrorClass("Cannot cross boundary\n");
211       mparam[0] = bparam[0]/step;
212       mparam[1] = bparam[1];
213       mparam[2] = bparam[2]/step;
214       mparam[3] = bparam[3]/step;
215       mparam[5] = bparam[5]/step;
216       mparam[4] = step;
217       mparam[0] = 0.;             // if crash of navigation take mean density 0
218       mparam[1] = 1000000;        // and infinite rad length
219       return bparam[0]/step;
220     }
221     mparam[6]+=1.;
222     step += snext;
223     bparam[1]    += snext/lparam[1];
224     bparam[2]    += snext*lparam[2];
225     bparam[3]    += snext*lparam[3];
226     bparam[5]    += snext*lparam[5];
227     bparam[0]    += snext*lparam[0];
228
229     if (snext>=length) break;
230     if (!currentnode) break;
231     length -= snext;
232     material = currentnode->GetVolume()->GetMedium()->GetMaterial();
233     lparam[0] = material->GetDensity();
234     lparam[1]  = material->GetRadLen();
235     lparam[2]  = material->GetA();
236     lparam[3]  = material->GetZ();
237     lparam[5]   = lparam[3]/lparam[2];
238     if (material->IsMixture()) {
239       TGeoMixture * mixture = (TGeoMixture*)material;
240       lparam[5]=0;
241       Double_t sum =0;
242       for (Int_t iel=0;iel<mixture->GetNelements();iel++){
243         sum+= mixture->GetWmixt()[iel];
244         lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
245       }
246       lparam[5]/=sum;
247     }
248     gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
249     snext = gGeoManager->GetStep();
250   }
251   mparam[0] = bparam[0]/step;
252   mparam[1] = bparam[1];
253   mparam[2] = bparam[2]/step;
254   mparam[3] = bparam[3]/step;
255   mparam[5] = bparam[5]/step;
256   return bparam[0]/step;
257 }
258
259
260 Bool_t 
261 AliTrackerBase::PropagateTrackTo(AliExternalTrackParam *track, Double_t xToGo, 
262                              Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Double_t sign){
263   //----------------------------------------------------------------
264   //
265   // Propagates the track to the plane X=xk (cm) using the magnetic field map 
266   // and correcting for the crossed material.
267   //
268   // mass     - mass used in propagation - used for energy loss correction
269   // maxStep  - maximal step for propagation
270   //
271   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
272   //
273   //----------------------------------------------------------------
274   const Double_t kEpsilon = 0.00001;
275   Double_t xpos     = track->GetX();
276   Double_t dir      = (xpos<xToGo) ? 1.:-1.;
277   //
278   while ( (xToGo-xpos)*dir > kEpsilon){
279     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
280     Double_t x    = xpos+step;
281     Double_t xyz0[3],xyz1[3],param[7];
282     track->GetXYZ(xyz0);   //starting global position
283
284     Double_t bz=GetBz(xyz0); // getting the local Bz
285
286     if (!track->GetXYZAt(x,bz,xyz1)) return kFALSE;   // no prolongation
287     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
288
289     if (TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return kFALSE;
290     if (!track->PropagateTo(x,bz))  return kFALSE;
291
292     MeanMaterialBudget(xyz0,xyz1,param);        
293     Double_t xrho=param[0]*param[4]*sign, xx0=param[1];
294
295     if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
296     if (rotateTo){
297       if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
298       track->GetXYZ(xyz0);   // global position
299       Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]); 
300       //
301       Double_t ca=TMath::Cos(alphan-track->GetAlpha()), 
302                sa=TMath::Sin(alphan-track->GetAlpha());
303       Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
304       Double_t sinNew =  sf*ca - cf*sa;
305       if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
306       if (!track->Rotate(alphan)) return kFALSE;
307     }
308     xpos = track->GetX();
309   }
310   return kTRUE;
311 }
312
313 Bool_t 
314 AliTrackerBase::PropagateTrackToBxByBz(AliExternalTrackParam *track,
315 Double_t xToGo, 
316                                    Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp,Double_t sign){
317   //----------------------------------------------------------------
318   //
319   // Propagates the track to the plane X=xk (cm)
320   // taking into account all the three components of the magnetic field 
321   // and correcting for the crossed material.
322   //
323   // mass     - mass used in propagation - used for energy loss correction
324   // maxStep  - maximal step for propagation
325   //
326   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
327   //
328   //----------------------------------------------------------------
329   const Double_t kEpsilon = 0.00001;
330   Double_t xpos     = track->GetX();
331   Double_t dir      = (xpos<xToGo) ? 1.:-1.;
332   //
333   while ( (xToGo-xpos)*dir > kEpsilon){
334     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
335     Double_t x    = xpos+step;
336     Double_t xyz0[3],xyz1[3],param[7];
337     track->GetXYZ(xyz0);   //starting global position
338
339     Double_t b[3]; GetBxByBz(xyz0,b); // getting the local Bx, By and Bz
340
341     if (!track->GetXYZAt(x,b[2],xyz1)) return kFALSE;   // no prolongation
342     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
343
344     if (TMath::Abs(track->GetSnpAt(x,b[2])) >= maxSnp) return kFALSE;
345     if (!track->PropagateToBxByBz(x,b))  return kFALSE;
346
347     MeanMaterialBudget(xyz0,xyz1,param);        
348     Double_t xrho=param[0]*param[4]*sign, xx0=param[1];
349
350     if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
351     if (rotateTo){
352       if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
353       track->GetXYZ(xyz0);   // global position
354       Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]); 
355       //
356       Double_t ca=TMath::Cos(alphan-track->GetAlpha()), 
357                sa=TMath::Sin(alphan-track->GetAlpha());
358       Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
359       Double_t sinNew =  sf*ca - cf*sa;
360       if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
361       if (!track->Rotate(alphan)) return kFALSE;
362     }
363     xpos = track->GetX();
364   }
365   return kTRUE;
366 }
367
368 Double_t AliTrackerBase::GetTrackPredictedChi2(AliExternalTrackParam *track,
369                                            Double_t mass, Double_t step,
370                                      const AliExternalTrackParam *backup) {
371   //
372   // This function brings the "track" with particle "mass" [GeV] 
373   // to the same local coord. system and the same reference plane as 
374   // of the "backup", doing it in "steps" [cm].
375   // Then, it calculates the 5D predicted Chi2 for these two tracks
376   //
377   Double_t chi2=kVeryBig;
378   Double_t alpha=backup->GetAlpha();
379   if (!track->Rotate(alpha)) return chi2;
380
381   Double_t xb=backup->GetX();
382   Double_t sign=(xb < track->GetX()) ? 1. : -1.;
383   if (!PropagateTrackTo(track,xb,mass,step,kFALSE,kAlmost1,sign)) return chi2;
384
385   chi2=track->GetPredictedChi2(backup);
386
387   return chi2;
388 }
389
390
391
392
393 Double_t AliTrackerBase::MakeC(Double_t x1,Double_t y1,
394                    Double_t x2,Double_t y2,
395                    Double_t x3,Double_t y3)
396 {
397   //-----------------------------------------------------------------
398   // Initial approximation of the track curvature
399   //-----------------------------------------------------------------
400   x3 -=x1;
401   x2 -=x1;
402   y3 -=y1;
403   y2 -=y1;
404   //  
405   Double_t det = x3*y2-x2*y3;
406   if (TMath::Abs(det)<1e-10) {
407     return 0;
408   }
409   //
410   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
411   Double_t x0 = x3*0.5-y3*u;
412   Double_t y0 = y3*0.5+x3*u;
413   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
414   if (det>0) c2*=-1;
415   return c2;
416 }
417
418 Double_t AliTrackerBase::MakeSnp(Double_t x1,Double_t y1,
419                    Double_t x2,Double_t y2,
420                    Double_t x3,Double_t y3)
421 {
422   //-----------------------------------------------------------------
423   // Initial approximation of the track snp
424   //-----------------------------------------------------------------
425   x3 -=x1;
426   x2 -=x1;
427   y3 -=y1;
428   y2 -=y1;
429   //  
430   Double_t det = x3*y2-x2*y3;
431   if (TMath::Abs(det)<1e-10) {
432     return 0;
433   }
434   //
435   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
436   Double_t x0 = x3*0.5-y3*u; 
437   Double_t y0 = y3*0.5+x3*u;
438   Double_t c2 = 1./TMath::Sqrt(x0*x0+y0*y0);
439   x0*=c2;
440   x0=TMath::Abs(x0);
441   if (y2*x2<0.) x0*=-1;
442   return x0;
443 }
444
445 Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
446                    Double_t x2,Double_t y2,
447                    Double_t z1,Double_t z2, Double_t c)
448 {
449   //-----------------------------------------------------------------
450   // Initial approximation of the tangent of the track dip angle
451   //-----------------------------------------------------------------
452   //
453   x2-=x1;
454   y2-=y1;
455   z2-=z1;
456   Double_t d  =  TMath::Sqrt(x2*x2+y2*y2);  // distance  straight line
457   if (TMath::Abs(d*c*0.5)>1) return 0;
458   Double_t   angle2    = TMath::ASin(d*c*0.5); 
459   angle2  = z2*TMath::Abs(c/(angle2*2.));
460   return angle2;
461 }
462
463
464 Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1, 
465                    Double_t x2,Double_t y2,
466                    Double_t z1,Double_t z2) 
467 {
468   //-----------------------------------------------------------------
469   // Initial approximation of the tangent of the track dip angle
470   //-----------------------------------------------------------------
471   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
472 }
473
474
475 AliExternalTrackParam * AliTrackerBase::MakeSeed( AliTrackPoint &point0, AliTrackPoint &point1, AliTrackPoint &point2){
476   //
477   // Make Seed  - AliExternalTrackParam from input 3 points   
478   // returning seed in local frame of point0
479   //
480   Double_t xyz0[3]={0,0,0};
481   Double_t xyz1[3]={0,0,0};
482   Double_t xyz2[3]={0,0,0};
483   Double_t alpha=point0.GetAngle();
484   Double_t xyz[3]={point0.GetX(),point0.GetY(),point0.GetZ()};
485   Double_t bxyz[3]; GetBxByBz(xyz,bxyz); 
486   Double_t bz = bxyz[2];
487   //
488   // get points in frame of point 0
489   //
490   AliTrackPoint p0r = point0.Rotate(alpha);
491   AliTrackPoint p1r = point1.Rotate(alpha);
492   AliTrackPoint p2r = point2.Rotate(alpha);
493   xyz0[0]=p0r.GetX();
494   xyz0[1]=p0r.GetY();
495   xyz0[2]=p0r.GetZ();
496   xyz1[0]=p1r.GetX();
497   xyz1[1]=p1r.GetY();
498   xyz1[2]=p1r.GetZ();
499   xyz2[0]=p2r.GetX();
500   xyz2[1]=p2r.GetY();
501   xyz2[2]=p2r.GetZ();
502   //
503   // make covariance estimate
504   //  
505   Double_t covar[15];
506   Double_t param[5]={0,0,0,0,0};
507   for (Int_t m=0; m<15; m++) covar[m]=0;
508   //
509   // calculate intitial param
510   param[0]=xyz0[1];              
511   param[1]=xyz0[2];
512   param[2]=MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
513   param[4]=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
514   param[3]=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2],param[4]);
515
516   //covariance matrix - only diagonal elements
517   //Double_t dist=p0r.GetX()-p2r.GetX();
518   Double_t deltaP=0;
519   covar[0]= p0r.GetCov()[3];
520   covar[2]= p0r.GetCov()[5];
521   //sigma snp
522   deltaP= (MakeSnp(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[2]);
523   covar[5]+= deltaP*deltaP;
524   deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[2]);
525   covar[5]+= deltaP*deltaP;
526   deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p1r.GetCov()[3]))-param[2]);
527   covar[5]+= deltaP*deltaP;
528   //sigma tgl
529   //
530   deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2]+TMath::Sqrt(p1r.GetCov()[5]),xyz1[2],param[4])-param[3];
531   covar[9]+= deltaP*deltaP;
532   deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2]+TMath::Sqrt(p1r.GetCov()[5]),param[4])-param[3];
533   covar[9]+= deltaP*deltaP;
534   //
535   
536   deltaP=MakeC(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[4];
537   covar[14]+= deltaP*deltaP;
538   deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[4];
539   covar[14]+= deltaP*deltaP;
540   deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p2r.GetCov()[3]))-param[4];
541   covar[14]+= deltaP*deltaP;
542   
543   covar[14]/=(bz*kB2C)*(bz*kB2C);
544   param[4]/=(bz*kB2C); // transform to 1/pt
545   AliExternalTrackParam * trackParam = new AliExternalTrackParam(xyz0[0],alpha,param, covar);
546   if (0) {
547     // consistency check  -to put warnings here 
548     // small disagrement once Track extrapolation used 
549     // nice agreement in seeds with MC track parameters - problem in extrapoloation - to be fixed
550     // to check later
551     Double_t y1,y2,z1,z2;
552     trackParam->GetYAt(xyz1[0],bz,y1);
553     trackParam->GetZAt(xyz1[0],bz,z1);
554     trackParam->GetYAt(xyz2[0],bz,y2);
555     trackParam->GetZAt(xyz2[0],bz,z2);
556     if (TMath::Abs(y1-xyz1[1])> TMath::Sqrt(p1r.GetCov()[3]*5)){
557       AliWarningClass("Seeding problem y1\n");
558     }
559     if (TMath::Abs(y2-xyz2[1])> TMath::Sqrt(p2r.GetCov()[3]*5)){
560       AliWarningClass("Seeding problem y2\n");
561     }
562     if (TMath::Abs(z1-xyz1[2])> TMath::Sqrt(p1r.GetCov()[5]*5)){
563       AliWarningClass("Seeding problem z1\n");
564     }
565   }
566   return trackParam;  
567
568
569 Double_t  AliTrackerBase::FitTrack(AliExternalTrackParam * trackParam, AliTrackPointArray *pointArray, Double_t mass, Double_t maxStep){
570   //
571   // refit the track  - trackParam using the points in point array  
572   //
573   const Double_t kMaxSnp=0.99;
574   if (!trackParam) return 0;
575   Int_t  npoints=pointArray->GetNPoints();
576   AliTrackPoint point,point2;
577   Double_t pointPos[2]={0,0};
578   Double_t pointCov[3]={0,0,0};
579   // choose coordinate frame
580   // in standard way the coordinate frame should be changed point by point
581   // Some problems with rotation observed
582   // rotate method of AliExternalTrackParam should be revisited
583   pointArray->GetPoint(point,0);
584   pointArray->GetPoint(point2,npoints-1);
585   Double_t alpha=TMath::ATan2(point.GetY()-point2.GetY(), point.GetX()-point2.GetX());
586   
587   for (Int_t ipoint=npoints-1; ipoint>0; ipoint-=1){
588     pointArray->GetPoint(point,ipoint);
589     AliTrackPoint pr = point.Rotate(alpha);
590     trackParam->Rotate(alpha);
591     Bool_t status = PropagateTrackTo(trackParam,pr.GetX(),mass,maxStep,kFALSE,kMaxSnp);
592     if(!status){
593       AliWarningClass("Problem to propagate\n");    
594       break;
595     }
596     if (TMath::Abs(trackParam->GetSnp())>kMaxSnp){ 
597       AliWarningClass("sin(phi) > kMaxSnp \n");
598       break;
599     }
600     pointPos[0]=pr.GetY();//local y
601     pointPos[1]=pr.GetZ();//local z
602     pointCov[0]=pr.GetCov()[3];//simay^2
603     pointCov[1]=pr.GetCov()[4];//sigmayz
604     pointCov[2]=pr.GetCov()[5];//sigmaz^2
605     trackParam->Update(pointPos,pointCov); 
606   }
607   return 0;
608 }
609
610
611
612 void AliTrackerBase::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
613   //
614   // Update track 1 with track 2
615   //
616   //
617   //
618   TMatrixD vecXk(5,1);    // X vector
619   TMatrixD covXk(5,5);    // X covariance 
620   TMatrixD matHk(5,5);    // vector to mesurement
621   TMatrixD measR(5,5);    // measurement error 
622   TMatrixD vecZk(5,1);    // measurement
623   //
624   TMatrixD vecYk(5,1);    // Innovation or measurement residual
625   TMatrixD matHkT(5,5);
626   TMatrixD matSk(5,5);    // Innovation (or residual) covariance
627   TMatrixD matKk(5,5);    // Optimal Kalman gain
628   TMatrixD mat1(5,5);     // update covariance matrix
629   TMatrixD covXk2(5,5);   // 
630   TMatrixD covOut(5,5);
631   //
632   Double_t *param1=(Double_t*) track1.GetParameter();
633   Double_t *covar1=(Double_t*) track1.GetCovariance();
634   Double_t *param2=(Double_t*) track2.GetParameter();
635   Double_t *covar2=(Double_t*) track2.GetCovariance();
636   //
637   // copy data to the matrix
638   for (Int_t ipar=0; ipar<5; ipar++){
639     for (Int_t jpar=0; jpar<5; jpar++){
640       covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
641       measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
642       matHk(ipar,jpar)=0;
643       mat1(ipar,jpar)=0;
644     }
645     vecXk(ipar,0) = param1[ipar];
646     vecZk(ipar,0) = param2[ipar];
647     matHk(ipar,ipar)=1;
648     mat1(ipar,ipar)=1;
649   }
650   //
651   //
652   //
653   //
654   //
655   vecYk = vecZk-matHk*vecXk;                 // Innovation or measurement residual
656   matHkT=matHk.T(); matHk.T();
657   matSk = (matHk*(covXk*matHkT))+measR;      // Innovation (or residual) covariance
658   matSk.Invert();
659   matKk = (covXk*matHkT)*matSk;              //  Optimal Kalman gain
660   vecXk += matKk*vecYk;                      //  updated vector 
661   covXk2 = (mat1-(matKk*matHk));
662   covOut =  covXk2*covXk; 
663   //
664   //
665   //
666   // copy from matrix to parameters
667   if (0) {
668     vecXk.Print();
669     vecZk.Print();
670     //
671     measR.Print();
672     covXk.Print();
673     covOut.Print();
674     //
675     track1.Print();
676     track2.Print();
677   }
678
679   for (Int_t ipar=0; ipar<5; ipar++){
680     param1[ipar]= vecXk(ipar,0) ;
681     for (Int_t jpar=0; jpar<5; jpar++){
682       covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
683     }
684   }
685 }
686
687