1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.27 2003/02/25 16:47:58 hristov
19 allow back propagation for more than 1 event (J.Chudoba)
21 Revision 1.26 2003/02/19 08:49:46 hristov
22 Track time measurement (S.Radomski)
24 Revision 1.25 2003/01/28 16:43:35 hristov
25 Additional protection: to be discussed with the Root team (M.Ivanov)
27 Revision 1.24 2002/11/19 16:13:24 hristov
28 stdlib.h included to declare exit() on HP
30 Revision 1.23 2002/11/19 11:50:08 hristov
31 Removing CONTAINERS (Yu.Belikov)
33 Revision 1.19 2002/07/19 07:31:40 kowal2
34 Improvement in tracking by J. Belikov
36 Revision 1.18 2002/05/13 07:33:52 kowal2
37 Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
38 in the case of defined region of interests.
40 Revision 1.17 2002/03/18 17:59:13 kowal2
41 Chnges in the pad geometry - 3 pad lengths introduced.
43 Revision 1.16 2001/11/08 16:39:03 hristov
44 Additional protection (M.Masera)
46 Revision 1.15 2001/11/08 16:36:33 hristov
47 Updated 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)
49 Revision 1.14 2001/10/21 19:04:55 hristov
50 Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
52 Revision 1.13 2001/07/23 12:01:30 hristov
55 Revision 1.12 2001/07/20 14:32:44 kowal2
56 Processing of many events possible now
58 Revision 1.11 2001/05/23 08:50:10 hristov
61 Revision 1.10 2001/05/16 14:57:25 alibrary
62 New files for folders and Stack
64 Revision 1.9 2001/05/11 07:16:56 hristov
65 Fix needed on Sun and Alpha
67 Revision 1.8 2001/05/08 15:00:15 hristov
68 Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
70 Revision 1.5 2000/12/20 07:51:59 kowal2
71 Changes suggested by Alessandra and Paolo to avoid overlapped
72 data fields in encapsulated classes.
74 Revision 1.4 2000/11/02 07:27:16 kowal2
77 Revision 1.2 2000/06/30 12:07:50 kowal2
78 Updated from the TPC-PreRelease branch
80 Revision 1.1.2.1 2000/06/25 08:53:55 kowal2
81 Splitted from AliTPCtracking
85 //-------------------------------------------------------
86 // Implementation of the TPC tracker
88 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
89 //-------------------------------------------------------
90 #include <TObjArray.h>
93 #include <Riostream.h>
95 #include "AliTPCtracker.h"
96 #include "AliTPCcluster.h"
97 #include "AliTPCParam.h"
98 #include "AliClusters.h"
100 #include "TVector2.h"
103 //_____________________________________________________________________________
104 AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
105 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
107 //---------------------------------------------------------------------
108 // The main TPC tracker constructor
109 //---------------------------------------------------------------------
110 fInnerSec=new AliTPCSector[fkNIS];
111 fOuterSec=new AliTPCSector[fkNOS];
114 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
115 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
117 fParam = (AliTPCParam*) par;
121 //_____________________________________________________________________________
122 AliTPCtracker::~AliTPCtracker() {
123 //------------------------------------------------------------------
124 // TPC tracker destructor
125 //------------------------------------------------------------------
134 //_____________________________________________________________________________
135 Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
138 // Parametrised error of the cluster reconstruction (pad direction)
141 const Float_t kArphi=0.41818e-2;
142 const Float_t kBrphi=0.17460e-4;
143 const Float_t kCrphi=0.30993e-2;
144 const Float_t kDrphi=0.41061e-3;
146 pt=TMath::Abs(pt)*1000.;
149 Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
150 if (s<0.4e-3) s=0.4e-3;
151 s*=1.3; //Iouri Belikov
156 //_____________________________________________________________________________
157 Double_t SigmaZ2(Double_t r, Double_t tgl)
160 // Parametrised error of the cluster reconstruction (drift direction)
163 const Float_t kAz=0.39614e-2;
164 const Float_t kBz=0.22443e-4;
165 const Float_t kCz=0.51504e-1;
169 Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
170 if (s<0.4e-3) s=0.4e-3;
171 s*=1.3; //Iouri Belikov
176 //_____________________________________________________________________________
177 Double_t f1(Double_t x1,Double_t y1,
178 Double_t x2,Double_t y2,
179 Double_t x3,Double_t y3)
181 //-----------------------------------------------------------------
182 // Initial approximation of the track curvature
183 //-----------------------------------------------------------------
184 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
185 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
186 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
187 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
188 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
190 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
192 return -xr*yr/sqrt(xr*xr+yr*yr);
196 //_____________________________________________________________________________
197 Double_t f2(Double_t x1,Double_t y1,
198 Double_t x2,Double_t y2,
199 Double_t x3,Double_t y3)
201 //-----------------------------------------------------------------
202 // Initial approximation of the track curvature times center of curvature
203 //-----------------------------------------------------------------
204 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
205 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
206 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
207 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
208 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
210 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
212 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
215 //_____________________________________________________________________________
216 Double_t f3(Double_t x1,Double_t y1,
217 Double_t x2,Double_t y2,
218 Double_t z1,Double_t z2)
220 //-----------------------------------------------------------------
221 // Initial approximation of the tangent of the track dip angle
222 //-----------------------------------------------------------------
223 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
226 //_____________________________________________________________________________
227 void AliTPCtracker::LoadClusters() {
228 //-----------------------------------------------------------------
229 // This function loads TPC clusters.
230 //-----------------------------------------------------------------
231 if (!gFile->IsOpen()) {
232 cerr<<"AliTPCtracker::LoadClusters : "<<
233 "file with clusters has not been open !\n";
238 sprintf(name,"TreeC_TPC_%d",GetEventNumber());
239 TTree *cTree=(TTree*)gFile->Get(name);
241 cerr<<"AliTPCtracker::LoadClusters : "<<
242 "can't get the tree with TPC clusters !\n";
246 TBranch *branch=cTree->GetBranch("Segment");
248 cerr<<"AliTPCtracker::LoadClusters : "<<
249 "can't get the segment branch !\n";
252 AliClusters carray, *addr=&carray;
253 carray.SetClass("AliTPCcluster");
255 branch->SetAddress(&addr);
257 Int_t nentr=(Int_t)cTree->GetEntries();
259 for (Int_t i=0; i<nentr; i++) {
262 Int_t ncl=carray.GetArray()->GetEntriesFast();
264 Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows();
265 Int_t id=carray.GetID();
266 if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) {
267 cerr<<"AliTPCtracker::LoadClusters : "<<
271 Int_t outindex = 2*fkNIS*nir;
274 Int_t row = id - sec*nir;
276 AliTPCRow &padrow=fInnerSec[sec][row];
278 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
279 padrow.InsertCluster(c,sec,row);
284 Int_t row = id - sec*nor;
286 AliTPCRow &padrow=fOuterSec[sec][row];
288 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
289 padrow.InsertCluster(c,sec+fkNIS,row);
293 carray.GetArray()->Clear();
298 //_____________________________________________________________________________
299 void AliTPCtracker::UnloadClusters() {
300 //-----------------------------------------------------------------
301 // This function unloads TPC clusters.
302 //-----------------------------------------------------------------
304 for (i=0; i<fkNIS; i++) {
305 Int_t nr=fInnerSec->GetNRows();
306 for (Int_t n=0; n<nr; n++) fInnerSec[i][n].ResetClusters();
308 for (i=0; i<fkNOS; i++) {
309 Int_t nr=fOuterSec->GetNRows();
310 for (Int_t n=0; n<nr; n++) fOuterSec[i][n].ResetClusters();
314 //_____________________________________________________________________________
315 Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
316 //-----------------------------------------------------------------
317 // This function tries to find a track prolongation.
318 //-----------------------------------------------------------------
319 Double_t xt=t.GetX();
320 const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
321 Int_t(0.5*fSectors->GetNRows());
322 Int_t tryAgain=kSKIP;
324 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
325 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
326 if (alpha < 0. ) alpha += 2.*TMath::Pi();
327 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
329 Int_t nrows=fSectors->GetRowNumber(xt)-1;
330 for (Int_t nr=nrows; nr>=rf; nr--) {
331 Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
332 if (!t.PropagateTo(x)) return 0;
336 Double_t maxchi2=kMaxCHI2;
337 const AliTPCRow &krow=fSectors[s][nr];
338 Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
339 Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
340 Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
341 Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
344 if (t.GetNumberOfClusters()>4)
345 cerr<<t.GetNumberOfClusters()
346 <<"FindProlongation warning: Too broad road !\n";
351 for (Int_t i=krow.Find(y-road); i<krow; i++) {
352 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
353 if (c->GetY() > y+road) break;
354 if (c->IsUsed()) continue;
355 if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
356 Double_t chi2=t.GetPredictedChi2(c);
357 if (chi2 > maxchi2) continue;
360 index=krow.GetIndex(i);
364 Float_t l=fSectors->GetPadPitchWidth();
365 Float_t corr=1.; if (nr>63) corr=0.67; // new (third) pad response !
366 t.SetSampledEdx(cl->GetQ()/l*corr,t.GetNumberOfClusters());
367 if (!t.Update(cl,maxchi2,index)) {
368 if (!tryAgain--) return 0;
369 } else tryAgain=kSKIP;
371 if (tryAgain==0) break;
374 if (!t.Rotate(fSectors->GetAlpha())) return 0;
375 } else if (y <-ymax) {
377 if (!t.Rotate(-fSectors->GetAlpha())) return 0;
385 //_____________________________________________________________________________
387 Int_t AliTPCtracker::FollowRefitInward(AliTPCseed *seed, AliTPCtrack *track) {
389 // This function propagates seed inward TPC using old clusters
392 // Sylwester Radomski, GSI
396 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
397 TVector2::Phi_0_2pi(alpha);
398 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
400 Int_t idx=-1, sec=-1, row=-1;
401 Int_t nc = seed->GetLabel(); // index to start with
403 idx=track->GetClusterIndex(nc);
404 sec=(idx&0xff000000)>>24;
405 row=(idx&0x00ff0000)>>16;
407 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
408 else { if (sec < fkNIS) row=-1; }
410 Int_t nr=fSectors->GetNRows();
411 for (Int_t i=0; i<nr; i++) {
413 Double_t x=fSectors->GetX(i);
414 if (!seed->PropagateTo(x)) return 0;
416 // Change sector and rotate seed if necessary
418 Double_t ymax=fSectors->GetMaxY(i);
419 Double_t y=seed->GetY();
423 if (!seed->Rotate(fSectors->GetAlpha())) return 0;
424 } else if (y <-ymax) {
426 if (!seed->Rotate(-fSectors->GetAlpha())) return 0;
429 // try to find a cluster
433 Double_t maxchi2 = kMaxCHI2;
435 //cout << i << " " << row << " " << nc << endl;
439 // accept already found cluster
440 AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
441 Double_t chi2 = seed->GetPredictedChi2(c);
447 if (++nc < track->GetNumberOfClusters() ) {
448 idx=track->GetClusterIndex(nc);
449 sec=(idx&0xff000000)>>24;
450 row=(idx&0x00ff0000)>>16;
453 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
454 else { if (sec < fkNIS) row=-1; }
460 //cout << "Assigned" << endl;
461 Float_t l=fSectors->GetPadPitchWidth();
462 Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
463 seed->SetSampledEdx(cl->GetQ()/l*corr,seed->GetNumberOfClusters());
464 seed->Update(cl,maxchi2,index);
472 //_____________________________________________________________________________
473 Int_t AliTPCtracker::FollowBackProlongation
474 (AliTPCseed& seed, const AliTPCtrack &track) {
475 //-----------------------------------------------------------------
476 // This function propagates tracks back through the TPC
477 //-----------------------------------------------------------------
478 Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
479 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
480 if (alpha < 0. ) alpha += 2.*TMath::Pi();
481 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
483 Int_t idx=-1, sec=-1, row=-1;
484 Int_t nc=seed.GetLabel(); //index of the cluster to start with
486 idx=track.GetClusterIndex(nc);
487 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
489 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
490 else { if (sec < fkNIS) row=-1; }
492 Int_t nr=fSectors->GetNRows();
493 for (Int_t i=0; i<nr; i++) {
494 Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
496 if (!seed.PropagateTo(x)) return 0;
498 Double_t y=seed.GetY();
501 if (!seed.Rotate(fSectors->GetAlpha())) return 0;
502 } else if (y <-ymax) {
504 if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
509 Double_t maxchi2=kMaxCHI2;
510 Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
511 Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
512 Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
513 Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
515 cerr<<seed.GetNumberOfClusters()
516 <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n";
521 Int_t accepted=seed.GetNumberOfClusters();
523 //try to accept already found cluster
524 AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
526 if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
527 index=idx; cl=c; maxchi2=chi2;
528 } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
531 idx=track.GetClusterIndex(nc);
532 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
534 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
535 else { if (sec < fkNIS) row=-1; }
539 //try to fill the gap
540 const AliTPCRow &krow=fSectors[s][i];
543 for (Int_t i=krow.Find(y-road); i<krow; i++) {
544 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
545 if (c->GetY() > y+road) break;
546 if (c->IsUsed()) continue;
547 if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
548 Double_t chi2=seed.GetPredictedChi2(c);
549 if (chi2 > maxchi2) continue;
552 index=krow.GetIndex(i);
558 Float_t l=fSectors->GetPadPitchWidth();
559 Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
560 seed.SetSampledEdx(cl->GetQ()/l*corr,seed.GetNumberOfClusters());
561 seed.Update(cl,maxchi2,index);
571 //_____________________________________________________________________________
572 void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
573 //-----------------------------------------------------------------
574 // This function creates track seeds.
575 //-----------------------------------------------------------------
576 Double_t x[5], c[15];
578 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
579 Double_t cs=cos(alpha), sn=sin(alpha);
581 Double_t x1 =fSectors->GetX(i1);
582 Double_t xx2=fSectors->GetX(i2);
584 for (Int_t ns=0; ns<fN; ns++) {
585 Int_t nl=fSectors[(ns-1+fN)%fN][i2];
586 Int_t nm=fSectors[ns][i2];
587 Int_t nu=fSectors[(ns+1)%fN][i2];
588 const AliTPCRow& kr1=fSectors[ns][i1];
589 for (Int_t is=0; is < kr1; is++) {
590 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
591 for (Int_t js=0; js < nl+nm+nu; js++) {
592 const AliTPCcluster *kcl;
594 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
597 const AliTPCRow& kr2=fSectors[(ns-1+fN)%fN][i2];
599 y2=kcl->GetY(); z2=kcl->GetZ();
604 const AliTPCRow& kr2=fSectors[ns][i2];
606 x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
608 const AliTPCRow& kr2=fSectors[(ns+1)%fN][i2];
610 y2=kcl->GetY(); z2=kcl->GetZ();
615 Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
616 if (TMath::Abs(zz-z2)>5.) continue;
618 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
619 if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
623 x[4]=f1(x1,y1,x2,y2,x3,y3);
624 if (TMath::Abs(x[4]) >= 0.0066) continue;
625 x[2]=f2(x1,y1,x2,y2,x3,y3);
626 //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
627 x[3]=f3(x1,y1,x2,y2,z1,z2);
628 if (TMath::Abs(x[3]) > 1.2) continue;
629 Double_t a=asin(x[2]);
630 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
631 if (TMath::Abs(zv-z3)>10.) continue;
633 Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
634 Double_t sy2=kcl->GetSigmaY2(), sz2=kcl->GetSigmaZ2();
635 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
636 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
638 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
639 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
640 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
641 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
642 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
643 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
644 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
645 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
646 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
647 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
651 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
652 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
653 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
654 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
655 c[13]=f30*sy1*f40+f32*sy2*f42;
656 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
658 UInt_t index=kr1.GetIndex(is);
659 AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
660 Float_t l=fSectors->GetPadPitchWidth();
661 track->SetSampledEdx(kr1[is]->GetQ()/l,0);
663 Int_t rc=FollowProlongation(*track, i2);
664 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
665 else fSeeds->AddLast(track);
671 //_____________________________________________________________________________
672 Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
673 //-----------------------------------------------------------------
674 // This function reades track seeds.
675 //-----------------------------------------------------------------
676 TDirectory *savedir=gDirectory;
678 TFile *in=(TFile*)inp;
680 cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
685 TTree *seedTree=(TTree*)in->Get("Seeds");
687 cerr<<"AliTPCtracker::ReadSeeds(): ";
688 cerr<<"can't get a tree with track seeds !\n";
691 AliTPCtrack *seed=new AliTPCtrack;
692 seedTree->SetBranchAddress("tracks",&seed);
694 if (fSeeds==0) fSeeds=new TObjArray(15000);
696 Int_t n=(Int_t)seedTree->GetEntries();
697 for (Int_t i=0; i<n; i++) {
698 seedTree->GetEvent(i);
699 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
704 delete seedTree; //Thanks to Mariana Bondila
711 //_____________________________________________________________________________
712 Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
713 //-----------------------------------------------------------------
714 // This is a track finder.
715 //-----------------------------------------------------------------
716 TDirectory *savedir=gDirectory;
719 TFile *in=(TFile*)inp;
721 cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
726 if (!out->IsOpen()) {
727 cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
736 sprintf(tname,"TreeT_TPC_%d",GetEventNumber());
737 TTree tracktree(tname,"Tree with TPC tracks");
738 AliTPCtrack *iotrack=0;
739 tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,3);
742 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
743 Int_t nrows=nlow+nup;
745 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
746 fSectors=fOuterSec; fN=fkNOS;
747 fSeeds=new TObjArray(15000);
748 MakeSeeds(nup-1, nup-1-gap);
749 MakeSeeds(nup-1-shift, nup-1-shift-gap);
754 Int_t nseed=fSeeds->GetEntriesFast();
755 for (Int_t i=0; i<nseed; i++) {
756 //tracking in the outer sectors
757 fSectors=fOuterSec; fN=fkNOS;
759 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
760 if (!FollowProlongation(t)) {
761 delete fSeeds->RemoveAt(i);
765 //tracking in the inner sectors
766 fSectors=fInnerSec; fN=fkNIS;
768 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
769 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
770 if (alpha < 0. ) alpha += 2.*TMath::Pi();
771 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
773 alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
775 if (t.Rotate(alpha)) {
776 if (FollowProlongation(t)) {
777 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
779 CookLabel(pt,0.1); //For comparison only
787 delete fSeeds->RemoveAt(i);
792 cerr<<"Number of found tracks : "<<found<<endl;
797 fSeeds->Clear(); delete fSeeds; fSeeds=0;
801 //_____________________________________________________________________________
803 Int_t AliTPCtracker::RefitInward(TFile *in, TFile *out) {
805 // The function propagates tracks throught TPC inward
806 // using already associated clusters.
808 // in - file with back propagated tracks
809 // outs - file with inward propagation
811 // Sylwester Radomski, GSI
817 cout << "Input File not open !\n" << endl;
821 if (!out->IsOpen()) {
822 cout << "Output File not open !\n" << endl;
826 TDirectory *savedir = gDirectory;
829 // Connect to input tree
831 TTree *inputTree = (TTree*)in->Get("seedsTPCtoTRD_0");
832 TBranch *inBranch = inputTree->GetBranch("tracks");
833 AliTPCtrack *inTrack = 0;
834 inBranch->SetAddress(&inTrack);
836 Int_t nTracks = (Int_t)inputTree->GetEntries();
838 // Create output tree
840 AliTPCtrack *outTrack = new AliTPCtrack();
841 TTree *outputTree = new TTree("tracksTPC_0","Refited TPC tracks");
842 outputTree->Branch("tracks", "AliTPCtrack", &outTrack, 32000, 3);
846 for(Int_t t=0; t < nTracks; t++) {
850 inputTree->GetEvent(t);
851 AliTPCseed *seed = new AliTPCseed(*inTrack, inTrack->GetAlpha());
853 seed->ResetCovariance();
856 fSectors = fInnerSec;
857 Int_t res = FollowRefitInward(seed, inTrack);
859 Int_t nc = seed->GetNumberOfClusters();
861 fSectors = fOuterSec;
862 res = FollowRefitInward(seed, inTrack);
863 UseClusters(seed, nc);
866 seed->PropagateTo(fParam->GetInnerRadiusLow());
867 outTrack = (AliTPCtrack*)seed;
868 outTrack->SetLabel(inTrack->GetLabel());
876 cout << "Refitted " << nRefited << " tracks." << endl;
880 if (inputTree) delete inputTree;
881 if (outputTree) delete outputTree;
889 //_____________________________________________________________________________
890 Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
891 //-----------------------------------------------------------------
892 // This function propagates tracks back through the TPC.
893 //-----------------------------------------------------------------
894 fSeeds=new TObjArray(15000);
896 TFile *in=(TFile*)inp;
897 TDirectory *savedir=gDirectory;
900 cerr<<"AliTPCtracker::PropagateBack(): ";
901 cerr<<"file with back propagated ITS tracks is not open !\n";
905 if (!out->IsOpen()) {
906 cerr<<"AliTPCtracker::PropagateBack(): ";
907 cerr<<"file for back propagated TPC tracks is not open !\n";
915 sprintf(tName,"TreeT_ITSb_%d",GetEventNumber());
916 TTree *bckTree=(TTree*)in->Get(tName);
918 cerr<<"AliTPCtracker::PropagateBack() ";
919 cerr<<"can't get a tree with back propagated ITS tracks !\n";
922 AliTPCtrack *bckTrack=new AliTPCtrack;
923 bckTree->SetBranchAddress("tracks",&bckTrack);
925 sprintf(tName,"TreeT_TPC_%d",GetEventNumber());
926 TTree *tpcTree=(TTree*)in->Get(tName);
928 cerr<<"AliTPCtracker::PropagateBack() ";
929 cerr<<"can't get a tree with TPC tracks !\n";
932 AliTPCtrack *tpcTrack=new AliTPCtrack;
933 tpcTree->SetBranchAddress("tracks",&tpcTrack);
935 //*** Prepare an array of tracks to be back propagated
936 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
937 Int_t nrows=nlow+nup;
940 // Match ITS tracks with old TPC tracks
941 // the tracks do not have to be sorted [SR, GSI, 18.02.2003]
943 // the algorithm is linear and uses LUT for sorting
944 // cost of the algorithm: 2 * number of tracks
947 TObjArray tracks(15000);
949 Int_t tpcN= (Int_t)tpcTree->GetEntries();
950 Int_t bckN= (bckTree)? (Int_t)bckTree->GetEntries() : 0;
953 const Int_t nLab = 20000; // limit on number of primaries (arbitrary)
955 for(Int_t i=0; i<nLab; i++) lut[i] = -1;
958 for(Int_t i=0; i<bckN; i++) {
959 bckTree->GetEvent(i);
960 Int_t lab = TMath::Abs(bckTrack->GetLabel());
961 if (lab < nLab) lut[lab] = i;
965 for (Int_t i=0; i<tpcN; i++) {
966 tpcTree->GetEvent(i);
967 Double_t alpha=tpcTrack->GetAlpha();
971 // No ITS - use TPC track only
973 fSeeds->AddLast(new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha()));
974 tracks.AddLast(new AliTPCtrack(*tpcTrack));
979 // discard not prolongated tracks (to be discussed)
981 Int_t lab = TMath::Abs(tpcTrack->GetLabel());
982 if (lab > nLab || lut[lab] < 0) continue;
984 bckTree->GetEvent(lut[lab]);
985 bckTrack->Rotate(alpha-bckTrack->GetAlpha());
987 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
988 tracks.AddLast(new AliTPCtrack(*tpcTrack));
995 TObjArray tracks(15000);
997 Int_t tpcN=(Int_t)tpcTree->GetEntries();
998 Int_t bckN=(Int_t)bckTree->GetEntries();
999 if (j<bckN) bckTree->GetEvent(j++);
1000 for (i=0; i<tpcN; i++) {
1001 tpcTree->GetEvent(i);
1002 Double_t alpha=tpcTrack->GetAlpha();
1004 TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
1005 if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
1006 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
1007 bckTree->GetEvent(j++);
1009 tpcTrack->ResetCovariance();
1010 fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
1012 tracks.AddLast(new AliTPCtrack(*tpcTrack));
1018 // tree name seedsTPCtoTRD as expected by TRD and as
1019 // discussed and decided in Strasbourg (May 2002)
1020 // [SR, GSI, 18.02.2003]
1022 sprintf(tName,"seedsTPCtoTRD_%d",GetEventNumber());
1023 TTree backTree(tName,"Tree with back propagated TPC tracks");
1025 AliTPCtrack *otrack=0;
1026 backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
1028 //*** Back propagation through inner sectors
1029 fSectors=fInnerSec; fN=fkNIS;
1030 Int_t nseed=fSeeds->GetEntriesFast();
1031 for (i=0; i<nseed; i++) {
1032 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
1033 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
1035 Int_t nc=t.GetNumberOfClusters();
1036 s.SetLabel(nc-1); //set number of the cluster to start with
1038 if (FollowBackProlongation(s,t)) {
1042 delete fSeeds->RemoveAt(i);
1045 //*** Back propagation through outer sectors
1046 fSectors=fOuterSec; fN=fkNOS;
1048 for (i=0; i<nseed; i++) {
1049 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
1051 Int_t nc=s.GetNumberOfClusters();
1052 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
1054 Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1055 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1056 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1057 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1059 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1060 alpha-=s.GetAlpha();
1062 if (s.Rotate(alpha)) {
1063 if (FollowBackProlongation(s,t)) {
1064 if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
1066 s.SetLabel(t.GetLabel());
1069 // Propagate to outer reference plane for comparison reasons
1070 // reason for keeping fParam object [SR, GSI, 18.02.2003]
1072 ps->PropagateTo(fParam->GetOuterRadiusUp());
1075 cerr<<found++<<'\r';
1080 delete fSeeds->RemoveAt(i);
1085 cerr<<"Number of seeds: "<<nseed<<endl;
1086 cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
1087 cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
1092 delete bckTree; //Thanks to Mariana Bondila
1093 delete tpcTree; //Thanks to Mariana Bondila
1100 //_________________________________________________________________________
1101 AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
1102 //--------------------------------------------------------------------
1103 // Return pointer to a given cluster
1104 //--------------------------------------------------------------------
1105 Int_t sec=(index&0xff000000)>>24;
1106 Int_t row=(index&0x00ff0000)>>16;
1107 Int_t ncl=(index&0x0000ffff)>>00;
1109 const AliTPCcluster *cl=0;
1112 cl=fInnerSec[sec][row].GetUnsortedCluster(ncl);
1115 cl=fOuterSec[sec][row].GetUnsortedCluster(ncl);
1118 return (AliCluster*)cl;
1121 //__________________________________________________________________________
1122 void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1123 //--------------------------------------------------------------------
1124 //This function "cooks" a track label. If label<0, this track is fake.
1125 //--------------------------------------------------------------------
1126 Int_t noc=t->GetNumberOfClusters();
1127 Int_t *lb=new Int_t[noc];
1128 Int_t *mx=new Int_t[noc];
1129 AliCluster **clusters=new AliCluster*[noc];
1132 for (i=0; i<noc; i++) {
1134 Int_t index=t->GetClusterIndex(i);
1135 clusters[i]=GetCluster(index);
1138 Int_t lab=123456789;
1139 for (i=0; i<noc; i++) {
1140 AliCluster *c=clusters[i];
1141 lab=TMath::Abs(c->GetLabel(0));
1143 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
1149 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
1151 for (i=0; i<noc; i++) {
1152 AliCluster *c=clusters[i];
1153 if (TMath::Abs(c->GetLabel(1)) == lab ||
1154 TMath::Abs(c->GetLabel(2)) == lab ) max++;
1157 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
1160 Int_t tail=Int_t(0.10*noc);
1162 for (i=1; i<=tail; i++) {
1163 AliCluster *c=clusters[noc-i];
1164 if (lab == TMath::Abs(c->GetLabel(0)) ||
1165 lab == TMath::Abs(c->GetLabel(1)) ||
1166 lab == TMath::Abs(c->GetLabel(2))) max++;
1168 if (max < Int_t(0.5*tail)) lab=-lab;
1178 //_________________________________________________________________________
1179 void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
1180 //-----------------------------------------------------------------------
1181 // Setup inner sector
1182 //-----------------------------------------------------------------------
1184 fAlpha=par->GetInnerAngle();
1185 fAlphaShift=par->GetInnerAngleShift();
1186 fPadPitchWidth=par->GetInnerPadPitchWidth();
1187 f1PadPitchLength=par->GetInnerPadPitchLength();
1188 f2PadPitchLength=f1PadPitchLength;
1190 fN=par->GetNRowLow();
1191 fRow=new AliTPCRow[fN];
1192 for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
1194 fAlpha=par->GetOuterAngle();
1195 fAlphaShift=par->GetOuterAngleShift();
1196 fPadPitchWidth=par->GetOuterPadPitchWidth();
1197 f1PadPitchLength=par->GetOuter1PadPitchLength();
1198 f2PadPitchLength=par->GetOuter2PadPitchLength();
1200 fN=par->GetNRowUp();
1201 fRow=new AliTPCRow[fN];
1202 for (Int_t i=0; i<fN; i++){
1203 fRow[i].SetX(par->GetPadRowRadiiUp(i));
1208 //_________________________________________________________________________
1209 void AliTPCtracker::
1210 AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) {
1211 //-----------------------------------------------------------------------
1212 // Insert a cluster into this pad row in accordence with its y-coordinate
1213 //-----------------------------------------------------------------------
1214 if (fN==kMaxClusterPerRow) {
1215 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
1218 Int_t index=(((sec<<8)+row)<<16)+fN;
1222 fClusterArray[0]=*c; fClusters[fN++]=fClusterArray;
1228 AliTPCcluster *buff=new AliTPCcluster[size];
1229 memcpy(buff,fClusterArray,fSize*sizeof(AliTPCcluster));
1230 for (Int_t i=0; i<fN; i++)
1231 fClusters[i]=buff+(fClusters[i]-fClusterArray);
1232 delete[] fClusterArray;
1237 Int_t i=Find(c->GetY());
1238 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
1239 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
1241 fClusters[i]=fClusterArray+fN; fClusterArray[fN++]=*c;
1244 //___________________________________________________________________
1245 Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
1246 //-----------------------------------------------------------------------
1247 // Return the index of the nearest cluster
1248 //-----------------------------------------------------------------------
1250 if (y <= fClusters[0]->GetY()) return 0;
1251 if (y > fClusters[fN-1]->GetY()) return fN;
1252 Int_t b=0, e=fN-1, m=(b+e)/2;
1253 for (; b<e; m=(b+e)/2) {
1254 if (y > fClusters[m]->GetY()) b=m+1;
1260 //_____________________________________________________________________________
1261 void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1262 //-----------------------------------------------------------------
1263 // This funtion calculates dE/dX within the "low" and "up" cuts.
1264 //-----------------------------------------------------------------
1266 Int_t nc=GetNumberOfClusters();
1268 Int_t swap;//stupid sorting
1271 for (i=0; i<nc-1; i++) {
1272 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
1273 Float_t tmp=fdEdxSample[i];
1274 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
1279 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
1281 for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
1286 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1288 Double_t log1=TMath::Log(p+0.45), log2=TMath::Log(p+0.12);
1290 if (dedx < 34 + 30/(p+0.45)/(p+0.45) + 24*log1) {SetMass(0.13957); return;}
1291 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.49368); return;}
1292 SetMass(0.93827); return;
1296 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.13957); return;}
1297 SetMass(0.93827); return;
1300 SetMass(0.13957); return;