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