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.30 2003/02/28 16:13:32 hristov
21 Revision 1.29 2003/02/28 15:18:16 hristov
22 Corrections suggested by J.Chudoba
24 Revision 1.28 2003/02/27 16:15:52 hristov
25 Code for inward refitting (S.Radomski)
27 Revision 1.27 2003/02/25 16:47:58 hristov
28 allow back propagation for more than 1 event (J.Chudoba)
30 Revision 1.26 2003/02/19 08:49:46 hristov
31 Track time measurement (S.Radomski)
33 Revision 1.25 2003/01/28 16:43:35 hristov
34 Additional protection: to be discussed with the Root team (M.Ivanov)
36 Revision 1.24 2002/11/19 16:13:24 hristov
37 stdlib.h included to declare exit() on HP
39 Revision 1.23 2002/11/19 11:50:08 hristov
40 Removing CONTAINERS (Yu.Belikov)
42 Revision 1.19 2002/07/19 07:31:40 kowal2
43 Improvement in tracking by J. Belikov
45 Revision 1.18 2002/05/13 07:33:52 kowal2
46 Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
47 in the case of defined region of interests.
49 Revision 1.17 2002/03/18 17:59:13 kowal2
50 Chnges in the pad geometry - 3 pad lengths introduced.
52 Revision 1.16 2001/11/08 16:39:03 hristov
53 Additional protection (M.Masera)
55 Revision 1.15 2001/11/08 16:36:33 hristov
56 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)
58 Revision 1.14 2001/10/21 19:04:55 hristov
59 Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
61 Revision 1.13 2001/07/23 12:01:30 hristov
64 Revision 1.12 2001/07/20 14:32:44 kowal2
65 Processing of many events possible now
67 Revision 1.11 2001/05/23 08:50:10 hristov
70 Revision 1.10 2001/05/16 14:57:25 alibrary
71 New files for folders and Stack
73 Revision 1.9 2001/05/11 07:16:56 hristov
74 Fix needed on Sun and Alpha
76 Revision 1.8 2001/05/08 15:00:15 hristov
77 Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
79 Revision 1.5 2000/12/20 07:51:59 kowal2
80 Changes suggested by Alessandra and Paolo to avoid overlapped
81 data fields in encapsulated classes.
83 Revision 1.4 2000/11/02 07:27:16 kowal2
86 Revision 1.2 2000/06/30 12:07:50 kowal2
87 Updated from the TPC-PreRelease branch
89 Revision 1.1.2.1 2000/06/25 08:53:55 kowal2
90 Splitted from AliTPCtracking
94 //-------------------------------------------------------
95 // Implementation of the TPC tracker
97 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
98 //-------------------------------------------------------
99 #include <TObjArray.h>
102 #include <Riostream.h>
104 #include "AliTPCtracker.h"
105 #include "AliTPCcluster.h"
106 #include "AliTPCParam.h"
107 #include "AliClusters.h"
109 #include "TVector2.h"
112 //_____________________________________________________________________________
113 AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
114 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
116 //---------------------------------------------------------------------
117 // The main TPC tracker constructor
118 //---------------------------------------------------------------------
119 fInnerSec=new AliTPCSector[fkNIS];
120 fOuterSec=new AliTPCSector[fkNOS];
123 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
124 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
126 fParam = (AliTPCParam*) par;
130 //_____________________________________________________________________________
131 AliTPCtracker::~AliTPCtracker() {
132 //------------------------------------------------------------------
133 // TPC tracker destructor
134 //------------------------------------------------------------------
143 //_____________________________________________________________________________
144 Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
147 // Parametrised error of the cluster reconstruction (pad direction)
150 const Float_t kArphi=0.41818e-2;
151 const Float_t kBrphi=0.17460e-4;
152 const Float_t kCrphi=0.30993e-2;
153 const Float_t kDrphi=0.41061e-3;
155 pt=TMath::Abs(pt)*1000.;
158 Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
159 if (s<0.4e-3) s=0.4e-3;
160 s*=1.3; //Iouri Belikov
165 //_____________________________________________________________________________
166 Double_t SigmaZ2(Double_t r, Double_t tgl)
169 // Parametrised error of the cluster reconstruction (drift direction)
172 const Float_t kAz=0.39614e-2;
173 const Float_t kBz=0.22443e-4;
174 const Float_t kCz=0.51504e-1;
178 Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
179 if (s<0.4e-3) s=0.4e-3;
180 s*=1.3; //Iouri Belikov
185 //_____________________________________________________________________________
186 Double_t f1(Double_t x1,Double_t y1,
187 Double_t x2,Double_t y2,
188 Double_t x3,Double_t y3)
190 //-----------------------------------------------------------------
191 // Initial approximation of the track curvature
192 //-----------------------------------------------------------------
193 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
194 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
195 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
196 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
197 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
199 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
201 return -xr*yr/sqrt(xr*xr+yr*yr);
205 //_____________________________________________________________________________
206 Double_t f2(Double_t x1,Double_t y1,
207 Double_t x2,Double_t y2,
208 Double_t x3,Double_t y3)
210 //-----------------------------------------------------------------
211 // Initial approximation of the track curvature times center of curvature
212 //-----------------------------------------------------------------
213 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
214 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
215 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
216 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
217 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
219 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
221 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
224 //_____________________________________________________________________________
225 Double_t f3(Double_t x1,Double_t y1,
226 Double_t x2,Double_t y2,
227 Double_t z1,Double_t z2)
229 //-----------------------------------------------------------------
230 // Initial approximation of the tangent of the track dip angle
231 //-----------------------------------------------------------------
232 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
235 //_____________________________________________________________________________
236 Int_t AliTPCtracker::LoadClusters() {
237 //-----------------------------------------------------------------
238 // This function loads TPC clusters.
239 //-----------------------------------------------------------------
240 if (!gFile->IsOpen()) {
241 cerr<<"AliTPCtracker::LoadClusters : "<<
242 "file with clusters has not been open !\n";
247 sprintf(name,"TreeC_TPC_%d",GetEventNumber());
248 TTree *cTree=(TTree*)gFile->Get(name);
250 cerr<<"AliTPCtracker::LoadClusters : "<<
251 "can't get the tree with TPC clusters !\n";
255 TBranch *branch=cTree->GetBranch("Segment");
257 cerr<<"AliTPCtracker::LoadClusters : "<<
258 "can't get the segment branch !\n";
261 // AliClusters carray, *addr=&carray;
262 AliClusters carray, *addr=&carray;
263 carray.SetClass("AliTPCcluster");
265 branch->SetAddress(&addr);
267 Int_t nentr=(Int_t)cTree->GetEntries();
269 for (Int_t i=0; i<nentr; i++) {
272 Int_t ncl=carray.GetArray()->GetEntriesFast();
274 Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows();
275 Int_t id=carray.GetID();
276 if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) {
277 cerr<<"AliTPCtracker::LoadClusters : "<<
281 Int_t outindex = 2*fkNIS*nir;
284 Int_t row = id - sec*nir;
286 AliTPCRow &padrow=fInnerSec[sec][row];
288 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
289 padrow.InsertCluster(c,sec,row);
294 Int_t row = id - sec*nor;
296 AliTPCRow &padrow=fOuterSec[sec][row];
298 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
299 padrow.InsertCluster(c,sec+fkNIS,row);
303 carray.GetArray()->Clear();
309 //_____________________________________________________________________________
310 void AliTPCtracker::UnloadClusters() {
311 //-----------------------------------------------------------------
312 // This function unloads TPC clusters.
313 //-----------------------------------------------------------------
315 for (i=0; i<fkNIS; i++) {
316 Int_t nr=fInnerSec->GetNRows();
317 for (Int_t n=0; n<nr; n++) fInnerSec[i][n].ResetClusters();
319 for (i=0; i<fkNOS; i++) {
320 Int_t nr=fOuterSec->GetNRows();
321 for (Int_t n=0; n<nr; n++) fOuterSec[i][n].ResetClusters();
325 //_____________________________________________________________________________
326 Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
327 //-----------------------------------------------------------------
328 // This function tries to find a track prolongation.
329 //-----------------------------------------------------------------
330 Double_t xt=t.GetX();
331 const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
332 Int_t(0.5*fSectors->GetNRows());
333 Int_t tryAgain=kSKIP;
335 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
336 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
337 if (alpha < 0. ) alpha += 2.*TMath::Pi();
338 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
340 Int_t nrows=fSectors->GetRowNumber(xt)-1;
341 for (Int_t nr=nrows; nr>=rf; nr--) {
342 Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
343 if (!t.PropagateTo(x)) return 0;
347 Double_t maxchi2=kMaxCHI2;
348 const AliTPCRow &krow=fSectors[s][nr];
349 Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
350 Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
351 Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
352 Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
355 if (t.GetNumberOfClusters()>4)
356 cerr<<t.GetNumberOfClusters()
357 <<"FindProlongation warning: Too broad road !\n";
362 for (Int_t i=krow.Find(y-road); i<krow; i++) {
363 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
364 if (c->GetY() > y+road) break;
365 if (c->IsUsed()) continue;
366 if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
367 Double_t chi2=t.GetPredictedChi2(c);
368 if (chi2 > maxchi2) continue;
371 index=krow.GetIndex(i);
375 Float_t l=fSectors->GetPadPitchWidth();
376 Float_t corr=1.; if (nr>63) corr=0.67; // new (third) pad response !
377 t.SetSampledEdx(cl->GetQ()/l*corr,t.GetNumberOfClusters());
378 if (!t.Update(cl,maxchi2,index)) {
379 if (!tryAgain--) return 0;
380 } else tryAgain=kSKIP;
382 if (tryAgain==0) break;
385 if (!t.Rotate(fSectors->GetAlpha())) return 0;
386 } else if (y <-ymax) {
388 if (!t.Rotate(-fSectors->GetAlpha())) return 0;
396 //_____________________________________________________________________________
398 Int_t AliTPCtracker::FollowRefitInward(AliTPCseed *seed, AliTPCtrack *track) {
400 // This function propagates seed inward TPC using old clusters
403 // Sylwester Radomski, GSI
407 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
408 TVector2::Phi_0_2pi(alpha);
409 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
411 Int_t idx=-1, sec=-1, row=-1;
412 Int_t nc = seed->GetLabel(); // index to start with
414 idx=track->GetClusterIndex(nc);
415 sec=(idx&0xff000000)>>24;
416 row=(idx&0x00ff0000)>>16;
418 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
419 else { if (sec < fkNIS) row=-1; }
421 Int_t nr=fSectors->GetNRows();
422 for (Int_t i=0; i<nr; i++) {
424 Double_t x=fSectors->GetX(i);
425 if (!seed->PropagateTo(x)) return 0;
427 // Change sector and rotate seed if necessary
429 Double_t ymax=fSectors->GetMaxY(i);
430 Double_t y=seed->GetY();
434 if (!seed->Rotate(fSectors->GetAlpha())) return 0;
435 } else if (y <-ymax) {
437 if (!seed->Rotate(-fSectors->GetAlpha())) return 0;
440 // try to find a cluster
444 Double_t maxchi2 = kMaxCHI2;
446 //cout << i << " " << row << " " << nc << endl;
450 // accept already found cluster
451 AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
452 Double_t chi2 = seed->GetPredictedChi2(c);
458 if (++nc < track->GetNumberOfClusters() ) {
459 idx=track->GetClusterIndex(nc);
460 sec=(idx&0xff000000)>>24;
461 row=(idx&0x00ff0000)>>16;
464 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
465 else { if (sec < fkNIS) row=-1; }
471 //cout << "Assigned" << endl;
472 Float_t l=fSectors->GetPadPitchWidth();
473 Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
474 seed->SetSampledEdx(cl->GetQ()/l*corr,seed->GetNumberOfClusters());
475 seed->Update(cl,maxchi2,index);
483 //_____________________________________________________________________________
484 Int_t AliTPCtracker::FollowBackProlongation
485 (AliTPCseed& seed, const AliTPCtrack &track) {
486 //-----------------------------------------------------------------
487 // This function propagates tracks back through the TPC
488 //-----------------------------------------------------------------
489 Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
490 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
491 if (alpha < 0. ) alpha += 2.*TMath::Pi();
492 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
494 Int_t idx=-1, sec=-1, row=-1;
495 Int_t nc=seed.GetLabel(); //index of the cluster to start with
497 idx=track.GetClusterIndex(nc);
498 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
500 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
501 else { if (sec < fkNIS) row=-1; }
503 Int_t nr=fSectors->GetNRows();
504 for (Int_t i=0; i<nr; i++) {
505 Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
507 if (!seed.PropagateTo(x)) return 0;
509 Double_t y=seed.GetY();
512 if (!seed.Rotate(fSectors->GetAlpha())) return 0;
513 } else if (y <-ymax) {
515 if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
520 Double_t maxchi2=kMaxCHI2;
521 Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
522 Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
523 Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
524 Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
526 cerr<<seed.GetNumberOfClusters()
527 <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n";
532 Int_t accepted=seed.GetNumberOfClusters();
534 //try to accept already found cluster
535 AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
537 if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
538 index=idx; cl=c; maxchi2=chi2;
539 } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
542 idx=track.GetClusterIndex(nc);
543 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
545 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
546 else { if (sec < fkNIS) row=-1; }
550 //try to fill the gap
551 const AliTPCRow &krow=fSectors[s][i];
554 for (Int_t i=krow.Find(y-road); i<krow; i++) {
555 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
556 if (c->GetY() > y+road) break;
557 if (c->IsUsed()) continue;
558 if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
559 Double_t chi2=seed.GetPredictedChi2(c);
560 if (chi2 > maxchi2) continue;
563 index=krow.GetIndex(i);
569 Float_t l=fSectors->GetPadPitchWidth();
570 Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
571 seed.SetSampledEdx(cl->GetQ()/l*corr,seed.GetNumberOfClusters());
572 seed.Update(cl,maxchi2,index);
582 //_____________________________________________________________________________
583 void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
584 //-----------------------------------------------------------------
585 // This function creates track seeds.
586 //-----------------------------------------------------------------
587 Double_t x[5], c[15];
589 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
590 Double_t cs=cos(alpha), sn=sin(alpha);
592 Double_t x1 =fSectors->GetX(i1);
593 Double_t xx2=fSectors->GetX(i2);
595 for (Int_t ns=0; ns<fN; ns++) {
596 Int_t nl=fSectors[(ns-1+fN)%fN][i2];
597 Int_t nm=fSectors[ns][i2];
598 Int_t nu=fSectors[(ns+1)%fN][i2];
599 const AliTPCRow& kr1=fSectors[ns][i1];
600 for (Int_t is=0; is < kr1; is++) {
601 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
602 for (Int_t js=0; js < nl+nm+nu; js++) {
603 const AliTPCcluster *kcl;
605 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
608 const AliTPCRow& kr2=fSectors[(ns-1+fN)%fN][i2];
610 y2=kcl->GetY(); z2=kcl->GetZ();
615 const AliTPCRow& kr2=fSectors[ns][i2];
617 x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
619 const AliTPCRow& kr2=fSectors[(ns+1)%fN][i2];
621 y2=kcl->GetY(); z2=kcl->GetZ();
626 Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
627 if (TMath::Abs(zz-z2)>5.) continue;
629 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
630 if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
634 x[4]=f1(x1,y1,x2,y2,x3,y3);
635 if (TMath::Abs(x[4]) >= 0.0066) continue;
636 x[2]=f2(x1,y1,x2,y2,x3,y3);
637 //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
638 x[3]=f3(x1,y1,x2,y2,z1,z2);
639 if (TMath::Abs(x[3]) > 1.2) continue;
640 Double_t a=asin(x[2]);
641 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
642 if (TMath::Abs(zv-z3)>10.) continue;
644 Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
645 Double_t sy2=kcl->GetSigmaY2(), sz2=kcl->GetSigmaZ2();
646 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
647 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
649 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
650 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
651 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
652 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
653 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
654 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
655 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
656 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
657 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
658 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
662 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
663 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
664 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
665 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
666 c[13]=f30*sy1*f40+f32*sy2*f42;
667 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
669 UInt_t index=kr1.GetIndex(is);
670 AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
671 Float_t l=fSectors->GetPadPitchWidth();
672 track->SetSampledEdx(kr1[is]->GetQ()/l,0);
674 Int_t rc=FollowProlongation(*track, i2);
675 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
676 else fSeeds->AddLast(track);
682 //_____________________________________________________________________________
683 Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
684 //-----------------------------------------------------------------
685 // This function reades track seeds.
686 //-----------------------------------------------------------------
687 TDirectory *savedir=gDirectory;
689 TFile *in=(TFile*)inp;
691 cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
696 TTree *seedTree=(TTree*)in->Get("Seeds");
698 cerr<<"AliTPCtracker::ReadSeeds(): ";
699 cerr<<"can't get a tree with track seeds !\n";
702 AliTPCtrack *seed=new AliTPCtrack;
703 seedTree->SetBranchAddress("tracks",&seed);
705 if (fSeeds==0) fSeeds=new TObjArray(15000);
707 Int_t n=(Int_t)seedTree->GetEntries();
708 for (Int_t i=0; i<n; i++) {
709 seedTree->GetEvent(i);
710 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
715 delete seedTree; //Thanks to Mariana Bondila
722 //_____________________________________________________________________________
723 Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
724 //-----------------------------------------------------------------
725 // This is a track finder.
726 //-----------------------------------------------------------------
727 TDirectory *savedir=gDirectory;
730 TFile *in=(TFile*)inp;
732 cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
737 if (!out->IsOpen()) {
738 cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
747 sprintf(tname,"TreeT_TPC_%d",GetEventNumber());
748 TTree tracktree(tname,"Tree with TPC tracks");
749 AliTPCtrack *iotrack=0;
750 tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,3);
753 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
754 Int_t nrows=nlow+nup;
756 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
757 fSectors=fOuterSec; fN=fkNOS;
758 fSeeds=new TObjArray(15000);
759 MakeSeeds(nup-1, nup-1-gap);
760 MakeSeeds(nup-1-shift, nup-1-shift-gap);
765 Int_t nseed=fSeeds->GetEntriesFast();
766 for (Int_t i=0; i<nseed; i++) {
767 //tracking in the outer sectors
768 fSectors=fOuterSec; fN=fkNOS;
770 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
771 if (!FollowProlongation(t)) {
772 delete fSeeds->RemoveAt(i);
776 //tracking in the inner sectors
777 fSectors=fInnerSec; fN=fkNIS;
779 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
780 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
781 if (alpha < 0. ) alpha += 2.*TMath::Pi();
782 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
784 alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
786 if (t.Rotate(alpha)) {
787 if (FollowProlongation(t)) {
788 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
790 CookLabel(pt,0.1); //For comparison only
791 pt->PropagateTo(fParam->GetInnerRadiusLow());
796 // cerr<<found<<'\r';
800 delete fSeeds->RemoveAt(i);
805 cerr<<"Number of found tracks : "<<found<<endl;
810 fSeeds->Clear(); delete fSeeds; fSeeds=0;
814 //_____________________________________________________________________________
816 Int_t AliTPCtracker::RefitInward(TFile *in, TFile *out) {
818 // The function propagates tracks throught TPC inward
819 // using already associated clusters.
821 // in - file with back propagated tracks
822 // outs - file with inward propagation
824 // Sylwester Radomski, GSI
830 cout << "Input File not open !\n" << endl;
834 if (!out->IsOpen()) {
835 cout << "Output File not open !\n" << endl;
839 TDirectory *savedir = gDirectory;
842 // Connect to input tree
844 TTree *inputTree = (TTree*)in->Get("seedsTPCtoTRD_0");
845 TBranch *inBranch = inputTree->GetBranch("tracks");
846 AliTPCtrack *inTrack = 0;
847 inBranch->SetAddress(&inTrack);
849 Int_t nTracks = (Int_t)inputTree->GetEntries();
851 // Create output tree
853 AliTPCtrack *outTrack = new AliTPCtrack();
854 TTree *outputTree = new TTree("tracksTPC_0","Refited TPC tracks");
855 outputTree->Branch("tracks", "AliTPCtrack", &outTrack, 32000, 3);
859 for(Int_t t=0; t < nTracks; t++) {
863 inputTree->GetEvent(t);
864 AliTPCseed *seed = new AliTPCseed(*inTrack, inTrack->GetAlpha());
866 seed->ResetCovariance();
869 fSectors = fInnerSec;
870 Int_t res = FollowRefitInward(seed, inTrack);
872 Int_t nc = seed->GetNumberOfClusters();
874 fSectors = fOuterSec;
875 res = FollowRefitInward(seed, inTrack);
876 UseClusters(seed, nc);
879 seed->PropagateTo(fParam->GetInnerRadiusLow());
880 outTrack = (AliTPCtrack*)seed;
881 outTrack->SetLabel(inTrack->GetLabel());
889 cout << "Refitted " << nRefited << " tracks." << endl;
893 if (inputTree) delete inputTree;
894 if (outputTree) delete outputTree;
902 //_____________________________________________________________________________
903 Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
904 //-----------------------------------------------------------------
905 // This function propagates tracks back through the TPC.
906 //-----------------------------------------------------------------
907 return PropagateBack(inp, NULL, out);
910 //_____________________________________________________________________________
911 Int_t AliTPCtracker::PropagateBack(const TFile *inp, const TFile *inp2, TFile *out) {
912 //-----------------------------------------------------------------
913 // This function propagates tracks back through the TPC.
914 //-----------------------------------------------------------------
915 fSeeds=new TObjArray(15000);
917 TFile *in=(TFile*)inp;
918 TFile *in2=(TFile*)inp2;
919 TDirectory *savedir=gDirectory;
922 cerr<<"AliTPCtracker::PropagateBack(): ";
923 cerr<<"file with TPC (or back propagated ITS) tracks is not open !\n";
927 if (!out->IsOpen()) {
928 cerr<<"AliTPCtracker::PropagateBack(): ";
929 cerr<<"file for back propagated TPC tracks is not open !\n";
937 sprintf(tName,"TreeT_ITSb_%d",GetEventNumber());
938 TTree *bckTree=(TTree*)in->Get(tName);
939 if (!bckTree && inp2) bckTree=(TTree*)in2->Get(tName);
941 cerr<<"AliTPCtracker::PropagateBack() ";
942 cerr<<"can't get a tree with back propagated ITS tracks !\n";
945 AliTPCtrack *bckTrack=new AliTPCtrack;
946 if (bckTree) bckTree->SetBranchAddress("tracks",&bckTrack);
948 sprintf(tName,"TreeT_TPC_%d",GetEventNumber());
949 TTree *tpcTree=(TTree*)in->Get(tName);
951 cerr<<"AliTPCtracker::PropagateBack() ";
952 cerr<<"can't get a tree with TPC tracks !\n";
955 AliTPCtrack *tpcTrack=new AliTPCtrack;
956 tpcTree->SetBranchAddress("tracks",&tpcTrack);
958 //*** Prepare an array of tracks to be back propagated
959 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
960 Int_t nrows=nlow+nup;
963 // Match ITS tracks with old TPC tracks
964 // the tracks do not have to be sorted [SR, GSI, 18.02.2003]
966 // the algorithm is linear and uses LUT for sorting
967 // cost of the algorithm: 2 * number of tracks
970 TObjArray tracks(15000);
972 Int_t tpcN= (Int_t)tpcTree->GetEntries();
973 Int_t bckN= (bckTree)? (Int_t)bckTree->GetEntries() : 0;
976 const Int_t nLab = 200000; // limit on number of primaries (arbitrary)
978 for(Int_t i=0; i<nLab; i++) lut[i] = -1;
981 for(Int_t i=0; i<bckN; i++) {
982 bckTree->GetEvent(i);
983 Int_t lab = TMath::Abs(bckTrack->GetLabel());
984 if (lab < nLab) lut[lab] = i;
986 cerr << "AliTPCtracker: The size of the LUT is too small\n";
991 for (Int_t i=0; i<tpcN; i++) {
992 tpcTree->GetEvent(i);
993 Double_t alpha=tpcTrack->GetAlpha();
997 // No ITS - use TPC track only
999 AliTPCseed * seed = new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha());
1000 seed->ResetCovariance();
1001 fSeeds->AddLast(seed);
1002 tracks.AddLast(new AliTPCtrack(*tpcTrack));
1007 // discard not prolongated tracks (to be discussed)
1009 Int_t lab = TMath::Abs(tpcTrack->GetLabel());
1010 if (lab > nLab || lut[lab] < 0) continue;
1012 bckTree->GetEvent(lut[lab]);
1013 bckTrack->Rotate(alpha-bckTrack->GetAlpha());
1015 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
1016 tracks.AddLast(new AliTPCtrack(*tpcTrack));
1023 TObjArray tracks(15000);
1025 Int_t tpcN=(Int_t)tpcTree->GetEntries();
1026 Int_t bckN=(Int_t)bckTree->GetEntries();
1027 if (j<bckN) bckTree->GetEvent(j++);
1028 for (i=0; i<tpcN; i++) {
1029 tpcTree->GetEvent(i);
1030 Double_t alpha=tpcTrack->GetAlpha();
1032 TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
1033 if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
1034 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
1035 bckTree->GetEvent(j++);
1037 tpcTrack->ResetCovariance();
1038 fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
1040 tracks.AddLast(new AliTPCtrack(*tpcTrack));
1046 // tree name seedsTPCtoTRD as expected by TRD and as
1047 // discussed and decided in Strasbourg (May 2002)
1048 // [SR, GSI, 18.02.2003]
1050 sprintf(tName,"seedsTPCtoTRD_%d",GetEventNumber());
1051 TTree backTree(tName,"Tree with back propagated TPC tracks");
1053 AliTPCtrack *otrack=0;
1054 backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
1056 //*** Back propagation through inner sectors
1057 fSectors=fInnerSec; fN=fkNIS;
1058 Int_t nseed=fSeeds->GetEntriesFast();
1059 for (i=0; i<nseed; i++) {
1060 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
1061 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
1063 Int_t nc=t.GetNumberOfClusters();
1064 s.SetLabel(nc-1); //set number of the cluster to start with
1066 if (FollowBackProlongation(s,t)) {
1070 delete fSeeds->RemoveAt(i);
1073 //*** Back propagation through outer sectors
1074 fSectors=fOuterSec; fN=fkNOS;
1076 for (i=0; i<nseed; i++) {
1077 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
1079 Int_t nc=s.GetNumberOfClusters();
1080 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
1082 Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1083 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1084 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1085 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1087 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1088 alpha-=s.GetAlpha();
1090 if (s.Rotate(alpha)) {
1091 if (FollowBackProlongation(s,t)) {
1092 if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
1094 s.SetLabel(t.GetLabel());
1097 // Propagate to outer reference plane for comparison reasons
1098 // reason for keeping fParam object [SR, GSI, 18.02.2003]
1100 ps->PropagateTo(fParam->GetOuterRadiusUp());
1104 // cerr<<found<<'\r';
1109 delete fSeeds->RemoveAt(i);
1114 cerr<<"Number of seeds: "<<nseed<<endl;
1115 cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
1116 cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
1121 if (bckTree) delete bckTree; //Thanks to Mariana Bondila
1122 delete tpcTree; //Thanks to Mariana Bondila
1129 //_________________________________________________________________________
1130 AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
1131 //--------------------------------------------------------------------
1132 // Return pointer to a given cluster
1133 //--------------------------------------------------------------------
1134 Int_t sec=(index&0xff000000)>>24;
1135 Int_t row=(index&0x00ff0000)>>16;
1136 Int_t ncl=(index&0x0000ffff)>>00;
1138 const AliTPCcluster *cl=0;
1141 cl=fInnerSec[sec][row].GetUnsortedCluster(ncl);
1144 cl=fOuterSec[sec][row].GetUnsortedCluster(ncl);
1147 return (AliCluster*)cl;
1150 //__________________________________________________________________________
1151 void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1152 //--------------------------------------------------------------------
1153 //This function "cooks" a track label. If label<0, this track is fake.
1154 //--------------------------------------------------------------------
1155 Int_t noc=t->GetNumberOfClusters();
1156 Int_t *lb=new Int_t[noc];
1157 Int_t *mx=new Int_t[noc];
1158 AliCluster **clusters=new AliCluster*[noc];
1161 for (i=0; i<noc; i++) {
1163 Int_t index=t->GetClusterIndex(i);
1164 clusters[i]=GetCluster(index);
1167 Int_t lab=123456789;
1168 for (i=0; i<noc; i++) {
1169 AliCluster *c=clusters[i];
1170 lab=TMath::Abs(c->GetLabel(0));
1172 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
1178 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
1180 for (i=0; i<noc; i++) {
1181 AliCluster *c=clusters[i];
1182 if (TMath::Abs(c->GetLabel(1)) == lab ||
1183 TMath::Abs(c->GetLabel(2)) == lab ) max++;
1186 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
1189 Int_t tail=Int_t(0.10*noc);
1191 for (i=1; i<=tail; i++) {
1192 AliCluster *c=clusters[noc-i];
1193 if (lab == TMath::Abs(c->GetLabel(0)) ||
1194 lab == TMath::Abs(c->GetLabel(1)) ||
1195 lab == TMath::Abs(c->GetLabel(2))) max++;
1197 if (max < Int_t(0.5*tail)) lab=-lab;
1207 //_________________________________________________________________________
1208 void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
1209 //-----------------------------------------------------------------------
1210 // Setup inner sector
1211 //-----------------------------------------------------------------------
1213 fAlpha=par->GetInnerAngle();
1214 fAlphaShift=par->GetInnerAngleShift();
1215 fPadPitchWidth=par->GetInnerPadPitchWidth();
1216 f1PadPitchLength=par->GetInnerPadPitchLength();
1217 f2PadPitchLength=f1PadPitchLength;
1219 fN=par->GetNRowLow();
1220 fRow=new AliTPCRow[fN];
1221 for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
1223 fAlpha=par->GetOuterAngle();
1224 fAlphaShift=par->GetOuterAngleShift();
1225 fPadPitchWidth=par->GetOuterPadPitchWidth();
1226 f1PadPitchLength=par->GetOuter1PadPitchLength();
1227 f2PadPitchLength=par->GetOuter2PadPitchLength();
1229 fN=par->GetNRowUp();
1230 fRow=new AliTPCRow[fN];
1231 for (Int_t i=0; i<fN; i++){
1232 fRow[i].SetX(par->GetPadRowRadiiUp(i));
1237 //_________________________________________________________________________
1238 void AliTPCtracker::
1239 AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) {
1240 //-----------------------------------------------------------------------
1241 // Insert a cluster into this pad row in accordence with its y-coordinate
1242 //-----------------------------------------------------------------------
1243 if (fN==kMaxClusterPerRow) {
1244 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
1247 Int_t index=(((sec<<8)+row)<<16)+fN;
1251 fClusterArray[0]=*c; fClusters[fN++]=fClusterArray;
1257 AliTPCcluster *buff=new AliTPCcluster[size];
1258 memcpy(buff,fClusterArray,fSize*sizeof(AliTPCcluster));
1259 for (Int_t i=0; i<fN; i++)
1260 fClusters[i]=buff+(fClusters[i]-fClusterArray);
1261 delete[] fClusterArray;
1266 Int_t i=Find(c->GetY());
1267 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
1268 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
1270 fClusters[i]=fClusterArray+fN; fClusterArray[fN++]=*c;
1273 //___________________________________________________________________
1274 Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
1275 //-----------------------------------------------------------------------
1276 // Return the index of the nearest cluster
1277 //-----------------------------------------------------------------------
1279 if (y <= fClusters[0]->GetY()) return 0;
1280 if (y > fClusters[fN-1]->GetY()) return fN;
1281 Int_t b=0, e=fN-1, m=(b+e)/2;
1282 for (; b<e; m=(b+e)/2) {
1283 if (y > fClusters[m]->GetY()) b=m+1;
1289 //_____________________________________________________________________________
1290 void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1291 //-----------------------------------------------------------------
1292 // This funtion calculates dE/dX within the "low" and "up" cuts.
1293 //-----------------------------------------------------------------
1295 Int_t nc=GetNumberOfClusters();
1297 Int_t swap;//stupid sorting
1300 for (i=0; i<nc-1; i++) {
1301 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
1302 Float_t tmp=fdEdxSample[i];
1303 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
1308 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
1310 for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
1315 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1317 Double_t log1=TMath::Log(p+0.45), log2=TMath::Log(p+0.12);
1319 if (dedx < 34 + 30/(p+0.45)/(p+0.45) + 24*log1) {SetMass(0.13957); return;}
1320 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.49368); return;}
1321 SetMass(0.93827); return;
1325 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.13957); return;}
1326 SetMass(0.93827); return;
1329 SetMass(0.13957); return;