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