]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCtracker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Corrections for tracking in arbitrary magnenetic field. Changes towards a concept...
[u/mrichter/AliRoot.git] / TPC / AliTPCtracker.cxx
CommitLineData
73042f01 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$Log$
be9c5115 18Revision 1.5 2000/12/20 07:51:59 kowal2
19Changes suggested by Alessandra and Paolo to avoid overlapped
20data fields in encapsulated classes.
21
d4cf1daa 22Revision 1.4 2000/11/02 07:27:16 kowal2
23code corrections
24
98ab53f4 25Revision 1.2 2000/06/30 12:07:50 kowal2
26Updated from the TPC-PreRelease branch
27
73042f01 28Revision 1.1.2.1 2000/06/25 08:53:55 kowal2
29Splitted from AliTPCtracking
30
31*/
32
33//-------------------------------------------------------
34// Implementation of the TPC tracker
35//
36// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
37//-------------------------------------------------------
38
73042f01 39#include <TObjArray.h>
40#include <TFile.h>
be5b9287 41#include <TTree.h>
42#include <iostream.h>
43
44#include "AliTPCtracker.h"
45#include "AliTPCcluster.h"
b9de75e1 46#include "AliTPCParam.h"
73042f01 47#include "AliTPCClustersRow.h"
48
b9de75e1 49//_____________________________________________________________________________
50AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
51fkNIS(par->GetNInnerSector()/2),
52fkNOS(par->GetNOuterSector()/2)
53{
54 //---------------------------------------------------------------------
55 // The main TPC tracker constructor
56 //---------------------------------------------------------------------
57 fInnerSec=new AliTPCSector[fkNIS];
58 fOuterSec=new AliTPCSector[fkNOS];
59
60 Int_t i;
61 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
62 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
63
64 fN=0; fSectors=0;
65
66 fClustersArray.Setup(par);
67 fClustersArray.SetClusterType("AliTPCcluster");
68 fClustersArray.ConnectTree("Segment Tree");
69
70 fSeeds=0;
71}
72
73//_____________________________________________________________________________
74AliTPCtracker::~AliTPCtracker() {
75 //------------------------------------------------------------------
76 // TPC tracker destructor
77 //------------------------------------------------------------------
78 delete[] fInnerSec;
79 delete[] fOuterSec;
80 delete fSeeds;
81}
73042f01 82
83//_____________________________________________________________________________
84Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
85{
86 //
87 // Parametrised error of the cluster reconstruction (pad direction)
88 //
89 // Sigma rphi
90 const Float_t kArphi=0.41818e-2;
91 const Float_t kBrphi=0.17460e-4;
92 const Float_t kCrphi=0.30993e-2;
93 const Float_t kDrphi=0.41061e-3;
94
95 pt=TMath::Abs(pt)*1000.;
96 Double_t x=r/pt;
97 tgl=TMath::Abs(tgl);
98 Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
99 if (s<0.4e-3) s=0.4e-3;
100 s*=1.3; //Iouri Belikov
101
102 return s;
103}
104
105//_____________________________________________________________________________
106Double_t SigmaZ2(Double_t r, Double_t tgl)
107{
108 //
109 // Parametrised error of the cluster reconstruction (drift direction)
110 //
111 // Sigma z
112 const Float_t kAz=0.39614e-2;
113 const Float_t kBz=0.22443e-4;
114 const Float_t kCz=0.51504e-1;
115
116
117 tgl=TMath::Abs(tgl);
118 Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
119 if (s<0.4e-3) s=0.4e-3;
120 s*=1.3; //Iouri Belikov
121
122 return s;
123}
124
125//_____________________________________________________________________________
126inline Double_t f1(Double_t x1,Double_t y1,
127 Double_t x2,Double_t y2,
128 Double_t x3,Double_t y3)
129{
130 //-----------------------------------------------------------------
131 // Initial approximation of the track curvature
132 //-----------------------------------------------------------------
133 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
134 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
135 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
136 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
137 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
138
139 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
140
141 return -xr*yr/sqrt(xr*xr+yr*yr);
142}
143
144
145//_____________________________________________________________________________
146inline Double_t f2(Double_t x1,Double_t y1,
147 Double_t x2,Double_t y2,
148 Double_t x3,Double_t y3)
149{
150 //-----------------------------------------------------------------
151 // Initial approximation of the track curvature times center of curvature
152 //-----------------------------------------------------------------
153 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
154 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
155 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
156 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
157 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
158
159 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
160
161 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
162}
163
164//_____________________________________________________________________________
165inline Double_t f3(Double_t x1,Double_t y1,
166 Double_t x2,Double_t y2,
167 Double_t z1,Double_t z2)
168{
169 //-----------------------------------------------------------------
170 // Initial approximation of the tangent of the track dip angle
171 //-----------------------------------------------------------------
172 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
173}
174
175//_____________________________________________________________________________
b9de75e1 176void AliTPCtracker::LoadOuterSectors() {
177 //-----------------------------------------------------------------
178 // This function fills outer TPC sectors with clusters.
179 //-----------------------------------------------------------------
180 UInt_t index;
181 Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
182 for (Int_t i=0; i<j; i++) {
183 AliSegmentID *s=fClustersArray.LoadEntry(i);
184 Int_t sec,row;
185 AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
186 par->AdjustSectorRow(s->GetID(),sec,row);
187 if (sec<fkNIS*2) continue;
188 AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
189 Int_t ncl=clrow->GetArray()->GetEntriesFast();
190 while (ncl--) {
191 AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
192 index=(((sec<<8)+row)<<16)+ncl;
193 fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
194 }
195 }
196
197 fN=fkNOS;
198 fSectors=fOuterSec;
199}
200
201//_____________________________________________________________________________
202void AliTPCtracker::UnloadOuterSectors() {
203 //-----------------------------------------------------------------
204 // This function clears outer TPC sectors.
205 //-----------------------------------------------------------------
206 Int_t nup=fOuterSec->GetNRows();
207 for (Int_t i=0; i<fkNOS; i++) {
208 for (Int_t j=0; j<nup; j++) {
209 if (fClustersArray.GetRow(i+fkNIS*2,j))
210 fClustersArray.ClearRow(i+fkNIS*2,j);
211 if (fClustersArray.GetRow(i+fkNIS*2+fkNOS,j))
212 fClustersArray.ClearRow(i+fkNIS*2+fkNOS,j);
213 }
214 }
215
216 fN=0;
217 fSectors=0;
218}
219
220//_____________________________________________________________________________
221void AliTPCtracker::LoadInnerSectors() {
222 //-----------------------------------------------------------------
223 // This function fills inner TPC sectors with clusters.
224 //-----------------------------------------------------------------
225 UInt_t index;
226 Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
227 for (Int_t i=0; i<j; i++) {
228 AliSegmentID *s=fClustersArray.LoadEntry(i);
229 Int_t sec,row;
230 AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
231 par->AdjustSectorRow(s->GetID(),sec,row);
232 if (sec>=fkNIS*2) continue;
233 AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
234 Int_t ncl=clrow->GetArray()->GetEntriesFast();
235 while (ncl--) {
236 AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
237 index=(((sec<<8)+row)<<16)+ncl;
238 fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
239 }
240 }
241
242 fN=fkNIS;
243 fSectors=fInnerSec;
244}
245
246//_____________________________________________________________________________
247void AliTPCtracker::UnloadInnerSectors() {
248 //-----------------------------------------------------------------
249 // This function clears inner TPC sectors.
250 //-----------------------------------------------------------------
251 Int_t nlow=fInnerSec->GetNRows();
252 for (Int_t i=0; i<fkNIS; i++) {
253 for (Int_t j=0; j<nlow; j++) {
254 if (fClustersArray.GetRow(i,j)) fClustersArray.ClearRow(i,j);
255 if (fClustersArray.GetRow(i+fkNIS,j)) fClustersArray.ClearRow(i+fkNIS,j);
256 }
257 }
258
259 fN=0;
260 fSectors=0;
261}
262
263//_____________________________________________________________________________
264Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
73042f01 265 //-----------------------------------------------------------------
266 // This function tries to find a track prolongation.
267 //-----------------------------------------------------------------
b9de75e1 268 Double_t xt=t.GetX();
269 const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
270 Int_t(0.5*fSectors->GetNRows());
73042f01 271 Int_t tryAgain=kSKIP;
73042f01 272
b9de75e1 273 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
274 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
275 if (alpha < 0. ) alpha += 2.*TMath::Pi();
276 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
277
278 for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) {
279 Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
73042f01 280 if (!t.PropagateTo(x)) return 0;
281
282 AliTPCcluster *cl=0;
283 UInt_t index=0;
b9de75e1 284 Double_t maxchi2=kMaxCHI2;
285 const AliTPCRow &krow=fSectors[s][nr];
286 Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
287 Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
73042f01 288 Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
289 Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
290
b9de75e1 291 if (road>kMaxROAD) {
73042f01 292 if (t.GetNumberOfClusters()>4)
293 cerr<<t.GetNumberOfClusters()
294 <<"FindProlongation warning: Too broad road !\n";
295 return 0;
296 }
297
298 if (krow) {
299 for (Int_t i=krow.Find(y-road); i<krow; i++) {
be5b9287 300 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
73042f01 301 if (c->GetY() > y+road) break;
302 if (c->IsUsed()) continue;
303 if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
304 Double_t chi2=t.GetPredictedChi2(c);
305 if (chi2 > maxchi2) continue;
306 maxchi2=chi2;
307 cl=c;
308 index=krow.GetIndex(i);
309 }
310 }
311 if (cl) {
b9de75e1 312 Float_t l=fSectors->GetPadPitchWidth();
73042f01 313 t.SetSampledEdx(cl->GetQ()/l,t.GetNumberOfClusters());
be9c5115 314 if (!t.Update(cl,maxchi2,index)) {
315 if (!tryAgain--) return 0;
316 } else tryAgain=kSKIP;
73042f01 317 } else {
318 if (tryAgain==0) break;
319 if (y > ymax) {
b9de75e1 320 s = (s+1) % fN;
321 if (!t.Rotate(fSectors->GetAlpha())) return 0;
73042f01 322 } else if (y <-ymax) {
b9de75e1 323 s = (s-1+fN) % fN;
324 if (!t.Rotate(-fSectors->GetAlpha())) return 0;
73042f01 325 }
326 tryAgain--;
327 }
328 }
329
330 return 1;
b9de75e1 331}
332
333//_____________________________________________________________________________
334Int_t AliTPCtracker::FollowBackProlongation
335(AliTPCseed& seed, const AliTPCtrack &track) {
336 //-----------------------------------------------------------------
337 // This function propagates tracks back through the TPC
338 //-----------------------------------------------------------------
339 Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
340 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
341 if (alpha < 0. ) alpha += 2.*TMath::Pi();
342 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
343
344 Int_t idx=-1, sec=-1, row=-1;
345 Int_t nc=seed.GetLabel(); //index of the cluster to start with
346 if (nc--) {
347 idx=track.GetClusterIndex(nc);
348 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
349 }
350 if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; }
351 else { if (sec < 2*fkNIS) row=-1; }
352
353 Int_t nr=fSectors->GetNRows();
354 for (Int_t i=0; i<nr; i++) {
355 Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
356
357 if (!seed.PropagateTo(x)) return 0;
358
359 Double_t y=seed.GetY();
360 if (y > ymax) {
361 s = (s+1) % fN;
362 if (!seed.Rotate(fSectors->GetAlpha())) return 0;
363 } else if (y <-ymax) {
364 s = (s-1+fN) % fN;
365 if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
366 }
367
368 AliTPCcluster *cl=0;
369 Int_t index=0;
370 Double_t maxchi2=kMaxCHI2;
371 Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
372 Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
373 Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
374 Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
375 if (road>kMaxROAD) {
376 cerr<<seed.GetNumberOfClusters()
377 <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n";
378 return 0;
379 }
380
381
382 Int_t accepted=seed.GetNumberOfClusters();
383 if (row==i) {
384 //try to accept already found cluster
385 AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
386 Double_t chi2;
387 if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
388 index=idx; cl=c; maxchi2=chi2;
389 } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
390
391 if (nc--) {
392 idx=track.GetClusterIndex(nc);
393 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
394 }
395 if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; }
396 else { if (sec < 2*fkNIS) row=-1; }
73042f01 397
b9de75e1 398 }
399 if (!cl) {
400 //try to fill the gap
401 const AliTPCRow &krow=fSectors[s][i];
402 if (accepted>27)
403 if (krow) {
404 for (Int_t i=krow.Find(y-road); i<krow; i++) {
405 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
406 if (c->GetY() > y+road) break;
407 if (c->IsUsed()) continue;
408 if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
409 Double_t chi2=seed.GetPredictedChi2(c);
410 if (chi2 > maxchi2) continue;
411 maxchi2=chi2;
412 cl=c;
413 index=krow.GetIndex(i);
414 }
415 }
416 }
417
418 if (cl) {
419 Float_t l=fSectors->GetPadPitchWidth();
420 seed.SetSampledEdx(cl->GetQ()/l,seed.GetNumberOfClusters());
421 seed.Update(cl,maxchi2,index);
422 }
423
424 }
425
426 seed.SetLabel(nc);
427
428 return 1;
73042f01 429}
430
431//_____________________________________________________________________________
b9de75e1 432void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
73042f01 433 //-----------------------------------------------------------------
434 // This function creates track seeds.
435 //-----------------------------------------------------------------
b9de75e1 436 if (fSeeds==0) fSeeds=new TObjArray(15000);
437
73042f01 438 Double_t x[5], c[15];
439
b9de75e1 440 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
73042f01 441 Double_t cs=cos(alpha), sn=sin(alpha);
442
b9de75e1 443 Double_t x1 =fOuterSec->GetX(i1);
444 Double_t xx2=fOuterSec->GetX(i2);
73042f01 445
b9de75e1 446 for (Int_t ns=0; ns<fkNOS; ns++) {
447 Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
448 Int_t nm=fOuterSec[ns][i2];
449 Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
450 const AliTPCRow& kr1=fOuterSec[ns][i1];
73042f01 451 for (Int_t is=0; is < kr1; is++) {
452 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
453 for (Int_t js=0; js < nl+nm+nu; js++) {
454 const AliTPCcluster *kcl;
455 Double_t x2, y2, z2;
456 Double_t x3=0.,y3=0.;
457
458 if (js<nl) {
b9de75e1 459 const AliTPCRow& kr2=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
73042f01 460 kcl=kr2[js];
461 y2=kcl->GetY(); z2=kcl->GetZ();
462 x2= xx2*cs+y2*sn;
463 y2=-xx2*sn+y2*cs;
464 } else
465 if (js<nl+nm) {
b9de75e1 466 const AliTPCRow& kr2=fOuterSec[ns][i2];
73042f01 467 kcl=kr2[js-nl];
468 x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
469 } else {
b9de75e1 470 const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];
73042f01 471 kcl=kr2[js-nl-nm];
472 y2=kcl->GetY(); z2=kcl->GetZ();
473 x2=xx2*cs-y2*sn;
474 y2=xx2*sn+y2*cs;
475 }
476
477 Double_t zz=z1 - z1/x1*(x1-x2);
478 if (TMath::Abs(zz-z2)>5.) continue;
479
480 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
481 if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
482
483 x[0]=y1;
484 x[1]=z1;
b9de75e1 485 x[4]=f1(x1,y1,x2,y2,x3,y3);
486 if (TMath::Abs(x[4]) >= 0.0066) continue;
be5b9287 487 x[2]=f2(x1,y1,x2,y2,x3,y3);
b9de75e1 488 //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
489 x[3]=f3(x1,y1,x2,y2,z1,z2);
490 if (TMath::Abs(x[3]) > 1.2) continue;
be5b9287 491 Double_t a=asin(x[2]);
b9de75e1 492 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
73042f01 493 if (TMath::Abs(zv)>10.) continue;
494
495 Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
496 Double_t sy2=kcl->GetSigmaY2(), sz2=kcl->GetSigmaZ2();
497 Double_t sy3=100*0.025, sy=0.1, sz=0.1;
498
b9de75e1 499 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
500 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
501 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
be5b9287 502 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
503 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
b9de75e1 504 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
505 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
506 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
507 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
508 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
73042f01 509
510 c[0]=sy1;
511 c[1]=0.; c[2]=sz1;
b9de75e1 512 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
513 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
514 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
515 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
516 c[13]=f30*sy1*f40+f32*sy2*f42;
517 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
73042f01 518
519 UInt_t index=kr1.GetIndex(is);
520 AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
b9de75e1 521 Float_t l=fOuterSec->GetPadPitchWidth();
73042f01 522 track->SetSampledEdx(kr1[is]->GetQ()/l,0);
523
b9de75e1 524 Int_t rc=FollowProlongation(*track, i2);
be9c5115 525 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
b9de75e1 526 else fSeeds->AddLast(track);
73042f01 527 }
528 }
529 }
530}
531
532//_____________________________________________________________________________
b9de75e1 533Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
73042f01 534 //-----------------------------------------------------------------
b9de75e1 535 // This function reades track seeds.
73042f01 536 //-----------------------------------------------------------------
537 TDirectory *savedir=gDirectory;
538
b9de75e1 539 TFile *in=(TFile*)inp;
540 if (!in->IsOpen()) {
541 cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
be5b9287 542 return 1;
73042f01 543 }
544
b9de75e1 545 in->cd();
546 TTree *seedTree=(TTree*)in->Get("Seeds");
547 if (!seedTree) {
548 cerr<<"AliTPCtracker::ReadSeeds(): ";
549 cerr<<"can't get a tree with track seeds !\n";
550 return 2;
551 }
552 AliTPCtrack *seed=new AliTPCtrack;
553 seedTree->SetBranchAddress("tracks",&seed);
554
555 if (fSeeds==0) fSeeds=new TObjArray(15000);
73042f01 556
b9de75e1 557 Int_t n=(Int_t)seedTree->GetEntries();
558 for (Int_t i=0; i<n; i++) {
559 seedTree->GetEvent(i);
560 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
561 }
562
563 delete seed;
564 savedir->cd();
73042f01 565
b9de75e1 566 return 0;
567}
73042f01 568
b9de75e1 569//_____________________________________________________________________________
570Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
571 //-----------------------------------------------------------------
572 // This is a track finder.
573 //-----------------------------------------------------------------
574 TDirectory *savedir=gDirectory;
73042f01 575
b9de75e1 576 if (inp) {
577 TFile *in=(TFile*)inp;
578 if (!in->IsOpen()) {
579 cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
580 return 1;
581 }
582 }
73042f01 583
b9de75e1 584 if (!out->IsOpen()) {
585 cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
586 return 2;
73042f01 587 }
588
b9de75e1 589 out->cd();
590 TTree tracktree("TPCf","Tree with TPC tracks");
591 AliTPCtrack *iotrack=0;
592 tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
593
594 LoadOuterSectors();
595
73042f01 596 //find track seeds
b9de75e1 597 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
73042f01 598 Int_t nrows=nlow+nup;
b9de75e1 599 if (fSeeds==0) {
600 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
601 MakeSeeds(nup-1, nup-1-gap);
602 MakeSeeds(nup-1-shift, nup-1-shift-gap);
603 }
604 fSeeds->Sort();
73042f01 605
606 //tracking in outer sectors
b9de75e1 607 Int_t nseed=fSeeds->GetEntriesFast();
608 Int_t i;
73042f01 609 for (i=0; i<nseed; i++) {
b9de75e1 610 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
611 if (FollowProlongation(t)) {
612 UseClusters(&t);
73042f01 613 continue;
614 }
b9de75e1 615 delete fSeeds->RemoveAt(i);
73042f01 616 }
b9de75e1 617 UnloadOuterSectors();
73042f01 618
619 //tracking in inner sectors
b9de75e1 620 LoadInnerSectors();
73042f01 621 Int_t found=0;
622 for (i=0; i<nseed; i++) {
b9de75e1 623 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
73042f01 624 if (!pt) continue;
625 Int_t nc=t.GetNumberOfClusters();
626
b9de75e1 627 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
73042f01 628 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
629 if (alpha < 0. ) alpha += 2.*TMath::Pi();
b9de75e1 630 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
73042f01 631
b9de75e1 632 alpha=ns*fInnerSec->GetAlpha() + fInnerSec->GetAlphaShift() - t.GetAlpha();
73042f01 633
634 if (t.Rotate(alpha)) {
b9de75e1 635 if (FollowProlongation(t)) {
73042f01 636 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
637 t.CookdEdx();
73042f01 638 iotrack=pt;
639 tracktree.Fill();
b9de75e1 640 UseClusters(&t,nc);
641 cerr<<found++<<'\r';
73042f01 642 }
643 }
644 }
b9de75e1 645 delete fSeeds->RemoveAt(i);
73042f01 646 }
b9de75e1 647 UnloadInnerSectors();
648
73042f01 649 tracktree.Write();
650
b9de75e1 651 cerr<<"Number of found tracks : "<<found<<endl;
652
653 savedir->cd();
654
655 return 0;
656}
657
658//_____________________________________________________________________________
659Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
660 //-----------------------------------------------------------------
661 // This function propagates tracks back through the TPC.
662 //-----------------------------------------------------------------
663 fSeeds=new TObjArray(15000);
664
665 TFile *in=(TFile*)inp;
666 TDirectory *savedir=gDirectory;
667
668 if (!in->IsOpen()) {
669 cerr<<"AliTPCtracker::PropagateBack(): ";
670 cerr<<"file with back propagated ITS tracks is not open !\n";
671 return 1;
73042f01 672 }
673
b9de75e1 674 if (!out->IsOpen()) {
675 cerr<<"AliTPCtracker::PropagateBack(): ";
676 cerr<<"file for back propagated TPC tracks is not open !\n";
677 return 2;
678 }
679
680 in->cd();
681 TTree *bckTree=(TTree*)in->Get("ITSb");
682 if (!bckTree) {
683 cerr<<"AliTPCtracker::PropagateBack() ";
684 cerr<<"can't get a tree with back propagated ITS tracks !\n";
685 return 3;
686 }
687 AliTPCtrack *bckTrack=new AliTPCtrack;
688 bckTree->SetBranchAddress("tracks",&bckTrack);
689
690 TTree *tpcTree=(TTree*)in->Get("TPCf");
691 if (!tpcTree) {
692 cerr<<"AliTPCtracker::PropagateBack() ";
693 cerr<<"can't get a tree with TPC tracks !\n";
694 return 4;
695 }
696 AliTPCtrack *tpcTrack=new AliTPCtrack;
697 tpcTree->SetBranchAddress("tracks",&tpcTrack);
698
699//*** Prepare an array of tracks to be back propagated
700 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
701 Int_t nrows=nlow+nup;
702
703 TObjArray tracks(15000);
704 Int_t i=0,j=0;
705 Int_t tpcN=(Int_t)tpcTree->GetEntries();
706 Int_t bckN=(Int_t)bckTree->GetEntries();
707 if (j<bckN) bckTree->GetEvent(j++);
708 for (i=0; i<tpcN; i++) {
709 tpcTree->GetEvent(i);
710 Double_t alpha=tpcTrack->GetAlpha();
711 if (j<bckN &&
712 TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
713 if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
714 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
715 bckTree->GetEvent(j++);
716 } else {
717 tpcTrack->ResetCovariance();
718 fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
719 }
720 tracks.AddLast(new AliTPCtrack(*tpcTrack));
721 }
722
723 out->cd();
724 TTree backTree("TPCb","Tree with back propagated TPC tracks");
725 AliTPCtrack *otrack=0;
726 backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
727
728//*** Back propagation through inner sectors
729 LoadInnerSectors();
730 Int_t nseed=fSeeds->GetEntriesFast();
731 for (i=0; i<nseed; i++) {
732 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
733 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
73042f01 734
b9de75e1 735 Int_t nc=t.GetNumberOfClusters();
736 s.SetLabel(nc-1); //set number of the cluster to start with
73042f01 737
b9de75e1 738 if (FollowBackProlongation(s,t)) {
739 UseClusters(&s);
740 continue;
741 }
742 delete fSeeds->RemoveAt(i);
743 }
744 UnloadInnerSectors();
745
746//*** Back propagation through outer sectors
747 LoadOuterSectors();
748 Int_t found=0;
749 for (i=0; i<nseed; i++) {
750 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
751 if (!ps) continue;
752 Int_t nc=s.GetNumberOfClusters();
753 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
754
755 Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
756 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
757 if (alpha < 0. ) alpha += 2.*TMath::Pi();
758 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
759
760 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
761 alpha-=s.GetAlpha();
762
763 if (s.Rotate(alpha)) {
764 if (FollowBackProlongation(s,t)) {
765 if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
766 s.CookdEdx();
767 s.SetLabel(t.GetLabel());
768 UseClusters(&s,nc);
769 otrack=ps;
770 backTree.Fill();
771 cerr<<found++<<'\r';
772 continue;
773 }
774 }
775 }
776 delete fSeeds->RemoveAt(i);
777 }
778 UnloadOuterSectors();
779
780 backTree.Write();
73042f01 781 savedir->cd();
b9de75e1 782 cerr<<"Number of seeds: "<<nseed<<endl;
783 cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
784 cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
785
786 delete bckTrack;
787 delete tpcTrack;
be5b9287 788
789 return 0;
73042f01 790}
791
b9de75e1 792//_________________________________________________________________________
793AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
794 //--------------------------------------------------------------------
795 // Return pointer to a given cluster
796 //--------------------------------------------------------------------
797 Int_t sec=(index&0xff000000)>>24;
798 Int_t row=(index&0x00ff0000)>>16;
799 Int_t ncl=(index&0x0000ffff)>>00;
800
801 AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
802 return (AliCluster*)(*clrow)[ncl];
803}
804
805//__________________________________________________________________________
806void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
807 //--------------------------------------------------------------------
808 //This function "cooks" a track label. If label<0, this track is fake.
809 //--------------------------------------------------------------------
810 Int_t noc=t->GetNumberOfClusters();
811 Int_t *lb=new Int_t[noc];
812 Int_t *mx=new Int_t[noc];
813 AliCluster **clusters=new AliCluster*[noc];
814
815 Int_t i;
816 for (i=0; i<noc; i++) {
817 lb[i]=mx[i]=0;
818 Int_t index=t->GetClusterIndex(i);
819 clusters[i]=GetCluster(index);
820 }
821
822 Int_t lab=123456789;
823 for (i=0; i<noc; i++) {
824 AliCluster *c=clusters[i];
825 lab=TMath::Abs(c->GetLabel(0));
826 Int_t j;
827 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
828 lb[j]=lab;
829 (mx[j])++;
830 }
831
832 Int_t max=0;
833 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
834
835 for (i=0; i<noc; i++) {
836 AliCluster *c=clusters[i];
837 if (TMath::Abs(c->GetLabel(1)) == lab ||
838 TMath::Abs(c->GetLabel(2)) == lab ) max++;
839 }
840
841 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
842
843 else {
844 Int_t tail=Int_t(0.10*noc);
845 max=0;
846 for (i=1; i<=tail; i++) {
847 AliCluster *c=clusters[noc-i];
848 if (lab == TMath::Abs(c->GetLabel(0)) ||
849 lab == TMath::Abs(c->GetLabel(1)) ||
850 lab == TMath::Abs(c->GetLabel(2))) max++;
851 }
852 if (max < Int_t(0.5*tail)) lab=-lab;
853 }
854
855 t->SetLabel(lab);
856
857 delete[] lb;
858 delete[] mx;
859 delete[] clusters;
860}
861
862//_________________________________________________________________________
863void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
864 //-----------------------------------------------------------------------
865 // Setup inner sector
866 //-----------------------------------------------------------------------
867 if (f==0) {
868 fAlpha=par->GetInnerAngle();
869 fAlphaShift=par->GetInnerAngleShift();
870 fPadPitchWidth=par->GetInnerPadPitchWidth();
871 fPadPitchLength=par->GetInnerPadPitchLength();
872
873 fN=par->GetNRowLow();
874 fRow=new AliTPCRow[fN];
875 for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
876 } else {
877 fAlpha=par->GetOuterAngle();
878 fAlphaShift=par->GetOuterAngleShift();
879 fPadPitchWidth=par->GetOuterPadPitchWidth();
880 fPadPitchLength=par->GetOuterPadPitchLength();
881
882 fN=par->GetNRowUp();
883 fRow=new AliTPCRow[fN];
884 for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiUp(i));
885 }
886}
887
73042f01 888//_________________________________________________________________________
889void
890AliTPCtracker::AliTPCRow::InsertCluster(const AliTPCcluster* c, UInt_t index) {
891 //-----------------------------------------------------------------------
892 // Insert a cluster into this pad row in accordence with its y-coordinate
893 //-----------------------------------------------------------------------
b9de75e1 894 if (fN==kMaxClusterPerRow) {
73042f01 895 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
896 }
897 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
898 Int_t i=Find(c->GetY());
899 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
900 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
901 fIndex[i]=index; fClusters[i]=c; fN++;
902}
903
904//___________________________________________________________________
905Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
906 //-----------------------------------------------------------------------
907 // Return the index of the nearest cluster
908 //-----------------------------------------------------------------------
909 if (y <= fClusters[0]->GetY()) return 0;
910 if (y > fClusters[fN-1]->GetY()) return fN;
911 Int_t b=0, e=fN-1, m=(b+e)/2;
912 for (; b<e; m=(b+e)/2) {
913 if (y > fClusters[m]->GetY()) b=m+1;
914 else e=m;
915 }
916 return m;
917}
918
919//_____________________________________________________________________________
920void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
921 //-----------------------------------------------------------------
922 // This funtion calculates dE/dX within the "low" and "up" cuts.
923 //-----------------------------------------------------------------
924 Int_t i;
be9c5115 925 Int_t nc=GetNumberOfClusters();
73042f01 926
927 Int_t swap;//stupid sorting
928 do {
929 swap=0;
be9c5115 930 for (i=0; i<nc-1; i++) {
d4cf1daa 931 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
932 Float_t tmp=fdEdxSample[i];
933 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
73042f01 934 swap++;
935 }
936 } while (swap);
937
be9c5115 938 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
73042f01 939 Float_t dedx=0;
d4cf1daa 940 for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
73042f01 941 dedx /= (nu-nl+1);
942 SetdEdx(dedx);
943}
944
73042f01 945
946