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