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