Bug fix - Teminate - as virtula function
[u/mrichter/AliRoot.git] / TPC / AliTPCTracklet.cxx
CommitLineData
69209e2c 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
69209e2c 24#include "AliTPCTracklet.h"
25#include "TObjArray.h"
9318a5b4 26#include "TLinearFitter.h"
69209e2c 27#include "AliTPCseed.h"
28#include "AliESDVertex.h"
9318a5b4 29#include "AliTracker.h"
30#include "TTreeStream.h"
31#include "TRandom3.h"
32#include "TDecompChol.h"
33
34#include <iostream>
35using namespace std;
69209e2c 36
37ClassImp(AliTPCTracklet)
38
9318a5b4 39const Double_t AliTPCTracklet::kB2C=0.299792458e-3;
40
69209e2c 41AliTPCTracklet::AliTPCTracklet()
9318a5b4 42 : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(-1),fOuter(0),
43 fInner(0),fPrimary(0) {
69209e2c 44 ////
45 // The default constructor. It is intended to be used for I/O only.
46 ////
47}
48
9318a5b4 49AliTPCTracklet::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) {
69209e2c 53 ////
54 // Contructor for a tracklet out of a track. Only clusters within a given
55 // sector are used.
56 ///
69209e2c 57
9318a5b4 58 //TODO: only kalman works
59
69209e2c 60 for (Int_t i=0;i<160;++i) {
9318a5b4 61 AliTPCclusterMI *c=track->GetClusterPointer(i);
62 if (c&&c->GetDetector()==sector)
63 ++fNClusters;
69209e2c 64 }
9318a5b4 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;
69209e2c 73 }
74 }
9318a5b4 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;
69209e2c 87 }
69209e2c 88
9318a5b4 89}
69209e2c 90
c94a79e1 91AliTPCTracklet::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) {
9318a5b4 95 //TODO: write it!
69209e2c 96}
97
98AliTPCTracklet::AliTPCTracklet(const AliTPCTracklet &t)
c94a79e1 99 : TObject(t),fNClusters(t.fNClusters),fNStoredClusters(t.fNStoredClusters),fClusters(0),
9318a5b4 100 fSector(t.fSector),fOuter(0),fInner(0),
69209e2c 101 fPrimary(0) {
102 ////
103 // The copy constructor. You can copy tracklets!
104 ////
105
9318a5b4 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 }
69209e2c 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
119AliTPCTracklet& AliTPCTracklet::operator=(const AliTPCTracklet &t) {
120 ////
121 // The assignment constructor. You can assign tracklets!
122 ////
69209e2c 123 if (this!=&t) {
9318a5b4 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;
69209e2c 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
171AliTPCTracklet::~AliTPCTracklet() {
172 //
173 // The destructor. Yes, you can even destruct tracklets.
174 //
9318a5b4 175 delete fClusters;
69209e2c 176 delete fOuter;
177 delete fInner;
178 delete fPrimary;
179}
180
9318a5b4 181void 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...
7eaa723e 242 if (fPrimary) if (!fPrimary->Rotate(fInner->GetAlpha())) {
9318a5b4 243 delete fPrimary;
244 fPrimary=0;
245 }
246 }
247 delete innerSeed;
248
249 delete t;
250}
251
252void 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;
c94a79e1 267 case kKalman:
268 case kRiemann:
269 break;
9318a5b4 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
312void 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
376void 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
c94a79e1 447Bool_t AliTPCTracklet::Riemann2Helix(Double_t *a,Double_t */*ca*/,
448 Double_t *b,Double_t */*cb*/,
9318a5b4 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
492TObjArray AliTPCTracklet::CreateTracklets(const AliTPCseed *track,
493 TrackType type,
494 Bool_t storeClusters,
69209e2c 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) {
9318a5b4 506 AliTPCclusterMI *c=track->GetClusterPointer(i);
69209e2c 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) {
9318a5b4 515 tracklets.Add(new AliTPCTracklet(track,indices[i],type,storeClusters));
69209e2c 516 }
517 return tracklets;
518}
9318a5b4 519
c94a79e1 520TObjArray AliTPCTracklet::CreateTracklets(const TObjArray &/*clusters*/,
521 TrackType /*type*/,
522 Bool_t /*storeClusters*/,
523 Int_t /*minClusters*/,
524 Int_t /*maxTracklets*/) {
525 // TODO!
526
71d2416b 527 TObjArray tracklets;
528 return tracklets;
9318a5b4 529}
530
531Bool_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
584double AliTPCTracklet::GetBz(Double_t *xyz) {
585 if (AliTracker::UniformField())
586 return AliTracker::GetBz();
587 else
588 return AliTracker::GetBz(xyz);
589}
590
591void 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
626TEllipse 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
643void 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 */
c94a79e1 653 TTreeSRedirector ds(filename);
9318a5b4 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}