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