PCONE QBEU had last z defined twice.
[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$
f03e3423 18Revision 1.16 2001/11/08 16:39:03 hristov
19Additional protection (M.Masera)
20
1c8d9aad 21Revision 1.15 2001/11/08 16:36:33 hristov
22Updated 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)
23
7f6ddf58 24Revision 1.14 2001/10/21 19:04:55 hristov
25Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
26
f38c8ae5 27Revision 1.13 2001/07/23 12:01:30 hristov
28Initialisation added
29
7afb86d1 30Revision 1.12 2001/07/20 14:32:44 kowal2
31Processing of many events possible now
32
afc42102 33Revision 1.11 2001/05/23 08:50:10 hristov
34Weird inline removed
35
bdbd0f7a 36Revision 1.10 2001/05/16 14:57:25 alibrary
37New files for folders and Stack
38
9e1a0ddb 39Revision 1.9 2001/05/11 07:16:56 hristov
40Fix needed on Sun and Alpha
41
3d39dd8b 42Revision 1.8 2001/05/08 15:00:15 hristov
43Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
44
be9c5115 45Revision 1.5 2000/12/20 07:51:59 kowal2
46Changes suggested by Alessandra and Paolo to avoid overlapped
47data fields in encapsulated classes.
48
d4cf1daa 49Revision 1.4 2000/11/02 07:27:16 kowal2
50code corrections
51
98ab53f4 52Revision 1.2 2000/06/30 12:07:50 kowal2
53Updated from the TPC-PreRelease branch
54
73042f01 55Revision 1.1.2.1 2000/06/25 08:53:55 kowal2
56Splitted from AliTPCtracking
57
58*/
59
60//-------------------------------------------------------
61// Implementation of the TPC tracker
62//
63// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
64//-------------------------------------------------------
65
73042f01 66#include <TObjArray.h>
67#include <TFile.h>
be5b9287 68#include <TTree.h>
69#include <iostream.h>
70
71#include "AliTPCtracker.h"
72#include "AliTPCcluster.h"
b9de75e1 73#include "AliTPCParam.h"
73042f01 74#include "AliTPCClustersRow.h"
75
b9de75e1 76//_____________________________________________________________________________
afc42102 77AliTPCtracker::AliTPCtracker(const AliTPCParam *par, Int_t eventn):
7f6ddf58 78AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
b9de75e1 79{
80 //---------------------------------------------------------------------
81 // The main TPC tracker constructor
82 //---------------------------------------------------------------------
83 fInnerSec=new AliTPCSector[fkNIS];
84 fOuterSec=new AliTPCSector[fkNOS];
85
86 Int_t i;
87 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
88 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
89
90 fN=0; fSectors=0;
91
92 fClustersArray.Setup(par);
93 fClustersArray.SetClusterType("AliTPCcluster");
b9de75e1 94
afc42102 95 char cname[100];
96 if (eventn==-1) {
97 sprintf(cname,"TreeC_TPC");
98 }
99 else {
100 sprintf(cname,"TreeC_TPC_%d",eventn);
101 }
102
103 fClustersArray.ConnectTree(cname);
104
105 fEventN = eventn;
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//_____________________________________________________________________________
b9de75e1 215void AliTPCtracker::LoadOuterSectors() {
216 //-----------------------------------------------------------------
217 // This function fills outer TPC sectors with clusters.
218 //-----------------------------------------------------------------
219 UInt_t index;
220 Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
221 for (Int_t i=0; i<j; i++) {
222 AliSegmentID *s=fClustersArray.LoadEntry(i);
223 Int_t sec,row;
224 AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
225 par->AdjustSectorRow(s->GetID(),sec,row);
226 if (sec<fkNIS*2) continue;
227 AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
228 Int_t ncl=clrow->GetArray()->GetEntriesFast();
229 while (ncl--) {
230 AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
231 index=(((sec<<8)+row)<<16)+ncl;
232 fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
233 }
234 }
235
236 fN=fkNOS;
237 fSectors=fOuterSec;
238}
239
240//_____________________________________________________________________________
241void AliTPCtracker::UnloadOuterSectors() {
242 //-----------------------------------------------------------------
243 // This function clears outer TPC sectors.
244 //-----------------------------------------------------------------
245 Int_t nup=fOuterSec->GetNRows();
246 for (Int_t i=0; i<fkNOS; i++) {
247 for (Int_t j=0; j<nup; j++) {
248 if (fClustersArray.GetRow(i+fkNIS*2,j))
249 fClustersArray.ClearRow(i+fkNIS*2,j);
250 if (fClustersArray.GetRow(i+fkNIS*2+fkNOS,j))
251 fClustersArray.ClearRow(i+fkNIS*2+fkNOS,j);
252 }
253 }
254
255 fN=0;
256 fSectors=0;
257}
258
259//_____________________________________________________________________________
260void AliTPCtracker::LoadInnerSectors() {
261 //-----------------------------------------------------------------
262 // This function fills inner TPC sectors with clusters.
263 //-----------------------------------------------------------------
264 UInt_t index;
265 Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
266 for (Int_t i=0; i<j; i++) {
267 AliSegmentID *s=fClustersArray.LoadEntry(i);
268 Int_t sec,row;
269 AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
270 par->AdjustSectorRow(s->GetID(),sec,row);
271 if (sec>=fkNIS*2) continue;
272 AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
273 Int_t ncl=clrow->GetArray()->GetEntriesFast();
274 while (ncl--) {
275 AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
276 index=(((sec<<8)+row)<<16)+ncl;
277 fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
278 }
279 }
280
281 fN=fkNIS;
282 fSectors=fInnerSec;
283}
284
285//_____________________________________________________________________________
286void AliTPCtracker::UnloadInnerSectors() {
287 //-----------------------------------------------------------------
288 // This function clears inner TPC sectors.
289 //-----------------------------------------------------------------
290 Int_t nlow=fInnerSec->GetNRows();
291 for (Int_t i=0; i<fkNIS; i++) {
292 for (Int_t j=0; j<nlow; j++) {
293 if (fClustersArray.GetRow(i,j)) fClustersArray.ClearRow(i,j);
294 if (fClustersArray.GetRow(i+fkNIS,j)) fClustersArray.ClearRow(i+fkNIS,j);
295 }
296 }
297
298 fN=0;
299 fSectors=0;
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();
73042f01 353 t.SetSampledEdx(cl->GetQ()/l,t.GetNumberOfClusters());
be9c5115 354 if (!t.Update(cl,maxchi2,index)) {
355 if (!tryAgain--) return 0;
356 } else tryAgain=kSKIP;
73042f01 357 } else {
358 if (tryAgain==0) break;
359 if (y > ymax) {
b9de75e1 360 s = (s+1) % fN;
361 if (!t.Rotate(fSectors->GetAlpha())) return 0;
73042f01 362 } else if (y <-ymax) {
b9de75e1 363 s = (s-1+fN) % fN;
364 if (!t.Rotate(-fSectors->GetAlpha())) return 0;
73042f01 365 }
366 tryAgain--;
367 }
368 }
369
370 return 1;
b9de75e1 371}
372
373//_____________________________________________________________________________
374Int_t AliTPCtracker::FollowBackProlongation
375(AliTPCseed& seed, const AliTPCtrack &track) {
376 //-----------------------------------------------------------------
377 // This function propagates tracks back through the TPC
378 //-----------------------------------------------------------------
379 Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
380 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
381 if (alpha < 0. ) alpha += 2.*TMath::Pi();
382 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
383
384 Int_t idx=-1, sec=-1, row=-1;
385 Int_t nc=seed.GetLabel(); //index of the cluster to start with
386 if (nc--) {
387 idx=track.GetClusterIndex(nc);
388 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
389 }
390 if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; }
391 else { if (sec < 2*fkNIS) row=-1; }
392
393 Int_t nr=fSectors->GetNRows();
394 for (Int_t i=0; i<nr; i++) {
395 Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
396
397 if (!seed.PropagateTo(x)) return 0;
398
399 Double_t y=seed.GetY();
400 if (y > ymax) {
401 s = (s+1) % fN;
402 if (!seed.Rotate(fSectors->GetAlpha())) return 0;
403 } else if (y <-ymax) {
404 s = (s-1+fN) % fN;
405 if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
406 }
407
408 AliTPCcluster *cl=0;
409 Int_t index=0;
410 Double_t maxchi2=kMaxCHI2;
411 Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
412 Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
413 Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
414 Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
415 if (road>kMaxROAD) {
416 cerr<<seed.GetNumberOfClusters()
417 <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n";
418 return 0;
419 }
420
421
422 Int_t accepted=seed.GetNumberOfClusters();
423 if (row==i) {
424 //try to accept already found cluster
425 AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
426 Double_t chi2;
427 if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
428 index=idx; cl=c; maxchi2=chi2;
429 } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
430
431 if (nc--) {
432 idx=track.GetClusterIndex(nc);
433 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
434 }
435 if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; }
436 else { if (sec < 2*fkNIS) row=-1; }
73042f01 437
b9de75e1 438 }
439 if (!cl) {
440 //try to fill the gap
441 const AliTPCRow &krow=fSectors[s][i];
442 if (accepted>27)
443 if (krow) {
444 for (Int_t i=krow.Find(y-road); i<krow; i++) {
445 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
446 if (c->GetY() > y+road) break;
447 if (c->IsUsed()) continue;
448 if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
449 Double_t chi2=seed.GetPredictedChi2(c);
450 if (chi2 > maxchi2) continue;
451 maxchi2=chi2;
452 cl=c;
453 index=krow.GetIndex(i);
454 }
455 }
456 }
457
458 if (cl) {
459 Float_t l=fSectors->GetPadPitchWidth();
460 seed.SetSampledEdx(cl->GetQ()/l,seed.GetNumberOfClusters());
461 seed.Update(cl,maxchi2,index);
462 }
463
464 }
465
466 seed.SetLabel(nc);
467
468 return 1;
73042f01 469}
470
471//_____________________________________________________________________________
b9de75e1 472void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
73042f01 473 //-----------------------------------------------------------------
474 // This function creates track seeds.
475 //-----------------------------------------------------------------
b9de75e1 476 if (fSeeds==0) fSeeds=new TObjArray(15000);
477
73042f01 478 Double_t x[5], c[15];
479
b9de75e1 480 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
73042f01 481 Double_t cs=cos(alpha), sn=sin(alpha);
482
b9de75e1 483 Double_t x1 =fOuterSec->GetX(i1);
484 Double_t xx2=fOuterSec->GetX(i2);
73042f01 485
b9de75e1 486 for (Int_t ns=0; ns<fkNOS; ns++) {
487 Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
488 Int_t nm=fOuterSec[ns][i2];
489 Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
490 const AliTPCRow& kr1=fOuterSec[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) {
b9de75e1 499 const AliTPCRow& kr2=fOuterSec[(ns-1+fkNOS)%fkNOS][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) {
b9de75e1 506 const AliTPCRow& kr2=fOuterSec[ns][i2];
73042f01 507 kcl=kr2[js-nl];
508 x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
509 } else {
b9de75e1 510 const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][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();
7f6ddf58 537 Double_t sy3=400*3./12., sy=0.1, sz=0.1;
73042f01 538
b9de75e1 539 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
540 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
541 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
be5b9287 542 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
543 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
b9de75e1 544 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
545 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
546 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
547 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
548 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
73042f01 549
550 c[0]=sy1;
551 c[1]=0.; c[2]=sz1;
b9de75e1 552 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
553 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
554 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
555 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
556 c[13]=f30*sy1*f40+f32*sy2*f42;
557 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
73042f01 558
559 UInt_t index=kr1.GetIndex(is);
560 AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
b9de75e1 561 Float_t l=fOuterSec->GetPadPitchWidth();
73042f01 562 track->SetSampledEdx(kr1[is]->GetQ()/l,0);
563
b9de75e1 564 Int_t rc=FollowProlongation(*track, i2);
be9c5115 565 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
b9de75e1 566 else fSeeds->AddLast(track);
73042f01 567 }
568 }
569 }
570}
571
572//_____________________________________________________________________________
b9de75e1 573Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
73042f01 574 //-----------------------------------------------------------------
b9de75e1 575 // This function reades track seeds.
73042f01 576 //-----------------------------------------------------------------
577 TDirectory *savedir=gDirectory;
578
b9de75e1 579 TFile *in=(TFile*)inp;
580 if (!in->IsOpen()) {
581 cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
be5b9287 582 return 1;
73042f01 583 }
584
b9de75e1 585 in->cd();
586 TTree *seedTree=(TTree*)in->Get("Seeds");
587 if (!seedTree) {
588 cerr<<"AliTPCtracker::ReadSeeds(): ";
589 cerr<<"can't get a tree with track seeds !\n";
590 return 2;
591 }
592 AliTPCtrack *seed=new AliTPCtrack;
593 seedTree->SetBranchAddress("tracks",&seed);
594
595 if (fSeeds==0) fSeeds=new TObjArray(15000);
73042f01 596
b9de75e1 597 Int_t n=(Int_t)seedTree->GetEntries();
598 for (Int_t i=0; i<n; i++) {
599 seedTree->GetEvent(i);
600 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
601 }
602
603 delete seed;
f38c8ae5 604
605 delete seedTree; //Thanks to Mariana Bondila
606
b9de75e1 607 savedir->cd();
73042f01 608
b9de75e1 609 return 0;
610}
73042f01 611
b9de75e1 612//_____________________________________________________________________________
613Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
614 //-----------------------------------------------------------------
615 // This is a track finder.
616 //-----------------------------------------------------------------
617 TDirectory *savedir=gDirectory;
73042f01 618
b9de75e1 619 if (inp) {
620 TFile *in=(TFile*)inp;
621 if (!in->IsOpen()) {
622 cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
623 return 1;
624 }
625 }
73042f01 626
b9de75e1 627 if (!out->IsOpen()) {
628 cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
629 return 2;
73042f01 630 }
631
b9de75e1 632 out->cd();
7f6ddf58 633
634 char tname[100];
635 sprintf(tname,"TreeT_TPC_%d",fEventN);
636 TTree tracktree(tname,"Tree with TPC tracks");
b9de75e1 637 AliTPCtrack *iotrack=0;
638 tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
639
640 LoadOuterSectors();
641
73042f01 642 //find track seeds
b9de75e1 643 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
73042f01 644 Int_t nrows=nlow+nup;
b9de75e1 645 if (fSeeds==0) {
646 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
647 MakeSeeds(nup-1, nup-1-gap);
648 MakeSeeds(nup-1-shift, nup-1-shift-gap);
649 }
650 fSeeds->Sort();
73042f01 651
652 //tracking in outer sectors
b9de75e1 653 Int_t nseed=fSeeds->GetEntriesFast();
654 Int_t i;
73042f01 655 for (i=0; i<nseed; i++) {
b9de75e1 656 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
657 if (FollowProlongation(t)) {
658 UseClusters(&t);
73042f01 659 continue;
660 }
b9de75e1 661 delete fSeeds->RemoveAt(i);
73042f01 662 }
7f6ddf58 663 //UnloadOuterSectors();
73042f01 664
665 //tracking in inner sectors
b9de75e1 666 LoadInnerSectors();
73042f01 667 Int_t found=0;
668 for (i=0; i<nseed; i++) {
b9de75e1 669 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
73042f01 670 if (!pt) continue;
671 Int_t nc=t.GetNumberOfClusters();
672
b9de75e1 673 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
73042f01 674 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
675 if (alpha < 0. ) alpha += 2.*TMath::Pi();
b9de75e1 676 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
73042f01 677
b9de75e1 678 alpha=ns*fInnerSec->GetAlpha() + fInnerSec->GetAlphaShift() - t.GetAlpha();
73042f01 679
680 if (t.Rotate(alpha)) {
b9de75e1 681 if (FollowProlongation(t)) {
73042f01 682 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
683 t.CookdEdx();
7f6ddf58 684 CookLabel(pt,0.1); //For comparison only
73042f01 685 iotrack=pt;
686 tracktree.Fill();
b9de75e1 687 UseClusters(&t,nc);
688 cerr<<found++<<'\r';
73042f01 689 }
690 }
691 }
b9de75e1 692 delete fSeeds->RemoveAt(i);
73042f01 693 }
b9de75e1 694 UnloadInnerSectors();
7f6ddf58 695 UnloadOuterSectors();
b9de75e1 696
7f6ddf58 697 tracktree.Write();
73042f01 698
b9de75e1 699 cerr<<"Number of found tracks : "<<found<<endl;
700
701 savedir->cd();
702
703 return 0;
704}
705
706//_____________________________________________________________________________
707Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
708 //-----------------------------------------------------------------
709 // This function propagates tracks back through the TPC.
710 //-----------------------------------------------------------------
711 fSeeds=new TObjArray(15000);
712
713 TFile *in=(TFile*)inp;
714 TDirectory *savedir=gDirectory;
715
716 if (!in->IsOpen()) {
717 cerr<<"AliTPCtracker::PropagateBack(): ";
718 cerr<<"file with back propagated ITS tracks is not open !\n";
719 return 1;
73042f01 720 }
721
b9de75e1 722 if (!out->IsOpen()) {
723 cerr<<"AliTPCtracker::PropagateBack(): ";
724 cerr<<"file for back propagated TPC tracks is not open !\n";
725 return 2;
726 }
727
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
777 LoadInnerSectors();
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 }
792 UnloadInnerSectors();
793
794//*** Back propagation through outer sectors
795 LoadOuterSectors();
796 Int_t found=0;
797 for (i=0; i<nseed; i++) {
798 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
799 if (!ps) continue;
800 Int_t nc=s.GetNumberOfClusters();
801 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
802
803 Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
804 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
805 if (alpha < 0. ) alpha += 2.*TMath::Pi();
806 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
807
808 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
809 alpha-=s.GetAlpha();
810
811 if (s.Rotate(alpha)) {
812 if (FollowBackProlongation(s,t)) {
813 if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
814 s.CookdEdx();
815 s.SetLabel(t.GetLabel());
816 UseClusters(&s,nc);
817 otrack=ps;
818 backTree.Fill();
819 cerr<<found++<<'\r';
820 continue;
821 }
822 }
823 }
824 delete fSeeds->RemoveAt(i);
825 }
826 UnloadOuterSectors();
827
828 backTree.Write();
73042f01 829 savedir->cd();
b9de75e1 830 cerr<<"Number of seeds: "<<nseed<<endl;
831 cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
832 cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
833
834 delete bckTrack;
835 delete tpcTrack;
be5b9287 836
f38c8ae5 837 delete bckTree; //Thanks to Mariana Bondila
838 delete tpcTree; //Thanks to Mariana Bondila
839
be5b9287 840 return 0;
73042f01 841}
842
843//_________________________________________________________________________
b9de75e1 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
9e1a0ddb 852 AliTPCClustersRow *clrow=((AliTPCtracker *) this)->fClustersArray.GetRow(sec,row);
b9de75e1 853 return (AliCluster*)(*clrow)[ncl];
854}
855
856//__________________________________________________________________________
857void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
858 //--------------------------------------------------------------------
859 //This function "cooks" a track label. If label<0, this track is fake.
860 //--------------------------------------------------------------------
861 Int_t noc=t->GetNumberOfClusters();
862 Int_t *lb=new Int_t[noc];
863 Int_t *mx=new Int_t[noc];
864 AliCluster **clusters=new AliCluster*[noc];
865
866 Int_t i;
867 for (i=0; i<noc; i++) {
868 lb[i]=mx[i]=0;
869 Int_t index=t->GetClusterIndex(i);
870 clusters[i]=GetCluster(index);
871 }
872
873 Int_t lab=123456789;
874 for (i=0; i<noc; i++) {
875 AliCluster *c=clusters[i];
876 lab=TMath::Abs(c->GetLabel(0));
877 Int_t j;
878 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
879 lb[j]=lab;
880 (mx[j])++;
881 }
882
883 Int_t max=0;
884 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
885
886 for (i=0; i<noc; i++) {
887 AliCluster *c=clusters[i];
888 if (TMath::Abs(c->GetLabel(1)) == lab ||
889 TMath::Abs(c->GetLabel(2)) == lab ) max++;
890 }
891
892 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
893
894 else {
895 Int_t tail=Int_t(0.10*noc);
896 max=0;
897 for (i=1; i<=tail; i++) {
898 AliCluster *c=clusters[noc-i];
899 if (lab == TMath::Abs(c->GetLabel(0)) ||
900 lab == TMath::Abs(c->GetLabel(1)) ||
901 lab == TMath::Abs(c->GetLabel(2))) max++;
902 }
903 if (max < Int_t(0.5*tail)) lab=-lab;
904 }
905
906 t->SetLabel(lab);
907
908 delete[] lb;
909 delete[] mx;
910 delete[] clusters;
911}
912
913//_________________________________________________________________________
914void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
915 //-----------------------------------------------------------------------
916 // Setup inner sector
917 //-----------------------------------------------------------------------
918 if (f==0) {
919 fAlpha=par->GetInnerAngle();
920 fAlphaShift=par->GetInnerAngleShift();
921 fPadPitchWidth=par->GetInnerPadPitchWidth();
922 fPadPitchLength=par->GetInnerPadPitchLength();
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
942//_________________________________________________________________________
73042f01 943void
944AliTPCtracker::AliTPCRow::InsertCluster(const AliTPCcluster* c, UInt_t index) {
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 }
951 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
952 Int_t i=Find(c->GetY());
953 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
954 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
955 fIndex[i]=index; fClusters[i]=c; fN++;
956}
957
958//___________________________________________________________________
959Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
960 //-----------------------------------------------------------------------
961 // Return the index of the nearest cluster
962 //-----------------------------------------------------------------------
963 if (y <= fClusters[0]->GetY()) return 0;
964 if (y > fClusters[fN-1]->GetY()) return fN;
965 Int_t b=0, e=fN-1, m=(b+e)/2;
966 for (; b<e; m=(b+e)/2) {
967 if (y > fClusters[m]->GetY()) b=m+1;
968 else e=m;
969 }
970 return m;
971}
972
973//_____________________________________________________________________________
974void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
975 //-----------------------------------------------------------------
976 // This funtion calculates dE/dX within the "low" and "up" cuts.
977 //-----------------------------------------------------------------
978 Int_t i;
be9c5115 979 Int_t nc=GetNumberOfClusters();
73042f01 980
981 Int_t swap;//stupid sorting
982 do {
983 swap=0;
be9c5115 984 for (i=0; i<nc-1; i++) {
d4cf1daa 985 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
986 Float_t tmp=fdEdxSample[i];
987 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
73042f01 988 swap++;
989 }
990 } while (swap);
991
be9c5115 992 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
73042f01 993 Float_t dedx=0;
d4cf1daa 994 for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
73042f01 995 dedx /= (nu-nl+1);
996 SetdEdx(dedx);
7f6ddf58 997
998 //Very rough PID
999 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1000
1001 if (p<0.6) {
1002 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
1003 if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
1004 SetMass(0.93827); return;
1005 }
1006
1007 if (p<1.2) {
1008 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
1009 SetMass(0.93827); return;
1010 }
1011
1012 SetMass(0.13957); return;
1013
73042f01 1014}
1015
73042f01 1016
1017