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