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.14 2001/11/14 10:50:46 cblume
19 Changes in digits IO. Add merging of summable digits
21 Revision 1.13 2001/05/30 12:17:47 hristov
22 Loop variables declared once
24 Revision 1.12 2001/05/28 17:07:58 hristov
25 Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
27 Revision 1.8 2000/12/20 13:00:44 cblume
28 Modifications for the HP-compiler
30 Revision 1.7 2000/12/08 16:07:02 cblume
31 Update of the tracking by Sergei
33 Revision 1.6 2000/11/30 17:38:08 cblume
34 Changes to get in line with new STEER and EVGEN
36 Revision 1.5 2000/11/14 14:40:27 cblume
37 Correction for the Sun compiler (kTRUE and kFALSE)
39 Revision 1.4 2000/11/10 14:57:52 cblume
40 Changes in the geometry constants for the DEC compiler
42 Revision 1.3 2000/10/15 23:40:01 cblume
45 Revision 1.2 2000/10/06 16:49:46 cblume
48 Revision 1.1.2.2 2000/10/04 16:34:58 cblume
49 Replace include files by forward declarations
51 Revision 1.1.2.1 2000/09/22 14:47:52 cblume
56 ////////////////////////////////////////////////////////////////////////////////
60 ////////////////////////////////////////////////////////////////////////////////
71 #include "AliTRDgeometry.h"
72 #include "AliTRDrecPoint.h"
73 #include "AliTRDcluster.h"
74 #include "AliTRDtrack.h"
75 #include "AliTRDtrackingSector.h"
76 #include "AliTRDtimeBin.h"
78 #include "AliTRDtracker.h"
80 ClassImp(AliTRDtracker)
82 const Float_t AliTRDtracker::fgkSeedDepth = 0.5;
83 const Float_t AliTRDtracker::fgkSeedStep = 0.05;
84 const Float_t AliTRDtracker::fgkSeedGap = 0.25;
86 const Float_t AliTRDtracker::fgkMaxSeedDeltaZ12 = 40.;
87 const Float_t AliTRDtracker::fgkMaxSeedDeltaZ = 25.;
88 const Float_t AliTRDtracker::fgkMaxSeedC = 0.0052;
89 const Float_t AliTRDtracker::fgkMaxSeedTan = 1.2;
90 const Float_t AliTRDtracker::fgkMaxSeedVertexZ = 150.;
92 const Double_t AliTRDtracker::fgkSeedErrorSY = 0.2;
93 const Double_t AliTRDtracker::fgkSeedErrorSY3 = 2.5;
94 const Double_t AliTRDtracker::fgkSeedErrorSZ = 0.1;
96 const Float_t AliTRDtracker::fgkMinClustersInSeed = 0.7;
98 const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5;
99 const Float_t AliTRDtracker::fgkSkipDepth = 0.05;
100 const Float_t AliTRDtracker::fgkLabelFraction = 0.5;
101 const Float_t AliTRDtracker::fgkWideRoad = 20.;
103 const Double_t AliTRDtracker::fgkMaxChi2 = 24.;
105 //____________________________________________________________________
106 AliTRDtracker::AliTRDtracker()
109 // Default constructor
125 //____________________________________________________________________
126 AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title)
130 // TRD tracker contructor
137 fClusters = new TObjArray(2000);
139 fSeeds = new TObjArray(20000);
141 fTracks = new TObjArray(10000);
146 //___________________________________________________________________
147 AliTRDtracker::~AliTRDtracker()
160 //___________________________________________________________________
161 void AliTRDtracker::Clusters2Tracks(TH1F *hs, TH1F *hd)
164 // Do the trackfinding
170 Int_t fTotalNofTimeBins = fGeom->GetTimeMax() * AliTRDgeometry::Nplan();
171 Int_t nSteps = (Int_t) (fgkSeedDepth / fgkSeedStep);
172 Int_t gap = (Int_t) (fTotalNofTimeBins * fgkSeedGap);
173 Int_t step = (Int_t) (fTotalNofTimeBins * fgkSeedStep);
178 if (!fClusters) return;
180 AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect];
181 SetUpSectors(fTrSec);
183 // make a first turn looking for seed ends in the same (n,n)
184 // and in the adjacent sectors (n,n+1)
186 for(i=0; i<nSteps; i++) {
187 printf("step %d out of %d \n", i+1, nSteps);
188 outer=fTotalNofTimeBins-1-i*step; inner=outer-gap;
189 MakeSeeds(inner,outer, fTrSec, 1, hs, hd);
190 FindTracks(fTrSec, hs, hd);
193 // make a second turn looking for seed ends in next-to-adjacent
196 for(i=0; i<nSteps; i++) {
197 printf("step %d out of %d \n", i+1, nSteps);
198 outer=fTotalNofTimeBins-1-i*step; inner=outer-gap;
199 MakeSeeds(inner, outer, fTrSec, 2, hs, hd);
200 FindTracks(fTrSec,hs,hd);
205 //_____________________________________________________________________
206 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) const
209 // Parametrised "expected" error of the cluster reconstruction in Y
212 Double_t s = 0.08 * 0.08;
217 //_____________________________________________________________________
218 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl) const
221 // Parametrised "expected" error of the cluster reconstruction in Z
224 Double_t s = 6 * 6 /12.;
229 //_____________________________________________________________________
230 Double_t f1trd(Double_t x1,Double_t y1,
231 Double_t x2,Double_t y2,
232 Double_t x3,Double_t y3)
235 // Initial approximation of the track curvature
236 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
239 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
240 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
241 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
242 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
243 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
245 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
247 return -xr*yr/sqrt(xr*xr+yr*yr);
251 //_____________________________________________________________________
252 Double_t f2trd(Double_t x1,Double_t y1,
253 Double_t x2,Double_t y2,
254 Double_t x3,Double_t y3)
257 // Initial approximation of the track curvature times center of curvature
258 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
261 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
262 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
263 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
264 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
265 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
267 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
269 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
273 //_____________________________________________________________________
274 Double_t f3trd(Double_t x1,Double_t y1,
275 Double_t x2,Double_t y2,
276 Double_t z1,Double_t z2)
279 // Initial approximation of the tangent of the track dip angle
280 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
283 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
288 //___________________________________________________________________
290 Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec,
291 Int_t s, Int_t rf, Int_t matchedIndex,
295 // Starting from current position on track=t this function tries
296 // to extrapolate the track up to timeBin=rf and to confirm prolongation
297 // if a close cluster is found. *sec is a pointer to allocated
298 // array of sectors, in which the initial sector has index=s.
304 // Bool_t good_match;
306 const Int_t kTimeBinsToSkip=Int_t(fgkSkipDepth*sec->GetNtimeBins());
307 Int_t tryAgain=kTimeBinsToSkip;
309 Double_t alpha=AliTRDgeometry::GetAlpha();
311 Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5);
315 for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) {
317 Double_t x=sec->GetX(nr);
318 Double_t ymax=x*TMath::Tan(0.5*alpha);
320 rho = 0.00295; x0 = 11.09; // TEC
321 if(sec->TECframe(nr,t.GetY(),t.GetZ())) {
322 rho = 1.7; x0 = 33.0; // G10 frame
324 if(TMath::Abs(x - t.GetX()) > 3) {
325 rho = 0.0559; x0 = 55.6; // radiator
327 if (!t.PropagateTo(x,x0,rho,0.139)) {
328 cerr<<"Can't propagate to x = "<<x<<endl;
332 AliTRDtimeBin& timeBin=sec[s][nr];
333 Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
334 Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
335 Double_t road=25.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
337 if (road>fgkWideRoad) {
338 if (t.GetNclusters() > 4) {
339 cerr<<t.GetNclusters()<<" FindProlongation: Road is too wide !\n";
348 Double_t maxChi2=fgkMaxChi2;
352 for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
353 AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
355 // good_match = kFALSE;
356 // if((c->GetTrackIndex(0) == matchedIndex) ||
357 // (c->GetTrackIndex(1) == matchedIndex) ||
358 // (c->GetTrackIndex(2) == matchedIndex)) good_match = kTRUE;
359 // if(hsame) hsame->Fill(TMath::Abs(c->GetY()-y)/road);
360 // if(hdiff) hdiff->Fill(road);
362 if (c->GetY() > y+road) break;
363 if (c->IsUsed() > 0) continue;
365 // if(good_match) hsame->Fill(TMath::Abs(c->GetZ()-z));
366 // else hdiff->Fill(TMath::Abs(c->GetZ()-z));
368 // if(!good_match) continue;
370 if((c->GetZ()-z)*(c->GetZ()-z) > 3 * 12 * sz2) continue;
372 Double_t chi2=t.GetPredictedChi2(c);
374 // if((c->GetTrackIndex(0) == matchedIndex) ||
375 // (c->GetTrackIndex(1) == matchedIndex) ||
376 // (c->GetTrackIndex(2) == matchedIndex))
377 // hdiff->Fill(chi2);
381 if (chi2 > maxChi2) continue;
384 index=timeBin.GetIndex(i);
389 for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
390 AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
392 if (c->GetY() > y+road) break;
393 if (c->IsUsed() > 0) continue;
394 if((c->GetZ()-z)*(c->GetZ()-z) > 3.25 * 12 * sz2) continue;
396 Double_t chi2=t.GetPredictedChi2(c);
400 if (chi2 > maxChi2) continue;
403 index=timeBin.GetIndex(i);
411 t.Update(cl,maxChi2,index);
413 tryAgain=kTimeBinsToSkip;
415 if (tryAgain==0) break;
417 // cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl;
419 if (!t.Rotate(alpha)) {
420 cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
423 } else if (y <-ymax) {
424 // cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl;
426 if (!t.Rotate(-alpha)) {
427 cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
431 if(!sec->TECframe(nr,y,z)) tryAgain--;
440 //_____________________________________________________________________________
441 void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile)
444 // Opens a ROOT-file with TRD-clusters and reads the cluster-tree in
447 ReadClusters(fClusters, clusterfile);
449 // get geometry from the file with hits
451 TFile *fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(hitfile);
453 printf("AliTRDtracker::Open -- ");
454 printf("Open the ALIROOT-file %s.\n",hitfile);
455 fInputFile = new TFile(hitfile,"UPDATE");
458 printf("AliTRDtracker::Open -- ");
459 printf("%s is already open.\n",hitfile);
462 // Get AliRun object from file or create it if not on file
464 gAlice = (AliRun*) fInputFile->Get("gAlice");
466 printf("AliTRDtracker::GetEvent -- ");
467 printf("AliRun object found on file.\n");
470 printf("AliTRDtracker::GetEvent -- ");
471 printf("Could not find AliRun object.\n");
475 // Import the Trees for the event nEvent in the file
476 Int_t nparticles = gAlice->GetEvent(fEvent);
477 cerr<<"nparticles = "<<nparticles<<endl;
479 if (nparticles <= 0) {
480 printf("AliTRDtracker::GetEvent -- ");
481 printf("No entries in the trees for event %d.\n",fEvent);
485 AliTRD *trd = (AliTRD*) gAlice->GetDetector("TRD");
486 fGeom = trd->GetGeometry();
491 //_____________________________________________________________________________
492 void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec)
495 // Fills clusters into TRD tracking_sectors
496 // Note that the numbering scheme for the TRD tracking_sectors
497 // differs from that of TRD sectors
500 for (Int_t i=0; i<AliTRDgeometry::Nsect(); i++) sec[i].SetUp();
502 // Sort clusters into AliTRDtimeBin's within AliTRDtrackSector's
504 cerr<<"SetUpSectors: sorting clusters"<<endl;
506 Int_t ncl=fClusters->GetEntriesFast();
509 printf("\r %d left ",ncl);
510 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
511 Int_t detector=c->GetDetector(), localTimeBin=c->GetLocalTimeBin();
512 Int_t sector=fGeom->GetSector(detector);
514 Int_t trackingSector = AliTRDgeometry::kNsect - sector - 1;
516 Int_t tb=sec[sector].GetTimeBin(detector,localTimeBin);
518 sec[trackingSector][tb].InsertCluster(c,index);
525 //_____________________________________________________________________________
526 void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer,
527 AliTRDtrackingSector* fTrSec, Int_t turn,
531 // Creates track seeds using clusters in timeBins=i1,i2
534 Int_t i2 = inner, i1 = outer;
536 Int_t nprim = 85600/2;
541 Bool_t match = false;
546 Double_t x[5], c[15];
547 Int_t maxSec=AliTRDgeometry::kNsect;
549 Double_t alpha=AliTRDgeometry::GetAlpha();
550 Double_t shift=AliTRDgeometry::GetAlpha()/2.;
551 Double_t cs=cos(alpha), sn=sin(alpha);
552 Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
554 Double_t x1 =fTrSec[0].GetX(i1);
555 Double_t xx2=fTrSec[0].GetX(i2);
560 if((turn != 1)&&(turn != 2)) {
561 printf("*** Error in MakeSeeds: unexpected turn = %d \n", turn);
566 for (Int_t ns=0; ns<maxSec; ns++) {
568 printf("\n MakeSeeds: sector %d \n", ns);
570 Int_t nl2=fTrSec[(ns-2+maxSec)%maxSec][i2];
571 Int_t nl=fTrSec[(ns-1+maxSec)%maxSec][i2];
572 Int_t nm=fTrSec[ns][i2];
573 Int_t nu=fTrSec[(ns+1)%maxSec][i2];
574 Int_t nu2=fTrSec[(ns+2)%maxSec][i2];
576 AliTRDtimeBin& r1=fTrSec[ns][i1];
578 for (Int_t is=0; is < r1; is++) {
579 Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
580 for(Int_t ii=0; ii<3; ii++) to[ii] = r1[is]->GetTrackIndex(ii);
582 for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
584 const AliTRDcluster *cl;
586 Double_t x3=0., y3=0.;
589 if(turn != 2) continue;
590 AliTRDtimeBin& r2=fTrSec[(ns-2+maxSec)%maxSec][i2];
592 y2=cl->GetY(); z2=cl->GetZ();
593 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
598 else if (js<nl2+nl) {
599 if(turn != 1) continue;
600 AliTRDtimeBin& r2=fTrSec[(ns-1+maxSec)%maxSec][i2];
602 y2=cl->GetY(); z2=cl->GetZ();
603 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
609 else if (js<nl2+nl+nm) {
610 if(turn != 1) continue;
611 AliTRDtimeBin& r2=fTrSec[ns][i2];
613 x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
614 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
616 else if (js<nl2+nl+nm+nu) {
617 if(turn != 1) continue;
618 AliTRDtimeBin& r2=fTrSec[(ns+1)%maxSec][i2];
620 y2=cl->GetY(); z2=cl->GetZ();
621 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
628 if(turn != 2) continue;
629 AliTRDtimeBin& r2=fTrSec[(ns+2)%maxSec][i2];
630 cl=r2[js-nl2-nl-nm-nu];
631 y2=cl->GetY(); z2=cl->GetZ();
632 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
641 for (Int_t ii=0; ii<3; ii++) {
642 // cerr<<"ti["<<ii<<"] = "<<ti[ii]<<"; to["<<ii<<"] = "<<to[ii]<<endl;
643 if(ti[ii] < 0) continue;
644 if(ti[ii] >= nprim) continue;
645 for (Int_t kk=0; kk<3; kk++) {
646 if(to[kk] < 0) continue;
647 if(to[kk] >= nprim) continue;
648 if(ti[ii] == to[kk]) {
649 //cerr<<"ti["<<ii<<"] = "<<ti[ii]<<" = "<<to[kk]<<" = to["<<kk<<"]"<<endl;
650 matchedIndex = ti[ii];
656 if(TMath::Abs(z1-z2) > fgkMaxSeedDeltaZ12) continue;
658 Double_t zz=z1 - z1/x1*(x1-x2);
660 if (TMath::Abs(zz-z2)>fgkMaxSeedDeltaZ) continue;
662 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
663 if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
667 x[2]=f1trd(x1,y1,x2,y2,x3,y3);
669 if (TMath::Abs(x[2]) > fgkMaxSeedC) continue;
671 x[3]=f2trd(x1,y1,x2,y2,x3,y3);
673 if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue;
675 x[4]=f3trd(x1,y1,x2,y2,z1,z2);
677 if (TMath::Abs(x[4]) > fgkMaxSeedTan) continue;
679 Double_t a=asin(x[3]);
680 Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3]));
682 if (TMath::Abs(zv)>fgkMaxSeedVertexZ) continue;
684 Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
685 Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2();
686 Double_t sy3=fgkSeedErrorSY3, sy=fgkSeedErrorSY, sz=fgkSeedErrorSZ;
688 Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
689 Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
690 Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
691 Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy;
692 Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy;
693 Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy;
694 Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy;
695 Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz;
696 Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy;
697 Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz;
701 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24;
702 c[6]=f30*sy1; c[7]=0.; c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24;
703 c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34;
704 c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22;
705 c[13]=f40*sy1*f30+f42*sy2*f32;
706 c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43;
708 UInt_t index=r1.GetIndex(is);
710 AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
712 Int_t rc=FindProlongation(*track,fTrSec,ns,i2,matchedIndex,hsame,hdiff);
714 // if (match) hsame->Fill((Float_t) track->GetNclusters());
715 // else hdiff->Fill((Float_t) track->GetNclusters());
720 (track->GetNclusters() < (i1-i2)*fgkMinClustersInSeed)) delete track;
722 fSeeds->AddLast(track); fNseeds++;
723 printf("\r found seed %d ", fNseeds);
732 //_____________________________________________________________________________
733 void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename)
736 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
737 // from the file. The names of the cluster tree and branches
738 // should match the ones used in AliTRDclusterizer::WriteClusters()
741 TDirectory *savedir=gDirectory;
743 TFile *file = TFile::Open(filename);
744 if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;}
747 sprintf(treeName,"TreeR%d_TRD",fEvent);
748 TTree *clusterTree = (TTree*) file->Get(treeName);
750 TObjArray *clusterArray = new TObjArray(400);
752 clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
754 Int_t nEntries = (Int_t) clusterTree->GetEntries();
755 printf("found %d entries in %s.\n",nEntries,clusterTree->GetName());
757 // Loop through all entries in the tree
759 AliTRDcluster *c = 0;
761 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
764 nbytes += clusterTree->GetEvent(iEntry);
766 // Get the number of points in the detector
767 Int_t nCluster = clusterArray->GetEntriesFast();
768 printf("Read %d clusters from entry %d \n", nCluster, iEntry);
770 // Loop through all TRD digits
771 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
772 c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
773 AliTRDcluster *co = new AliTRDcluster(*c);
774 co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
776 delete clusterArray->RemoveAt(iCluster);
786 //___________________________________________________________________
787 void AliTRDtracker::FindTracks(AliTRDtrackingSector* fTrSec, TH1F *hs, TH1F *hd)
790 // Finds tracks in TRD
796 Int_t numOfTimeBins = fTrSec[0].GetNtimeBins();
797 Int_t nseed=fSeeds->GetEntriesFast();
800 for (Int_t i=0; i<nseed; i++) {
801 cerr<<"FindTracks: seed "<<i+1<<" out of "<<nseed<<endl;
803 AliTRDtrack& t=*((AliTRDtrack*)fSeeds->UncheckedAt(i));
805 nSeedClusters = t.GetNclusters();
806 Double_t alpha=t.GetAlpha();
808 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
809 if (alpha < 0. ) alpha += 2.*TMath::Pi();
810 Int_t ns=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
812 Int_t label = GetTrackLabel(t);
814 if (FindProlongation(t,fTrSec,ns,0,label,hsame,hdiff)) {
815 cerr<<"No of clusters in the track = "<<t.GetNclusters()<<endl;
816 if (t.GetNclusters() >= Int_t(fgkMinClustersInTrack*numOfTimeBins)) {
817 Int_t label = GetTrackLabel(t);
822 AliTRDtrack *pt = new AliTRDtrack(t);
823 fTracks->AddLast(pt); fNtracks++;
825 cerr<<"found track "<<fNtracks<<endl;
828 delete fSeeds->RemoveAt(i);
833 //__________________________________________________________________
834 void AliTRDtracker::UseClusters(AliTRDtrack t)
840 Int_t ncl=t.GetNclusters();
841 for (Int_t i=0; i<ncl; i++) {
842 Int_t index = t.GetClusterIndex(i);
843 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
849 //__________________________________________________________________
850 Int_t AliTRDtracker::GetTrackLabel(AliTRDtrack t)
856 Int_t label=123456789, index, i, j;
857 Int_t ncl=t.GetNclusters();
858 const Int_t kRange = AliTRDgeometry::kNplan * fGeom->GetTimeMax();
861 // Int_t s[kRange][2];
862 Int_t **s = new Int_t* [kRange];
863 for (i=0; i<kRange; i++) {
866 for (i=0; i<kRange; i++) {
872 for (i=0; i<ncl; i++) {
873 index=t.GetClusterIndex(i);
874 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
875 t0=c->GetTrackIndex(0);
876 t1=c->GetTrackIndex(1);
877 t2=c->GetTrackIndex(2);
880 for (i=0; i<ncl; i++) {
881 index=t.GetClusterIndex(i);
882 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
883 for (Int_t k=0; k<3; k++) {
884 label=c->GetTrackIndex(k);
885 labelAdded=kFALSE; j=0;
887 while ( (!labelAdded) && ( j < kRange ) ) {
888 if (s[j][0]==label || s[j][1]==0) {
902 for (i=0; i<kRange; i++) {
904 max=s[i][1]; label=s[i][0];
908 if(max > ncl*fgkLabelFraction) return label;
912 //___________________________________________________________________
913 Int_t AliTRDtracker::WriteTracks(const Char_t *filename)
916 // Write the tracks to the output file
919 TDirectory *savedir=gDirectory;
921 //TFile *out=TFile::Open(filename,"RECREATE");
922 TFile *out = (TFile*) gROOT->GetListOfFiles()->FindObject(filename);
924 printf("AliTRDtracker::Open -- ");
925 printf("Open the ALIROOT-file %s.\n",filename);
926 out = new TFile(filename,"RECREATE");
929 printf("AliTRDtracker::Open -- ");
930 printf("%s is already open.\n",filename);
934 sprintf(treeName,"TreeT%d_TRD",fEvent);
935 TTree tracktree(treeName,"Tree with TRD tracks");
937 AliTRDtrack *iotrack=0;
938 tracktree.Branch("tracks","AliTRDtrack",&iotrack,32000,0);
940 Int_t ntracks=fTracks->GetEntriesFast();
942 for (Int_t i=0; i<ntracks; i++) {
943 AliTRDtrack *pt=(AliTRDtrack*)fTracks->UncheckedAt(i);
946 cerr<<"WriteTracks: put track "<<i<<" in the tree"<<endl;
954 cerr<<"WriteTracks: done"<<endl;