Line data Source code
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 "AliTPCreco.h"
29 : #include "AliESDVertex.h"
30 : #include "AliTracker.h"
31 : #include "TTreeStream.h"
32 : #include "TRandom3.h"
33 : #include "TDecompChol.h"
34 :
35 : #include <iostream>
36 : using namespace std;
37 :
38 16 : ClassImp(AliTPCTracklet)
39 :
40 : const Double_t AliTPCTracklet::kB2C=0.299792458e-3;
41 : Float_t AliTPCTracklet::fgEdgeCutY=3;
42 : Float_t AliTPCTracklet::fgEdgeCutX=0;
43 :
44 0 : AliTPCTracklet::AliTPCTracklet()
45 0 : : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(-1),fOuter(0),
46 0 : fInner(0),fPrimary(0) {
47 : ////
48 : // The default constructor. It is intended to be used for I/O only.
49 : ////
50 0 : }
51 :
52 0 : AliTPCTracklet::AliTPCTracklet(const AliTPCseed *track,Int_t sector,
53 : TrackType type,Bool_t storeClusters)
54 0 : : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(sector),fOuter(0),
55 0 : fInner(0),fPrimary(0) {
56 : ////
57 : // Contructor for a tracklet out of a track. Only clusters within a given
58 : // sector are used.
59 : ///
60 :
61 : //TODO: only kalman works
62 :
63 0 : for (Int_t i=0;i<kMaxRow;++i) {
64 0 : AliTPCclusterMI *c=track->GetClusterPointer(i);
65 0 : if (c && RejectCluster(c)) continue;
66 0 : if (c&&c->GetDetector()==sector)
67 0 : ++fNClusters;
68 0 : }
69 :
70 0 : if (storeClusters) {
71 0 : fClusters=new AliTPCclusterMI[fNClusters];
72 0 : for (Int_t i=0;i<kMaxRow;++i) {
73 0 : AliTPCclusterMI *c=track->GetClusterPointer(i);
74 0 : if (c && RejectCluster(c)) continue;
75 0 : if (c&&c->GetDetector()==sector)
76 0 : fClusters[fNStoredClusters]=*c;
77 0 : ++fNStoredClusters;
78 0 : }
79 0 : }
80 :
81 0 : switch (type) {
82 : case kKalman:
83 0 : FitKalman(track,sector);
84 : break;
85 : case kLinear:
86 : case kQuadratic:
87 0 : FitLinear(track,sector,type);
88 : break;
89 : case kRiemann:
90 0 : FitRiemann(track,sector);
91 : break;
92 : }
93 :
94 0 : }
95 :
96 0 : AliTPCTracklet::AliTPCTracklet(const TObjArray &/*clusters*/,Int_t sector,
97 : TrackType /*type*/,Bool_t /*storeClusters*/)
98 0 : : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(sector),fOuter(0),
99 0 : fInner(0),fPrimary(0) {
100 : //TODO: write it!
101 0 : }
102 :
103 : AliTPCTracklet::AliTPCTracklet(const AliTPCTracklet &t)
104 0 : : TObject(t),fNClusters(t.fNClusters),fNStoredClusters(t.fNStoredClusters),fClusters(0),
105 0 : fSector(t.fSector),fOuter(0),fInner(0),
106 0 : fPrimary(0) {
107 : ////
108 : // The copy constructor. You can copy tracklets!
109 : ////
110 :
111 0 : if (t.fClusters) {
112 0 : fClusters=new AliTPCclusterMI[t.fNStoredClusters];
113 0 : for (int i=0;i<t.fNStoredClusters;++i)
114 0 : fClusters[i]=t.fClusters[i];
115 0 : }
116 0 : if (t.fOuter)
117 0 : fOuter=new AliExternalTrackParam(*t.fOuter);
118 0 : if (t.fInner)
119 0 : fInner=new AliExternalTrackParam(*t.fInner);
120 0 : if (t.fPrimary)
121 0 : fPrimary=new AliExternalTrackParam(*t.fPrimary);
122 0 : }
123 :
124 : AliTPCTracklet& AliTPCTracklet::operator=(const AliTPCTracklet &t) {
125 : ////
126 : // The assignment constructor. You can assign tracklets!
127 : ////
128 0 : if (this!=&t) {
129 0 : fNClusters=t.fNClusters;
130 0 : delete fClusters;
131 0 : if (t.fClusters) {
132 0 : fClusters=new AliTPCclusterMI[t.fNStoredClusters];
133 0 : for (int i=0;i<t.fNStoredClusters;++i)
134 0 : fClusters[i]=t.fClusters[i];
135 0 : }
136 : else
137 0 : fClusters=0;
138 0 : fSector=t.fSector;
139 0 : if (t.fOuter) {
140 0 : if (fOuter)
141 0 : *fOuter=*t.fOuter;
142 : else
143 0 : fOuter=new AliExternalTrackParam(*t.fOuter);
144 : }
145 : else {
146 0 : delete fOuter;
147 0 : fOuter=0;
148 : }
149 :
150 0 : if (t.fInner) {
151 0 : if (fInner)
152 0 : *fInner=*t.fInner;
153 : else
154 0 : fInner=new AliExternalTrackParam(*t.fInner);
155 : }
156 : else {
157 0 : delete fInner;
158 0 : fInner=0;
159 : }
160 :
161 0 : if (t.fPrimary) {
162 0 : if (fPrimary)
163 0 : *fPrimary=*t.fPrimary;
164 : else
165 0 : fPrimary=new AliExternalTrackParam(*t.fPrimary);
166 : }
167 : else {
168 0 : delete fPrimary;
169 0 : fPrimary=0;
170 : }
171 : }
172 0 : return *this;
173 0 : }
174 :
175 0 : AliTPCTracklet::~AliTPCTracklet() {
176 : //
177 : // The destructor. Yes, you can even destruct tracklets.
178 : //
179 0 : delete fClusters;
180 0 : delete fOuter;
181 0 : delete fInner;
182 0 : delete fPrimary;
183 0 : }
184 :
185 :
186 :
187 :
188 :
189 : void AliTPCTracklet::FitKalman(const AliTPCseed *seed,Int_t sector) {
190 : //
191 : // Fit using Kalman filter
192 : //
193 0 : AliTPCseed *track=new AliTPCseed(*seed);
194 0 : if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
195 0 : delete track;
196 0 : return;
197 : }
198 : // fit from inner to outer row
199 0 : Double_t covar[15];
200 0 : for (Int_t i=0;i<15;i++) covar[i]=0;
201 0 : covar[0]=1.*1.;
202 0 : covar[2]=1.*1.;
203 0 : covar[5]=1.*1./(64.*64.);
204 0 : covar[9]=1.*1./(64.*64.);
205 0 : covar[14]=0; // keep pt
206 : Float_t xmin=1000, xmax=-10000;
207 : Int_t imin=158, imax=0;
208 0 : for (Int_t i=0;i<kMaxRow;i++) {
209 0 : AliTPCclusterMI *c=track->GetClusterPointer(i);
210 0 : if (!c) continue;
211 0 : if (c->GetDetector()!=sector) continue;
212 0 : if (c->GetX()<xmin) xmin=c->GetX();
213 0 : if (c->GetX()>xmax) xmax=c->GetX();
214 0 : if (i<imin) imin=i;
215 0 : if (i>imax) imax=i;
216 0 : }
217 0 : if(imax-imin<10) {
218 0 : delete track;
219 0 : return;
220 : }
221 :
222 0 : for (Float_t x=track->GetX(); x<xmin; x++) track->PropagateTo(x);
223 0 : track->AddCovariance(covar);
224 : //
225 0 : AliExternalTrackParam paramIn;
226 0 : AliExternalTrackParam paramOut;
227 : Bool_t isOK=kTRUE;
228 : //
229 : //
230 : //
231 0 : for (Int_t i=imin; i<=imax; i++){
232 0 : AliTPCclusterMI *c=track->GetClusterPointer(i);
233 0 : if (!c) continue;
234 0 : Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
235 0 : Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
236 0 : AliTPCseed::GetError(c, track,cov[0],cov[2]);
237 0 : cov[0]*=cov[0];
238 0 : cov[2]*=cov[2];
239 0 : if (!track->PropagateTo(r[0])) {
240 : isOK=kFALSE;
241 0 : break;
242 : }
243 0 : if (RejectCluster(c,track)) continue;
244 0 : if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
245 0 : }
246 0 : if (!isOK) { delete track; return;}
247 0 : track->AddCovariance(covar);
248 : //
249 : //
250 : //
251 0 : for (Int_t i=imax; i>=imin; i--){
252 0 : AliTPCclusterMI *c=track->GetClusterPointer(i);
253 0 : if (!c) continue;
254 0 : Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
255 0 : Double_t cov[3]={0.01,0.,0.01};
256 0 : AliTPCseed::GetError(c, track,cov[0],cov[2]);
257 0 : cov[0]*=cov[0];
258 0 : cov[2]*=cov[2];
259 0 : if (!track->PropagateTo(r[0])) {
260 : isOK=kFALSE;
261 0 : break;
262 : }
263 0 : if (RejectCluster(c,track)) continue;
264 0 : if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
265 0 : }
266 0 : if (!isOK) { delete track; return;}
267 0 : paramIn = *track;
268 0 : track->AddCovariance(covar);
269 : //
270 : //
271 0 : for (Int_t i=imin; i<=imax; i++){
272 0 : AliTPCclusterMI *c=track->GetClusterPointer(i);
273 0 : if (!c) continue;
274 0 : Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
275 0 : Double_t cov[3]={0.01,0.,0.01};
276 0 : AliTPCseed::GetError(c, track,cov[0],cov[2]);
277 0 : cov[0]*=cov[0];
278 0 : cov[2]*=cov[2];
279 0 : if (!track->PropagateTo(r[0])) {
280 : isOK=kFALSE;
281 0 : break;
282 : }
283 0 : if (RejectCluster(c,track)) continue;
284 0 : if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
285 0 : }
286 0 : if (!isOK) { delete track; return;}
287 0 : paramOut=*track;
288 : //
289 : //
290 : //
291 0 : fOuter=new AliExternalTrackParam(paramOut);
292 0 : fInner=new AliExternalTrackParam(paramIn);
293 : //
294 0 : delete track;
295 0 : }
296 :
297 :
298 :
299 :
300 : void AliTPCTracklet::FitLinear(const AliTPCseed *track,Int_t sector,
301 : TrackType type) {
302 0 : TLinearFitter fy(1);
303 0 : TLinearFitter fz(1);
304 0 : fy.StoreData(kFALSE);
305 0 : fz.StoreData(kFALSE);
306 0 : switch (type) {
307 : case kLinear:
308 0 : fy.SetFormula("1 ++ x");
309 0 : fz.SetFormula("1 ++ x");
310 : break;
311 : case kQuadratic:
312 0 : fy.SetFormula("1 ++ x ++ x*x");
313 0 : 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 0 : for (Int_t i=0;i<kMaxRow;++i) {
322 0 : AliTPCclusterMI *c=track->GetClusterPointer(i);
323 0 : if (c && RejectCluster(c)) continue;
324 0 : if (c&&c->GetDetector()==sector) {
325 0 : Double_t x=c->GetX();
326 0 : fy.AddPoint(&x,c->GetY());
327 0 : fz.AddPoint(&x,c->GetZ());
328 0 : xmax=TMath::Max(xmax,x);
329 0 : xmin=TMath::Min(xmin,x);
330 0 : }
331 0 : }
332 0 : fy.Eval();
333 0 : fz.Eval();
334 0 : Double_t a[3]={fy.GetParameter(0),
335 0 : fy.GetParameter(1),
336 0 : type==kQuadratic?fy.GetParameter(2):0.};
337 0 : Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
338 0 : fy.GetCovarianceMatrixElement(1,0),
339 0 : fy.GetCovarianceMatrixElement(1,1),
340 0 : type==kQuadratic?fy.GetCovarianceMatrixElement(2,0):0.,
341 0 : type==kQuadratic?fy.GetCovarianceMatrixElement(2,1):0.,
342 0 : type==kQuadratic?fy.GetCovarianceMatrixElement(2,2):0.};
343 0 : for (int i=0;i<6;++i) ca[i]*=fy.GetChisquare()/fNClusters;
344 0 : Double_t b[2]={fz.GetParameter(0),
345 0 : fz.GetParameter(1)};
346 0 : Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
347 0 : fz.GetCovarianceMatrixElement(1,0),
348 0 : fz.GetCovarianceMatrixElement(1,1)};
349 0 : for (int i=0;i<3;++i) cb[i]*=fz.GetChisquare()/fNClusters;
350 0 : Double_t p[5];
351 0 : Double_t c[15];
352 0 : Double_t alpha=track->GetAlpha();
353 0 : Quadratic2Helix(a,ca,b,cb,0.,p,c);
354 0 : fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
355 0 : Quadratic2Helix(a,ca,b,cb,xmin,p,c);
356 0 : fInner=new AliExternalTrackParam(xmin,alpha,p,c);
357 0 : Quadratic2Helix(a,ca,b,cb,xmax,p,c);
358 0 : fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
359 0 : }
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 0 : Double_t a0=x0*x0*a[2] + x0*a[1] + a[0];
371 0 : Double_t a1=2.*x0*a[2] + a[1];
372 : Double_t a2= a[2];
373 0 : Double_t ca00=ca[0]+x0*(2.*ca[1]+x0*(ca[2]+2.*ca[3]+x0*(2.*ca[4]+x0*ca[5])));
374 0 : Double_t ca10=ca[1]+x0*(ca[2]+2.*ca[3]+x0*(3.*ca[4]+x0*2.*ca[5]));
375 0 : Double_t ca11=ca[2]+x0*4.*(ca[4]+x0*ca[5]);
376 0 : Double_t ca20=ca[3]+x0*(ca[4]+x0*ca[5]);
377 0 : Double_t ca21=ca[3]+x0*2.*ca[5];
378 : Double_t ca22=ca[5];
379 :
380 0 : Double_t b0=x0*b[1] + b[0];
381 : Double_t b1= b[1];
382 0 : Double_t cb00=cb[0]+x0*(2.*cb[1]+x0*cb[2]);
383 0 : Double_t cb10=cb[1]+x0*cb[2];
384 : Double_t cb11=cb[2];
385 :
386 : // transform to helix parameters
387 0 : Double_t f =1.+a1*a1;
388 0 : Double_t f2 =f*f;
389 0 : Double_t fi =1./f;
390 0 : Double_t fi12=TMath::Sqrt(fi);
391 0 : Double_t fi32=fi*fi12;
392 0 : Double_t fi2 =fi*fi;
393 0 : Double_t fi52=fi2*fi12;
394 0 : Double_t fi3 =fi2*fi;
395 0 : Double_t fi5 =fi2*fi3;
396 :
397 0 : Double_t xyz[3]={0.}; // TODO...
398 0 : Double_t fc=1./(GetBz(xyz)*kB2C);
399 :
400 0 : p[0]=a0; // y0
401 0 : p[1]=b0; // z0
402 0 : p[2]=a1*fi12; // snp
403 0 : p[3]=b1; // tgl
404 0 : p[4]=2.*a2*fi32*fc; // 1/pt
405 :
406 0 : c[0] =ca00; // y0-y0
407 0 : c[1] =0.; // z0-y0
408 0 : c[2] =cb00; // z0-z0
409 0 : c[3] =ca10*fi32; // snp-y0
410 0 : c[4] =0.; // snp-z0
411 0 : c[5] =ca11*fi3; // snp-snp
412 0 : c[6] =0.; // tgl-y0
413 0 : c[7] =cb10; // tgl-z0
414 0 : c[8] =0.; // tgl-snp
415 0 : c[9] =cb11; // tgl-tgl
416 0 : c[10]=2.*(-3.*a1*a2*ca10+f*ca20)*fi3*fc; // 1/pt-y0
417 0 : c[11]=0.; // 1/pt-z0
418 0 : c[12]=2.*(-3.*a1*a2*ca11+f*ca21)*fi52*fc; // 1/pt-snp
419 0 : c[13]=0.; // 1/pt-tgl
420 0 : c[14]=(-12.*a1*a2*(-3.*a1*a2*ca11+2.*f*ca21)+4.*f2*ca22)*fi5
421 0 : *fc*fc; // 1/pt-1/pt
422 0 : }
423 :
424 :
425 : void AliTPCTracklet::FitRiemann(const AliTPCseed *track,Int_t sector) {
426 0 : TLinearFitter fy(2);
427 0 : fy.StoreData(kFALSE);
428 0 : fy.SetFormula("hyp2");
429 : Double_t xmax=-1.;
430 : Double_t xmin=1000.;
431 0 : for (Int_t i=0;i<kMaxRow;++i) {
432 0 : AliTPCclusterMI *c=track->GetClusterPointer(i);
433 0 : if (c && RejectCluster(c)) continue;
434 0 : if (c&&c->GetDetector()==sector) {
435 0 : Double_t x=c->GetX();
436 0 : Double_t y=c->GetY();
437 0 : Double_t xy[2]={x,y};
438 0 : Double_t r=x*x+y*y;
439 : Double_t errx=1.,erry=1.;//TODO!
440 0 : Double_t err=TMath::Sqrt(4.*x*x*errx+4.*y*y*erry);
441 : err=1.;
442 0 : fy.AddPoint(xy,r,err);
443 0 : xmax=TMath::Max(xmax,x);
444 0 : xmin=TMath::Min(xmin,x);
445 0 : }
446 0 : }
447 0 : fy.Eval();
448 0 : Double_t a[3]={fy.GetParameter(0),
449 0 : fy.GetParameter(1),
450 0 : fy.GetParameter(2)};
451 0 : Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
452 0 : fy.GetCovarianceMatrixElement(1,0),
453 0 : fy.GetCovarianceMatrixElement(1,1),
454 0 : fy.GetCovarianceMatrixElement(2,0),
455 0 : fy.GetCovarianceMatrixElement(2,1),
456 0 : fy.GetCovarianceMatrixElement(2,2)};
457 :
458 0 : TLinearFitter fz(1);
459 0 : fz.StoreData(kFALSE);
460 0 : fz.SetFormula("hyp1");
461 0 : 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 0 : Double_t phi=0.;
465 0 : for (Int_t i=0;i<kMaxRow;++i) {
466 0 : AliTPCclusterMI *c=track->GetClusterPointer(i);
467 0 : if (c && RejectCluster(c)) continue;
468 0 : if (c&&c->GetDetector()==sector) {
469 0 : Double_t x=c->GetX();
470 0 : Double_t y=c->GetY();
471 0 : Double_t dx=x-oldx;
472 0 : Double_t dy=y-oldy;
473 0 : phi+=2.*TMath::Abs(TMath::ATan2(.5*TMath::Sqrt(dx*dx+dy*dy),R));
474 : Double_t err=1.;
475 0 : fz.AddPoint(&phi,c->GetZ(),err);
476 : oldx=x;
477 : oldy=y;
478 0 : }
479 0 : }
480 0 : fz.Eval();
481 0 : Double_t b[2]={fz.GetParameter(0),
482 0 : fz.GetParameter(1)};
483 0 : Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
484 0 : fz.GetCovarianceMatrixElement(1,0),
485 0 : fz.GetCovarianceMatrixElement(1,1)};
486 :
487 0 : Double_t p[5];
488 0 : Double_t c[15];
489 0 : Double_t alpha=track->GetAlpha();
490 0 : if (Riemann2Helix(a,ca,b,cb,0.,p,c))
491 0 : fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
492 0 : if (Riemann2Helix(a,ca,b,cb,xmin,p,c))
493 0 : fInner=new AliExternalTrackParam(xmin,alpha,p,c);
494 0 : if (Riemann2Helix(a,ca,b,cb,xmax,p,c))
495 0 : fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
496 0 : }
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 0 : Double_t xr0=.5*a[1];
505 0 : Double_t yr0=.5*a[2];
506 0 : Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
507 0 : Double_t dx=x0-xr0;
508 0 : if (dx*dx>=R*R) return kFALSE;
509 0 : Double_t dy=TMath::Sqrt((R-dx)*(R+dx)); //sign!!
510 0 : if (TMath::Abs(yr0+dy)>TMath::Abs(yr0-dy))
511 0 : dy=-dy;
512 0 : Double_t y0=yr0+dy;
513 0 : Double_t tgp=-dx/dy; //TODO: dy!=0
514 0 : Double_t z0=b[0]+TMath::ATan(tgp)*b[1];
515 0 : Double_t xyz[3]={x0,y0,z0};
516 0 : Double_t fc=1./(GetBz(xyz)*kB2C);
517 : fc=1;
518 0 : p[0]=y0; // y0
519 0 : p[1]=z0; // z0
520 0 : p[2]=tgp/TMath::Sqrt(1.+tgp*tgp); // snp
521 0 : p[3]=b[1]; // tgl
522 0 : p[4]=1./R*fc; // 1/pt
523 :
524 0 : c[0] =0.; // y0-y0
525 0 : c[1] =0.; // z0-y0
526 0 : c[2] =0.; // z0-z0
527 0 : c[3] =0.; // snp-y0
528 0 : c[4] =0.; // snp-z0
529 0 : c[5] =0.; // snp-snp
530 0 : c[6] =0.; // tgl-y0
531 0 : c[7] =0.; // tgl-z0
532 0 : c[8] =0.; // tgl-snp
533 0 : c[9] =0.; // tgl-tgl
534 0 : c[10]=0.; // 1/pt-y0
535 0 : c[11]=0.; // 1/pt-z0
536 0 : c[12]=0.; // 1/pt-snp
537 0 : c[13]=0.; // 1/pt-tgl
538 0 : c[14]=0.; // 1/pt-1/pt
539 :
540 : return kTRUE;
541 0 : }
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 0 : Int_t sectors[72]={0};
556 0 : for (Int_t i=0;i<kMaxRow;++i) {
557 0 : AliTPCclusterMI *c=track->GetClusterPointer(i);
558 0 : if (c && RejectCluster(c)) continue;
559 0 : if (c)
560 0 : ++sectors[c->GetDetector()];
561 0 : }
562 0 : Int_t indices[72];
563 0 : TMath::Sort(72,sectors,indices);
564 0 : TObjArray tracklets;
565 0 : if (maxTracklets>72) maxTracklets=72; // just to protect against "users".
566 0 : for (Int_t i=0;i<maxTracklets&§ors[indices[i]]>=minClusters;++i) {
567 0 : tracklets.Add(new AliTPCTracklet(track,indices[i],type,storeClusters));
568 : }
569 : return tracklets;
570 0 : }
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 0 : TObjArray tracklets;
580 : return tracklets;
581 0 : }
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 0 : if (t1.GetInner() && t1.GetOuter() &&
601 0 : t2.GetInner() && t2.GetOuter()) {
602 0 : if (t1.GetOuter()->GetX()<t2.GetInner()->GetX()) {
603 0 : t1m=new AliExternalTrackParam(*t1.GetOuter());
604 0 : t2m=new AliExternalTrackParam(*t2.GetInner());
605 0 : }
606 : else {
607 0 : t1m=new AliExternalTrackParam(*t1.GetInner());
608 0 : t2m=new AliExternalTrackParam(*t2.GetOuter());
609 : }
610 0 : Double_t mx=.5*(t1m->GetX()+t2m->GetX());
611 : //Double_t b1,b2;
612 0 : Double_t xyz[3];
613 0 : t1m->GetXYZ(xyz);
614 : //b1=GetBz(xyz);
615 0 : Double_t b1[3]; AliTracker::GetBxByBz(xyz,b1);
616 0 : t2m->GetXYZ(xyz);
617 : //b2=GetBz(xyz);
618 0 : Double_t b2[3]; AliTracker::GetBxByBz(xyz,b2);
619 0 : if (t1m->Rotate(t2m->GetAlpha())
620 : //&& t1m->PropagateTo(mx,b1)
621 : //&& t2m->PropagateTo(mx,b2));
622 0 : && t1m->PropagateToBxByBz(mx,b1)
623 0 : && t2m->PropagateToBxByBz(mx,b2));
624 : else
625 0 : if (t2m->Rotate(t1m->GetAlpha())
626 : //&& t1m->PropagateTo(mx,b1)
627 : //&& t2m->PropagateTo(mx,b2));
628 0 : && t1m->PropagateToBxByBz(mx,b1)
629 0 : && t2m->PropagateToBxByBz(mx,b2));
630 : else {
631 0 : delete t1m;
632 0 : delete t2m;
633 0 : t1m=t2m=0;
634 : }
635 0 : }
636 : else {
637 0 : t1m=t2m=0;
638 : }
639 0 : return t1m&&t2m;
640 0 : }
641 :
642 : double AliTPCTracklet::GetBz(Double_t *xyz)
643 : {
644 0 : 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 0 : static TRandom3 random;
662 0 : Double_t *c2= new Double_t[ndim*ndim];
663 : Int_t k=0;
664 0 : for (Int_t i=0;i<ndim;++i)
665 0 : for (Int_t j=0;j<=i;++j)
666 0 : c2[i*ndim+j]=c2[j*ndim+i]=c[k++];
667 0 : TMatrixDSym cm(ndim,c2);
668 0 : delete[] c2;
669 0 : TDecompChol chol(cm);
670 0 : chol.Decompose();
671 0 : const TVectorD pv(ndim);
672 0 : const_cast<TVectorD*>(&pv)->Use(ndim,const_cast<Double_t*>(p));
673 0 : TVectorD xv(ndim);
674 0 : xv.Use(ndim,x);
675 0 : for (Int_t i=0;i<ndim;++i)
676 0 : xv[i]=random.Gaus();
677 0 : TMatrixD L=chol.GetU();
678 0 : L.T();
679 0 : xv=L*xv+pv;
680 0 : }
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 0 : Double_t det1=1./(sx*sy-sxy*sxy);
688 0 : Double_t a=sy*det1;
689 0 : Double_t b=-sxy*det1;
690 0 : Double_t c=sx*det1;
691 0 : Double_t d=c-a;
692 0 : Double_t s=TMath::Sqrt(d*d+4.*b*b);
693 0 : Double_t r1=TMath::Sqrt(.5*(a+c-s));
694 0 : Double_t r2=TMath::Sqrt(.5*(a+c+s));
695 0 : Double_t alpha=.5*TMath::ATan2(2.*b,d);
696 0 : return TEllipse(x,y,r1,r2,0.,360.,alpha*TMath::RadToDeg());
697 0 : }
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 0 : TTreeSRedirector ds(filename);
710 0 : Double_t p[5]={0.};
711 0 : 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 0 : for (Int_t i=0;i<10000;++i) {
717 0 : Double_t x[5];
718 0 : RandomND(5,p,c,x);
719 0 : ds<<"AliTPCTrackletDebug"
720 0 : <<"p0="<<x[0]
721 0 : <<"p1="<<x[1]
722 0 : <<"p2="<<x[2]
723 0 : <<"p3="<<x[3]
724 0 : <<"p4="<<x[4]
725 0 : <<"\n";
726 0 : }
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 0 : }
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 0 : Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
768 0 : Float_t dist = edgeY - TMath::Abs(cl->GetY());
769 0 : if (param) dist = edgeY - TMath::Abs(param->GetY());
770 0 : if (dist<fgEdgeCutY) isReject=kTRUE;
771 0 : if (cl->GetType()<0) isReject=kTRUE;
772 0 : return isReject;
773 : }
|
