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.8 2000/12/20 13:00:44 cblume
19 Modifications for the HP-compiler
21 Revision 1.7 2000/12/08 16:07:02 cblume
22 Update of the tracking by Sergei
24 Revision 1.6 2000/11/30 17:38:08 cblume
25 Changes to get in line with new STEER and EVGEN
27 Revision 1.5 2000/11/14 14:40:27 cblume
28 Correction for the Sun compiler (kTRUE and kFALSE)
30 Revision 1.4 2000/11/10 14:57:52 cblume
31 Changes in the geometry constants for the DEC compiler
33 Revision 1.3 2000/10/15 23:40:01 cblume
36 Revision 1.2 2000/10/06 16:49:46 cblume
39 Revision 1.1.2.2 2000/10/04 16:34:58 cblume
40 Replace include files by forward declarations
42 Revision 1.1.2.1 2000/09/22 14:47:52 cblume
56 #include "AliTRDgeometry.h"
57 #include "AliTRDrecPoint.h"
58 #include "AliTRDcluster.h"
59 #include "AliTRDtrack.h"
60 #include "AliTRDtrackingSector.h"
61 #include "AliTRDtimeBin.h"
63 #include "AliTRDtracker.h"
65 ClassImp(AliTRDtracker)
67 const Int_t AliTRDtracker::fSeedGap = 35;
68 const Int_t AliTRDtracker::fSeedStep = 5;
71 const Float_t AliTRDtracker::fMinClustersInTrack = 0.5;
72 const Float_t AliTRDtracker::fMinClustersInSeed = 0.5;
73 const Float_t AliTRDtracker::fSeedDepth = 0.5;
74 const Float_t AliTRDtracker::fSkipDepth = 0.2;
75 const Float_t AliTRDtracker::fMaxSeedDeltaZ = 30.;
76 const Float_t AliTRDtracker::fMaxSeedC = 0.01;
77 const Float_t AliTRDtracker::fMaxSeedTan = 1.2;
78 const Float_t AliTRDtracker::fMaxSeedVertexZ = 200.;
79 const Float_t AliTRDtracker::fLabelFraction = 0.5;
80 const Float_t AliTRDtracker::fWideRoad = 30.;
82 const Double_t AliTRDtracker::fMaxChi2 = 12.;
83 const Double_t AliTRDtracker::fSeedErrorSY = 0.1;
84 const Double_t AliTRDtracker::fSeedErrorSY3 = 2.5;
85 const Double_t AliTRDtracker::fSeedErrorSZ = 0.1;
87 //____________________________________________________________________
88 AliTRDtracker::AliTRDtracker()
91 // Default constructor
106 //____________________________________________________________________
107 AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title)
114 fClusters = new TObjArray(2000);
116 fSeeds = new TObjArray(20000);
118 fTracks = new TObjArray(10000);
122 //___________________________________________________________________
123 AliTRDtracker::~AliTRDtracker()
131 //___________________________________________________________________
132 void AliTRDtracker::Clusters2Tracks()
135 Int_t fTotalNofTimeBins = fGeom->GetTimeMax() * AliTRDgeometry::Nplan();
136 Int_t nSteps = (Int_t) (fTotalNofTimeBins * fSeedDepth)/fSeedStep;
138 for(Int_t i=0; i<nSteps; i++) {
139 printf("step %d out of %d \n", i+1, nSteps);
140 outer=fTotalNofTimeBins-1-i*fSeedStep; inner=outer-fSeedGap;
141 MakeSeeds(inner,outer);
146 //_____________________________________________________________________
147 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt)
149 // Parametrised "expected" error of the cluster reconstruction in Y
155 //_____________________________________________________________________
156 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl)
158 // Parametrised "expected" error of the cluster reconstruction in Z
160 // Take an example pad size for the time being, check back
162 Double_t s, pad = fGeom->GetRowPadSize(0,2,0);
167 //_____________________________________________________________________
168 inline Double_t f1trd(Double_t x1,Double_t y1,
169 Double_t x2,Double_t y2,
170 Double_t x3,Double_t y3)
172 // Initial approximation of the track curvature
173 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
175 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
176 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
177 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
178 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
179 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
181 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
183 return -xr*yr/sqrt(xr*xr+yr*yr);
186 //_____________________________________________________________________
187 inline Double_t f2trd(Double_t x1,Double_t y1,
188 Double_t x2,Double_t y2,
189 Double_t x3,Double_t y3)
191 // Initial approximation of the track curvature times center of curvature
192 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
194 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
195 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
196 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
197 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
198 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
200 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
202 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
205 //_____________________________________________________________________
206 inline Double_t f3trd(Double_t x1,Double_t y1,
207 Double_t x2,Double_t y2,
208 Double_t z1,Double_t z2)
210 // Initial approximation of the tangent of the track dip angle
211 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
213 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
217 //___________________________________________________________________
219 Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec,
222 // Starting from current position on track=t this function tries
223 // to extrapolate the track up to timeBin=rf and to confirm prolongation
224 // if a close cluster is found. *sec is a pointer to allocated
225 // array of sectors, in which the initial sector has index=s.
227 const Int_t TIME_BINS_TO_SKIP=Int_t(fSkipDepth*sec->GetNtimeBins());
228 Int_t try_again=TIME_BINS_TO_SKIP;
230 Double_t alpha=AliTRDgeometry::GetAlpha();
232 Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5);
234 for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) {
235 Double_t x=sec->GetX(nr), ymax=sec->GetMaxY(nr);
236 if (!t.PropagateTo(x)) {
237 cerr<<"Can't propagate to x = "<<x<<endl;
244 Double_t max_chi2=fMaxChi2;
246 AliTRDtimeBin& time_bin=sec[s][nr];
247 Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
248 Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
249 Double_t road=5.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
251 if (road>fWideRoad) {
252 if (t.GetNclusters() > 4) {
253 cerr<<t.GetNclusters()<<" FindProlongation: Road is too wide !\n";
259 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
260 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
261 if (c->GetY() > y+road) break;
262 if (c->IsUsed() > 0) continue;
264 if((c->GetZ()-z)*(c->GetZ()-z) > 25. + sz2) continue;
266 Double_t chi2=t.GetPredictedChi2(c);
268 if (chi2 > max_chi2) continue;
271 index=time_bin.GetIndex(i);
276 // Float_t l=sec->GetPitch();
277 // t.SetSampledEdx(cl->fQ/l,Int_t(t));
279 t.Update(cl,max_chi2,index);
281 try_again=TIME_BINS_TO_SKIP;
283 if (try_again==0) break;
285 cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl;
287 if (!t.Rotate(alpha)) {
288 cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
291 } else if (y <-ymax) {
292 cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl;
294 if (!t.Rotate(-alpha)) {
295 cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
308 //_____________________________________________________________________________
309 void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile)
311 // Opens a ROOT-file with TRD-clusters and reads in the cluster-tree
313 ReadClusters(fClusters, clusterfile);
315 // get geometry from the file with hits
317 TFile *fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(hitfile);
319 printf("AliTRDtracker::Open -- ");
320 printf("Open the ALIROOT-file %s.\n",hitfile);
321 fInputFile = new TFile(hitfile);
324 printf("AliTRDtracker::Open -- ");
325 printf("%s is already open.\n",hitfile);
328 // Get AliRun object from file or create it if not on file
330 gAlice = (AliRun*) fInputFile->Get("gAlice");
332 printf("AliTRDtracker::GetEvent -- ");
333 printf("AliRun object found on file.\n");
336 printf("AliTRDtracker::GetEvent -- ");
337 printf("Could not find AliRun object.\n");
340 // Import the Trees for the event nEvent in the file
341 Int_t nparticles = gAlice->GetEvent(fEvent);
342 cerr<<"nparticles = "<<nparticles<<endl;
344 if (nparticles <= 0) {
345 printf("AliTRDtracker::GetEvent -- ");
346 printf("No entries in the trees for event %d.\n",fEvent);
349 AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
350 fGeom = TRD->GetGeometry();
355 //_____________________________________________________________________________
356 void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec)
358 // Fills clusters into TRD tracking_sectors
359 // Note that the numbering scheme for the TRD tracking_sectors
360 // differs from that of TRD sectors
362 for (Int_t i=0; i<AliTRDgeometry::Nsect(); i++) sec[i].SetUp();
364 // Sort clusters into AliTRDtimeBin's within AliTRDtrackSector's
366 cerr<<"MakeSeeds: sorting clusters"<<endl;
368 Int_t ncl=fClusters->GetEntriesFast();
371 printf("\r %d left",ncl);
372 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
373 Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin();
374 Int_t sector=fGeom->GetSector(detector);
376 Int_t tracking_sector=sector;
377 if(sector > 0) tracking_sector = AliTRDgeometry::kNsect - sector;
379 Int_t tb=sec[sector].GetTimeBin(detector,local_time_bin);
381 sec[tracking_sector][tb].InsertCluster(c,index);
387 //_____________________________________________________________________________
388 void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer)
390 // Creates track seeds using clusters in timeBins=i1,i2
392 Int_t i2 = inner, i1 = outer;
394 if (!fClusters) return;
396 AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect];
397 SetUpSectors(fTrSec);
401 Double_t x[5], c[15];
402 Int_t max_sec=AliTRDgeometry::kNsect;
404 Double_t alpha=AliTRDgeometry::GetAlpha();
405 Double_t shift=AliTRDgeometry::GetAlpha()/2.;
406 Double_t cs=cos(alpha), sn=sin(alpha);
408 Double_t x1 =fTrSec[0].GetX(i1);
409 Double_t xx2=fTrSec[0].GetX(i2);
411 for (Int_t ns=0; ns<max_sec; ns++) {
413 Int_t nl=fTrSec[(ns-1+max_sec)%max_sec][i2];
414 Int_t nm=fTrSec[ns][i2];
415 Int_t nu=fTrSec[(ns+1)%max_sec][i2];
417 AliTRDtimeBin& r1=fTrSec[ns][i1];
419 for (Int_t is=0; is < r1; is++) {
420 Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
422 for (Int_t js=0; js < nl+nm+nu; js++) {
423 const AliTRDcluster *cl;
425 Double_t x3=0.,y3=0.;
428 AliTRDtimeBin& r2=fTrSec[(ns-1+max_sec)%max_sec][i2];
430 y2=cl->GetY(); z2=cl->GetZ();
437 AliTRDtimeBin& r2=fTrSec[ns][i2];
439 x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
441 AliTRDtimeBin& r2=fTrSec[(ns+1)%max_sec][i2];
443 y2=cl->GetY(); z2=cl->GetZ();
450 Double_t zz=z1 - z1/x1*(x1-x2);
452 if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue;
454 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
455 if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
459 x[2]=f1trd(x1,y1,x2,y2,x3,y3);
461 if (TMath::Abs(x[2]) >= fMaxSeedC) continue;
463 x[3]=f2trd(x1,y1,x2,y2,x3,y3);
465 if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue;
467 x[4]=f3trd(x1,y1,x2,y2,z1,z2);
469 if (TMath::Abs(x[4]) > fMaxSeedTan) continue;
471 Double_t a=asin(x[3]);
472 Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3]));
473 if (TMath::Abs(zv)>fMaxSeedVertexZ) continue;
475 Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
476 Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2();
477 Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ;
479 Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
480 Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
481 Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
482 Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy;
483 Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy;
484 Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy;
485 Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy;
486 Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz;
487 Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy;
488 Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz;
492 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24;
493 c[6]=f30*sy1; c[7]=0.; c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24;
494 c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34;
495 c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22;
496 c[13]=f40*sy1*f30+f42*sy2*f32;
497 c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43;
499 UInt_t index=r1.GetIndex(is);
500 AliTRDtrack *track=new AliTRDtrack(index, x, c, x1, ns*alpha+shift);
502 Int_t rc=FindProlongation(*track,fTrSec,ns,i2);
505 (track->GetNclusters() < (i1-i2)*fMinClustersInSeed)) delete track;
507 fSeeds->AddLast(track); fNseeds++;
508 cerr<<"found seed "<<fNseeds<<endl;
517 //___________________________________________________________________
518 void AliTRDtracker::FindTracks()
520 if (!fClusters) return;
522 AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect];
523 SetUpSectors(fTrSec);
527 Int_t num_of_time_bins = fTrSec[0].GetNtimeBins();
528 Int_t nseed=fSeeds->GetEntriesFast();
531 for (Int_t i=0; i<nseed; i++) {
532 cerr<<"FindTracks: seed "<<i+1<<" out of "<<nseed<<endl;
534 AliTRDtrack& t=*((AliTRDtrack*)fSeeds->UncheckedAt(i));
536 nSeedClusters = t.GetNclusters();
537 Double_t alpha=t.GetAlpha();
539 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
540 if (alpha < 0. ) alpha += 2.*TMath::Pi();
541 Int_t ns=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
543 if (FindProlongation(t,fTrSec,ns)) {
544 cerr<<"No of clusters in the track = "<<t.GetNclusters()<<endl;
545 if (t.GetNclusters() >= Int_t(fMinClustersInTrack*num_of_time_bins)) {
546 Int_t label = GetTrackLabel(t);
550 AliTRDtrack *pt = new AliTRDtrack(t);
551 fTracks->AddLast(pt); fNtracks++;
553 cerr<<"found track "<<fNtracks<<endl;
556 delete fSeeds->RemoveAt(i);
560 //__________________________________________________________________
561 void AliTRDtracker::UseClusters(AliTRDtrack t) {
562 Int_t ncl=t.GetNclusters();
563 for (Int_t i=0; i<ncl; i++) {
564 Int_t index = t.GetClusterIndex(i);
565 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
570 //__________________________________________________________________
571 Int_t AliTRDtracker::GetTrackLabel(AliTRDtrack t) {
573 Int_t label=123456789, index, i, j;
574 Int_t ncl=t.GetNclusters();
575 const Int_t range = AliTRDgeometry::kNplan * fGeom->GetTimeMax();
578 // Int_t s[range][2];
579 Int_t **s = new Int_t* [range];
580 for (i=0; i<range; i++) {
583 for (i=0; i<range; i++) {
589 for (i=0; i<ncl; i++) {
590 index=t.GetClusterIndex(i);
591 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
592 t0=c->GetTrackIndex(0);
593 t1=c->GetTrackIndex(1);
594 t2=c->GetTrackIndex(2);
597 for (i=0; i<ncl; i++) {
598 index=t.GetClusterIndex(i);
599 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
600 for (Int_t k=0; k<3; k++) {
601 label=c->GetTrackIndex(k);
602 label_added=kFALSE; j=0;
604 while ( (!label_added) && ( j < range ) ) {
605 if (s[j][0]==label || s[j][1]==0) {
619 for (i=0; i<range; i++) {
621 max=s[i][1]; label=s[i][0];
625 if(max > ncl*fLabelFraction) return label;
629 //___________________________________________________________________
631 Int_t AliTRDtracker::WriteTracks(const Char_t *filename) {
633 TDirectory *savedir=gDirectory;
635 TFile *out=TFile::Open(filename,"RECREATE");
637 TTree tracktree("TreeT","Tree with TRD tracks");
639 AliTRDtrack *iotrack=0;
640 tracktree.Branch("tracks","AliTRDtrack",&iotrack,32000,0);
642 Int_t ntracks=fTracks->GetEntriesFast();
644 for (Int_t i=0; i<ntracks; i++) {
645 AliTRDtrack *pt=(AliTRDtrack*)fTracks->UncheckedAt(i);
648 cerr<<"WriteTracks: put track "<<i<<" in the tree"<<endl;
656 cerr<<"WriteTracks: done"<<endl;
660 //_____________________________________________________________________________
661 void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename,
665 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
666 // from the file. The names of the cluster tree and branches
667 // should match the ones used in AliTRDclusterizer::WriteClusters()
670 TDirectory *savedir=gDirectory;
672 TFile *file = TFile::Open(filename);
673 if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;}
675 TTree *tree = (TTree*)file->Get("ClusterTree");
676 Int_t nentr = (Int_t) tree->GetEntries();
677 printf("found %d entries in %s.\n",nentr,tree->GetName());
680 TObjArray *ioArray = new TObjArray(400);
683 branch = tree->GetBranch("RecPoints");
685 for (Int_t i=0; i<nentr; i++) {
686 branch->SetAddress(&ioArray);
688 Int_t npoints = ioArray->GetEntriesFast();
689 printf("Read %d rec. points from entry %d \n", npoints, i);
691 for(Int_t j=0; j<npoints; j++) {
692 AliTRDrecPoint *p=(AliTRDrecPoint*)ioArray->UncheckedAt(j);
694 ioArray->RemoveAt(j);
699 branch = tree->GetBranch("Clusters");
701 for (Int_t i=0; i<nentr; i++) {
702 branch->SetAddress(&ioArray);
704 Int_t npoints = ioArray->GetEntriesFast();
705 printf("Read %d clusters from entry %d \n", npoints, i);
707 for(Int_t j=0; j<npoints; j++) {
708 AliTRDcluster *c=(AliTRDcluster*)ioArray->UncheckedAt(j);
710 ioArray->RemoveAt(j);