]> 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
Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
[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$
11c7b548 18Revision 1.17 2002/03/18 17:59:13 kowal2
19Chnges in the pad geometry - 3 pad lengths introduced.
20
f03e3423 21Revision 1.16 2001/11/08 16:39:03 hristov
22Additional protection (M.Masera)
23
1c8d9aad 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
7f6ddf58 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
f38c8ae5 30Revision 1.13 2001/07/23 12:01:30 hristov
31Initialisation added
32
7afb86d1 33Revision 1.12 2001/07/20 14:32:44 kowal2
34Processing of many events possible now
35
afc42102 36Revision 1.11 2001/05/23 08:50:10 hristov
37Weird inline removed
38
bdbd0f7a 39Revision 1.10 2001/05/16 14:57:25 alibrary
40New files for folders and Stack
41
9e1a0ddb 42Revision 1.9 2001/05/11 07:16:56 hristov
43Fix needed on Sun and Alpha
44
3d39dd8b 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
be9c5115 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
d4cf1daa 52Revision 1.4 2000/11/02 07:27:16 kowal2
53code corrections
54
98ab53f4 55Revision 1.2 2000/06/30 12:07:50 kowal2
56Updated from the TPC-PreRelease branch
57
73042f01 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
73042f01 69#include <TObjArray.h>
70#include <TFile.h>
be5b9287 71#include <TTree.h>
72#include <iostream.h>
73
74#include "AliTPCtracker.h"
75#include "AliTPCcluster.h"
b9de75e1 76#include "AliTPCParam.h"
73042f01 77#include "AliTPCClustersRow.h"
78
b9de75e1 79//_____________________________________________________________________________
afc42102 80AliTPCtracker::AliTPCtracker(const AliTPCParam *par, Int_t eventn):
7f6ddf58 81AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
b9de75e1 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");
b9de75e1 97
afc42102 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;
b9de75e1 109 fSeeds=0;
110}
111
112//_____________________________________________________________________________
113AliTPCtracker::~AliTPCtracker() {
114 //------------------------------------------------------------------
115 // TPC tracker destructor
116 //------------------------------------------------------------------
117 delete[] fInnerSec;
118 delete[] fOuterSec;
1c8d9aad 119 if (fSeeds) {
120 fSeeds->Delete();
121 delete fSeeds;
122 }
b9de75e1 123}
73042f01 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//_____________________________________________________________________________
bdbd0f7a 168Double_t f1(Double_t x1,Double_t y1,
73042f01 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//_____________________________________________________________________________
bdbd0f7a 188Double_t f2(Double_t x1,Double_t y1,
73042f01 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//_____________________________________________________________________________
bdbd0f7a 207Double_t f3(Double_t x1,Double_t y1,
73042f01 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//_____________________________________________________________________________
b9de75e1 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) {
73042f01 307 //-----------------------------------------------------------------
308 // This function tries to find a track prolongation.
309 //-----------------------------------------------------------------
b9de75e1 310 Double_t xt=t.GetX();
311 const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
312 Int_t(0.5*fSectors->GetNRows());
73042f01 313 Int_t tryAgain=kSKIP;
73042f01 314
b9de75e1 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
f03e3423 320 Int_t nrows=fSectors->GetRowNumber(xt)-1;
321 for (Int_t nr=nrows; nr>=rf; nr--) {
b9de75e1 322 Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
73042f01 323 if (!t.PropagateTo(x)) return 0;
324
325 AliTPCcluster *cl=0;
326 UInt_t index=0;
b9de75e1 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);
73042f01 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
b9de75e1 334 if (road>kMaxROAD) {
73042f01 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++) {
be5b9287 343 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
73042f01 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) {
b9de75e1 355 Float_t l=fSectors->GetPadPitchWidth();
73042f01 356 t.SetSampledEdx(cl->GetQ()/l,t.GetNumberOfClusters());
be9c5115 357 if (!t.Update(cl,maxchi2,index)) {
358 if (!tryAgain--) return 0;
359 } else tryAgain=kSKIP;
73042f01 360 } else {
361 if (tryAgain==0) break;
362 if (y > ymax) {
b9de75e1 363 s = (s+1) % fN;
364 if (!t.Rotate(fSectors->GetAlpha())) return 0;
73042f01 365 } else if (y <-ymax) {
b9de75e1 366 s = (s-1+fN) % fN;
367 if (!t.Rotate(-fSectors->GetAlpha())) return 0;
73042f01 368 }
369 tryAgain--;
370 }
371 }
372
373 return 1;
b9de75e1 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; }
73042f01 440
b9de75e1 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;
73042f01 472}
473
474//_____________________________________________________________________________
b9de75e1 475void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
73042f01 476 //-----------------------------------------------------------------
477 // This function creates track seeds.
478 //-----------------------------------------------------------------
b9de75e1 479 if (fSeeds==0) fSeeds=new TObjArray(15000);
480
73042f01 481 Double_t x[5], c[15];
482
b9de75e1 483 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
73042f01 484 Double_t cs=cos(alpha), sn=sin(alpha);
485
b9de75e1 486 Double_t x1 =fOuterSec->GetX(i1);
487 Double_t xx2=fOuterSec->GetX(i2);
73042f01 488
b9de75e1 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];
73042f01 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;
7f6ddf58 499 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
73042f01 500
501 if (js<nl) {
b9de75e1 502 const AliTPCRow& kr2=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
73042f01 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) {
b9de75e1 509 const AliTPCRow& kr2=fOuterSec[ns][i2];
73042f01 510 kcl=kr2[js-nl];
511 x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
512 } else {
b9de75e1 513 const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];
73042f01 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
7f6ddf58 520 Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
73042f01 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;
b9de75e1 528 x[4]=f1(x1,y1,x2,y2,x3,y3);
529 if (TMath::Abs(x[4]) >= 0.0066) continue;
be5b9287 530 x[2]=f2(x1,y1,x2,y2,x3,y3);
b9de75e1 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;
be5b9287 534 Double_t a=asin(x[2]);
b9de75e1 535 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
7f6ddf58 536 if (TMath::Abs(zv-z3)>10.) continue;
73042f01 537
538 Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
539 Double_t sy2=kcl->GetSigmaY2(), sz2=kcl->GetSigmaZ2();
7f6ddf58 540 Double_t sy3=400*3./12., sy=0.1, sz=0.1;
73042f01 541
b9de75e1 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;
be5b9287 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;
b9de75e1 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;
73042f01 552
553 c[0]=sy1;
554 c[1]=0.; c[2]=sz1;
b9de75e1 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;
73042f01 561
562 UInt_t index=kr1.GetIndex(is);
563 AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
b9de75e1 564 Float_t l=fOuterSec->GetPadPitchWidth();
73042f01 565 track->SetSampledEdx(kr1[is]->GetQ()/l,0);
566
b9de75e1 567 Int_t rc=FollowProlongation(*track, i2);
be9c5115 568 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
b9de75e1 569 else fSeeds->AddLast(track);
73042f01 570 }
571 }
572 }
573}
574
575//_____________________________________________________________________________
b9de75e1 576Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
73042f01 577 //-----------------------------------------------------------------
b9de75e1 578 // This function reades track seeds.
73042f01 579 //-----------------------------------------------------------------
580 TDirectory *savedir=gDirectory;
581
b9de75e1 582 TFile *in=(TFile*)inp;
583 if (!in->IsOpen()) {
584 cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
be5b9287 585 return 1;
73042f01 586 }
587
b9de75e1 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);
73042f01 599
b9de75e1 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;
f38c8ae5 607
608 delete seedTree; //Thanks to Mariana Bondila
609
b9de75e1 610 savedir->cd();
73042f01 611
b9de75e1 612 return 0;
613}
73042f01 614
b9de75e1 615//_____________________________________________________________________________
616Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
617 //-----------------------------------------------------------------
618 // This is a track finder.
619 //-----------------------------------------------------------------
620 TDirectory *savedir=gDirectory;
73042f01 621
b9de75e1 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 }
73042f01 629
b9de75e1 630 if (!out->IsOpen()) {
631 cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
632 return 2;
73042f01 633 }
634
b9de75e1 635 out->cd();
7f6ddf58 636
637 char tname[100];
638 sprintf(tname,"TreeT_TPC_%d",fEventN);
639 TTree tracktree(tname,"Tree with TPC tracks");
b9de75e1 640 AliTPCtrack *iotrack=0;
641 tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
642
643 LoadOuterSectors();
644
73042f01 645 //find track seeds
b9de75e1 646 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
73042f01 647 Int_t nrows=nlow+nup;
b9de75e1 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();
73042f01 654
655 //tracking in outer sectors
b9de75e1 656 Int_t nseed=fSeeds->GetEntriesFast();
657 Int_t i;
73042f01 658 for (i=0; i<nseed; i++) {
b9de75e1 659 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
660 if (FollowProlongation(t)) {
661 UseClusters(&t);
73042f01 662 continue;
663 }
b9de75e1 664 delete fSeeds->RemoveAt(i);
73042f01 665 }
7f6ddf58 666 //UnloadOuterSectors();
73042f01 667
668 //tracking in inner sectors
b9de75e1 669 LoadInnerSectors();
73042f01 670 Int_t found=0;
671 for (i=0; i<nseed; i++) {
b9de75e1 672 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
73042f01 673 if (!pt) continue;
674 Int_t nc=t.GetNumberOfClusters();
675
b9de75e1 676 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
73042f01 677 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
678 if (alpha < 0. ) alpha += 2.*TMath::Pi();
b9de75e1 679 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
73042f01 680
b9de75e1 681 alpha=ns*fInnerSec->GetAlpha() + fInnerSec->GetAlphaShift() - t.GetAlpha();
73042f01 682
683 if (t.Rotate(alpha)) {
b9de75e1 684 if (FollowProlongation(t)) {
73042f01 685 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
686 t.CookdEdx();
7f6ddf58 687 CookLabel(pt,0.1); //For comparison only
73042f01 688 iotrack=pt;
689 tracktree.Fill();
b9de75e1 690 UseClusters(&t,nc);
691 cerr<<found++<<'\r';
73042f01 692 }
693 }
694 }
b9de75e1 695 delete fSeeds->RemoveAt(i);
73042f01 696 }
b9de75e1 697 UnloadInnerSectors();
7f6ddf58 698 UnloadOuterSectors();
b9de75e1 699
7f6ddf58 700 tracktree.Write();
73042f01 701
b9de75e1 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;
73042f01 723 }
724
b9de75e1 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();
f38c8ae5 732 TTree *bckTree=(TTree*)in->Get("TreeT_ITSb_0");
b9de75e1 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
f38c8ae5 741 TTree *tpcTree=(TTree*)in->Get("TreeT_TPC_0");
b9de75e1 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();
f38c8ae5 775 TTree backTree("TreeT_TPCb_0","Tree with back propagated TPC tracks");
b9de75e1 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;
73042f01 785
b9de75e1 786 Int_t nc=t.GetNumberOfClusters();
787 s.SetLabel(nc-1); //set number of the cluster to start with
73042f01 788
b9de75e1 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();
73042f01 832 savedir->cd();
b9de75e1 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;
be5b9287 839
f38c8ae5 840 delete bckTree; //Thanks to Mariana Bondila
841 delete tpcTree; //Thanks to Mariana Bondila
842
be5b9287 843 return 0;
73042f01 844}
845
b9de75e1 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
9e1a0ddb 855 AliTPCClustersRow *clrow=((AliTPCtracker *) this)->fClustersArray.GetRow(sec,row);
b9de75e1 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();
f03e3423 934 f1PadPitchLength=par->GetOuter1PadPitchLength();
935 f2PadPitchLength=par->GetOuter2PadPitchLength();
b9de75e1 936
937 fN=par->GetNRowUp();
938 fRow=new AliTPCRow[fN];
f03e3423 939 for (Int_t i=0; i<fN; i++){
940 fRow[i].SetX(par->GetPadRowRadiiUp(i));
941 }
b9de75e1 942 }
943}
944
73042f01 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 //-----------------------------------------------------------------------
b9de75e1 951 if (fN==kMaxClusterPerRow) {
73042f01 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 //-----------------------------------------------------------------------
11c7b548 966 if(fN<=0) return 0;
73042f01 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;
be9c5115 983 Int_t nc=GetNumberOfClusters();
73042f01 984
985 Int_t swap;//stupid sorting
986 do {
987 swap=0;
be9c5115 988 for (i=0; i<nc-1; i++) {
d4cf1daa 989 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
990 Float_t tmp=fdEdxSample[i];
991 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
73042f01 992 swap++;
993 }
994 } while (swap);
995
be9c5115 996 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
73042f01 997 Float_t dedx=0;
d4cf1daa 998 for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
73042f01 999 dedx /= (nu-nl+1);
1000 SetdEdx(dedx);
7f6ddf58 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
73042f01 1018}
1019
73042f01 1020
1021