]> 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
Possibility to stop tracking at a given layer (default is 1). Added new method AliITS...
[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$
e24ea474 18Revision 1.30 2003/02/28 16:13:32 hristov
19Typos corrected
20
2afb3f88 21Revision 1.29 2003/02/28 15:18:16 hristov
22Corrections suggested by J.Chudoba
23
fd3694a3 24Revision 1.28 2003/02/27 16:15:52 hristov
25Code for inward refitting (S.Radomski)
26
08dba16b 27Revision 1.27 2003/02/25 16:47:58 hristov
28allow back propagation for more than 1 event (J.Chudoba)
29
438cd838 30Revision 1.26 2003/02/19 08:49:46 hristov
31Track time measurement (S.Radomski)
32
e3f7105e 33Revision 1.25 2003/01/28 16:43:35 hristov
34Additional protection: to be discussed with the Root team (M.Ivanov)
35
646a38dc 36Revision 1.24 2002/11/19 16:13:24 hristov
37stdlib.h included to declare exit() on HP
38
1f9a65c4 39Revision 1.23 2002/11/19 11:50:08 hristov
40Removing CONTAINERS (Yu.Belikov)
41
b9d0a01d 42Revision 1.19 2002/07/19 07:31:40 kowal2
43Improvement in tracking by J. Belikov
44
024a7fe9 45Revision 1.18 2002/05/13 07:33:52 kowal2
46Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
47in the case of defined region of interests.
48
11c7b548 49Revision 1.17 2002/03/18 17:59:13 kowal2
50Chnges in the pad geometry - 3 pad lengths introduced.
51
f03e3423 52Revision 1.16 2001/11/08 16:39:03 hristov
53Additional protection (M.Masera)
54
1c8d9aad 55Revision 1.15 2001/11/08 16:36:33 hristov
56Updated 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)
57
7f6ddf58 58Revision 1.14 2001/10/21 19:04:55 hristov
59Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
60
f38c8ae5 61Revision 1.13 2001/07/23 12:01:30 hristov
62Initialisation added
63
7afb86d1 64Revision 1.12 2001/07/20 14:32:44 kowal2
65Processing of many events possible now
66
afc42102 67Revision 1.11 2001/05/23 08:50:10 hristov
68Weird inline removed
69
bdbd0f7a 70Revision 1.10 2001/05/16 14:57:25 alibrary
71New files for folders and Stack
72
9e1a0ddb 73Revision 1.9 2001/05/11 07:16:56 hristov
74Fix needed on Sun and Alpha
75
3d39dd8b 76Revision 1.8 2001/05/08 15:00:15 hristov
77Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
78
be9c5115 79Revision 1.5 2000/12/20 07:51:59 kowal2
80Changes suggested by Alessandra and Paolo to avoid overlapped
81data fields in encapsulated classes.
82
d4cf1daa 83Revision 1.4 2000/11/02 07:27:16 kowal2
84code corrections
85
98ab53f4 86Revision 1.2 2000/06/30 12:07:50 kowal2
87Updated from the TPC-PreRelease branch
88
73042f01 89Revision 1.1.2.1 2000/06/25 08:53:55 kowal2
90Splitted from AliTPCtracking
91
92*/
93
94//-------------------------------------------------------
95// Implementation of the TPC tracker
96//
97// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
98//-------------------------------------------------------
73042f01 99#include <TObjArray.h>
100#include <TFile.h>
be5b9287 101#include <TTree.h>
19364939 102#include <Riostream.h>
be5b9287 103
104#include "AliTPCtracker.h"
105#include "AliTPCcluster.h"
b9de75e1 106#include "AliTPCParam.h"
bcc04e2a 107#include "AliClusters.h"
73042f01 108
08dba16b 109#include "TVector2.h"
110
1f9a65c4 111#include <stdlib.h>
b9de75e1 112//_____________________________________________________________________________
bcc04e2a 113AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
7f6ddf58 114AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
b9de75e1 115{
116 //---------------------------------------------------------------------
117 // The main TPC tracker constructor
118 //---------------------------------------------------------------------
119 fInnerSec=new AliTPCSector[fkNIS];
120 fOuterSec=new AliTPCSector[fkNOS];
121
122 Int_t i;
123 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
124 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
125
e3f7105e 126 fParam = (AliTPCParam*) par;
b9de75e1 127 fSeeds=0;
128}
129
130//_____________________________________________________________________________
131AliTPCtracker::~AliTPCtracker() {
132 //------------------------------------------------------------------
133 // TPC tracker destructor
134 //------------------------------------------------------------------
135 delete[] fInnerSec;
136 delete[] fOuterSec;
1c8d9aad 137 if (fSeeds) {
138 fSeeds->Delete();
139 delete fSeeds;
140 }
b9de75e1 141}
73042f01 142
143//_____________________________________________________________________________
144Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
145{
146 //
147 // Parametrised error of the cluster reconstruction (pad direction)
148 //
149 // Sigma rphi
150 const Float_t kArphi=0.41818e-2;
151 const Float_t kBrphi=0.17460e-4;
152 const Float_t kCrphi=0.30993e-2;
153 const Float_t kDrphi=0.41061e-3;
154
155 pt=TMath::Abs(pt)*1000.;
156 Double_t x=r/pt;
157 tgl=TMath::Abs(tgl);
158 Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
159 if (s<0.4e-3) s=0.4e-3;
160 s*=1.3; //Iouri Belikov
161
162 return s;
163}
164
165//_____________________________________________________________________________
166Double_t SigmaZ2(Double_t r, Double_t tgl)
167{
168 //
169 // Parametrised error of the cluster reconstruction (drift direction)
170 //
171 // Sigma z
172 const Float_t kAz=0.39614e-2;
173 const Float_t kBz=0.22443e-4;
174 const Float_t kCz=0.51504e-1;
175
176
177 tgl=TMath::Abs(tgl);
178 Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
179 if (s<0.4e-3) s=0.4e-3;
180 s*=1.3; //Iouri Belikov
181
182 return s;
183}
184
185//_____________________________________________________________________________
bdbd0f7a 186Double_t f1(Double_t x1,Double_t y1,
73042f01 187 Double_t x2,Double_t y2,
188 Double_t x3,Double_t y3)
189{
190 //-----------------------------------------------------------------
191 // Initial approximation of the track curvature
192 //-----------------------------------------------------------------
193 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
194 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
195 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
196 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
197 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
198
199 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
200
201 return -xr*yr/sqrt(xr*xr+yr*yr);
202}
203
204
205//_____________________________________________________________________________
bdbd0f7a 206Double_t f2(Double_t x1,Double_t y1,
73042f01 207 Double_t x2,Double_t y2,
208 Double_t x3,Double_t y3)
209{
210 //-----------------------------------------------------------------
211 // Initial approximation of the track curvature times center of curvature
212 //-----------------------------------------------------------------
213 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
214 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
215 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
216 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
217 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
218
219 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
220
221 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
222}
223
224//_____________________________________________________________________________
bdbd0f7a 225Double_t f3(Double_t x1,Double_t y1,
73042f01 226 Double_t x2,Double_t y2,
227 Double_t z1,Double_t z2)
228{
229 //-----------------------------------------------------------------
230 // Initial approximation of the tangent of the track dip angle
231 //-----------------------------------------------------------------
232 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
233}
234
235//_____________________________________________________________________________
e24ea474 236Int_t AliTPCtracker::LoadClusters() {
b9de75e1 237 //-----------------------------------------------------------------
bcc04e2a 238 // This function loads TPC clusters.
b9de75e1 239 //-----------------------------------------------------------------
bcc04e2a 240 if (!gFile->IsOpen()) {
241 cerr<<"AliTPCtracker::LoadClusters : "<<
242 "file with clusters has not been open !\n";
e24ea474 243 return 1;
b9de75e1 244 }
245
bcc04e2a 246 Char_t name[100];
247 sprintf(name,"TreeC_TPC_%d",GetEventNumber());
248 TTree *cTree=(TTree*)gFile->Get(name);
249 if (!cTree) {
250 cerr<<"AliTPCtracker::LoadClusters : "<<
251 "can't get the tree with TPC clusters !\n";
e24ea474 252 return 2;
b9de75e1 253 }
254
bcc04e2a 255 TBranch *branch=cTree->GetBranch("Segment");
256 if (!branch) {
257 cerr<<"AliTPCtracker::LoadClusters : "<<
258 "can't get the segment branch !\n";
e24ea474 259 return 3;
b9de75e1 260 }
e24ea474 261// AliClusters carray, *addr=&carray;
646a38dc 262 AliClusters carray, *addr=&carray;
263 carray.SetClass("AliTPCcluster");
264 carray.SetArray(0);
bcc04e2a 265 branch->SetAddress(&addr);
266
267 Int_t nentr=(Int_t)cTree->GetEntries();
268
269 for (Int_t i=0; i<nentr; i++) {
270 cTree->GetEvent(i);
271
272 Int_t ncl=carray.GetArray()->GetEntriesFast();
273
274 Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows();
275 Int_t id=carray.GetID();
276 if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) {
277 cerr<<"AliTPCtracker::LoadClusters : "<<
278 "wrong index !\n";
279 exit(1);
280 }
281 Int_t outindex = 2*fkNIS*nir;
282 if (id<outindex) {
283 Int_t sec = id/nir;
284 Int_t row = id - sec*nir;
285 sec %= fkNIS;
286 AliTPCRow &padrow=fInnerSec[sec][row];
287 while (ncl--) {
288 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
289 padrow.InsertCluster(c,sec,row);
290 }
291 } else {
292 id -= outindex;
293 Int_t sec = id/nor;
294 Int_t row = id - sec*nor;
295 sec %= fkNOS;
296 AliTPCRow &padrow=fOuterSec[sec][row];
297 while (ncl--) {
298 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
299 padrow.InsertCluster(c,sec+fkNIS,row);
300 }
301 }
b9de75e1 302
bcc04e2a 303 carray.GetArray()->Clear();
304 }
305 delete cTree;
e24ea474 306 return 0;
b9de75e1 307}
308
309//_____________________________________________________________________________
bcc04e2a 310void AliTPCtracker::UnloadClusters() {
b9de75e1 311 //-----------------------------------------------------------------
bcc04e2a 312 // This function unloads TPC clusters.
b9de75e1 313 //-----------------------------------------------------------------
bcc04e2a 314 Int_t i;
315 for (i=0; i<fkNIS; i++) {
316 Int_t nr=fInnerSec->GetNRows();
317 for (Int_t n=0; n<nr; n++) fInnerSec[i][n].ResetClusters();
318 }
319 for (i=0; i<fkNOS; i++) {
320 Int_t nr=fOuterSec->GetNRows();
321 for (Int_t n=0; n<nr; n++) fOuterSec[i][n].ResetClusters();
b9de75e1 322 }
b9de75e1 323}
324
325//_____________________________________________________________________________
326Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
73042f01 327 //-----------------------------------------------------------------
328 // This function tries to find a track prolongation.
329 //-----------------------------------------------------------------
b9de75e1 330 Double_t xt=t.GetX();
331 const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
332 Int_t(0.5*fSectors->GetNRows());
73042f01 333 Int_t tryAgain=kSKIP;
73042f01 334
b9de75e1 335 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
336 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
337 if (alpha < 0. ) alpha += 2.*TMath::Pi();
338 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
339
f03e3423 340 Int_t nrows=fSectors->GetRowNumber(xt)-1;
341 for (Int_t nr=nrows; nr>=rf; nr--) {
b9de75e1 342 Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
73042f01 343 if (!t.PropagateTo(x)) return 0;
344
345 AliTPCcluster *cl=0;
346 UInt_t index=0;
b9de75e1 347 Double_t maxchi2=kMaxCHI2;
348 const AliTPCRow &krow=fSectors[s][nr];
349 Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
350 Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
73042f01 351 Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
352 Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
353
b9de75e1 354 if (road>kMaxROAD) {
73042f01 355 if (t.GetNumberOfClusters()>4)
356 cerr<<t.GetNumberOfClusters()
357 <<"FindProlongation warning: Too broad road !\n";
358 return 0;
359 }
360
361 if (krow) {
362 for (Int_t i=krow.Find(y-road); i<krow; i++) {
be5b9287 363 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
73042f01 364 if (c->GetY() > y+road) break;
365 if (c->IsUsed()) continue;
366 if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
367 Double_t chi2=t.GetPredictedChi2(c);
368 if (chi2 > maxchi2) continue;
369 maxchi2=chi2;
370 cl=c;
371 index=krow.GetIndex(i);
372 }
373 }
374 if (cl) {
b9de75e1 375 Float_t l=fSectors->GetPadPitchWidth();
024a7fe9 376 Float_t corr=1.; if (nr>63) corr=0.67; // new (third) pad response !
377 t.SetSampledEdx(cl->GetQ()/l*corr,t.GetNumberOfClusters());
be9c5115 378 if (!t.Update(cl,maxchi2,index)) {
379 if (!tryAgain--) return 0;
380 } else tryAgain=kSKIP;
73042f01 381 } else {
382 if (tryAgain==0) break;
383 if (y > ymax) {
b9de75e1 384 s = (s+1) % fN;
385 if (!t.Rotate(fSectors->GetAlpha())) return 0;
73042f01 386 } else if (y <-ymax) {
b9de75e1 387 s = (s-1+fN) % fN;
388 if (!t.Rotate(-fSectors->GetAlpha())) return 0;
73042f01 389 }
390 tryAgain--;
391 }
392 }
393
394 return 1;
b9de75e1 395}
08dba16b 396//_____________________________________________________________________________
397
398Int_t AliTPCtracker::FollowRefitInward(AliTPCseed *seed, AliTPCtrack *track) {
399 //
400 // This function propagates seed inward TPC using old clusters
401 // from track.
402 //
403 // Sylwester Radomski, GSI
404 // 26.02.2003
405 //
406
407 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
408 TVector2::Phi_0_2pi(alpha);
409 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
410
411 Int_t idx=-1, sec=-1, row=-1;
412 Int_t nc = seed->GetLabel(); // index to start with
413
414 idx=track->GetClusterIndex(nc);
415 sec=(idx&0xff000000)>>24;
416 row=(idx&0x00ff0000)>>16;
417
418 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
419 else { if (sec < fkNIS) row=-1; }
420
421 Int_t nr=fSectors->GetNRows();
422 for (Int_t i=0; i<nr; i++) {
423
424 Double_t x=fSectors->GetX(i);
425 if (!seed->PropagateTo(x)) return 0;
426
427 // Change sector and rotate seed if necessary
428
429 Double_t ymax=fSectors->GetMaxY(i);
430 Double_t y=seed->GetY();
431
432 if (y > ymax) {
433 s = (s+1) % fN;
434 if (!seed->Rotate(fSectors->GetAlpha())) return 0;
435 } else if (y <-ymax) {
436 s = (s-1+fN) % fN;
437 if (!seed->Rotate(-fSectors->GetAlpha())) return 0;
438 }
439
440 // try to find a cluster
441
442 AliTPCcluster *cl=0;
443 Int_t index = 0;
444 Double_t maxchi2 = kMaxCHI2;
445
446 //cout << i << " " << row << " " << nc << endl;
447
448 if (row==i) {
449
450 // accept already found cluster
451 AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
452 Double_t chi2 = seed->GetPredictedChi2(c);
453
454 index=idx;
455 cl=c;
456 maxchi2=chi2;
457
458 if (++nc < track->GetNumberOfClusters() ) {
459 idx=track->GetClusterIndex(nc);
460 sec=(idx&0xff000000)>>24;
461 row=(idx&0x00ff0000)>>16;
462 }
463
464 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
465 else { if (sec < fkNIS) row=-1; }
466
467 }
468
469 if (cl) {
470
471 //cout << "Assigned" << endl;
472 Float_t l=fSectors->GetPadPitchWidth();
473 Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
474 seed->SetSampledEdx(cl->GetQ()/l*corr,seed->GetNumberOfClusters());
475 seed->Update(cl,maxchi2,index);
476 }
477 }
478
479 seed->SetLabel(nc);
480 return 1;
481}
b9de75e1 482
483//_____________________________________________________________________________
484Int_t AliTPCtracker::FollowBackProlongation
485(AliTPCseed& seed, const AliTPCtrack &track) {
486 //-----------------------------------------------------------------
487 // This function propagates tracks back through the TPC
488 //-----------------------------------------------------------------
489 Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
490 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
491 if (alpha < 0. ) alpha += 2.*TMath::Pi();
492 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
493
494 Int_t idx=-1, sec=-1, row=-1;
495 Int_t nc=seed.GetLabel(); //index of the cluster to start with
496 if (nc--) {
497 idx=track.GetClusterIndex(nc);
498 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
499 }
bcc04e2a 500 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
501 else { if (sec < fkNIS) row=-1; }
b9de75e1 502
503 Int_t nr=fSectors->GetNRows();
504 for (Int_t i=0; i<nr; i++) {
505 Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
506
507 if (!seed.PropagateTo(x)) return 0;
508
509 Double_t y=seed.GetY();
510 if (y > ymax) {
511 s = (s+1) % fN;
512 if (!seed.Rotate(fSectors->GetAlpha())) return 0;
513 } else if (y <-ymax) {
514 s = (s-1+fN) % fN;
515 if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
516 }
517
518 AliTPCcluster *cl=0;
519 Int_t index=0;
520 Double_t maxchi2=kMaxCHI2;
521 Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
522 Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
523 Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
524 Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
525 if (road>kMaxROAD) {
526 cerr<<seed.GetNumberOfClusters()
527 <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n";
528 return 0;
529 }
530
531
532 Int_t accepted=seed.GetNumberOfClusters();
533 if (row==i) {
534 //try to accept already found cluster
535 AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
536 Double_t chi2;
537 if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
538 index=idx; cl=c; maxchi2=chi2;
539 } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
540
541 if (nc--) {
542 idx=track.GetClusterIndex(nc);
543 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
544 }
bcc04e2a 545 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
546 else { if (sec < fkNIS) row=-1; }
73042f01 547
b9de75e1 548 }
549 if (!cl) {
550 //try to fill the gap
551 const AliTPCRow &krow=fSectors[s][i];
552 if (accepted>27)
553 if (krow) {
554 for (Int_t i=krow.Find(y-road); i<krow; i++) {
555 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
556 if (c->GetY() > y+road) break;
557 if (c->IsUsed()) continue;
558 if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
559 Double_t chi2=seed.GetPredictedChi2(c);
560 if (chi2 > maxchi2) continue;
561 maxchi2=chi2;
562 cl=c;
563 index=krow.GetIndex(i);
564 }
565 }
566 }
567
568 if (cl) {
569 Float_t l=fSectors->GetPadPitchWidth();
024a7fe9 570 Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
571 seed.SetSampledEdx(cl->GetQ()/l*corr,seed.GetNumberOfClusters());
b9de75e1 572 seed.Update(cl,maxchi2,index);
573 }
574
575 }
576
577 seed.SetLabel(nc);
578
579 return 1;
73042f01 580}
581
582//_____________________________________________________________________________
b9de75e1 583void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
73042f01 584 //-----------------------------------------------------------------
585 // This function creates track seeds.
586 //-----------------------------------------------------------------
587 Double_t x[5], c[15];
588
bcc04e2a 589 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
73042f01 590 Double_t cs=cos(alpha), sn=sin(alpha);
591
bcc04e2a 592 Double_t x1 =fSectors->GetX(i1);
593 Double_t xx2=fSectors->GetX(i2);
73042f01 594
bcc04e2a 595 for (Int_t ns=0; ns<fN; ns++) {
596 Int_t nl=fSectors[(ns-1+fN)%fN][i2];
597 Int_t nm=fSectors[ns][i2];
598 Int_t nu=fSectors[(ns+1)%fN][i2];
599 const AliTPCRow& kr1=fSectors[ns][i1];
73042f01 600 for (Int_t is=0; is < kr1; is++) {
601 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
602 for (Int_t js=0; js < nl+nm+nu; js++) {
603 const AliTPCcluster *kcl;
604 Double_t x2, y2, z2;
7f6ddf58 605 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
73042f01 606
607 if (js<nl) {
bcc04e2a 608 const AliTPCRow& kr2=fSectors[(ns-1+fN)%fN][i2];
73042f01 609 kcl=kr2[js];
610 y2=kcl->GetY(); z2=kcl->GetZ();
611 x2= xx2*cs+y2*sn;
612 y2=-xx2*sn+y2*cs;
613 } else
614 if (js<nl+nm) {
bcc04e2a 615 const AliTPCRow& kr2=fSectors[ns][i2];
73042f01 616 kcl=kr2[js-nl];
617 x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
618 } else {
bcc04e2a 619 const AliTPCRow& kr2=fSectors[(ns+1)%fN][i2];
73042f01 620 kcl=kr2[js-nl-nm];
621 y2=kcl->GetY(); z2=kcl->GetZ();
622 x2=xx2*cs-y2*sn;
623 y2=xx2*sn+y2*cs;
624 }
625
7f6ddf58 626 Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
73042f01 627 if (TMath::Abs(zz-z2)>5.) continue;
628
629 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
630 if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
631
632 x[0]=y1;
633 x[1]=z1;
b9de75e1 634 x[4]=f1(x1,y1,x2,y2,x3,y3);
635 if (TMath::Abs(x[4]) >= 0.0066) continue;
be5b9287 636 x[2]=f2(x1,y1,x2,y2,x3,y3);
b9de75e1 637 //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
638 x[3]=f3(x1,y1,x2,y2,z1,z2);
639 if (TMath::Abs(x[3]) > 1.2) continue;
be5b9287 640 Double_t a=asin(x[2]);
b9de75e1 641 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
7f6ddf58 642 if (TMath::Abs(zv-z3)>10.) continue;
73042f01 643
644 Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
645 Double_t sy2=kcl->GetSigmaY2(), sz2=kcl->GetSigmaZ2();
024a7fe9 646 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
647 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
73042f01 648
b9de75e1 649 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
650 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
651 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
be5b9287 652 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
653 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
b9de75e1 654 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
655 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
656 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
657 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
658 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
73042f01 659
660 c[0]=sy1;
661 c[1]=0.; c[2]=sz1;
b9de75e1 662 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
663 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
664 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
665 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
666 c[13]=f30*sy1*f40+f32*sy2*f42;
667 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
73042f01 668
669 UInt_t index=kr1.GetIndex(is);
670 AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
bcc04e2a 671 Float_t l=fSectors->GetPadPitchWidth();
73042f01 672 track->SetSampledEdx(kr1[is]->GetQ()/l,0);
673
b9de75e1 674 Int_t rc=FollowProlongation(*track, i2);
be9c5115 675 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
b9de75e1 676 else fSeeds->AddLast(track);
73042f01 677 }
678 }
679 }
680}
681
682//_____________________________________________________________________________
b9de75e1 683Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
73042f01 684 //-----------------------------------------------------------------
b9de75e1 685 // This function reades track seeds.
73042f01 686 //-----------------------------------------------------------------
687 TDirectory *savedir=gDirectory;
688
b9de75e1 689 TFile *in=(TFile*)inp;
690 if (!in->IsOpen()) {
691 cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
be5b9287 692 return 1;
73042f01 693 }
694
b9de75e1 695 in->cd();
696 TTree *seedTree=(TTree*)in->Get("Seeds");
697 if (!seedTree) {
698 cerr<<"AliTPCtracker::ReadSeeds(): ";
699 cerr<<"can't get a tree with track seeds !\n";
700 return 2;
701 }
702 AliTPCtrack *seed=new AliTPCtrack;
703 seedTree->SetBranchAddress("tracks",&seed);
704
705 if (fSeeds==0) fSeeds=new TObjArray(15000);
73042f01 706
b9de75e1 707 Int_t n=(Int_t)seedTree->GetEntries();
708 for (Int_t i=0; i<n; i++) {
709 seedTree->GetEvent(i);
710 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
711 }
712
713 delete seed;
f38c8ae5 714
715 delete seedTree; //Thanks to Mariana Bondila
716
b9de75e1 717 savedir->cd();
73042f01 718
b9de75e1 719 return 0;
720}
73042f01 721
b9de75e1 722//_____________________________________________________________________________
723Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
724 //-----------------------------------------------------------------
725 // This is a track finder.
726 //-----------------------------------------------------------------
727 TDirectory *savedir=gDirectory;
73042f01 728
b9de75e1 729 if (inp) {
730 TFile *in=(TFile*)inp;
731 if (!in->IsOpen()) {
732 cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
733 return 1;
734 }
735 }
73042f01 736
b9de75e1 737 if (!out->IsOpen()) {
738 cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
739 return 2;
73042f01 740 }
741
bcc04e2a 742 LoadClusters();
743
b9de75e1 744 out->cd();
7f6ddf58 745
746 char tname[100];
bcc04e2a 747 sprintf(tname,"TreeT_TPC_%d",GetEventNumber());
7f6ddf58 748 TTree tracktree(tname,"Tree with TPC tracks");
b9de75e1 749 AliTPCtrack *iotrack=0;
08dba16b 750 tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,3);
b9de75e1 751
73042f01 752 //find track seeds
b9de75e1 753 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
73042f01 754 Int_t nrows=nlow+nup;
b9de75e1 755 if (fSeeds==0) {
756 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
bcc04e2a 757 fSectors=fOuterSec; fN=fkNOS;
758 fSeeds=new TObjArray(15000);
b9de75e1 759 MakeSeeds(nup-1, nup-1-gap);
760 MakeSeeds(nup-1-shift, nup-1-shift-gap);
bcc04e2a 761 }
b9de75e1 762 fSeeds->Sort();
73042f01 763
bcc04e2a 764 Int_t found=0;
b9de75e1 765 Int_t nseed=fSeeds->GetEntriesFast();
bcc04e2a 766 for (Int_t i=0; i<nseed; i++) {
767 //tracking in the outer sectors
768 fSectors=fOuterSec; fN=fkNOS;
769
b9de75e1 770 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
bcc04e2a 771 if (!FollowProlongation(t)) {
772 delete fSeeds->RemoveAt(i);
73042f01 773 continue;
774 }
73042f01 775
bcc04e2a 776 //tracking in the inner sectors
777 fSectors=fInnerSec; fN=fkNIS;
73042f01 778
b9de75e1 779 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
73042f01 780 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
781 if (alpha < 0. ) alpha += 2.*TMath::Pi();
b9de75e1 782 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
73042f01 783
bcc04e2a 784 alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
73042f01 785
786 if (t.Rotate(alpha)) {
e3f7105e 787 if (FollowProlongation(t)) {
788 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
789 t.CookdEdx();
790 CookLabel(pt,0.1); //For comparison only
fd3694a3 791 pt->PropagateTo(fParam->GetInnerRadiusLow());
e3f7105e 792 iotrack=pt;
793 tracktree.Fill();
794 UseClusters(&t);
e24ea474 795 found++;
796// cerr<<found<<'\r';
e3f7105e 797 }
798 }
73042f01 799 }
b9de75e1 800 delete fSeeds->RemoveAt(i);
bcc04e2a 801 }
e3f7105e 802
7f6ddf58 803 tracktree.Write();
73042f01 804
b9de75e1 805 cerr<<"Number of found tracks : "<<found<<endl;
806
807 savedir->cd();
808
bcc04e2a 809 UnloadClusters();
810 fSeeds->Clear(); delete fSeeds; fSeeds=0;
811
b9de75e1 812 return 0;
813}
08dba16b 814//_____________________________________________________________________________
815
816Int_t AliTPCtracker::RefitInward(TFile *in, TFile *out) {
817 //
818 // The function propagates tracks throught TPC inward
819 // using already associated clusters.
820 //
821 // in - file with back propagated tracks
822 // outs - file with inward propagation
823 //
824 // Sylwester Radomski, GSI
825 // 26.02.2003
826 //
827
828
829 if (!in->IsOpen()) {
830 cout << "Input File not open !\n" << endl;
831 return 1;
832 }
833
834 if (!out->IsOpen()) {
835 cout << "Output File not open !\n" << endl;
836 return 2;
837 }
838
839 TDirectory *savedir = gDirectory;
840 LoadClusters();
841
842 // Connect to input tree
843 in->cd();
844 TTree *inputTree = (TTree*)in->Get("seedsTPCtoTRD_0");
845 TBranch *inBranch = inputTree->GetBranch("tracks");
846 AliTPCtrack *inTrack = 0;
847 inBranch->SetAddress(&inTrack);
848
849 Int_t nTracks = (Int_t)inputTree->GetEntries();
850
851 // Create output tree
852 out->cd();
853 AliTPCtrack *outTrack = new AliTPCtrack();
854 TTree *outputTree = new TTree("tracksTPC_0","Refited TPC tracks");
855 outputTree->Branch("tracks", "AliTPCtrack", &outTrack, 32000, 3);
856
857 Int_t nRefited = 0;
858
859 for(Int_t t=0; t < nTracks; t++) {
860
861 cout << t << "\r";
862
863 inputTree->GetEvent(t);
864 AliTPCseed *seed = new AliTPCseed(*inTrack, inTrack->GetAlpha());
865
866 seed->ResetCovariance();
867
868 seed->SetLabel(0);
869 fSectors = fInnerSec;
870 Int_t res = FollowRefitInward(seed, inTrack);
871 UseClusters(seed);
872 Int_t nc = seed->GetNumberOfClusters();
873
874 fSectors = fOuterSec;
875 res = FollowRefitInward(seed, inTrack);
876 UseClusters(seed, nc);
877
878 if (res) {
879 seed->PropagateTo(fParam->GetInnerRadiusLow());
880 outTrack = (AliTPCtrack*)seed;
881 outTrack->SetLabel(inTrack->GetLabel());
882 outputTree->Fill();
883 nRefited++;
884 }
885
886 delete seed;
887 }
888
889 cout << "Refitted " << nRefited << " tracks." << endl;
890
891 outputTree->Write();
892
893 if (inputTree) delete inputTree;
894 if (outputTree) delete outputTree;
895
896 savedir->cd();
897 UnloadClusters();
898
899 return 0;
900}
b9de75e1 901
902//_____________________________________________________________________________
903Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
e24ea474 904 //-----------------------------------------------------------------
905 // This function propagates tracks back through the TPC.
906 //-----------------------------------------------------------------
907 return PropagateBack(inp, NULL, out);
908}
909
910//_____________________________________________________________________________
911Int_t AliTPCtracker::PropagateBack(const TFile *inp, const TFile *inp2, TFile *out) {
b9de75e1 912 //-----------------------------------------------------------------
913 // This function propagates tracks back through the TPC.
914 //-----------------------------------------------------------------
915 fSeeds=new TObjArray(15000);
916
917 TFile *in=(TFile*)inp;
e24ea474 918 TFile *in2=(TFile*)inp2;
b9de75e1 919 TDirectory *savedir=gDirectory;
920
921 if (!in->IsOpen()) {
922 cerr<<"AliTPCtracker::PropagateBack(): ";
e24ea474 923 cerr<<"file with TPC (or back propagated ITS) tracks is not open !\n";
924 return 1;
73042f01 925 }
926
b9de75e1 927 if (!out->IsOpen()) {
928 cerr<<"AliTPCtracker::PropagateBack(): ";
929 cerr<<"file for back propagated TPC tracks is not open !\n";
930 return 2;
931 }
932
bcc04e2a 933 LoadClusters();
934
b9de75e1 935 in->cd();
438cd838 936 char tName[100];
937 sprintf(tName,"TreeT_ITSb_%d",GetEventNumber());
938 TTree *bckTree=(TTree*)in->Get(tName);
e24ea474 939 if (!bckTree && inp2) bckTree=(TTree*)in2->Get(tName);
b9de75e1 940 if (!bckTree) {
941 cerr<<"AliTPCtracker::PropagateBack() ";
942 cerr<<"can't get a tree with back propagated ITS tracks !\n";
fd3694a3 943 //return 3;
b9de75e1 944 }
945 AliTPCtrack *bckTrack=new AliTPCtrack;
fd3694a3 946 if (bckTree) bckTree->SetBranchAddress("tracks",&bckTrack);
b9de75e1 947
438cd838 948 sprintf(tName,"TreeT_TPC_%d",GetEventNumber());
949 TTree *tpcTree=(TTree*)in->Get(tName);
b9de75e1 950 if (!tpcTree) {
951 cerr<<"AliTPCtracker::PropagateBack() ";
952 cerr<<"can't get a tree with TPC tracks !\n";
953 return 4;
954 }
955 AliTPCtrack *tpcTrack=new AliTPCtrack;
956 tpcTree->SetBranchAddress("tracks",&tpcTrack);
957
e3f7105e 958 //*** Prepare an array of tracks to be back propagated
b9de75e1 959 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
960 Int_t nrows=nlow+nup;
961
e3f7105e 962 //
963 // Match ITS tracks with old TPC tracks
964 // the tracks do not have to be sorted [SR, GSI, 18.02.2003]
965 //
966 // the algorithm is linear and uses LUT for sorting
967 // cost of the algorithm: 2 * number of tracks
968 //
969
970 TObjArray tracks(15000);
971 Int_t i=0;
972 Int_t tpcN= (Int_t)tpcTree->GetEntries();
973 Int_t bckN= (bckTree)? (Int_t)bckTree->GetEntries() : 0;
974
975 // look up table
fd3694a3 976 const Int_t nLab = 200000; // limit on number of primaries (arbitrary)
e3f7105e 977 Int_t lut[nLab];
978 for(Int_t i=0; i<nLab; i++) lut[i] = -1;
979
980 if (bckTree) {
981 for(Int_t i=0; i<bckN; i++) {
982 bckTree->GetEvent(i);
983 Int_t lab = TMath::Abs(bckTrack->GetLabel());
984 if (lab < nLab) lut[lab] = i;
fd3694a3 985 else {
986 cerr << "AliTPCtracker: The size of the LUT is too small\n";
2afb3f88 987 }
e3f7105e 988 }
989 }
990
991 for (Int_t i=0; i<tpcN; i++) {
992 tpcTree->GetEvent(i);
993 Double_t alpha=tpcTrack->GetAlpha();
994
995 if (!bckTree) {
996
997 // No ITS - use TPC track only
998
2afb3f88 999 AliTPCseed * seed = new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha());
fd3694a3 1000 seed->ResetCovariance();
1001 fSeeds->AddLast(seed);
e3f7105e 1002 tracks.AddLast(new AliTPCtrack(*tpcTrack));
1003
1004 } else {
1005
1006 // with ITS
1007 // discard not prolongated tracks (to be discussed)
1008
1009 Int_t lab = TMath::Abs(tpcTrack->GetLabel());
1010 if (lab > nLab || lut[lab] < 0) continue;
1011
1012 bckTree->GetEvent(lut[lab]);
1013 bckTrack->Rotate(alpha-bckTrack->GetAlpha());
1014
1015 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
1016 tracks.AddLast(new AliTPCtrack(*tpcTrack));
1017 }
1018 }
1019
1020 // end of matching
1021
1022 /*
b9de75e1 1023 TObjArray tracks(15000);
e3f7105e 1024 Int_t i=0, j=0;
b9de75e1 1025 Int_t tpcN=(Int_t)tpcTree->GetEntries();
1026 Int_t bckN=(Int_t)bckTree->GetEntries();
1027 if (j<bckN) bckTree->GetEvent(j++);
1028 for (i=0; i<tpcN; i++) {
1029 tpcTree->GetEvent(i);
1030 Double_t alpha=tpcTrack->GetAlpha();
1031 if (j<bckN &&
1032 TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
1033 if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
1034 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
1035 bckTree->GetEvent(j++);
1036 } else {
1037 tpcTrack->ResetCovariance();
1038 fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
1039 }
1040 tracks.AddLast(new AliTPCtrack(*tpcTrack));
1041 }
e3f7105e 1042 */
1043
b9de75e1 1044 out->cd();
e3f7105e 1045
1046 // tree name seedsTPCtoTRD as expected by TRD and as
1047 // discussed and decided in Strasbourg (May 2002)
1048 // [SR, GSI, 18.02.2003]
1049
438cd838 1050 sprintf(tName,"seedsTPCtoTRD_%d",GetEventNumber());
1051 TTree backTree(tName,"Tree with back propagated TPC tracks");
08dba16b 1052
b9de75e1 1053 AliTPCtrack *otrack=0;
1054 backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
1055
08dba16b 1056 //*** Back propagation through inner sectors
bcc04e2a 1057 fSectors=fInnerSec; fN=fkNIS;
b9de75e1 1058 Int_t nseed=fSeeds->GetEntriesFast();
1059 for (i=0; i<nseed; i++) {
1060 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
1061 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
73042f01 1062
b9de75e1 1063 Int_t nc=t.GetNumberOfClusters();
1064 s.SetLabel(nc-1); //set number of the cluster to start with
73042f01 1065
b9de75e1 1066 if (FollowBackProlongation(s,t)) {
1067 UseClusters(&s);
1068 continue;
1069 }
1070 delete fSeeds->RemoveAt(i);
1071 }
b9de75e1 1072
1073//*** Back propagation through outer sectors
bcc04e2a 1074 fSectors=fOuterSec; fN=fkNOS;
b9de75e1 1075 Int_t found=0;
1076 for (i=0; i<nseed; i++) {
1077 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
1078 if (!ps) continue;
1079 Int_t nc=s.GetNumberOfClusters();
1080 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
1081
1082 Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1083 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1084 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1085 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1086
1087 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1088 alpha-=s.GetAlpha();
1089
1090 if (s.Rotate(alpha)) {
e3f7105e 1091 if (FollowBackProlongation(s,t)) {
1092 if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
1093 s.CookdEdx();
1094 s.SetLabel(t.GetLabel());
1095 UseClusters(&s,nc);
1096
1097 // Propagate to outer reference plane for comparison reasons
1098 // reason for keeping fParam object [SR, GSI, 18.02.2003]
1099
1100 ps->PropagateTo(fParam->GetOuterRadiusUp());
1101 otrack=ps;
1102 backTree.Fill();
e24ea474 1103 found++;
1104// cerr<<found<<'\r';
e3f7105e 1105 continue;
1106 }
1107 }
b9de75e1 1108 }
1109 delete fSeeds->RemoveAt(i);
1110 }
b9de75e1 1111
1112 backTree.Write();
73042f01 1113 savedir->cd();
b9de75e1 1114 cerr<<"Number of seeds: "<<nseed<<endl;
1115 cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
1116 cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
1117
1118 delete bckTrack;
1119 delete tpcTrack;
be5b9287 1120
e24ea474 1121 if (bckTree) delete bckTree; //Thanks to Mariana Bondila
f38c8ae5 1122 delete tpcTree; //Thanks to Mariana Bondila
1123
bcc04e2a 1124 UnloadClusters();
1125
be5b9287 1126 return 0;
73042f01 1127}
1128
b9de75e1 1129//_________________________________________________________________________
1130AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
1131 //--------------------------------------------------------------------
1132 // Return pointer to a given cluster
1133 //--------------------------------------------------------------------
1134 Int_t sec=(index&0xff000000)>>24;
1135 Int_t row=(index&0x00ff0000)>>16;
1136 Int_t ncl=(index&0x0000ffff)>>00;
1137
bcc04e2a 1138 const AliTPCcluster *cl=0;
1139
1140 if (sec<fkNIS) {
1141 cl=fInnerSec[sec][row].GetUnsortedCluster(ncl);
1142 } else {
1143 sec-=fkNIS;
1144 cl=fOuterSec[sec][row].GetUnsortedCluster(ncl);
1145 }
1146
1147 return (AliCluster*)cl;
b9de75e1 1148}
1149
1150//__________________________________________________________________________
1151void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1152 //--------------------------------------------------------------------
1153 //This function "cooks" a track label. If label<0, this track is fake.
1154 //--------------------------------------------------------------------
1155 Int_t noc=t->GetNumberOfClusters();
1156 Int_t *lb=new Int_t[noc];
1157 Int_t *mx=new Int_t[noc];
1158 AliCluster **clusters=new AliCluster*[noc];
1159
1160 Int_t i;
1161 for (i=0; i<noc; i++) {
1162 lb[i]=mx[i]=0;
1163 Int_t index=t->GetClusterIndex(i);
1164 clusters[i]=GetCluster(index);
1165 }
1166
1167 Int_t lab=123456789;
1168 for (i=0; i<noc; i++) {
1169 AliCluster *c=clusters[i];
1170 lab=TMath::Abs(c->GetLabel(0));
1171 Int_t j;
1172 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
1173 lb[j]=lab;
1174 (mx[j])++;
1175 }
1176
1177 Int_t max=0;
1178 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
1179
1180 for (i=0; i<noc; i++) {
1181 AliCluster *c=clusters[i];
1182 if (TMath::Abs(c->GetLabel(1)) == lab ||
1183 TMath::Abs(c->GetLabel(2)) == lab ) max++;
1184 }
1185
1186 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
1187
1188 else {
1189 Int_t tail=Int_t(0.10*noc);
1190 max=0;
1191 for (i=1; i<=tail; i++) {
1192 AliCluster *c=clusters[noc-i];
1193 if (lab == TMath::Abs(c->GetLabel(0)) ||
1194 lab == TMath::Abs(c->GetLabel(1)) ||
1195 lab == TMath::Abs(c->GetLabel(2))) max++;
1196 }
1197 if (max < Int_t(0.5*tail)) lab=-lab;
1198 }
1199
1200 t->SetLabel(lab);
1201
1202 delete[] lb;
1203 delete[] mx;
1204 delete[] clusters;
1205}
1206
1207//_________________________________________________________________________
1208void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
1209 //-----------------------------------------------------------------------
1210 // Setup inner sector
1211 //-----------------------------------------------------------------------
1212 if (f==0) {
1213 fAlpha=par->GetInnerAngle();
1214 fAlphaShift=par->GetInnerAngleShift();
1215 fPadPitchWidth=par->GetInnerPadPitchWidth();
024a7fe9 1216 f1PadPitchLength=par->GetInnerPadPitchLength();
1217 f2PadPitchLength=f1PadPitchLength;
b9de75e1 1218
1219 fN=par->GetNRowLow();
1220 fRow=new AliTPCRow[fN];
1221 for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
1222 } else {
1223 fAlpha=par->GetOuterAngle();
1224 fAlphaShift=par->GetOuterAngleShift();
1225 fPadPitchWidth=par->GetOuterPadPitchWidth();
f03e3423 1226 f1PadPitchLength=par->GetOuter1PadPitchLength();
1227 f2PadPitchLength=par->GetOuter2PadPitchLength();
b9de75e1 1228
1229 fN=par->GetNRowUp();
1230 fRow=new AliTPCRow[fN];
f03e3423 1231 for (Int_t i=0; i<fN; i++){
1232 fRow[i].SetX(par->GetPadRowRadiiUp(i));
1233 }
b9de75e1 1234 }
1235}
1236
73042f01 1237//_________________________________________________________________________
bcc04e2a 1238void AliTPCtracker::
1239AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) {
73042f01 1240 //-----------------------------------------------------------------------
1241 // Insert a cluster into this pad row in accordence with its y-coordinate
1242 //-----------------------------------------------------------------------
b9de75e1 1243 if (fN==kMaxClusterPerRow) {
73042f01 1244 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
1245 }
bcc04e2a 1246
1247 Int_t index=(((sec<<8)+row)<<16)+fN;
1248
1249 if (fN==0) {
1250 fIndex[0]=index;
1251 fClusterArray[0]=*c; fClusters[fN++]=fClusterArray;
1252 return;
1253 }
1254
1255 if (fN==fSize) {
1256 Int_t size=fSize*2;
1257 AliTPCcluster *buff=new AliTPCcluster[size];
1258 memcpy(buff,fClusterArray,fSize*sizeof(AliTPCcluster));
1259 for (Int_t i=0; i<fN; i++)
1260 fClusters[i]=buff+(fClusters[i]-fClusterArray);
1261 delete[] fClusterArray;
1262 fClusterArray=buff;
1263 fSize=size;
1264 }
1265
73042f01 1266 Int_t i=Find(c->GetY());
1267 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
1268 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
bcc04e2a 1269 fIndex[i]=index;
1270 fClusters[i]=fClusterArray+fN; fClusterArray[fN++]=*c;
73042f01 1271}
1272
1273//___________________________________________________________________
1274Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
1275 //-----------------------------------------------------------------------
1276 // Return the index of the nearest cluster
1277 //-----------------------------------------------------------------------
11c7b548 1278 if(fN<=0) return 0;
73042f01 1279 if (y <= fClusters[0]->GetY()) return 0;
1280 if (y > fClusters[fN-1]->GetY()) return fN;
1281 Int_t b=0, e=fN-1, m=(b+e)/2;
1282 for (; b<e; m=(b+e)/2) {
1283 if (y > fClusters[m]->GetY()) b=m+1;
1284 else e=m;
1285 }
1286 return m;
1287}
1288
1289//_____________________________________________________________________________
1290void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1291 //-----------------------------------------------------------------
1292 // This funtion calculates dE/dX within the "low" and "up" cuts.
1293 //-----------------------------------------------------------------
1294 Int_t i;
be9c5115 1295 Int_t nc=GetNumberOfClusters();
bcc04e2a 1296
1297 Int_t swap;//stupid sorting
1298 do {
1299 swap=0;
1300 for (i=0; i<nc-1; i++) {
1301 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
1302 Float_t tmp=fdEdxSample[i];
1303 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
1304 swap++;
1305 }
1306 } while (swap);
73042f01 1307
be9c5115 1308 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
73042f01 1309 Float_t dedx=0;
bcc04e2a 1310 for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
73042f01 1311 dedx /= (nu-nl+1);
1312 SetdEdx(dedx);
7f6ddf58 1313
1314 //Very rough PID
1315 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1316
024a7fe9 1317 Double_t log1=TMath::Log(p+0.45), log2=TMath::Log(p+0.12);
7f6ddf58 1318 if (p<0.6) {
024a7fe9 1319 if (dedx < 34 + 30/(p+0.45)/(p+0.45) + 24*log1) {SetMass(0.13957); return;}
1320 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.49368); return;}
7f6ddf58 1321 SetMass(0.93827); return;
1322 }
1323
1324 if (p<1.2) {
024a7fe9 1325 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.13957); return;}
7f6ddf58 1326 SetMass(0.93827); return;
1327 }
1328
1329 SetMass(0.13957); return;
1330
73042f01 1331}
1332
73042f01 1333
1334