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.35.2.3 2003/07/15 09:58:03 hristov
19 Corrected back-propagation (Yu.Belikov)
21 Revision 1.35.2.2 2003/07/14 09:19:33 hristov
22 TOF included in the combined PID (Yu.Belikov)
24 Revision 1.35.2.1 2003/07/11 10:53:01 hristov
25 Inward refit for TPC and TRD in the ESD schema (T.Kuhr)
27 Revision 1.35 2003/05/23 10:08:51 hristov
28 SetLabel replaced by SetNumber (Yu.Belikov)
30 Revision 1.34 2003/05/22 13:57:48 hristov
31 First implementation of ESD classes (Yu.Belikov)
33 Revision 1.32 2003/04/10 10:36:54 hristov
34 Code for unified TPC/TRD tracking (S.Radomski)
36 Revision 1.31 2003/03/19 17:14:11 hristov
37 Load/UnloadClusters added to the base class and the derived classes changed correspondingly. Possibility to give 2 input files for ITS and TPC tracks in PropagateBack. TRD tracker uses fEventN from the base class (T.Kuhr)
39 Revision 1.30 2003/02/28 16:13:32 hristov
42 Revision 1.29 2003/02/28 15:18:16 hristov
43 Corrections suggested by J.Chudoba
45 Revision 1.28 2003/02/27 16:15:52 hristov
46 Code for inward refitting (S.Radomski)
48 Revision 1.27 2003/02/25 16:47:58 hristov
49 allow back propagation for more than 1 event (J.Chudoba)
51 Revision 1.26 2003/02/19 08:49:46 hristov
52 Track time measurement (S.Radomski)
54 Revision 1.25 2003/01/28 16:43:35 hristov
55 Additional protection: to be discussed with the Root team (M.Ivanov)
57 Revision 1.24 2002/11/19 16:13:24 hristov
58 stdlib.h included to declare exit() on HP
60 Revision 1.23 2002/11/19 11:50:08 hristov
61 Removing CONTAINERS (Yu.Belikov)
63 Revision 1.19 2002/07/19 07:31:40 kowal2
64 Improvement in tracking by J. Belikov
66 Revision 1.18 2002/05/13 07:33:52 kowal2
67 Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
68 in the case of defined region of interests.
70 Revision 1.17 2002/03/18 17:59:13 kowal2
71 Chnges in the pad geometry - 3 pad lengths introduced.
73 Revision 1.16 2001/11/08 16:39:03 hristov
74 Additional protection (M.Masera)
76 Revision 1.15 2001/11/08 16:36:33 hristov
77 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)
79 Revision 1.14 2001/10/21 19:04:55 hristov
80 Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
82 Revision 1.13 2001/07/23 12:01:30 hristov
85 Revision 1.12 2001/07/20 14:32:44 kowal2
86 Processing of many events possible now
88 Revision 1.11 2001/05/23 08:50:10 hristov
91 Revision 1.10 2001/05/16 14:57:25 alibrary
92 New files for folders and Stack
94 Revision 1.9 2001/05/11 07:16:56 hristov
95 Fix needed on Sun and Alpha
97 Revision 1.8 2001/05/08 15:00:15 hristov
98 Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
100 Revision 1.5 2000/12/20 07:51:59 kowal2
101 Changes suggested by Alessandra and Paolo to avoid overlapped
102 data fields in encapsulated classes.
104 Revision 1.4 2000/11/02 07:27:16 kowal2
107 Revision 1.2 2000/06/30 12:07:50 kowal2
108 Updated from the TPC-PreRelease branch
110 Revision 1.1.2.1 2000/06/25 08:53:55 kowal2
111 Splitted from AliTPCtracking
115 //-------------------------------------------------------
116 // Implementation of the TPC tracker
118 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
119 //-------------------------------------------------------
120 #include <TObjArray.h>
128 #include "AliTPCtracker.h"
129 #include "AliTPCcluster.h"
130 #include "AliTPCParam.h"
131 #include "AliClusters.h"
133 ClassImp(AliTPCtracker)
135 //_____________________________________________________________________________
136 AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
137 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
139 //---------------------------------------------------------------------
140 // The main TPC tracker constructor
141 //---------------------------------------------------------------------
142 fInnerSec=new AliTPCSector[fkNIS];
143 fOuterSec=new AliTPCSector[fkNOS];
146 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
147 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
149 fParam = (AliTPCParam*) par;
161 //_____________________________________________________________________________
162 AliTPCtracker::~AliTPCtracker() {
163 //------------------------------------------------------------------
164 // TPC tracker destructor
165 //------------------------------------------------------------------
175 fBarrelFile->Close();
180 //_____________________________________________________________________________
182 void AliTPCtracker::SetBarrelTree(const char *mode) {
184 // Creates a tree for BarrelTracks
185 // mode = "back" or "refit"
190 if (!IsStoringBarrel()) return;
192 TDirectory *sav = gDirectory;
193 if (!fBarrelFile) fBarrelFile = new TFile("AliBarrelTracks.root", "UPDATE");
196 sprintf(buff, "BarrelTPC_%d_%s", GetEventNumber(), mode);
199 fBarrelTree = new TTree(buff, "Barrel TPC tracks");
201 if (!fBarrelArray) fBarrelArray = new TClonesArray("AliBarrelTrack", 4);
202 for(Int_t i=0; i<4; i++) new((*fBarrelArray)[i]) AliBarrelTrack();
204 fBarrelTree->Branch("tracks", &fBarrelArray);
208 //_____________________________________________________________________________
210 void AliTPCtracker::StoreBarrelTrack(AliTPCtrack *ps, Int_t refPlane, Int_t isIn) {
212 // Stores Track at a given reference plane
215 // isIn: 1 - backward, 2 - refit
218 if (!IsStoringBarrel()) return;
219 if (refPlane < 0 || refPlane > 4) return;
220 if (isIn > 2) return;
222 static Int_t nClusters;
224 static Double_t chi2;
227 Int_t newClusters, newWrong;
230 if ( (refPlane == 1 && isIn == kTrackBack) ||
231 (refPlane == 4 && isIn == kTrackRefit) ) {
233 fBarrelArray->Clear();
234 nClusters = nWrong = 0;
241 if (refPlane == 1) refX = fParam->GetInnerRadiusLow();
242 if (refPlane == 2) refX = fParam->GetInnerRadiusUp();
243 if (refPlane == 3) refX = fParam->GetOuterRadiusLow();
244 if (refPlane == 4) refX = fParam->GetOuterRadiusUp();
246 ps->PropagateTo(refX);
248 fBarrelTrack = (AliBarrelTrack*)(*fBarrelArray)[index++];
249 ps->GetBarrelTrack(fBarrelTrack);
251 newClusters = ps->GetNumberOfClusters() - nClusters;
252 newWrong = ps->GetNWrong() - nWrong;
253 newChi2 = ps->GetChi2() - chi2;
255 nClusters = ps->GetNumberOfClusters();
256 nWrong = ps->GetNWrong();
257 chi2 = ps->GetChi2();
259 fBarrelTrack->SetNClusters(newClusters, newChi2);
260 fBarrelTrack->SetNWrongClusters(newWrong);
261 fBarrelTrack->SetRefPlane(refPlane, isIn);
264 //_____________________________________________________________________________
265 Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
268 // Parametrised error of the cluster reconstruction (pad direction)
271 const Float_t kArphi=0.41818e-2;
272 const Float_t kBrphi=0.17460e-4;
273 const Float_t kCrphi=0.30993e-2;
274 const Float_t kDrphi=0.41061e-3;
276 pt=TMath::Abs(pt)*1000.;
279 Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
280 if (s<0.4e-3) s=0.4e-3;
281 s*=1.3; //Iouri Belikov
286 //_____________________________________________________________________________
287 Double_t SigmaZ2(Double_t r, Double_t tgl)
290 // Parametrised error of the cluster reconstruction (drift direction)
293 const Float_t kAz=0.39614e-2;
294 const Float_t kBz=0.22443e-4;
295 const Float_t kCz=0.51504e-1;
299 Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
300 if (s<0.4e-3) s=0.4e-3;
301 s*=1.3; //Iouri Belikov
306 //_____________________________________________________________________________
307 Double_t f1(Double_t x1,Double_t y1,
308 Double_t x2,Double_t y2,
309 Double_t x3,Double_t y3)
311 //-----------------------------------------------------------------
312 // Initial approximation of the track curvature
313 //-----------------------------------------------------------------
314 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
315 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
316 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
317 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
318 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
320 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
322 return -xr*yr/sqrt(xr*xr+yr*yr);
326 //_____________________________________________________________________________
327 Double_t f2(Double_t x1,Double_t y1,
328 Double_t x2,Double_t y2,
329 Double_t x3,Double_t y3)
331 //-----------------------------------------------------------------
332 // Initial approximation of the track curvature times center of curvature
333 //-----------------------------------------------------------------
334 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
335 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
336 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
337 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
338 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
340 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
342 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
345 //_____________________________________________________________________________
346 Double_t f3(Double_t x1,Double_t y1,
347 Double_t x2,Double_t y2,
348 Double_t z1,Double_t z2)
350 //-----------------------------------------------------------------
351 // Initial approximation of the tangent of the track dip angle
352 //-----------------------------------------------------------------
353 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
356 //_____________________________________________________________________________
357 Int_t AliTPCtracker::LoadClusters(TTree *cTree) {
358 //-----------------------------------------------------------------
359 // This function loads TPC clusters.
360 //-----------------------------------------------------------------
361 TBranch *branch=cTree->GetBranch("Segment");
363 Error("LoadClusters","Can't get the branch !");
367 AliClusters carray, *addr=&carray;
368 carray.SetClass("AliTPCcluster");
370 branch->SetAddress(&addr);
372 Int_t nentr=(Int_t)cTree->GetEntries();
374 for (Int_t i=0; i<nentr; i++) {
377 Int_t ncl=carray.GetArray()->GetEntriesFast();
379 Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows();
380 Int_t id=carray.GetID();
381 if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) {
382 Error("LoadClusters","Wrong index !");
385 Int_t outindex = 2*fkNIS*nir;
388 Int_t row = id - sec*nir;
390 AliTPCRow &padrow=fInnerSec[sec][row];
392 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
393 padrow.InsertCluster(c,sec,row);
398 Int_t row = id - sec*nor;
400 AliTPCRow &padrow=fOuterSec[sec][row];
402 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
403 padrow.InsertCluster(c,sec+fkNIS,row);
406 carray.GetArray()->Clear();
412 //_____________________________________________________________________________
413 void AliTPCtracker::UnloadClusters() {
414 //-----------------------------------------------------------------
415 // This function unloads TPC clusters.
416 //-----------------------------------------------------------------
418 for (i=0; i<fkNIS; i++) {
419 Int_t nr=fInnerSec->GetNRows();
420 for (Int_t n=0; n<nr; n++) fInnerSec[i][n].ResetClusters();
422 for (i=0; i<fkNOS; i++) {
423 Int_t nr=fOuterSec->GetNRows();
424 for (Int_t n=0; n<nr; n++) fOuterSec[i][n].ResetClusters();
428 //_____________________________________________________________________________
429 Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
430 //-----------------------------------------------------------------
431 // This function tries to find a track prolongation.
432 //-----------------------------------------------------------------
433 Double_t xt=t.GetX();
434 const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
435 Int_t(0.5*fSectors->GetNRows());
436 Int_t tryAgain=kSKIP;
438 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
439 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
440 if (alpha < 0. ) alpha += 2.*TMath::Pi();
441 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
443 Int_t nrows=fSectors->GetRowNumber(xt)-1;
444 for (Int_t nr=nrows; nr>=rf; nr--) {
445 Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
446 if (!t.PropagateTo(x)) return 0;
450 Double_t maxchi2=kMaxCHI2;
451 const AliTPCRow &krow=fSectors[s][nr];
452 Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
453 Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
454 Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
455 Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
458 if (t.GetNumberOfClusters()>4)
459 Warning("FindProlongation","Too broad road !");
464 for (Int_t i=krow.Find(y-road); i<krow; i++) {
465 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
466 if (c->GetY() > y+road) break;
467 if (c->IsUsed()) continue;
468 if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
469 Double_t chi2=t.GetPredictedChi2(c);
470 if (chi2 > maxchi2) continue;
473 index=krow.GetIndex(i);
477 Float_t l=fSectors->GetPadPitchWidth();
478 Float_t corr=1.; if (nr>63) corr=0.67; // new (third) pad response !
479 t.SetSampledEdx(cl->GetQ()/l*corr,t.GetNumberOfClusters());
480 if (!t.Update(cl,maxchi2,index)) {
481 if (!tryAgain--) return 0;
482 } else tryAgain=kSKIP;
484 if (tryAgain==0) break;
487 if (!t.Rotate(fSectors->GetAlpha())) return 0;
488 } else if (y <-ymax) {
490 if (!t.Rotate(-fSectors->GetAlpha())) return 0;
498 //_____________________________________________________________________________
500 Int_t AliTPCtracker::FollowRefitInward(AliTPCseed *seed, AliTPCtrack *track) {
502 // This function propagates seed inward TPC using old clusters
505 // Sylwester Radomski, GSI
511 Int_t nRows = fSectors->GetNRows();
512 for (Int_t iRow = nRows; iRow >= 0; iRow--) {
514 Double_t x = fSectors->GetX(iRow);
515 if (!seed->PropagateTo(x)) return 0;
517 // try to find an assigned cluster in this row
519 AliTPCcluster* cluster = NULL;
522 for (Int_t iCluster = 0; iCluster < track->GetNumberOfClusters(); iCluster++) {
523 idx = track->GetClusterIndex(iCluster);
524 sec = (idx&0xff000000)>>24;
525 Int_t row = (idx&0x00ff0000)>>16;
526 if (((fSectors == fInnerSec) && (sec >= fkNIS)) ||
527 ((fSectors == fOuterSec) && (sec < fkNIS))) continue;
529 cluster = (AliTPCcluster*) GetCluster(idx);
534 // update the track seed with the found cluster
537 Double_t dAlpha = fParam->GetAngle(sec) - seed->GetAlpha();
538 if (TMath::Abs(dAlpha) > 0.0001) {
539 if (!seed->Rotate(dAlpha)) return 0;
540 if (!seed->PropagateTo(x)) return 0;
543 seed->Update(cluster, seed->GetPredictedChi2(cluster), idx);
550 //_____________________________________________________________________________
551 Int_t AliTPCtracker::FollowBackProlongation
552 (AliTPCseed& seed, const AliTPCtrack &track) {
553 //-----------------------------------------------------------------
554 // This function propagates tracks back through the TPC
555 //-----------------------------------------------------------------
556 Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
557 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
558 if (alpha < 0. ) alpha += 2.*TMath::Pi();
559 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
561 Int_t idx=-1, sec=-1, row=-1;
562 Int_t nc=seed.GetNumber();
565 idx=track.GetClusterIndex(nc);
566 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
568 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
569 else { if (sec < fkNIS) row=-1; }
571 Int_t nr=fSectors->GetNRows();
572 for (Int_t i=0; i<nr; i++) {
573 Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
574 Double_t y=seed.GetYat(x);
578 if (!seed.Rotate(fSectors->GetAlpha())) return 0;
579 } else if (y <-ymax) {
581 if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
584 if (!seed.PropagateTo(x)) return 0;
588 Double_t maxchi2=kMaxCHI2;
589 Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
590 Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
591 Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
592 Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
594 Warning("FollowBackProlongation","Too broad road !");
598 Int_t accepted=seed.GetNumberOfClusters();
600 //try to accept already found cluster
601 AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
603 if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
604 index=idx; cl=c; maxchi2=chi2;
605 } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
608 idx=track.GetClusterIndex(nc);
609 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
611 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
612 else { if (sec < fkNIS) row=-1; }
616 //try to fill the gap
617 const AliTPCRow &krow=fSectors[s][i];
620 for (Int_t i=krow.Find(y-road); i<krow; i++) {
621 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
622 if (c->GetY() > y+road) break;
623 if (c->IsUsed()) continue;
624 if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
625 Double_t chi2=seed.GetPredictedChi2(c);
626 if (chi2 > maxchi2) continue;
629 index=krow.GetIndex(i);
635 Float_t l=fSectors->GetPadPitchWidth();
636 Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
637 seed.SetSampledEdx(cl->GetQ()/l*corr,seed.GetNumberOfClusters());
638 seed.Update(cl,maxchi2,index);
648 //_____________________________________________________________________________
649 void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
650 //-----------------------------------------------------------------
651 // This function creates track seeds.
652 //-----------------------------------------------------------------
653 Double_t x[5], c[15];
655 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
656 Double_t cs=cos(alpha), sn=sin(alpha);
658 Double_t x1 =fSectors->GetX(i1);
659 Double_t xx2=fSectors->GetX(i2);
661 for (Int_t ns=0; ns<fN; ns++) {
662 Int_t nl=fSectors[(ns-1+fN)%fN][i2];
663 Int_t nm=fSectors[ns][i2];
664 Int_t nu=fSectors[(ns+1)%fN][i2];
665 const AliTPCRow& kr1=fSectors[ns][i1];
666 for (Int_t is=0; is < kr1; is++) {
667 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
668 for (Int_t js=0; js < nl+nm+nu; js++) {
669 const AliTPCcluster *kcl;
671 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
674 const AliTPCRow& kr2=fSectors[(ns-1+fN)%fN][i2];
676 y2=kcl->GetY(); z2=kcl->GetZ();
681 const AliTPCRow& kr2=fSectors[ns][i2];
683 x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
685 const AliTPCRow& kr2=fSectors[(ns+1)%fN][i2];
687 y2=kcl->GetY(); z2=kcl->GetZ();
692 Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
693 if (TMath::Abs(zz-z2)>5.) continue;
695 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
697 Warning("MakeSeeds","Straight seed !");
702 x[4]=f1(x1,y1,x2,y2,x3,y3);
703 if (TMath::Abs(x[4]) >= 0.0066) continue;
704 x[2]=f2(x1,y1,x2,y2,x3,y3);
705 //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
706 x[3]=f3(x1,y1,x2,y2,z1,z2);
707 if (TMath::Abs(x[3]) > 1.2) continue;
708 Double_t a=asin(x[2]);
709 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
710 if (TMath::Abs(zv-z3)>10.) continue;
712 Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
713 Double_t sy2=kcl->GetSigmaY2(), sz2=kcl->GetSigmaZ2();
714 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
715 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
717 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
718 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
719 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
720 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
721 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
722 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
723 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
724 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
725 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
726 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
730 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
731 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
732 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
733 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
734 c[13]=f30*sy1*f40+f32*sy2*f42;
735 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
737 UInt_t index=kr1.GetIndex(is);
738 AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
739 Float_t l=fSectors->GetPadPitchWidth();
740 track->SetSampledEdx(kr1[is]->GetQ()/l,0);
742 Int_t rc=FollowProlongation(*track, i2);
743 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
744 else fSeeds->AddLast(track);
750 //_____________________________________________________________________________
751 Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
752 //-----------------------------------------------------------------
753 // This function reades track seeds.
754 //-----------------------------------------------------------------
755 TDirectory *savedir=gDirectory;
757 TFile *in=(TFile*)inp;
759 Error("ReadSeeds","Input file has not been open !");
764 TTree *seedTree=(TTree*)in->Get("Seeds");
766 Error("ReadSeeds","Can't get a tree with track seeds !");
769 AliTPCtrack *seed=new AliTPCtrack;
770 seedTree->SetBranchAddress("tracks",&seed);
772 if (fSeeds==0) fSeeds=new TObjArray(15000);
774 Int_t n=(Int_t)seedTree->GetEntries();
775 for (Int_t i=0; i<n; i++) {
776 seedTree->GetEvent(i);
777 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
782 delete seedTree; //Thanks to Mariana Bondila
789 //_____________________________________________________________________________
790 Int_t AliTPCtracker::Clusters2Tracks(AliESD *event) {
791 //-----------------------------------------------------------------
792 // This is a track finder.
793 // The clusters must be already loaded !
794 //-----------------------------------------------------------------
797 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
798 Int_t nrows=nlow+nup;
800 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
801 fSectors=fOuterSec; fN=fkNOS;
802 fSeeds=new TObjArray(15000);
803 MakeSeeds(nup-1, nup-1-gap);
804 MakeSeeds(nup-1-shift, nup-1-shift-gap);
808 Int_t nseed=fSeeds->GetEntriesFast();
809 for (Int_t i=0; i<nseed; i++) {
810 //tracking in the outer sectors
811 fSectors=fOuterSec; fN=fkNOS;
813 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
814 if (!FollowProlongation(t)) {
815 delete fSeeds->RemoveAt(i);
819 //tracking in the inner sectors
820 fSectors=fInnerSec; fN=fkNIS;
822 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
823 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
824 if (alpha < 0. ) alpha += 2.*TMath::Pi();
825 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
827 alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
829 if (t.Rotate(alpha)) {
830 if (FollowProlongation(t)) {
831 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
833 CookLabel(pt,0.1); //For comparison only
834 pt->PropagateTo(fParam->GetInnerRadiusLow());
836 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
838 event->AddTrack(&iotrack);
844 delete fSeeds->RemoveAt(i);
847 Info("Clusters2Tracks","Number of found tracks : %d",
848 event->GetNumberOfTracks());
850 fSeeds->Clear(); delete fSeeds; fSeeds=0;
855 //_____________________________________________________________________________
856 Int_t AliTPCtracker::Clusters2Tracks(TTree *cTree, TTree *tTree) {
857 //-----------------------------------------------------------------
858 // This is a track finder.
859 //-----------------------------------------------------------------
860 if (LoadClusters(cTree)) {
861 Error("Clusters2Tracks","Problem with loading the clusters...");
865 AliTPCtrack *iotrack=0;
866 TBranch *branch=tTree->GetBranch("tracks");
867 if (!branch) tTree->Branch("tracks","AliTPCtrack",&iotrack,32000,3);
868 else branch->SetAddress(&iotrack);
871 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
872 Int_t nrows=nlow+nup;
874 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
875 fSectors=fOuterSec; fN=fkNOS;
876 fSeeds=new TObjArray(15000);
877 MakeSeeds(nup-1, nup-1-gap);
878 MakeSeeds(nup-1-shift, nup-1-shift-gap);
883 Int_t nseed=fSeeds->GetEntriesFast();
885 for (Int_t i=0; i<nseed; i++) {
886 //tracking in the outer sectors
887 fSectors=fOuterSec; fN=fkNOS;
889 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
890 if (!FollowProlongation(t)) {
891 delete fSeeds->RemoveAt(i);
895 //tracking in the inner sectors
896 fSectors=fInnerSec; fN=fkNIS;
898 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
899 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
900 if (alpha < 0. ) alpha += 2.*TMath::Pi();
901 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
903 alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
905 if (t.Rotate(alpha)) {
906 if (FollowProlongation(t)) {
907 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
909 CookLabel(pt,0.1); //For comparison only
910 pt->PropagateTo(fParam->GetInnerRadiusLow());
918 delete fSeeds->RemoveAt(i);
921 Info("Clusters2Tracks","Number of found tracks : %d",found);
925 fSeeds->Clear(); delete fSeeds; fSeeds=0;
930 //_____________________________________________________________________________
931 Int_t AliTPCtracker::RefitInward(AliESD* event) {
933 // The function propagates tracks throught TPC inward
934 // using already associated clusters.
935 // The clusters must be already loaded !
938 Int_t nTracks = event->GetNumberOfTracks();
941 for (Int_t i = 0; i < nTracks; i++) {
942 AliESDtrack* track = event->GetTrack(i);
943 ULong_t status = track->GetStatus();
945 if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
946 if ( (status & AliESDtrack::kTPCrefit) != 0 ) continue;
948 AliTPCtrack* tpcTrack = new AliTPCtrack(*track);
949 AliTPCseed* seed = new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha());
951 fSectors = fOuterSec;
953 Int_t res = FollowRefitInward(seed, tpcTrack);
955 Int_t nc = seed->GetNumberOfClusters();
957 fSectors = fInnerSec;
959 res = FollowRefitInward(seed, tpcTrack);
960 UseClusters(seed, nc);
963 seed->PropagateTo(fParam->GetInnerRadiusLow());
964 seed->SetLabel(tpcTrack->GetLabel());
965 seed->SetdEdx(tpcTrack->GetdEdx());
966 track->UpdateTrackParams(seed, AliESDtrack::kTPCrefit);
974 Info("RefitInward","Number of refitted tracks : %d",nRefited);
979 //_____________________________________________________________________________
980 Int_t AliTPCtracker::RefitInward(TTree *in, TTree *out) {
982 // The function propagates tracks throught TPC inward
983 // using already associated clusters.
985 Error("RefitInward","This method is not converted to NewIO yet\n");
988 if (!inSeeds->IsOpen()) {
989 cout << "Input File with seeds not open !\n" << endl;
993 if (!inTracks->IsOpen()) {
994 cout << "Input File not open !\n" << endl;
998 if (!outTracks->IsOpen()) {
999 cout << "Output File not open !\n" << endl;
1003 TDirectory *savedir = gDirectory;
1006 // Connect to input seeds tree
1009 sprintf(tname, "seedTRDtoTPC_%d", GetEventNumber());
1010 TTree *seedsTree = (TTree*)inSeeds->Get(tname);
1011 TBranch *inSeedBranch = seedsTree->GetBranch("tracks");
1012 AliTPCtrack *inSeedTrack = 0;
1013 inSeedBranch->SetAddress(&inSeedTrack);
1015 Int_t nSeeds = (Int_t)seedsTree->GetEntries();
1017 // Connect to input tree
1019 sprintf(tname,"TreeT_TPC_%d",GetEventNumber());
1020 // sprintf(tname,"seedsTPCtoTRD_%d",GetEventNumber());
1021 TTree *inputTree = (TTree*)inTracks->Get(tname);
1022 TBranch *inBranch = inputTree->GetBranch("tracks");
1023 AliTPCtrack *inTrack = 0;
1024 inBranch->SetAddress(&inTrack);
1026 Int_t nTracks = (Int_t)inputTree->GetEntries();
1028 // Create output tree
1030 AliTPCtrack *outTrack = new AliTPCtrack();
1031 sprintf(tname, "tracksTPC_%d", GetEventNumber());
1032 TTree *outputTree = new TTree(tname,"Refited TPC tracks");
1033 outputTree->Branch("tracks", "AliTPCtrack", &outTrack, 32000, 3);
1037 // create table of track labels
1038 Int_t* inLab = new Int_t[nTracks];
1039 for (Int_t i = 0; i < nTracks; i++) inLab[i] = -1;
1041 // [SR, 01.04.2003] - Barrel tracks
1042 if (IsStoringBarrel()) SetBarrelTree("refit");
1044 for(Int_t t=0; t < nSeeds; t++) {
1046 seedsTree->GetEntry(t);
1047 // find TPC track for seed
1048 Int_t lab = TMath::Abs(inSeedTrack->GetLabel());
1049 for(Int_t i=0; i < nTracks; i++) {
1050 if (inLab[i] == lab) {
1051 inputTree->GetEntry(i);
1053 } else if (inLab[i] < 0) {
1054 inputTree->GetEntry(i);
1055 inLab[i] = TMath::Abs(inTrack->GetLabel());
1056 if (inLab[i] == lab) break;
1059 if (TMath::Abs(inTrack->GetLabel()) != lab) continue;
1060 AliTPCseed *seed = new AliTPCseed(*inSeedTrack, inTrack->GetAlpha());
1061 if (IsStoringBarrel()) StoreBarrelTrack(seed, 4, 2);
1063 fSectors = fOuterSec;
1065 Int_t res = FollowRefitInward(seed, inTrack);
1067 Int_t nc = seed->GetNumberOfClusters();
1069 if (IsStoringBarrel()) StoreBarrelTrack(seed, 3, 2);
1070 if (IsStoringBarrel()) StoreBarrelTrack(seed, 2, 2);
1072 fSectors = fInnerSec;
1074 res = FollowRefitInward(seed, inTrack);
1075 UseClusters(seed, nc);
1077 if (IsStoringBarrel()) StoreBarrelTrack(seed, 1, 2);
1080 seed->PropagateTo(fParam->GetInnerRadiusLow());
1081 outTrack = (AliTPCtrack*)seed;
1082 outTrack->SetLabel(inTrack->GetLabel());
1083 outTrack->SetdEdx(inTrack->GetdEdx());
1088 if (IsStoringBarrel()) fBarrelTree->Fill();
1092 cout << "Refitted " << nRefited << " tracks." << endl;
1095 outputTree->Write();
1100 if (IsStoringBarrel()) {
1102 fBarrelTree->Write();
1103 fBarrelFile->Flush();
1105 delete fBarrelArray;
1109 if (inputTree) delete inputTree;
1110 if (outputTree) delete outputTree;
1119 Int_t AliTPCtracker::PropagateBack(AliESD *event) {
1120 //-----------------------------------------------------------------
1121 // This function propagates tracks back through the TPC.
1122 // The clusters must be already loaded !
1123 //-----------------------------------------------------------------
1124 Int_t nentr=event->GetNumberOfTracks();
1125 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
1128 for (Int_t i=0; i<nentr; i++) {
1129 AliESDtrack *esd=event->GetTrack(i);
1130 ULong_t status=esd->GetStatus();
1132 if ( (status & AliESDtrack::kTPCin ) == 0 ) continue;
1133 if ( (status & AliESDtrack::kTPCout) != 0 ) continue;
1135 const AliTPCtrack t(*esd);
1136 AliTPCseed s(t,t.GetAlpha());
1138 if ( (status & AliESDtrack::kITSout) == 0 ) s.ResetCovariance();
1143 Int_t nc=t.GetNumberOfClusters();
1144 s.SetNumber(nc); //set number of the cluster to start with
1147 fSectors=fInnerSec; fN=fkNIS;
1149 Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1150 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1151 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1152 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1153 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1154 alpha-=s.GetAlpha();
1156 if (!s.Rotate(alpha)) continue;
1157 if (!FollowBackProlongation(s,t)) continue;
1162 fSectors=fOuterSec; fN=fkNOS;
1164 nc=s.GetNumberOfClusters();
1166 alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1167 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1168 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1169 ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1171 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1172 alpha-=s.GetAlpha();
1174 if (!s.Rotate(alpha)) continue;
1175 if (!FollowBackProlongation(s,t)) continue;
1177 Int_t nrows=fOuterSec->GetNRows()+fInnerSec->GetNRows();
1178 if (s.GetNumberOfClusters() < Int_t(0.4*nrows)) continue;
1180 s.PropagateTo(fParam->GetOuterRadiusUp());
1182 CookLabel(&s,0.1); //For comparison only
1184 esd->UpdateTrackParams(&s,AliESDtrack::kTPCout);
1187 Info("PropagateBack","Number of back propagated tracks: %d",ntrk);
1192 //_____________________________________________________________________________
1193 Int_t AliTPCtracker::PropagateBack(TTree *in, TTree *out) {
1194 //-----------------------------------------------------------------
1195 // This function propagates tracks back through the TPC.
1196 //-----------------------------------------------------------------
1197 Error("PropagateBack","This method is not converted to NewIO yet\n");
1200 fSeeds=new TObjArray(15000);
1202 TFile *in=(TFile*)inp;
1203 TFile *in2=(TFile*)inp2;
1204 TDirectory *savedir=gDirectory;
1206 if (!in->IsOpen()) {
1207 cerr<<"AliTPCtracker::PropagateBack(): ";
1208 cerr<<"file with TPC (or back propagated ITS) tracks is not open !\n";
1212 if (!out->IsOpen()) {
1213 cerr<<"AliTPCtracker::PropagateBack(): ";
1214 cerr<<"file for back propagated TPC tracks is not open !\n";
1222 sprintf(tName,"TreeT_ITSb_%d",GetEventNumber());
1223 TTree *bckTree=(TTree*)in->Get(tName);
1224 if (!bckTree && inp2) bckTree=(TTree*)in2->Get(tName);
1226 cerr<<"AliTPCtracker::PropagateBack() ";
1227 cerr<<"can't get a tree with back propagated ITS tracks !\n";
1230 AliTPCtrack *bckTrack=new AliTPCtrack;
1231 if (bckTree) bckTree->SetBranchAddress("tracks",&bckTrack);
1233 sprintf(tName,"TreeT_TPC_%d",GetEventNumber());
1234 TTree *tpcTree=(TTree*)in->Get(tName);
1236 cerr<<"AliTPCtracker::PropagateBack() ";
1237 cerr<<"can't get a tree with TPC tracks !\n";
1240 AliTPCtrack *tpcTrack=new AliTPCtrack;
1241 tpcTree->SetBranchAddress("tracks",&tpcTrack);
1243 // [SR, 01.04.2003] - Barrel tracks
1244 if (IsStoringBarrel()) SetBarrelTree("back");
1246 // Prepare an array of tracks to be back propagated
1247 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1248 Int_t nrows=nlow+nup;
1251 // Match ITS tracks with old TPC tracks
1252 // the tracks do not have to be sorted [SR, GSI, 18.02.2003]
1254 // the algorithm is linear and uses LUT for sorting
1255 // cost of the algorithm: 2 * number of tracks
1258 TObjArray tracks(15000);
1260 Int_t tpcN= (Int_t)tpcTree->GetEntries();
1261 Int_t bckN= (bckTree)? (Int_t)bckTree->GetEntries() : 0;
1264 const Int_t nLab = 200000; // limit on number of primaries (arbitrary)
1266 for(Int_t i=0; i<nLab; i++) lut[i] = -1;
1269 for(Int_t i=0; i<bckN; i++) {
1270 bckTree->GetEvent(i);
1271 Int_t lab = TMath::Abs(bckTrack->GetLabel());
1272 if (lab < nLab) lut[lab] = i;
1274 cerr << "AliTPCtracker: The size of the LUT is too small\n";
1279 for (Int_t i=0; i<tpcN; i++) {
1280 tpcTree->GetEvent(i);
1281 Double_t alpha=tpcTrack->GetAlpha();
1285 // No ITS - use TPC track only
1287 tpcTrack->ResetCovariance();
1288 AliTPCseed * seed = new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha());
1290 fSeeds->AddLast(seed);
1291 tracks.AddLast(new AliTPCtrack(*tpcTrack));
1296 // discard not prolongated tracks (to be discussed)
1298 Int_t lab = TMath::Abs(tpcTrack->GetLabel());
1299 if (lab > nLab || lut[lab] < 0) continue;
1301 bckTree->GetEvent(lut[lab]);
1302 bckTrack->Rotate(alpha-bckTrack->GetAlpha());
1304 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
1305 tracks.AddLast(new AliTPCtrack(*tpcTrack));
1313 // tree name seedsTPCtoTRD as expected by TRD and as
1314 // discussed and decided in Strasbourg (May 2002)
1315 // [SR, GSI, 18.02.2003]
1317 sprintf(tName,"seedsTPCtoTRD_%d",GetEventNumber());
1318 TTree backTree(tName,"Tree with back propagated TPC tracks");
1320 AliTPCtrack *otrack=0;
1321 backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
1326 Int_t nseed=fSeeds->GetEntriesFast();
1328 // loop changed [SR, 01.04.2003]
1330 for (i=0; i<nseed; i++) {
1333 if (IsStoringBarrel()) fBarrelArray->Clear();
1335 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
1336 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
1339 ps->ResetNRotation();
1341 if (IsStoringBarrel()) StoreBarrelTrack(ps, 1, 1);
1343 // Load outer sectors
1347 Int_t nc=t.GetNumberOfClusters();
1348 //s.SetLabel(nc-1); //set number of the cluster to start with
1351 if (FollowBackProlongation(s,t)) UseClusters(&s);
1353 fSeeds->RemoveAt(i);
1357 if (IsStoringBarrel()) StoreBarrelTrack(ps, 2, 1);
1358 if (IsStoringBarrel()) StoreBarrelTrack(ps, 3, 1);
1360 // Load outer sectors
1364 nc=s.GetNumberOfClusters();
1366 Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1367 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1368 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1369 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1371 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1372 alpha-=s.GetAlpha();
1374 if (s.Rotate(alpha)) {
1375 if (FollowBackProlongation(s,t)) {
1376 if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
1378 s.SetLabel(t.GetLabel());
1381 // Propagate to outer reference plane for comparison reasons
1382 // reason for keeping fParam object [SR, GSI, 18.02.2003]
1384 ps->PropagateTo(fParam->GetOuterRadiusUp());
1388 if (IsStoringBarrel()) StoreBarrelTrack(ps, 4, 1);
1389 cerr<<found++<<'\r';
1394 if (IsStoringBarrel()) fBarrelTree->Fill();
1395 delete fSeeds->RemoveAt(i);
1403 if (IsStoringBarrel()) {
1405 fBarrelTree->Write();
1406 fBarrelFile->Flush();
1408 delete fBarrelArray;
1413 cerr<<"Number of seeds: "<<nseed<<endl;
1414 cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
1415 cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
1420 if (bckTree) delete bckTree; //Thanks to Mariana Bondila
1421 delete tpcTree; //Thanks to Mariana Bondila
1429 //_________________________________________________________________________
1430 AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
1431 //--------------------------------------------------------------------
1432 // Return pointer to a given cluster
1433 //--------------------------------------------------------------------
1434 Int_t sec=(index&0xff000000)>>24;
1435 Int_t row=(index&0x00ff0000)>>16;
1436 Int_t ncl=(index&0x0000ffff)>>00;
1438 const AliTPCcluster *cl=0;
1441 cl=fInnerSec[sec][row].GetUnsortedCluster(ncl);
1444 cl=fOuterSec[sec][row].GetUnsortedCluster(ncl);
1447 return (AliCluster*)cl;
1450 //__________________________________________________________________________
1451 void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1452 //--------------------------------------------------------------------
1453 //This function "cooks" a track label. If label<0, this track is fake.
1454 //--------------------------------------------------------------------
1455 Int_t noc=t->GetNumberOfClusters();
1456 Int_t *lb=new Int_t[noc];
1457 Int_t *mx=new Int_t[noc];
1458 AliCluster **clusters=new AliCluster*[noc];
1461 for (i=0; i<noc; i++) {
1463 Int_t index=t->GetClusterIndex(i);
1464 clusters[i]=GetCluster(index);
1467 Int_t lab=123456789;
1468 for (i=0; i<noc; i++) {
1469 AliCluster *c=clusters[i];
1470 lab=TMath::Abs(c->GetLabel(0));
1472 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
1478 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
1480 for (i=0; i<noc; i++) {
1481 AliCluster *c=clusters[i];
1482 if (TMath::Abs(c->GetLabel(1)) == lab ||
1483 TMath::Abs(c->GetLabel(2)) == lab ) max++;
1486 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
1489 Int_t tail=Int_t(0.10*noc);
1491 for (i=1; i<=tail; i++) {
1492 AliCluster *c=clusters[noc-i];
1493 if (lab == TMath::Abs(c->GetLabel(0)) ||
1494 lab == TMath::Abs(c->GetLabel(1)) ||
1495 lab == TMath::Abs(c->GetLabel(2))) max++;
1497 if (max < Int_t(0.5*tail)) lab=-lab;
1507 //_________________________________________________________________________
1508 void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
1509 //-----------------------------------------------------------------------
1510 // Setup inner sector
1511 //-----------------------------------------------------------------------
1513 fAlpha=par->GetInnerAngle();
1514 fAlphaShift=par->GetInnerAngleShift();
1515 fPadPitchWidth=par->GetInnerPadPitchWidth();
1516 f1PadPitchLength=par->GetInnerPadPitchLength();
1517 f2PadPitchLength=f1PadPitchLength;
1519 fN=par->GetNRowLow();
1520 fRow=new AliTPCRow[fN];
1521 for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
1523 fAlpha=par->GetOuterAngle();
1524 fAlphaShift=par->GetOuterAngleShift();
1525 fPadPitchWidth=par->GetOuterPadPitchWidth();
1526 f1PadPitchLength=par->GetOuter1PadPitchLength();
1527 f2PadPitchLength=par->GetOuter2PadPitchLength();
1529 fN=par->GetNRowUp();
1530 fRow=new AliTPCRow[fN];
1531 for (Int_t i=0; i<fN; i++){
1532 fRow[i].SetX(par->GetPadRowRadiiUp(i));
1537 //_________________________________________________________________________
1538 void AliTPCtracker::
1539 AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) {
1540 //-----------------------------------------------------------------------
1541 // Insert a cluster into this pad row in accordence with its y-coordinate
1542 //-----------------------------------------------------------------------
1543 if (fN==kMaxClusterPerRow) {
1544 ::Error("InsertCluster","Too many clusters !");
1548 Int_t index=(((sec<<8)+row)<<16)+fN;
1551 fSize=kMaxClusterPerRow/8;
1552 fClusterArray=new AliTPCcluster[fSize];
1554 fClusterArray[0]=*c; fClusters[fN++]=fClusterArray;
1560 AliTPCcluster *buff=new AliTPCcluster[size];
1561 memcpy(buff,fClusterArray,fSize*sizeof(AliTPCcluster));
1562 for (Int_t i=0; i<fN; i++)
1563 fClusters[i]=buff+(fClusters[i]-fClusterArray);
1564 delete[] fClusterArray;
1569 Int_t i=Find(c->GetY());
1570 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
1571 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
1573 fClusters[i]=fClusterArray+fN; fClusterArray[fN++]=*c;
1576 //___________________________________________________________________
1577 Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
1578 //-----------------------------------------------------------------------
1579 // Return the index of the nearest cluster
1580 //-----------------------------------------------------------------------
1582 if (y <= fClusters[0]->GetY()) return 0;
1583 if (y > fClusters[fN-1]->GetY()) return fN;
1584 Int_t b=0, e=fN-1, m=(b+e)/2;
1585 for (; b<e; m=(b+e)/2) {
1586 if (y > fClusters[m]->GetY()) b=m+1;
1592 //_____________________________________________________________________________
1593 void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1594 //-----------------------------------------------------------------
1595 // This funtion calculates dE/dX within the "low" and "up" cuts.
1596 //-----------------------------------------------------------------
1598 Int_t nc=GetNumberOfClusters();
1600 Int_t swap;//stupid sorting
1603 for (i=0; i<nc-1; i++) {
1604 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
1605 Float_t tmp=fdEdxSample[i];
1606 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
1611 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
1613 for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
1618 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1620 Double_t log1=TMath::Log(p+0.45), log2=TMath::Log(p+0.12);
1622 if (dedx < 34 + 30/(p+0.45)/(p+0.45) + 24*log1) {SetMass(0.13957); return;}
1623 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.49368); return;}
1624 SetMass(0.93827); return;
1628 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.13957); return;}
1629 SetMass(0.93827); return;
1632 SetMass(0.13957); return;