Updated macro for ACORDE raw-data visualisation. We use now the ACORDE sym.names...
[u/mrichter/AliRoot.git] / TPC / AliTPCTracklet.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 ////
17 // This class stores a tracklet (a track that lives only in a single TPC
18 // sector). Its objects can be constructed out of TPCseeds, that are
19 // holding the necessary cluster information.
20 ////
21 ////
22 //// 
23
24 #include "AliTPCTracklet.h"
25 #include "TObjArray.h"
26 #include "TLinearFitter.h"
27 #include "AliTPCseed.h"
28 #include "AliESDVertex.h"
29 #include "AliTracker.h"
30 #include "TTreeStream.h"
31 #include "TRandom3.h"
32 #include "TDecompChol.h"
33
34 #include <iostream>
35 using namespace std;
36
37 ClassImp(AliTPCTracklet)
38
39 const Double_t AliTPCTracklet::kB2C=0.299792458e-3;
40
41 AliTPCTracklet::AliTPCTracklet() 
42   : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(-1),fOuter(0),
43     fInner(0),fPrimary(0) {
44   ////
45   // The default constructor. It is intended to be used for I/O only.
46   ////
47 }
48
49 AliTPCTracklet::AliTPCTracklet(const AliTPCseed *track,Int_t sector,
50                                TrackType type,Bool_t storeClusters)
51   : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(sector),fOuter(0),
52     fInner(0),fPrimary(0) {
53   ////
54   // Contructor for a tracklet out of a track. Only clusters within a given 
55   // sector are used.
56   ///
57
58   //TODO: only kalman works
59   
60   for (Int_t i=0;i<160;++i) {
61     AliTPCclusterMI *c=track->GetClusterPointer(i);
62     if (c&&c->GetDetector()==sector)
63       ++fNClusters;
64   }
65
66   if (storeClusters) {
67     fClusters=new AliTPCclusterMI[fNClusters];
68     for (Int_t i=0;i<160;++i) {
69       AliTPCclusterMI *c=track->GetClusterPointer(i);
70       if (c&&c->GetDetector()==sector)
71         fClusters[fNStoredClusters]=c;
72       ++fNStoredClusters;
73     }
74   }
75
76   switch (type) {
77   case kKalman:
78     FitKalman(track,sector);
79     break;
80   case kLinear:
81   case kQuadratic:
82     FitLinear(track,sector,type);
83     break;
84   case kRiemann:
85     FitRiemann(track,sector);
86     break;
87   }
88
89 }
90
91 AliTPCTracklet::AliTPCTracklet(const TObjArray &/*clusters*/,Int_t sector,
92                                TrackType /*type*/,Bool_t /*storeClusters*/)
93   : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(sector),fOuter(0),
94     fInner(0),fPrimary(0) {
95   //TODO: write it!
96 }
97
98 AliTPCTracklet::AliTPCTracklet(const AliTPCTracklet &t)
99   : TObject(t),fNClusters(t.fNClusters),fNStoredClusters(t.fNStoredClusters),fClusters(0),
100     fSector(t.fSector),fOuter(0),fInner(0),
101     fPrimary(0) {
102   ////
103   // The copy constructor. You can copy tracklets! 
104   ////
105
106   if (t.fClusters) {
107     fClusters=new AliTPCclusterMI[t.fNStoredClusters];
108     for (int i=0;i<t.fNStoredClusters;++i)
109       fClusters[i]=t.fClusters[i];
110   }
111   if (t.fOuter)
112     fOuter=new AliExternalTrackParam(*t.fOuter);
113   if (t.fInner)
114     fInner=new AliExternalTrackParam(*t.fInner);
115   if (t.fPrimary)
116     fPrimary=new AliExternalTrackParam(*t.fPrimary);
117 }
118
119 AliTPCTracklet& AliTPCTracklet::operator=(const AliTPCTracklet &t) {
120   ////
121   // The assignment constructor. You can assign tracklets!
122   ////
123   if (this!=&t) {
124     fNClusters=t.fNClusters;
125     fNStoredClusters=fNStoredClusters;
126     delete fClusters;
127     if (t.fClusters) {
128       fClusters=new AliTPCclusterMI[t.fNStoredClusters];
129       for (int i=0;i<t.fNStoredClusters;++i)
130         fClusters[i]=t.fClusters[i];
131     }
132     else
133       fClusters=0;
134     fSector=t.fSector;
135     if (t.fOuter) {
136       if (fOuter)
137         *fOuter=*t.fOuter;
138       else
139         fOuter=new AliExternalTrackParam(*t.fOuter);
140     }
141     else {
142       delete fOuter;
143       fOuter=0;
144     }
145
146     if (t.fInner) {
147       if (fInner)
148         *fInner=*t.fInner;
149       else
150         fInner=new AliExternalTrackParam(*t.fInner);
151     }
152     else {
153       delete fInner;
154       fInner=0;
155     }
156
157     if (t.fPrimary) {
158       if (fPrimary)
159         *fPrimary=*t.fPrimary;
160       else
161         fPrimary=new AliExternalTrackParam(*t.fPrimary);
162     }
163     else {
164       delete fPrimary;
165       fPrimary=0;
166     }
167   }
168   return *this;
169 }
170
171 AliTPCTracklet::~AliTPCTracklet() {
172   //
173   // The destructor. Yes, you can even destruct tracklets.
174   //
175   delete fClusters;
176   delete fOuter;
177   delete fInner;
178   delete fPrimary;
179 }
180
181 void AliTPCTracklet::FitKalman(const AliTPCseed *track,Int_t sector) {
182   AliTPCseed *t=new AliTPCseed(*track);
183   if (!t->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-t->GetAlpha())) {
184     delete t;
185     return;
186   }
187   // fit from inner to outer row
188   AliTPCseed *outerSeed=new AliTPCseed(*t);
189   Int_t n=0;
190   for (Int_t i=0;i<160;++i) {
191     AliTPCclusterMI *c=t->GetClusterPointer(i);
192     if (c&&c->GetDetector()==sector) {
193       if (n==1) {
194         outerSeed->ResetCovariance(100.);
195       }
196       ++n;
197       Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
198       Double_t cov[3]={0.1,0.,0.1}; //TODO: correct error parametrisation
199       if (!outerSeed->PropagateTo(r[0]) ||
200           !static_cast<AliExternalTrackParam*>(outerSeed)->Update(&r[1],cov)) {
201         delete outerSeed;
202         outerSeed=0;
203         break;
204       }
205     }
206   }
207   if (outerSeed)
208     fOuter=new AliExternalTrackParam(*outerSeed);
209   delete outerSeed;
210   // fit from outer to inner rows
211   AliTPCseed *innerSeed=new AliTPCseed(*t);
212   n=0;
213   for (Int_t i=159;i>=0;--i) {
214     AliTPCclusterMI *c=t->GetClusterPointer(i);
215     if (c&&c->GetDetector()==sector) {
216       if (n==1) {
217         innerSeed->ResetCovariance(100.);
218       }
219       ++n;
220       Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
221       Double_t cov[3]={0.1,0.,0.1};
222       if (!innerSeed->PropagateTo(r[0]) ||
223           !static_cast<AliExternalTrackParam*>(innerSeed)->Update(&r[1],cov)) {
224         delete innerSeed;
225         innerSeed=0;
226         break;
227       }
228     }
229   }
230   if (innerSeed)
231     fInner=new AliExternalTrackParam(*innerSeed);
232   // propagate to the primary vertex
233   if (innerSeed) {
234     AliTPCseed *primarySeed=new AliTPCseed(*innerSeed);
235     Double_t pos[]={0.,0.,0.};
236     Double_t sigma[]={.1,.1,.1}; //TODO: is this correct?
237     AliESDVertex vertex(pos,sigma);
238     if (primarySeed->PropagateToVertex(&vertex))
239       fPrimary=new AliExternalTrackParam(*primarySeed);
240     delete primarySeed;
241     // for better comparison one does not want to have alpha changed...
242     if (!fPrimary->Rotate(fInner->GetAlpha())) {
243       delete fPrimary;
244       fPrimary=0;
245     }
246   }
247   delete innerSeed;
248
249   delete t;
250 }
251
252 void AliTPCTracklet::FitLinear(const AliTPCseed *track,Int_t sector,
253                                TrackType type) {
254   TLinearFitter fy(1);
255   TLinearFitter fz(1);
256   fy.StoreData(kFALSE);
257   fz.StoreData(kFALSE);
258   switch (type) {
259   case kLinear:
260     fy.SetFormula("1 ++ x");
261     fz.SetFormula("1 ++ x");
262     break;
263   case kQuadratic:
264     fy.SetFormula("1 ++ x ++ x*x");
265     fz.SetFormula("1 ++ x");
266     break;
267   case kKalman:
268   case kRiemann:
269     break;
270   }
271   Double_t xmax=-1.;
272   Double_t xmin=1000.;
273   for (Int_t i=0;i<160;++i) {
274     AliTPCclusterMI *c=track->GetClusterPointer(i);
275     if (c&&c->GetDetector()==sector) {
276       Double_t x=c->GetX();
277       fy.AddPoint(&x,c->GetY());
278       fz.AddPoint(&x,c->GetZ());
279       xmax=TMath::Max(xmax,x);
280       xmin=TMath::Min(xmin,x);
281     }
282   }
283   fy.Eval();
284   fz.Eval();
285   Double_t a[3]={fy.GetParameter(0),
286                  fy.GetParameter(1),
287                  type==kQuadratic?fy.GetParameter(2):0.};
288   Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
289                   fy.GetCovarianceMatrixElement(1,0),
290                   fy.GetCovarianceMatrixElement(1,1),
291                   type==kQuadratic?fy.GetCovarianceMatrixElement(2,0):0.,
292                   type==kQuadratic?fy.GetCovarianceMatrixElement(2,1):0.,
293                   type==kQuadratic?fy.GetCovarianceMatrixElement(2,2):0.};
294   for (int i=0;i<6;++i) ca[i]*=fy.GetChisquare()/fNClusters;
295   Double_t b[2]={fz.GetParameter(0),
296                  fz.GetParameter(1)};
297   Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
298                   fz.GetCovarianceMatrixElement(1,0),
299                   fz.GetCovarianceMatrixElement(1,1)};
300   for (int i=0;i<3;++i) cb[i]*=fz.GetChisquare()/fNClusters;
301   Double_t p[5];
302   Double_t c[15];
303   Double_t alpha=track->GetAlpha();
304   Quadratic2Helix(a,ca,b,cb,0.,p,c);
305   fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
306   Quadratic2Helix(a,ca,b,cb,xmin,p,c);
307   fInner=new AliExternalTrackParam(xmin,alpha,p,c);
308   Quadratic2Helix(a,ca,b,cb,xmax,p,c);
309   fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
310 }
311   
312 void AliTPCTracklet::Quadratic2Helix(Double_t *a,Double_t *ca,
313                                      Double_t *b,Double_t *cb,
314                                      Double_t x0,
315                                      Double_t *p,Double_t *c) {
316   // y(x)=a[0]+a[1]*x+a[2]*x^2
317   // z(x)=b[0]+b[1]*x
318   // parametrises the corosponding helix at x0
319
320   // get the polynoms at x0
321   Double_t a0=x0*x0*a[2] + x0*a[1] + a[0];
322   Double_t a1=2.*x0*a[2] +     a[1];
323   Double_t a2=      a[2];
324   Double_t ca00=ca[0]+x0*(2.*ca[1]+x0*(ca[2]+2.*ca[3]+x0*(2.*ca[4]+x0*ca[5])));
325   Double_t ca10=ca[1]+x0*(ca[2]+2.*ca[3]+x0*(3.*ca[4]+x0*2.*ca[5]));
326   Double_t ca11=ca[2]+x0*4.*(ca[4]+x0*ca[5]);
327   Double_t ca20=ca[3]+x0*(ca[4]+x0*ca[5]);
328   Double_t ca21=ca[3]+x0*2.*ca[5];
329   Double_t ca22=ca[5];
330
331   Double_t b0=x0*b[1] + b[0];
332   Double_t b1=   b[1];
333   Double_t cb00=cb[0]+x0*(2.*cb[1]+x0*cb[2]);
334   Double_t cb10=cb[1]+x0*cb[2];
335   Double_t cb11=cb[2];
336
337   // transform to helix parameters
338   Double_t f   =1.+a1*a1;
339   Double_t f2  =f*f;
340   Double_t fi  =1./f; 
341   Double_t fi12=TMath::Sqrt(fi);
342   Double_t fi32=fi*fi12;
343   Double_t fi2 =fi*fi;
344   Double_t fi52=fi2*fi12;
345   Double_t fi3 =fi2*fi;
346   Double_t fi5 =fi2*fi3;
347   
348   Double_t xyz[3]={0.}; // TODO...
349   Double_t fc=1./(GetBz(xyz)*kB2C);
350
351   p[0]=a0;            // y0
352   p[1]=b0;            // z0
353   p[2]=a1*fi12;       // snp
354   p[3]=b1;            // tgl
355   p[4]=2.*a2*fi32*fc; // 1/pt
356
357   c[0] =ca00;      //  y0-y0
358   c[1] =0.;        //  z0-y0
359   c[2] =cb00;      //  z0-z0
360   c[3] =ca10*fi32; // snp-y0
361   c[4] =0.;        // snp-z0
362   c[5] =ca11*fi3;  // snp-snp
363   c[6] =0.;        // tgl-y0
364   c[7] =cb10;      // tgl-z0
365   c[8] =0.;        // tgl-snp
366   c[9] =cb11;      // tgl-tgl
367   c[10]=2.*(-3.*a1*a2*ca10+f*ca20)*fi3*fc;  // 1/pt-y0
368   c[11]=0.;                                 // 1/pt-z0
369   c[12]=2.*(-3.*a1*a2*ca11+f*ca21)*fi52*fc; // 1/pt-snp
370   c[13]=0.;                                 // 1/pt-tgl
371   c[14]=(-12.*a1*a2*(-3.*a1*a2*ca11+2.*f*ca21)+4.*f2*ca22)*fi5
372     *fc*fc;        // 1/pt-1/pt
373 }
374
375
376 void AliTPCTracklet::FitRiemann(const AliTPCseed *track,Int_t sector) {
377   TLinearFitter fy(2);
378   fy.StoreData(kFALSE);
379   fy.SetFormula("hyp2");
380   Double_t xmax=-1.;
381   Double_t xmin=1000.;
382   for (Int_t i=0;i<160;++i) {
383     AliTPCclusterMI *c=track->GetClusterPointer(i);
384     if (c&&c->GetDetector()==sector) {
385       Double_t x=c->GetX();
386       Double_t y=c->GetY();
387       Double_t xy[2]={x,y};
388       Double_t r=x*x+y*y;
389       Double_t errx=1.,erry=1.;//TODO!
390       Double_t err=TMath::Sqrt(4.*x*x*errx+4.*y*y*erry);
391       err=1.;
392       fy.AddPoint(xy,r,err);
393       xmax=TMath::Max(xmax,x);
394       xmin=TMath::Min(xmin,x);
395     }
396   }
397   fy.Eval();
398   Double_t a[3]={fy.GetParameter(0),
399                  fy.GetParameter(1),
400                  fy.GetParameter(2)};
401   Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
402                   fy.GetCovarianceMatrixElement(1,0),
403                   fy.GetCovarianceMatrixElement(1,1),
404                   fy.GetCovarianceMatrixElement(2,0),
405                   fy.GetCovarianceMatrixElement(2,1),
406                   fy.GetCovarianceMatrixElement(2,2)};
407
408   TLinearFitter fz(1);
409   fz.StoreData(kFALSE);
410   fz.SetFormula("hyp1");
411   Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
412   Double_t oldx=0.;
413   Double_t oldy=R;
414   Double_t phi=0.;
415   for (Int_t i=0;i<160;++i) {
416     AliTPCclusterMI *c=track->GetClusterPointer(i);
417     if (c&&c->GetDetector()==sector) {
418       Double_t x=c->GetX();
419       Double_t y=c->GetY();
420       Double_t dx=x-oldx;
421       Double_t dy=y-oldy;
422       phi+=2.*TMath::Abs(TMath::ATan2(.5*TMath::Sqrt(dx*dx+dy*dy),R));
423       Double_t err=1.;
424       fz.AddPoint(&phi,c->GetZ(),err);
425       oldx=x;
426       oldy=y;
427     }
428   }
429   fz.Eval();
430   Double_t b[2]={fz.GetParameter(0),
431                  fz.GetParameter(1)};
432   Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
433                   fz.GetCovarianceMatrixElement(1,0),
434                   fz.GetCovarianceMatrixElement(1,1)};
435
436   Double_t p[5];
437   Double_t c[15];
438   Double_t alpha=track->GetAlpha();
439   if (Riemann2Helix(a,ca,b,cb,0.,p,c))
440     fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
441   if (Riemann2Helix(a,ca,b,cb,xmin,p,c))
442     fInner=new AliExternalTrackParam(xmin,alpha,p,c);
443   if (Riemann2Helix(a,ca,b,cb,xmax,p,c))
444     fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
445 }
446
447 Bool_t AliTPCTracklet::Riemann2Helix(Double_t *a,Double_t */*ca*/,
448                                      Double_t *b,Double_t */*cb*/,
449                                      Double_t x0,
450                                      Double_t *p,Double_t *c) {
451   //TODO: signs!
452
453   Double_t xr0=.5*a[1];
454   Double_t yr0=.5*a[2];
455   Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
456   Double_t dx=x0-xr0;
457   if (dx*dx>=R*R) return kFALSE;
458   Double_t dy=TMath::Sqrt(R*R-dx*dx); //sign!!
459   if (TMath::Abs(yr0+dy)>TMath::Abs(yr0-dy))
460     dy=-dy;
461   Double_t y0=yr0+dy; 
462   Double_t tgp=-dx/dy; //TODO: dy!=0
463   Double_t z0=b[0]+TMath::ATan(tgp)*b[1];
464   Double_t xyz[3]={x0,y0,z0};
465   Double_t fc=1./(GetBz(xyz)*kB2C);
466   fc=1;
467   p[0]=y0;  // y0
468   p[1]=z0; // z0
469   p[2]=tgp/TMath::Sqrt(1.+tgp*tgp); // snp
470   p[3]=b[1];       // tgl
471   p[4]=1./R*fc;    // 1/pt
472
473   c[0] =0.;      //  y0-y0
474   c[1] =0.;        //  z0-y0
475   c[2] =0.;      //  z0-z0
476   c[3] =0.; // snp-y0
477   c[4] =0.;        // snp-z0
478   c[5] =0.;  // snp-snp
479   c[6] =0.;        // tgl-y0
480   c[7] =0.;      // tgl-z0
481   c[8] =0.;        // tgl-snp
482   c[9] =0.;      // tgl-tgl
483   c[10]=0.;  // 1/pt-y0
484   c[11]=0.;                                 // 1/pt-z0
485   c[12]=0.; // 1/pt-snp
486   c[13]=0.;                                 // 1/pt-tgl
487   c[14]=0.;        // 1/pt-1/pt
488
489   return kTRUE;
490 }
491
492 TObjArray AliTPCTracklet::CreateTracklets(const AliTPCseed *track,
493                                           TrackType type,
494                                           Bool_t storeClusters,
495                                           Int_t minClusters,
496                                           Int_t maxTracklets) {
497 // The tracklet factory: It creates several tracklets out of a track. They
498 // are created for sectors that fullfill the constraint of having enough
499 // clusters inside. Futhermore you can specify the maximum amount of
500 // tracklets that are to be created.
501 // The tracklets appear in a sorted fashion, beginning with those having the
502 // most clusters.
503
504   Int_t sectors[72]={0};
505   for (Int_t i=0;i<160;++i) {
506     AliTPCclusterMI *c=track->GetClusterPointer(i);
507     if (c)
508       ++sectors[c->GetDetector()];
509   }
510   Int_t indices[72];
511   TMath::Sort(72,sectors,indices);
512   TObjArray tracklets;
513   if (maxTracklets>72) maxTracklets=72; // just to protect against "users".
514   for (Int_t i=0;i<maxTracklets&&sectors[indices[i]]>=minClusters;++i) {
515     tracklets.Add(new AliTPCTracklet(track,indices[i],type,storeClusters));
516   }
517   return tracklets;
518 }
519
520 TObjArray AliTPCTracklet::CreateTracklets(const TObjArray &/*clusters*/,
521                                           TrackType /*type*/,
522                                           Bool_t /*storeClusters*/,
523                                           Int_t /*minClusters*/,
524                                           Int_t /*maxTracklets*/) {
525   // TODO!
526
527   TObjArray tracklets;
528   return tracklets;
529 }
530
531 Bool_t AliTPCTracklet::PropagateToMeanX(const AliTPCTracklet &t1,
532                                         const AliTPCTracklet &t2,
533                                         AliExternalTrackParam *&t1m,
534                                         AliExternalTrackParam *&t2m) {
535   // This function propagates two Tracklets to a common x-coordinate. This
536   // x is dermined as the one that is in the middle of the two tracklets (they
537   // are assumed to live on two distinct x-intervalls).
538   // The inner parametrisation of the outer Tracklet and the outer 
539   // parametrisation of the inner Tracklet are used and propagated to this
540   // common x. This result is saved not inside the Tracklets but two new
541   // ExternalTrackParams are created (that means you might want to delete
542   // them).
543   // In the case that the alpha angles of the Tracklets differ both angles
544   // are tried out for this propagation.
545   // In case of any failure kFALSE is returned, no AliExternalTrackParam
546   // is created und the pointers are set to 0.
547
548   if (t1.GetInner() && t1.GetOuter() && 
549       t2.GetInner() && t2.GetOuter()) {
550     if (t1.GetOuter()->GetX()<t2.GetInner()->GetX()) {
551       t1m=new AliExternalTrackParam(*t1.GetOuter());
552       t2m=new AliExternalTrackParam(*t2.GetInner());
553     }
554     else {
555       t1m=new AliExternalTrackParam(*t1.GetInner());
556       t2m=new AliExternalTrackParam(*t2.GetOuter());
557     }
558     Double_t mx=.5*(t1m->GetX()+t2m->GetX());
559     Double_t b1,b2;
560     Double_t xyz[3];
561     t1m->GetXYZ(xyz);
562     b1=GetBz(xyz);
563     t2m->GetXYZ(xyz);
564     b2=GetBz(xyz);
565     if (t1m->Rotate(t2m->GetAlpha()) 
566         && t1m->PropagateTo(mx,b1) 
567         && t2m->PropagateTo(mx,b2));
568     else
569       if (t2m->Rotate(t1m->GetAlpha())
570           && t1m->PropagateTo(mx,b1) 
571           && t2m->PropagateTo(mx,b2));
572       else {
573         delete t1m;
574         delete t2m;
575         t1m=t2m=0;
576       }
577   }
578   else {
579     t1m=t2m=0;
580   }
581   return t1m&&t2m;
582 }
583
584 double AliTPCTracklet::GetBz(Double_t *xyz) {
585   if (AliTracker::UniformField())
586     return AliTracker::GetBz();
587   else
588     return AliTracker::GetBz(xyz);
589 }
590
591 void AliTPCTracklet::RandomND(Int_t ndim,const Double_t *p,const Double_t *c,
592                               Double_t *x) {
593   // This function generates a n-dimensional random variable x with mean
594   // p and covariance c.
595   // That is done using the cholesky decomposition of the covariance matrix,
596   // Begin_Latex C=U^{t} U End_Latex, with Begin_Latex U End_Latex being an
597   // upper triangular matrix. Given a vector v of iid gausian random variables
598   // with variance 1 one obtains the asked result as: Begin_Latex x=U^t v 
599   // End_Latex.
600   // c is expected to be in a lower triangular format:
601   // c[0]
602   // c[1] c[2]
603   // c[3] c[4] c[5]
604   // etc.
605   static TRandom3 random;
606   Double_t *c2= new Double_t[ndim*ndim];
607   Int_t k=0;
608   for (Int_t i=0;i<ndim;++i)
609     for (Int_t j=0;j<=i;++j)
610       c2[i*ndim+j]=c2[j*ndim+i]=c[k++];
611   TMatrixDSym cm(ndim,c2);
612   delete[] c2;
613   TDecompChol chol(cm);
614   chol.Decompose();
615   const TVectorD pv(ndim);
616   const_cast<TVectorD*>(&pv)->Use(ndim,const_cast<Double_t*>(p));
617   TVectorD xv(ndim);
618   xv.Use(ndim,x);
619   for (Int_t i=0;i<ndim;++i)
620     xv[i]=random.Gaus();
621   TMatrixD L=chol.GetU();
622   L.T();
623   xv=L*xv+pv;
624 }
625
626 TEllipse AliTPCTracklet::ErrorEllipse(Double_t x,Double_t y,
627                                       Double_t sx,Double_t sy,Double_t sxy) {
628   /* Begin_Latex
629      r_{1,2}=1/2 (a+c#pm#sqrt{(a-c)^{2}+(2b)^{2}})
630   End_Latex */
631   Double_t det1=1./(sx*sy-sxy*sxy);
632   Double_t a=sy*det1;
633   Double_t b=-sxy*det1;
634   Double_t c=sx*det1;
635   Double_t d=c-a;
636   Double_t s=TMath::Sqrt(d*d+4.*b*b);
637   Double_t r1=TMath::Sqrt(.5*(a+c-s));
638   Double_t r2=TMath::Sqrt(.5*(a+c+s));
639   Double_t alpha=.5*TMath::ATan2(2.*b,d);
640   return TEllipse(x,y,r1,r2,0.,360.,alpha*TMath::RadToDeg());
641 }
642
643 void AliTPCTracklet::Test(const char* filename) {
644   /*
645     aliroot
646     AliTPCTracklet::Test("");
647     TFile f("AliTPCTrackletDebug.root");
648     TTree *t=f.Get("AliTPCTrackletDebug");
649     t->Draw("p0:p4");
650     TEllipse e=AliTPCTracklet::ErrorEllipse(0.,0.,4.,1.,1.8);
651     e.Draw();
652  */
653   TTreeSRedirector ds(filename);
654   Double_t p[5]={0.};
655   Double_t c[15]={4.,
656                   0.,4.,
657                   0.,0.,9.,
658                   0.,0.,0.,16.,
659                   1.8,0.,0.,0.,1.};
660   for (Int_t i=0;i<10000;++i) {
661     Double_t x[5];
662     RandomND(5,p,c,x);
663     ds<<"AliTPCTrackletDebug"
664       <<"p0="<<x[0]
665       <<"p1="<<x[1]
666       <<"p2="<<x[2]
667       <<"p3="<<x[3]
668       <<"p4="<<x[4]
669       <<"\n";
670   }
671
672   /*
673   Double_t b;
674   Double_t x=0.;
675   Double_t alpha=0.;
676   Double_t param[5]={0.};
677   Double_t covar[15]={1.,
678                       0.,1.,
679                       0.,0.,1.,
680                       0.,0.,0.,1.,
681                       0.,0.,0.,0.,1.};
682   AliExternalTrackParam track(x,alpha,param,covar);
683
684   
685
686   for (Int_t i=0;i<points.GetNPoints();++i) {
687     Double_t x=0.;
688     Double_t alpha=0.;
689     Double_t param[5]={0.};
690     Double_t covar[15]={1.,
691                         0.,1.,
692                         0.,0.,1.,
693                         0.,0.,0.,1.,
694                         0.,0.,0.,0.,1.};
695     AliExternalTrackParam track(x,alpha,param,covar);
696     for (x=90.;x<250.;x+=1.) {
697       track.PropagateTo(x,b);
698       AliTPCclusterMI c();
699     }
700   }
701   */
702 }