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.12 2001/05/28 17:07:58 hristov
19 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)
21 Revision 1.8 2000/12/20 13:00:44 cblume
22 Modifications for the HP-compiler
24 Revision 1.7 2000/12/08 16:07:02 cblume
25 Update of the tracking by Sergei
27 Revision 1.6 2000/11/30 17:38:08 cblume
28 Changes to get in line with new STEER and EVGEN
30 Revision 1.5 2000/11/14 14:40:27 cblume
31 Correction for the Sun compiler (kTRUE and kFALSE)
33 Revision 1.4 2000/11/10 14:57:52 cblume
34 Changes in the geometry constants for the DEC compiler
36 Revision 1.3 2000/10/15 23:40:01 cblume
39 Revision 1.2 2000/10/06 16:49:46 cblume
42 Revision 1.1.2.2 2000/10/04 16:34:58 cblume
43 Replace include files by forward declarations
45 Revision 1.1.2.1 2000/09/22 14:47:52 cblume
59 #include "AliTRDgeometry.h"
60 #include "AliTRDrecPoint.h"
61 #include "AliTRDcluster.h"
62 #include "AliTRDtrack.h"
63 #include "AliTRDtrackingSector.h"
64 #include "AliTRDtimeBin.h"
66 #include "AliTRDtracker.h"
68 ClassImp(AliTRDtracker)
70 const Float_t AliTRDtracker::fSeedDepth = 0.5;
71 const Float_t AliTRDtracker::fSeedStep = 0.05;
72 const Float_t AliTRDtracker::fSeedGap = 0.25;
74 const Float_t AliTRDtracker::fMaxSeedDeltaZ12 = 40.;
75 const Float_t AliTRDtracker::fMaxSeedDeltaZ = 25.;
76 const Float_t AliTRDtracker::fMaxSeedC = 0.0052;
77 const Float_t AliTRDtracker::fMaxSeedTan = 1.2;
78 const Float_t AliTRDtracker::fMaxSeedVertexZ = 150.;
80 const Double_t AliTRDtracker::fSeedErrorSY = 0.2;
81 const Double_t AliTRDtracker::fSeedErrorSY3 = 2.5;
82 const Double_t AliTRDtracker::fSeedErrorSZ = 0.1;
84 const Float_t AliTRDtracker::fMinClustersInSeed = 0.7;
86 const Float_t AliTRDtracker::fMinClustersInTrack = 0.5;
87 const Float_t AliTRDtracker::fSkipDepth = 0.05;
88 const Float_t AliTRDtracker::fLabelFraction = 0.5;
89 const Float_t AliTRDtracker::fWideRoad = 20.;
91 const Double_t AliTRDtracker::fMaxChi2 = 24.;
93 //____________________________________________________________________
94 AliTRDtracker::AliTRDtracker()
97 // Default constructor
113 //____________________________________________________________________
114 AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title)
121 fClusters = new TObjArray(2000);
123 fSeeds = new TObjArray(20000);
125 fTracks = new TObjArray(10000);
130 //___________________________________________________________________
131 AliTRDtracker::~AliTRDtracker()
139 //___________________________________________________________________
140 void AliTRDtracker::Clusters2Tracks(TH1F *hs, TH1F *hd)
145 Int_t fTotalNofTimeBins = fGeom->GetTimeMax() * AliTRDgeometry::Nplan();
146 Int_t nSteps = (Int_t) (fSeedDepth / fSeedStep);
147 Int_t gap = (Int_t) (fTotalNofTimeBins * fSeedGap);
148 Int_t step = (Int_t) (fTotalNofTimeBins * fSeedStep);
153 if (!fClusters) return;
155 AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect];
156 SetUpSectors(fTrSec);
158 // make a first turn looking for seed ends in the same (n,n)
159 // and in the adjacent sectors (n,n+1)
161 for(i=0; i<nSteps; i++) {
162 printf("step %d out of %d \n", i+1, nSteps);
163 outer=fTotalNofTimeBins-1-i*step; inner=outer-gap;
164 MakeSeeds(inner,outer, fTrSec, 1, hs, hd);
165 FindTracks(fTrSec, hs, hd);
168 // make a second turn looking for seed ends in next-to-adjacent
171 for(i=0; i<nSteps; i++) {
172 printf("step %d out of %d \n", i+1, nSteps);
173 outer=fTotalNofTimeBins-1-i*step; inner=outer-gap;
174 MakeSeeds(inner, outer, fTrSec, 2, hs, hd);
175 FindTracks(fTrSec,hs,hd);
180 //_____________________________________________________________________
181 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt)
183 // Parametrised "expected" error of the cluster reconstruction in Y
185 Double_t s = 0.08 * 0.08;
189 //_____________________________________________________________________
190 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl)
192 // Parametrised "expected" error of the cluster reconstruction in Z
194 Double_t s = 6 * 6 /12.;
198 //_____________________________________________________________________
199 Double_t f1trd(Double_t x1,Double_t y1,
200 Double_t x2,Double_t y2,
201 Double_t x3,Double_t y3)
203 // Initial approximation of the track curvature
204 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
206 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
207 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
208 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
209 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
210 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
212 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
214 return -xr*yr/sqrt(xr*xr+yr*yr);
217 //_____________________________________________________________________
218 Double_t f2trd(Double_t x1,Double_t y1,
219 Double_t x2,Double_t y2,
220 Double_t x3,Double_t y3)
222 // Initial approximation of the track curvature times center of curvature
223 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
225 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
226 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
227 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
228 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
229 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
231 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
233 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
236 //_____________________________________________________________________
237 Double_t f3trd(Double_t x1,Double_t y1,
238 Double_t x2,Double_t y2,
239 Double_t z1,Double_t z2)
241 // Initial approximation of the tangent of the track dip angle
242 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
244 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
248 //___________________________________________________________________
250 Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec,
251 Int_t s, Int_t rf, Int_t matched_index,
254 // Starting from current position on track=t this function tries
255 // to extrapolate the track up to timeBin=rf and to confirm prolongation
256 // if a close cluster is found. *sec is a pointer to allocated
257 // array of sectors, in which the initial sector has index=s.
263 // Bool_t good_match;
265 const Int_t TIME_BINS_TO_SKIP=Int_t(fSkipDepth*sec->GetNtimeBins());
266 Int_t try_again=TIME_BINS_TO_SKIP;
268 Double_t alpha=AliTRDgeometry::GetAlpha();
270 Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5);
274 for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) {
276 Double_t x=sec->GetX(nr);
277 Double_t ymax=x*TMath::Tan(0.5*alpha);
279 rho = 0.00295; x0 = 11.09; // TEC
280 if(sec->TECframe(nr,t.GetY(),t.GetZ())) {
281 rho = 1.7; x0 = 33.0; // G10 frame
283 if(TMath::Abs(x - t.GetX()) > 3) {
284 rho = 0.0559; x0 = 55.6; // radiator
286 if (!t.PropagateTo(x,x0,rho,0.139)) {
287 cerr<<"Can't propagate to x = "<<x<<endl;
291 AliTRDtimeBin& time_bin=sec[s][nr];
292 Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
293 Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
294 Double_t road=25.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
296 if (road>fWideRoad) {
297 if (t.GetNclusters() > 4) {
298 cerr<<t.GetNclusters()<<" FindProlongation: Road is too wide !\n";
307 Double_t max_chi2=fMaxChi2;
311 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
312 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
314 // good_match = kFALSE;
315 // if((c->GetTrackIndex(0) == matched_index) ||
316 // (c->GetTrackIndex(1) == matched_index) ||
317 // (c->GetTrackIndex(2) == matched_index)) good_match = kTRUE;
318 // if(hsame) hsame->Fill(TMath::Abs(c->GetY()-y)/road);
319 // if(hdiff) hdiff->Fill(road);
321 if (c->GetY() > y+road) break;
322 if (c->IsUsed() > 0) continue;
324 // if(good_match) hsame->Fill(TMath::Abs(c->GetZ()-z));
325 // else hdiff->Fill(TMath::Abs(c->GetZ()-z));
327 // if(!good_match) continue;
329 if((c->GetZ()-z)*(c->GetZ()-z) > 3 * 12 * sz2) continue;
331 Double_t chi2=t.GetPredictedChi2(c);
333 // if((c->GetTrackIndex(0) == matched_index) ||
334 // (c->GetTrackIndex(1) == matched_index) ||
335 // (c->GetTrackIndex(2) == matched_index))
336 // hdiff->Fill(chi2);
340 if (chi2 > max_chi2) continue;
343 index=time_bin.GetIndex(i);
348 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
349 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
351 if (c->GetY() > y+road) break;
352 if (c->IsUsed() > 0) continue;
353 if((c->GetZ()-z)*(c->GetZ()-z) > 3.25 * 12 * sz2) continue;
355 Double_t chi2=t.GetPredictedChi2(c);
359 if (chi2 > max_chi2) continue;
362 index=time_bin.GetIndex(i);
370 t.Update(cl,max_chi2,index);
372 try_again=TIME_BINS_TO_SKIP;
374 if (try_again==0) break;
376 // cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl;
378 if (!t.Rotate(alpha)) {
379 cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
382 } else if (y <-ymax) {
383 // cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl;
385 if (!t.Rotate(-alpha)) {
386 cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
390 if(!sec->TECframe(nr,y,z)) try_again--;
399 //_____________________________________________________________________________
400 void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile)
402 // Opens a ROOT-file with TRD-clusters and reads the cluster-tree in
404 ReadClusters(fClusters, clusterfile);
406 // get geometry from the file with hits
408 TFile *fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(hitfile);
410 printf("AliTRDtracker::Open -- ");
411 printf("Open the ALIROOT-file %s.\n",hitfile);
412 fInputFile = new TFile(hitfile,"UPDATE");
415 printf("AliTRDtracker::Open -- ");
416 printf("%s is already open.\n",hitfile);
419 // Get AliRun object from file or create it if not on file
421 gAlice = (AliRun*) fInputFile->Get("gAlice");
423 printf("AliTRDtracker::GetEvent -- ");
424 printf("AliRun object found on file.\n");
427 printf("AliTRDtracker::GetEvent -- ");
428 printf("Could not find AliRun object.\n");
432 // Import the Trees for the event nEvent in the file
433 Int_t nparticles = gAlice->GetEvent(fEvent);
434 cerr<<"nparticles = "<<nparticles<<endl;
436 if (nparticles <= 0) {
437 printf("AliTRDtracker::GetEvent -- ");
438 printf("No entries in the trees for event %d.\n",fEvent);
442 AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
443 fGeom = TRD->GetGeometry();
448 //_____________________________________________________________________________
449 void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec)
451 // Fills clusters into TRD tracking_sectors
452 // Note that the numbering scheme for the TRD tracking_sectors
453 // differs from that of TRD sectors
455 for (Int_t i=0; i<AliTRDgeometry::Nsect(); i++) sec[i].SetUp();
457 // Sort clusters into AliTRDtimeBin's within AliTRDtrackSector's
459 cerr<<"SetUpSectors: sorting clusters"<<endl;
461 Int_t ncl=fClusters->GetEntriesFast();
464 printf("\r %d left ",ncl);
465 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
466 Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin();
467 Int_t sector=fGeom->GetSector(detector);
469 Int_t tracking_sector = AliTRDgeometry::kNsect - sector - 1;
471 Int_t tb=sec[sector].GetTimeBin(detector,local_time_bin);
473 sec[tracking_sector][tb].InsertCluster(c,index);
480 //_____________________________________________________________________________
481 void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer,
482 AliTRDtrackingSector* fTrSec, Int_t turn,
485 // Creates track seeds using clusters in timeBins=i1,i2
487 Int_t i2 = inner, i1 = outer;
489 Int_t nprim = 85600/2;
494 Bool_t match = false;
499 Double_t x[5], c[15];
500 Int_t max_sec=AliTRDgeometry::kNsect;
502 Double_t alpha=AliTRDgeometry::GetAlpha();
503 Double_t shift=AliTRDgeometry::GetAlpha()/2.;
504 Double_t cs=cos(alpha), sn=sin(alpha);
505 Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
507 Double_t x1 =fTrSec[0].GetX(i1);
508 Double_t xx2=fTrSec[0].GetX(i2);
513 if((turn != 1)&&(turn != 2)) {
514 printf("*** Error in MakeSeeds: unexpected turn = %d \n", turn);
519 for (Int_t ns=0; ns<max_sec; ns++) {
521 printf("\n MakeSeeds: sector %d \n", ns);
523 Int_t nl2=fTrSec[(ns-2+max_sec)%max_sec][i2];
524 Int_t nl=fTrSec[(ns-1+max_sec)%max_sec][i2];
525 Int_t nm=fTrSec[ns][i2];
526 Int_t nu=fTrSec[(ns+1)%max_sec][i2];
527 Int_t nu2=fTrSec[(ns+2)%max_sec][i2];
529 AliTRDtimeBin& r1=fTrSec[ns][i1];
531 for (Int_t is=0; is < r1; is++) {
532 Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
533 for(Int_t ii=0; ii<3; ii++) to[ii] = r1[is]->GetTrackIndex(ii);
535 for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
537 const AliTRDcluster *cl;
539 Double_t x3=0., y3=0.;
542 if(turn != 2) continue;
543 AliTRDtimeBin& r2=fTrSec[(ns-2+max_sec)%max_sec][i2];
545 y2=cl->GetY(); z2=cl->GetZ();
546 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
551 else if (js<nl2+nl) {
552 if(turn != 1) continue;
553 AliTRDtimeBin& r2=fTrSec[(ns-1+max_sec)%max_sec][i2];
555 y2=cl->GetY(); z2=cl->GetZ();
556 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
562 else if (js<nl2+nl+nm) {
563 if(turn != 1) continue;
564 AliTRDtimeBin& r2=fTrSec[ns][i2];
566 x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
567 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
569 else if (js<nl2+nl+nm+nu) {
570 if(turn != 1) continue;
571 AliTRDtimeBin& r2=fTrSec[(ns+1)%max_sec][i2];
573 y2=cl->GetY(); z2=cl->GetZ();
574 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
581 if(turn != 2) continue;
582 AliTRDtimeBin& r2=fTrSec[(ns+2)%max_sec][i2];
583 cl=r2[js-nl2-nl-nm-nu];
584 y2=cl->GetY(); z2=cl->GetZ();
585 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
594 for (Int_t ii=0; ii<3; ii++) {
595 // cerr<<"ti["<<ii<<"] = "<<ti[ii]<<"; to["<<ii<<"] = "<<to[ii]<<endl;
596 if(ti[ii] < 0) continue;
597 if(ti[ii] >= nprim) continue;
598 for (Int_t kk=0; kk<3; kk++) {
599 if(to[kk] < 0) continue;
600 if(to[kk] >= nprim) continue;
601 if(ti[ii] == to[kk]) {
602 //cerr<<"ti["<<ii<<"] = "<<ti[ii]<<" = "<<to[kk]<<" = to["<<kk<<"]"<<endl;
603 matched_index = ti[ii];
609 if(TMath::Abs(z1-z2) > fMaxSeedDeltaZ12) continue;
611 Double_t zz=z1 - z1/x1*(x1-x2);
613 if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue;
615 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
616 if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
620 x[2]=f1trd(x1,y1,x2,y2,x3,y3);
622 if (TMath::Abs(x[2]) > fMaxSeedC) continue;
624 x[3]=f2trd(x1,y1,x2,y2,x3,y3);
626 if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue;
628 x[4]=f3trd(x1,y1,x2,y2,z1,z2);
630 if (TMath::Abs(x[4]) > fMaxSeedTan) continue;
632 Double_t a=asin(x[3]);
633 Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3]));
635 if (TMath::Abs(zv)>fMaxSeedVertexZ) continue;
637 Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
638 Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2();
639 Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ;
641 Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
642 Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
643 Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
644 Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy;
645 Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy;
646 Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy;
647 Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy;
648 Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz;
649 Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy;
650 Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz;
654 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24;
655 c[6]=f30*sy1; c[7]=0.; c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24;
656 c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34;
657 c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22;
658 c[13]=f40*sy1*f30+f42*sy2*f32;
659 c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43;
661 UInt_t index=r1.GetIndex(is);
663 AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
665 Int_t rc=FindProlongation(*track,fTrSec,ns,i2,matched_index,hsame,hdiff);
667 // if (match) hsame->Fill((Float_t) track->GetNclusters());
668 // else hdiff->Fill((Float_t) track->GetNclusters());
673 (track->GetNclusters() < (i1-i2)*fMinClustersInSeed)) delete track;
675 fSeeds->AddLast(track); fNseeds++;
676 printf("\r found seed %d ", fNseeds);
685 //_____________________________________________________________________________
686 void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename)
689 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
690 // from the file. The names of the cluster tree and branches
691 // should match the ones used in AliTRDclusterizer::WriteClusters()
694 TDirectory *savedir=gDirectory;
696 TFile *file = TFile::Open(filename);
697 if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;}
699 TTree *ClusterTree = (TTree*)file->Get("ClusterTree");
701 TObjArray *ClusterArray = new TObjArray(400);
703 ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray);
705 Int_t nEntries = (Int_t) ClusterTree->GetEntries();
706 printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
708 // Loop through all entries in the tree
710 AliTRDcluster *c = 0;
712 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
715 nbytes += ClusterTree->GetEvent(iEntry);
717 // Get the number of points in the detector
718 Int_t nCluster = ClusterArray->GetEntriesFast();
719 printf("Read %d clusters from entry %d \n", nCluster, iEntry);
721 // Loop through all TRD digits
722 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
723 c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
724 AliTRDcluster *co = new AliTRDcluster(*c);
725 co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
727 delete ClusterArray->RemoveAt(iCluster);
737 //___________________________________________________________________
738 void AliTRDtracker::FindTracks(AliTRDtrackingSector* fTrSec, TH1F *hs, TH1F *hd)
741 // Finds tracks in TRD
747 Int_t num_of_time_bins = fTrSec[0].GetNtimeBins();
748 Int_t nseed=fSeeds->GetEntriesFast();
751 for (Int_t i=0; i<nseed; i++) {
752 cerr<<"FindTracks: seed "<<i+1<<" out of "<<nseed<<endl;
754 AliTRDtrack& t=*((AliTRDtrack*)fSeeds->UncheckedAt(i));
756 nSeedClusters = t.GetNclusters();
757 Double_t alpha=t.GetAlpha();
759 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
760 if (alpha < 0. ) alpha += 2.*TMath::Pi();
761 Int_t ns=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
763 Int_t label = GetTrackLabel(t);
765 if (FindProlongation(t,fTrSec,ns,0,label,hsame,hdiff)) {
766 cerr<<"No of clusters in the track = "<<t.GetNclusters()<<endl;
767 if (t.GetNclusters() >= Int_t(fMinClustersInTrack*num_of_time_bins)) {
768 Int_t label = GetTrackLabel(t);
773 AliTRDtrack *pt = new AliTRDtrack(t);
774 fTracks->AddLast(pt); fNtracks++;
776 cerr<<"found track "<<fNtracks<<endl;
779 delete fSeeds->RemoveAt(i);
784 //__________________________________________________________________
785 void AliTRDtracker::UseClusters(AliTRDtrack t) {
786 Int_t ncl=t.GetNclusters();
787 for (Int_t i=0; i<ncl; i++) {
788 Int_t index = t.GetClusterIndex(i);
789 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
794 //__________________________________________________________________
795 Int_t AliTRDtracker::GetTrackLabel(AliTRDtrack t) {
797 Int_t label=123456789, index, i, j;
798 Int_t ncl=t.GetNclusters();
799 const Int_t range = AliTRDgeometry::kNplan * fGeom->GetTimeMax();
802 // Int_t s[range][2];
803 Int_t **s = new Int_t* [range];
804 for (i=0; i<range; i++) {
807 for (i=0; i<range; i++) {
813 for (i=0; i<ncl; i++) {
814 index=t.GetClusterIndex(i);
815 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
816 t0=c->GetTrackIndex(0);
817 t1=c->GetTrackIndex(1);
818 t2=c->GetTrackIndex(2);
821 for (i=0; i<ncl; i++) {
822 index=t.GetClusterIndex(i);
823 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
824 for (Int_t k=0; k<3; k++) {
825 label=c->GetTrackIndex(k);
826 label_added=kFALSE; j=0;
828 while ( (!label_added) && ( j < range ) ) {
829 if (s[j][0]==label || s[j][1]==0) {
843 for (i=0; i<range; i++) {
845 max=s[i][1]; label=s[i][0];
849 if(max > ncl*fLabelFraction) return label;
853 //___________________________________________________________________
855 Int_t AliTRDtracker::WriteTracks(const Char_t *filename) {
857 TDirectory *savedir=gDirectory;
859 //TFile *out=TFile::Open(filename,"RECREATE");
860 TFile *out = (TFile*) gROOT->GetListOfFiles()->FindObject(filename);
862 printf("AliTRDtracker::Open -- ");
863 printf("Open the ALIROOT-file %s.\n",filename);
864 out = new TFile(filename,"RECREATE");
867 printf("AliTRDtracker::Open -- ");
868 printf("%s is already open.\n",filename);
871 TTree tracktree("TreeT","Tree with TRD tracks");
873 AliTRDtrack *iotrack=0;
874 tracktree.Branch("tracks","AliTRDtrack",&iotrack,32000,0);
876 Int_t ntracks=fTracks->GetEntriesFast();
878 for (Int_t i=0; i<ntracks; i++) {
879 AliTRDtrack *pt=(AliTRDtrack*)fTracks->UncheckedAt(i);
882 cerr<<"WriteTracks: put track "<<i<<" in the tree"<<endl;
890 cerr<<"WriteTracks: done"<<endl;