]> 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
stdlib.h included to declare exit() on HP
[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$
1f9a65c4 18Revision 1.23 2002/11/19 11:50:08 hristov
19Removing CONTAINERS (Yu.Belikov)
20
b9d0a01d 21Revision 1.19 2002/07/19 07:31:40 kowal2
22Improvement in tracking by J. Belikov
23
024a7fe9 24Revision 1.18 2002/05/13 07:33:52 kowal2
25Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
26in the case of defined region of interests.
27
11c7b548 28Revision 1.17 2002/03/18 17:59:13 kowal2
29Chnges in the pad geometry - 3 pad lengths introduced.
30
f03e3423 31Revision 1.16 2001/11/08 16:39:03 hristov
32Additional protection (M.Masera)
33
1c8d9aad 34Revision 1.15 2001/11/08 16:36:33 hristov
35Updated 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)
36
7f6ddf58 37Revision 1.14 2001/10/21 19:04:55 hristov
38Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
39
f38c8ae5 40Revision 1.13 2001/07/23 12:01:30 hristov
41Initialisation added
42
7afb86d1 43Revision 1.12 2001/07/20 14:32:44 kowal2
44Processing of many events possible now
45
afc42102 46Revision 1.11 2001/05/23 08:50:10 hristov
47Weird inline removed
48
bdbd0f7a 49Revision 1.10 2001/05/16 14:57:25 alibrary
50New files for folders and Stack
51
9e1a0ddb 52Revision 1.9 2001/05/11 07:16:56 hristov
53Fix needed on Sun and Alpha
54
3d39dd8b 55Revision 1.8 2001/05/08 15:00:15 hristov
56Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
57
be9c5115 58Revision 1.5 2000/12/20 07:51:59 kowal2
59Changes suggested by Alessandra and Paolo to avoid overlapped
60data fields in encapsulated classes.
61
d4cf1daa 62Revision 1.4 2000/11/02 07:27:16 kowal2
63code corrections
64
98ab53f4 65Revision 1.2 2000/06/30 12:07:50 kowal2
66Updated from the TPC-PreRelease branch
67
73042f01 68Revision 1.1.2.1 2000/06/25 08:53:55 kowal2
69Splitted from AliTPCtracking
70
71*/
72
73//-------------------------------------------------------
74// Implementation of the TPC tracker
75//
76// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
77//-------------------------------------------------------
73042f01 78#include <TObjArray.h>
79#include <TFile.h>
be5b9287 80#include <TTree.h>
19364939 81#include <Riostream.h>
be5b9287 82
83#include "AliTPCtracker.h"
84#include "AliTPCcluster.h"
b9de75e1 85#include "AliTPCParam.h"
bcc04e2a 86#include "AliClusters.h"
73042f01 87
1f9a65c4 88#include <stdlib.h>
b9de75e1 89//_____________________________________________________________________________
bcc04e2a 90AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
7f6ddf58 91AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
b9de75e1 92{
93 //---------------------------------------------------------------------
94 // The main TPC tracker constructor
95 //---------------------------------------------------------------------
96 fInnerSec=new AliTPCSector[fkNIS];
97 fOuterSec=new AliTPCSector[fkNOS];
98
99 Int_t i;
100 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
101 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
102
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//_____________________________________________________________________________
bcc04e2a 212void AliTPCtracker::LoadClusters() {
b9de75e1 213 //-----------------------------------------------------------------
bcc04e2a 214 // This function loads TPC clusters.
b9de75e1 215 //-----------------------------------------------------------------
bcc04e2a 216 if (!gFile->IsOpen()) {
217 cerr<<"AliTPCtracker::LoadClusters : "<<
218 "file with clusters has not been open !\n";
219 return;
b9de75e1 220 }
221
bcc04e2a 222 Char_t name[100];
223 sprintf(name,"TreeC_TPC_%d",GetEventNumber());
224 TTree *cTree=(TTree*)gFile->Get(name);
225 if (!cTree) {
226 cerr<<"AliTPCtracker::LoadClusters : "<<
227 "can't get the tree with TPC clusters !\n";
228 return;
b9de75e1 229 }
230
bcc04e2a 231 TBranch *branch=cTree->GetBranch("Segment");
232 if (!branch) {
233 cerr<<"AliTPCtracker::LoadClusters : "<<
234 "can't get the segment branch !\n";
235 return;
b9de75e1 236 }
bcc04e2a 237 AliClusters carray, *addr=&carray; carray.SetClass("AliTPCcluster");
238 branch->SetAddress(&addr);
239
240 Int_t nentr=(Int_t)cTree->GetEntries();
241
242 for (Int_t i=0; i<nentr; i++) {
243 cTree->GetEvent(i);
244
245 Int_t ncl=carray.GetArray()->GetEntriesFast();
246
247 Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows();
248 Int_t id=carray.GetID();
249 if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) {
250 cerr<<"AliTPCtracker::LoadClusters : "<<
251 "wrong index !\n";
252 exit(1);
253 }
254 Int_t outindex = 2*fkNIS*nir;
255 if (id<outindex) {
256 Int_t sec = id/nir;
257 Int_t row = id - sec*nir;
258 sec %= fkNIS;
259 AliTPCRow &padrow=fInnerSec[sec][row];
260 while (ncl--) {
261 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
262 padrow.InsertCluster(c,sec,row);
263 }
264 } else {
265 id -= outindex;
266 Int_t sec = id/nor;
267 Int_t row = id - sec*nor;
268 sec %= fkNOS;
269 AliTPCRow &padrow=fOuterSec[sec][row];
270 while (ncl--) {
271 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
272 padrow.InsertCluster(c,sec+fkNIS,row);
273 }
274 }
b9de75e1 275
bcc04e2a 276 carray.GetArray()->Clear();
277 }
278 delete cTree;
b9de75e1 279}
280
281//_____________________________________________________________________________
bcc04e2a 282void AliTPCtracker::UnloadClusters() {
b9de75e1 283 //-----------------------------------------------------------------
bcc04e2a 284 // This function unloads TPC clusters.
b9de75e1 285 //-----------------------------------------------------------------
bcc04e2a 286 Int_t i;
287 for (i=0; i<fkNIS; i++) {
288 Int_t nr=fInnerSec->GetNRows();
289 for (Int_t n=0; n<nr; n++) fInnerSec[i][n].ResetClusters();
290 }
291 for (i=0; i<fkNOS; i++) {
292 Int_t nr=fOuterSec->GetNRows();
293 for (Int_t n=0; n<nr; n++) fOuterSec[i][n].ResetClusters();
b9de75e1 294 }
b9de75e1 295}
296
297//_____________________________________________________________________________
298Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
73042f01 299 //-----------------------------------------------------------------
300 // This function tries to find a track prolongation.
301 //-----------------------------------------------------------------
b9de75e1 302 Double_t xt=t.GetX();
303 const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
304 Int_t(0.5*fSectors->GetNRows());
73042f01 305 Int_t tryAgain=kSKIP;
73042f01 306
b9de75e1 307 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
308 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
309 if (alpha < 0. ) alpha += 2.*TMath::Pi();
310 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
311
f03e3423 312 Int_t nrows=fSectors->GetRowNumber(xt)-1;
313 for (Int_t nr=nrows; nr>=rf; nr--) {
b9de75e1 314 Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
73042f01 315 if (!t.PropagateTo(x)) return 0;
316
317 AliTPCcluster *cl=0;
318 UInt_t index=0;
b9de75e1 319 Double_t maxchi2=kMaxCHI2;
320 const AliTPCRow &krow=fSectors[s][nr];
321 Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
322 Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
73042f01 323 Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
324 Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
325
b9de75e1 326 if (road>kMaxROAD) {
73042f01 327 if (t.GetNumberOfClusters()>4)
328 cerr<<t.GetNumberOfClusters()
329 <<"FindProlongation warning: Too broad road !\n";
330 return 0;
331 }
332
333 if (krow) {
334 for (Int_t i=krow.Find(y-road); i<krow; i++) {
be5b9287 335 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
73042f01 336 if (c->GetY() > y+road) break;
337 if (c->IsUsed()) continue;
338 if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
339 Double_t chi2=t.GetPredictedChi2(c);
340 if (chi2 > maxchi2) continue;
341 maxchi2=chi2;
342 cl=c;
343 index=krow.GetIndex(i);
344 }
345 }
346 if (cl) {
b9de75e1 347 Float_t l=fSectors->GetPadPitchWidth();
024a7fe9 348 Float_t corr=1.; if (nr>63) corr=0.67; // new (third) pad response !
349 t.SetSampledEdx(cl->GetQ()/l*corr,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 }
bcc04e2a 386 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
387 else { if (sec < fkNIS) row=-1; }
b9de75e1 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 }
bcc04e2a 431 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
432 else { if (sec < 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();
024a7fe9 456 Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
457 seed.SetSampledEdx(cl->GetQ()/l*corr,seed.GetNumberOfClusters());
b9de75e1 458 seed.Update(cl,maxchi2,index);
459 }
460
461 }
462
463 seed.SetLabel(nc);
464
465 return 1;
73042f01 466}
467
468//_____________________________________________________________________________
b9de75e1 469void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
73042f01 470 //-----------------------------------------------------------------
471 // This function creates track seeds.
472 //-----------------------------------------------------------------
473 Double_t x[5], c[15];
474
bcc04e2a 475 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
73042f01 476 Double_t cs=cos(alpha), sn=sin(alpha);
477
bcc04e2a 478 Double_t x1 =fSectors->GetX(i1);
479 Double_t xx2=fSectors->GetX(i2);
73042f01 480
bcc04e2a 481 for (Int_t ns=0; ns<fN; ns++) {
482 Int_t nl=fSectors[(ns-1+fN)%fN][i2];
483 Int_t nm=fSectors[ns][i2];
484 Int_t nu=fSectors[(ns+1)%fN][i2];
485 const AliTPCRow& kr1=fSectors[ns][i1];
73042f01 486 for (Int_t is=0; is < kr1; is++) {
487 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
488 for (Int_t js=0; js < nl+nm+nu; js++) {
489 const AliTPCcluster *kcl;
490 Double_t x2, y2, z2;
7f6ddf58 491 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
73042f01 492
493 if (js<nl) {
bcc04e2a 494 const AliTPCRow& kr2=fSectors[(ns-1+fN)%fN][i2];
73042f01 495 kcl=kr2[js];
496 y2=kcl->GetY(); z2=kcl->GetZ();
497 x2= xx2*cs+y2*sn;
498 y2=-xx2*sn+y2*cs;
499 } else
500 if (js<nl+nm) {
bcc04e2a 501 const AliTPCRow& kr2=fSectors[ns][i2];
73042f01 502 kcl=kr2[js-nl];
503 x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
504 } else {
bcc04e2a 505 const AliTPCRow& kr2=fSectors[(ns+1)%fN][i2];
73042f01 506 kcl=kr2[js-nl-nm];
507 y2=kcl->GetY(); z2=kcl->GetZ();
508 x2=xx2*cs-y2*sn;
509 y2=xx2*sn+y2*cs;
510 }
511
7f6ddf58 512 Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
73042f01 513 if (TMath::Abs(zz-z2)>5.) continue;
514
515 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
516 if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
517
518 x[0]=y1;
519 x[1]=z1;
b9de75e1 520 x[4]=f1(x1,y1,x2,y2,x3,y3);
521 if (TMath::Abs(x[4]) >= 0.0066) continue;
be5b9287 522 x[2]=f2(x1,y1,x2,y2,x3,y3);
b9de75e1 523 //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
524 x[3]=f3(x1,y1,x2,y2,z1,z2);
525 if (TMath::Abs(x[3]) > 1.2) continue;
be5b9287 526 Double_t a=asin(x[2]);
b9de75e1 527 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
7f6ddf58 528 if (TMath::Abs(zv-z3)>10.) continue;
73042f01 529
530 Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
531 Double_t sy2=kcl->GetSigmaY2(), sz2=kcl->GetSigmaZ2();
024a7fe9 532 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
533 Double_t sy3=25000*x[4]*x[4]+0.1, 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);
bcc04e2a 557 Float_t l=fSectors->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
bcc04e2a 628 LoadClusters();
629
b9de75e1 630 out->cd();
7f6ddf58 631
632 char tname[100];
bcc04e2a 633 sprintf(tname,"TreeT_TPC_%d",GetEventNumber());
7f6ddf58 634 TTree tracktree(tname,"Tree with TPC tracks");
b9de75e1 635 AliTPCtrack *iotrack=0;
636 tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
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);
bcc04e2a 643 fSectors=fOuterSec; fN=fkNOS;
644 fSeeds=new TObjArray(15000);
b9de75e1 645 MakeSeeds(nup-1, nup-1-gap);
646 MakeSeeds(nup-1-shift, nup-1-shift-gap);
bcc04e2a 647 }
b9de75e1 648 fSeeds->Sort();
73042f01 649
bcc04e2a 650 Int_t found=0;
b9de75e1 651 Int_t nseed=fSeeds->GetEntriesFast();
bcc04e2a 652 for (Int_t i=0; i<nseed; i++) {
653 //tracking in the outer sectors
654 fSectors=fOuterSec; fN=fkNOS;
655
b9de75e1 656 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
bcc04e2a 657 if (!FollowProlongation(t)) {
658 delete fSeeds->RemoveAt(i);
73042f01 659 continue;
660 }
73042f01 661
bcc04e2a 662 //tracking in the inner sectors
663 fSectors=fInnerSec; fN=fkNIS;
73042f01 664
b9de75e1 665 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
73042f01 666 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
667 if (alpha < 0. ) alpha += 2.*TMath::Pi();
b9de75e1 668 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
73042f01 669
bcc04e2a 670 alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
73042f01 671
672 if (t.Rotate(alpha)) {
b9de75e1 673 if (FollowProlongation(t)) {
73042f01 674 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
675 t.CookdEdx();
7f6ddf58 676 CookLabel(pt,0.1); //For comparison only
73042f01 677 iotrack=pt;
678 tracktree.Fill();
bcc04e2a 679 UseClusters(&t);
b9de75e1 680 cerr<<found++<<'\r';
73042f01 681 }
682 }
683 }
b9de75e1 684 delete fSeeds->RemoveAt(i);
bcc04e2a 685 }
b9de75e1 686
7f6ddf58 687 tracktree.Write();
73042f01 688
b9de75e1 689 cerr<<"Number of found tracks : "<<found<<endl;
690
691 savedir->cd();
692
bcc04e2a 693 UnloadClusters();
694 fSeeds->Clear(); delete fSeeds; fSeeds=0;
695
b9de75e1 696 return 0;
697}
698
699//_____________________________________________________________________________
700Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
701 //-----------------------------------------------------------------
702 // This function propagates tracks back through the TPC.
703 //-----------------------------------------------------------------
704 fSeeds=new TObjArray(15000);
705
706 TFile *in=(TFile*)inp;
707 TDirectory *savedir=gDirectory;
708
709 if (!in->IsOpen()) {
710 cerr<<"AliTPCtracker::PropagateBack(): ";
711 cerr<<"file with back propagated ITS tracks is not open !\n";
712 return 1;
73042f01 713 }
714
b9de75e1 715 if (!out->IsOpen()) {
716 cerr<<"AliTPCtracker::PropagateBack(): ";
717 cerr<<"file for back propagated TPC tracks is not open !\n";
718 return 2;
719 }
720
bcc04e2a 721 LoadClusters();
722
b9de75e1 723 in->cd();
f38c8ae5 724 TTree *bckTree=(TTree*)in->Get("TreeT_ITSb_0");
b9de75e1 725 if (!bckTree) {
726 cerr<<"AliTPCtracker::PropagateBack() ";
727 cerr<<"can't get a tree with back propagated ITS tracks !\n";
728 return 3;
729 }
730 AliTPCtrack *bckTrack=new AliTPCtrack;
731 bckTree->SetBranchAddress("tracks",&bckTrack);
732
f38c8ae5 733 TTree *tpcTree=(TTree*)in->Get("TreeT_TPC_0");
b9de75e1 734 if (!tpcTree) {
735 cerr<<"AliTPCtracker::PropagateBack() ";
736 cerr<<"can't get a tree with TPC tracks !\n";
737 return 4;
738 }
739 AliTPCtrack *tpcTrack=new AliTPCtrack;
740 tpcTree->SetBranchAddress("tracks",&tpcTrack);
741
742//*** Prepare an array of tracks to be back propagated
743 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
744 Int_t nrows=nlow+nup;
745
746 TObjArray tracks(15000);
747 Int_t i=0,j=0;
748 Int_t tpcN=(Int_t)tpcTree->GetEntries();
749 Int_t bckN=(Int_t)bckTree->GetEntries();
750 if (j<bckN) bckTree->GetEvent(j++);
751 for (i=0; i<tpcN; i++) {
752 tpcTree->GetEvent(i);
753 Double_t alpha=tpcTrack->GetAlpha();
754 if (j<bckN &&
755 TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
756 if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
757 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
758 bckTree->GetEvent(j++);
759 } else {
760 tpcTrack->ResetCovariance();
761 fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
762 }
763 tracks.AddLast(new AliTPCtrack(*tpcTrack));
764 }
765
766 out->cd();
f38c8ae5 767 TTree backTree("TreeT_TPCb_0","Tree with back propagated TPC tracks");
b9de75e1 768 AliTPCtrack *otrack=0;
769 backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
770
771//*** Back propagation through inner sectors
bcc04e2a 772 fSectors=fInnerSec; fN=fkNIS;
b9de75e1 773 Int_t nseed=fSeeds->GetEntriesFast();
774 for (i=0; i<nseed; i++) {
775 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
776 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
73042f01 777
b9de75e1 778 Int_t nc=t.GetNumberOfClusters();
779 s.SetLabel(nc-1); //set number of the cluster to start with
73042f01 780
b9de75e1 781 if (FollowBackProlongation(s,t)) {
782 UseClusters(&s);
783 continue;
784 }
785 delete fSeeds->RemoveAt(i);
786 }
b9de75e1 787
788//*** Back propagation through outer sectors
bcc04e2a 789 fSectors=fOuterSec; fN=fkNOS;
b9de75e1 790 Int_t found=0;
791 for (i=0; i<nseed; i++) {
792 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
793 if (!ps) continue;
794 Int_t nc=s.GetNumberOfClusters();
795 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
796
797 Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
798 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
799 if (alpha < 0. ) alpha += 2.*TMath::Pi();
800 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
801
802 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
803 alpha-=s.GetAlpha();
804
805 if (s.Rotate(alpha)) {
806 if (FollowBackProlongation(s,t)) {
807 if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
808 s.CookdEdx();
809 s.SetLabel(t.GetLabel());
810 UseClusters(&s,nc);
811 otrack=ps;
812 backTree.Fill();
813 cerr<<found++<<'\r';
814 continue;
815 }
816 }
817 }
818 delete fSeeds->RemoveAt(i);
819 }
b9de75e1 820
821 backTree.Write();
73042f01 822 savedir->cd();
b9de75e1 823 cerr<<"Number of seeds: "<<nseed<<endl;
824 cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
825 cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
826
827 delete bckTrack;
828 delete tpcTrack;
be5b9287 829
f38c8ae5 830 delete bckTree; //Thanks to Mariana Bondila
831 delete tpcTree; //Thanks to Mariana Bondila
832
bcc04e2a 833 UnloadClusters();
834
be5b9287 835 return 0;
73042f01 836}
837
b9de75e1 838//_________________________________________________________________________
839AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
840 //--------------------------------------------------------------------
841 // Return pointer to a given cluster
842 //--------------------------------------------------------------------
843 Int_t sec=(index&0xff000000)>>24;
844 Int_t row=(index&0x00ff0000)>>16;
845 Int_t ncl=(index&0x0000ffff)>>00;
846
bcc04e2a 847 const AliTPCcluster *cl=0;
848
849 if (sec<fkNIS) {
850 cl=fInnerSec[sec][row].GetUnsortedCluster(ncl);
851 } else {
852 sec-=fkNIS;
853 cl=fOuterSec[sec][row].GetUnsortedCluster(ncl);
854 }
855
856 return (AliCluster*)cl;
b9de75e1 857}
858
859//__________________________________________________________________________
860void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
861 //--------------------------------------------------------------------
862 //This function "cooks" a track label. If label<0, this track is fake.
863 //--------------------------------------------------------------------
864 Int_t noc=t->GetNumberOfClusters();
865 Int_t *lb=new Int_t[noc];
866 Int_t *mx=new Int_t[noc];
867 AliCluster **clusters=new AliCluster*[noc];
868
869 Int_t i;
870 for (i=0; i<noc; i++) {
871 lb[i]=mx[i]=0;
872 Int_t index=t->GetClusterIndex(i);
873 clusters[i]=GetCluster(index);
874 }
875
876 Int_t lab=123456789;
877 for (i=0; i<noc; i++) {
878 AliCluster *c=clusters[i];
879 lab=TMath::Abs(c->GetLabel(0));
880 Int_t j;
881 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
882 lb[j]=lab;
883 (mx[j])++;
884 }
885
886 Int_t max=0;
887 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
888
889 for (i=0; i<noc; i++) {
890 AliCluster *c=clusters[i];
891 if (TMath::Abs(c->GetLabel(1)) == lab ||
892 TMath::Abs(c->GetLabel(2)) == lab ) max++;
893 }
894
895 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
896
897 else {
898 Int_t tail=Int_t(0.10*noc);
899 max=0;
900 for (i=1; i<=tail; i++) {
901 AliCluster *c=clusters[noc-i];
902 if (lab == TMath::Abs(c->GetLabel(0)) ||
903 lab == TMath::Abs(c->GetLabel(1)) ||
904 lab == TMath::Abs(c->GetLabel(2))) max++;
905 }
906 if (max < Int_t(0.5*tail)) lab=-lab;
907 }
908
909 t->SetLabel(lab);
910
911 delete[] lb;
912 delete[] mx;
913 delete[] clusters;
914}
915
916//_________________________________________________________________________
917void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
918 //-----------------------------------------------------------------------
919 // Setup inner sector
920 //-----------------------------------------------------------------------
921 if (f==0) {
922 fAlpha=par->GetInnerAngle();
923 fAlphaShift=par->GetInnerAngleShift();
924 fPadPitchWidth=par->GetInnerPadPitchWidth();
024a7fe9 925 f1PadPitchLength=par->GetInnerPadPitchLength();
926 f2PadPitchLength=f1PadPitchLength;
b9de75e1 927
928 fN=par->GetNRowLow();
929 fRow=new AliTPCRow[fN];
930 for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
931 } else {
932 fAlpha=par->GetOuterAngle();
933 fAlphaShift=par->GetOuterAngleShift();
934 fPadPitchWidth=par->GetOuterPadPitchWidth();
f03e3423 935 f1PadPitchLength=par->GetOuter1PadPitchLength();
936 f2PadPitchLength=par->GetOuter2PadPitchLength();
b9de75e1 937
938 fN=par->GetNRowUp();
939 fRow=new AliTPCRow[fN];
f03e3423 940 for (Int_t i=0; i<fN; i++){
941 fRow[i].SetX(par->GetPadRowRadiiUp(i));
942 }
b9de75e1 943 }
944}
945
73042f01 946//_________________________________________________________________________
bcc04e2a 947void AliTPCtracker::
948AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) {
73042f01 949 //-----------------------------------------------------------------------
950 // Insert a cluster into this pad row in accordence with its y-coordinate
951 //-----------------------------------------------------------------------
b9de75e1 952 if (fN==kMaxClusterPerRow) {
73042f01 953 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
954 }
bcc04e2a 955
956 Int_t index=(((sec<<8)+row)<<16)+fN;
957
958 if (fN==0) {
959 fIndex[0]=index;
960 fClusterArray[0]=*c; fClusters[fN++]=fClusterArray;
961 return;
962 }
963
964 if (fN==fSize) {
965 Int_t size=fSize*2;
966 AliTPCcluster *buff=new AliTPCcluster[size];
967 memcpy(buff,fClusterArray,fSize*sizeof(AliTPCcluster));
968 for (Int_t i=0; i<fN; i++)
969 fClusters[i]=buff+(fClusters[i]-fClusterArray);
970 delete[] fClusterArray;
971 fClusterArray=buff;
972 fSize=size;
973 }
974
73042f01 975 Int_t i=Find(c->GetY());
976 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
977 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
bcc04e2a 978 fIndex[i]=index;
979 fClusters[i]=fClusterArray+fN; fClusterArray[fN++]=*c;
73042f01 980}
981
982//___________________________________________________________________
983Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
984 //-----------------------------------------------------------------------
985 // Return the index of the nearest cluster
986 //-----------------------------------------------------------------------
11c7b548 987 if(fN<=0) return 0;
73042f01 988 if (y <= fClusters[0]->GetY()) return 0;
989 if (y > fClusters[fN-1]->GetY()) return fN;
990 Int_t b=0, e=fN-1, m=(b+e)/2;
991 for (; b<e; m=(b+e)/2) {
992 if (y > fClusters[m]->GetY()) b=m+1;
993 else e=m;
994 }
995 return m;
996}
997
998//_____________________________________________________________________________
999void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1000 //-----------------------------------------------------------------
1001 // This funtion calculates dE/dX within the "low" and "up" cuts.
1002 //-----------------------------------------------------------------
1003 Int_t i;
be9c5115 1004 Int_t nc=GetNumberOfClusters();
bcc04e2a 1005
1006 Int_t swap;//stupid sorting
1007 do {
1008 swap=0;
1009 for (i=0; i<nc-1; i++) {
1010 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
1011 Float_t tmp=fdEdxSample[i];
1012 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
1013 swap++;
1014 }
1015 } while (swap);
73042f01 1016
be9c5115 1017 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
73042f01 1018 Float_t dedx=0;
bcc04e2a 1019 for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
73042f01 1020 dedx /= (nu-nl+1);
1021 SetdEdx(dedx);
7f6ddf58 1022
1023 //Very rough PID
1024 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1025
024a7fe9 1026 Double_t log1=TMath::Log(p+0.45), log2=TMath::Log(p+0.12);
7f6ddf58 1027 if (p<0.6) {
024a7fe9 1028 if (dedx < 34 + 30/(p+0.45)/(p+0.45) + 24*log1) {SetMass(0.13957); return;}
1029 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.49368); return;}
7f6ddf58 1030 SetMass(0.93827); return;
1031 }
1032
1033 if (p<1.2) {
024a7fe9 1034 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.13957); return;}
7f6ddf58 1035 SetMass(0.93827); return;
1036 }
1037
1038 SetMass(0.13957); return;
1039
73042f01 1040}
1041
73042f01 1042
1043