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