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