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