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