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