Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of...
[u/mrichter/AliRoot.git] / TRD / AliTRDtracker.cxx
CommitLineData
46d29e70 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16/*
17$Log$
71d9fa7b 18Revision 1.8 2000/12/20 13:00:44 cblume
19Modifications for the HP-compiler
20
d1b06c24 21Revision 1.7 2000/12/08 16:07:02 cblume
22Update of the tracking by Sergei
23
bbf92647 24Revision 1.6 2000/11/30 17:38:08 cblume
25Changes to get in line with new STEER and EVGEN
26
1819f4bb 27Revision 1.5 2000/11/14 14:40:27 cblume
28Correction for the Sun compiler (kTRUE and kFALSE)
29
57527628 30Revision 1.4 2000/11/10 14:57:52 cblume
31Changes in the geometry constants for the DEC compiler
32
dd56b762 33Revision 1.3 2000/10/15 23:40:01 cblume
34Remove AliTRDconst
35
0e9c2ad5 36Revision 1.2 2000/10/06 16:49:46 cblume
37Made Getters const
38
46d29e70 39Revision 1.1.2.2 2000/10/04 16:34:58 cblume
40Replace include files by forward declarations
41
42Revision 1.1.2.1 2000/09/22 14:47:52 cblume
43Add the tracking code
44
45*/
bbf92647 46
1819f4bb 47#include <iostream.h>
46d29e70 48
49#include <TFile.h>
50#include <TROOT.h>
51#include <TBranch.h>
52#include <TTree.h>
53
54#include "AliRun.h"
46d29e70 55#include "AliTRD.h"
46d29e70 56#include "AliTRDgeometry.h"
57#include "AliTRDrecPoint.h"
58#include "AliTRDcluster.h"
59#include "AliTRDtrack.h"
60#include "AliTRDtrackingSector.h"
61#include "AliTRDtimeBin.h"
62
63#include "AliTRDtracker.h"
64
65ClassImp(AliTRDtracker)
66
a819a5f7 67 const Float_t AliTRDtracker::fSeedDepth = 0.5;
68 const Float_t AliTRDtracker::fSeedStep = 0.05;
69 const Float_t AliTRDtracker::fSeedGap = 0.25;
70
71 const Float_t AliTRDtracker::fMaxSeedDeltaZ12 = 40.;
72 const Float_t AliTRDtracker::fMaxSeedDeltaZ = 25.;
73 const Float_t AliTRDtracker::fMaxSeedC = 0.0052;
74 const Float_t AliTRDtracker::fMaxSeedTan = 1.2;
75 const Float_t AliTRDtracker::fMaxSeedVertexZ = 150.;
bbf92647 76
a819a5f7 77 const Double_t AliTRDtracker::fSeedErrorSY = 0.2;
78 const Double_t AliTRDtracker::fSeedErrorSY3 = 2.5;
79 const Double_t AliTRDtracker::fSeedErrorSZ = 0.1;
bbf92647 80
a819a5f7 81 const Float_t AliTRDtracker::fMinClustersInSeed = 0.7;
bbf92647 82
a819a5f7 83 const Float_t AliTRDtracker::fMinClustersInTrack = 0.5;
84 const Float_t AliTRDtracker::fSkipDepth = 0.05;
85 const Float_t AliTRDtracker::fLabelFraction = 0.5;
86 const Float_t AliTRDtracker::fWideRoad = 20.;
87
88 const Double_t AliTRDtracker::fMaxChi2 = 24.;
bbf92647 89
46d29e70 90//____________________________________________________________________
91AliTRDtracker::AliTRDtracker()
92{
93 //
94 // Default constructor
95 //
96
46d29e70 97 fEvent = 0;
46d29e70 98 fGeom = NULL;
bbf92647 99
46d29e70 100 fNclusters = 0;
101 fClusters = NULL;
102 fNseeds = 0;
103 fSeeds = NULL;
104 fNtracks = 0;
105 fTracks = NULL;
106
a819a5f7 107 fSY2corr = 0.025;
46d29e70 108}
109
110//____________________________________________________________________
111AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title)
112 :TNamed(name, title)
113{
46d29e70 114 fEvent = 0;
46d29e70 115 fGeom = NULL;
bbf92647 116
46d29e70 117 fNclusters = 0;
118 fClusters = new TObjArray(2000);
119 fNseeds = 0;
120 fSeeds = new TObjArray(20000);
121 fNtracks = 0;
122 fTracks = new TObjArray(10000);
123
a819a5f7 124 fSY2corr = 0.025;
46d29e70 125}
126
46d29e70 127//___________________________________________________________________
128AliTRDtracker::~AliTRDtracker()
129{
46d29e70 130 delete fClusters;
131 delete fTracks;
132 delete fSeeds;
133 delete fGeom;
134}
135
bbf92647 136//___________________________________________________________________
a819a5f7 137void AliTRDtracker::Clusters2Tracks(TH1F *hs, TH1F *hd)
bbf92647 138{
139 Int_t inner, outer;
140 Int_t fTotalNofTimeBins = fGeom->GetTimeMax() * AliTRDgeometry::Nplan();
a819a5f7 141 Int_t nSteps = (Int_t) (fSeedDepth / fSeedStep);
142 Int_t gap = (Int_t) (fTotalNofTimeBins * fSeedGap);
143 Int_t step = (Int_t) (fTotalNofTimeBins * fSeedStep);
144
145
146 // nSteps = 1;
147
148 if (!fClusters) return;
149
150 AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect];
151 SetUpSectors(fTrSec);
152
153 // make a first turn looking for seed ends in the same (n,n)
154 // and in the adjacent sectors (n,n+1)
bbf92647 155
156 for(Int_t i=0; i<nSteps; i++) {
157 printf("step %d out of %d \n", i+1, nSteps);
a819a5f7 158 outer=fTotalNofTimeBins-1-i*step; inner=outer-gap;
159 MakeSeeds(inner,outer, fTrSec, 1, hs, hd);
160 FindTracks(fTrSec, hs, hd);
bbf92647 161 }
a819a5f7 162
163 // make a second turn looking for seed ends in next-to-adjacent
164 // sectors (n,n+2)
165
166 for(Int_t i=0; i<nSteps; i++) {
167 printf("step %d out of %d \n", i+1, nSteps);
168 outer=fTotalNofTimeBins-1-i*step; inner=outer-gap;
169 MakeSeeds(inner, outer, fTrSec, 2, hs, hd);
170 FindTracks(fTrSec,hs,hd);
171 }
172
bbf92647 173}
46d29e70 174
175//_____________________________________________________________________
bbf92647 176Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt)
46d29e70 177{
178 // Parametrised "expected" error of the cluster reconstruction in Y
179
a819a5f7 180 Double_t s = 0.08 * 0.08;
46d29e70 181 return s;
182}
183
184//_____________________________________________________________________
bbf92647 185Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl)
46d29e70 186{
187 // Parametrised "expected" error of the cluster reconstruction in Z
188
a819a5f7 189 Double_t s = 6 * 6 /12.;
46d29e70 190 return s;
191}
192
193//_____________________________________________________________________
bdbd0f7a 194Double_t f1trd(Double_t x1,Double_t y1,
a819a5f7 195 Double_t x2,Double_t y2,
196 Double_t x3,Double_t y3)
46d29e70 197{
198 // Initial approximation of the track curvature
199 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
200
201 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
202 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
203 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
204 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
205 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
206
207 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
208
209 return -xr*yr/sqrt(xr*xr+yr*yr);
210}
211
212//_____________________________________________________________________
bdbd0f7a 213Double_t f2trd(Double_t x1,Double_t y1,
a819a5f7 214 Double_t x2,Double_t y2,
215 Double_t x3,Double_t y3)
46d29e70 216{
217 // Initial approximation of the track curvature times center of curvature
218 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
219
220 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
221 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
222 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
223 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
224 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
225
226 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
227
228 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
229}
230
231//_____________________________________________________________________
bdbd0f7a 232Double_t f3trd(Double_t x1,Double_t y1,
a819a5f7 233 Double_t x2,Double_t y2,
234 Double_t z1,Double_t z2)
46d29e70 235{
236 // Initial approximation of the tangent of the track dip angle
237 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
238
239 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
240}
241
242
243//___________________________________________________________________
244
bbf92647 245Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec,
a819a5f7 246 Int_t s, Int_t rf, Int_t matched_index,
247 TH1F *hs, TH1F *hd)
46d29e70 248{
249 // Starting from current position on track=t this function tries
250 // to extrapolate the track up to timeBin=rf and to confirm prolongation
251 // if a close cluster is found. *sec is a pointer to allocated
252 // array of sectors, in which the initial sector has index=s.
253
a819a5f7 254
255 // TH1F *hsame = hs;
256 // TH1F *hdiff = hd;
257
258 // Bool_t good_match;
259
bbf92647 260 const Int_t TIME_BINS_TO_SKIP=Int_t(fSkipDepth*sec->GetNtimeBins());
46d29e70 261 Int_t try_again=TIME_BINS_TO_SKIP;
262
0e9c2ad5 263 Double_t alpha=AliTRDgeometry::GetAlpha();
46d29e70 264
265 Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5);
266
a819a5f7 267 Double_t x0, rho;
268
46d29e70 269 for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) {
a819a5f7 270
271 Double_t x=sec->GetX(nr);
272 Double_t ymax=x*TMath::Tan(0.5*alpha);
273
274 rho = 0.00295; x0 = 11.09; // TEC
275 if(sec->TECframe(nr,t.GetY(),t.GetZ())) {
276 rho = 1.7; x0 = 33.0; // G10 frame
277 }
278 if(TMath::Abs(x - t.GetX()) > 3) {
279 rho = 0.0559; x0 = 55.6; // radiator
280 }
281 if (!t.PropagateTo(x,x0,rho,0.139)) {
46d29e70 282 cerr<<"Can't propagate to x = "<<x<<endl;
283 return 0;
284 }
285
46d29e70 286 AliTRDtimeBin& time_bin=sec[s][nr];
bbf92647 287 Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
288 Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
a819a5f7 289 Double_t road=25.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
46d29e70 290
bbf92647 291 if (road>fWideRoad) {
46d29e70 292 if (t.GetNclusters() > 4) {
293 cerr<<t.GetNclusters()<<" FindProlongation: Road is too wide !\n";
294 }
295 return 0;
296 }
297
a819a5f7 298 AliTRDcluster *cl=0;
299 UInt_t index=0;
300 // Int_t ncl = 0;
301
302 Double_t max_chi2=fMaxChi2;
303
46d29e70 304 if (time_bin) {
a819a5f7 305
46d29e70 306 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
307 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
a819a5f7 308
309 // good_match = kFALSE;
310 // if((c->GetTrackIndex(0) == matched_index) ||
311 // (c->GetTrackIndex(1) == matched_index) ||
312 // (c->GetTrackIndex(2) == matched_index)) good_match = kTRUE;
313 // if(hsame) hsame->Fill(TMath::Abs(c->GetY()-y)/road);
314 // if(hdiff) hdiff->Fill(road);
315
46d29e70 316 if (c->GetY() > y+road) break;
317 if (c->IsUsed() > 0) continue;
318
a819a5f7 319 // if(good_match) hsame->Fill(TMath::Abs(c->GetZ()-z));
320 // else hdiff->Fill(TMath::Abs(c->GetZ()-z));
321
322 // if(!good_match) continue;
323
324 if((c->GetZ()-z)*(c->GetZ()-z) > 3 * 12 * sz2) continue;
46d29e70 325
326 Double_t chi2=t.GetPredictedChi2(c);
327
a819a5f7 328 // if((c->GetTrackIndex(0) == matched_index) ||
329 // (c->GetTrackIndex(1) == matched_index) ||
330 // (c->GetTrackIndex(2) == matched_index))
331 // hdiff->Fill(chi2);
332
333 // ncl++;
334
46d29e70 335 if (chi2 > max_chi2) continue;
336 max_chi2=chi2;
337 cl=c;
338 index=time_bin.GetIndex(i);
339 }
a819a5f7 340
341 if(!cl) {
342
343 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
344 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
345
346 if (c->GetY() > y+road) break;
347 if (c->IsUsed() > 0) continue;
348 if((c->GetZ()-z)*(c->GetZ()-z) > 3.25 * 12 * sz2) continue;
349
350 Double_t chi2=t.GetPredictedChi2(c);
351
352 // ncl++;
353
354 if (chi2 > max_chi2) continue;
355 max_chi2=chi2;
356 cl=c;
357 index=time_bin.GetIndex(i);
358 }
359 }
360
46d29e70 361 }
a819a5f7 362
46d29e70 363 if (cl) {
364
46d29e70 365 t.Update(cl,max_chi2,index);
366
367 try_again=TIME_BINS_TO_SKIP;
368 } else {
369 if (try_again==0) break;
370 if (y > ymax) {
a819a5f7 371 // cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl;
46d29e70 372 s = (s+1) % ns;
373 if (!t.Rotate(alpha)) {
374 cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
375 return 0;
376 }
377 } else if (y <-ymax) {
a819a5f7 378 // cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl;
46d29e70 379 s = (s-1+ns) % ns;
380 if (!t.Rotate(-alpha)) {
381 cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
382 return 0;
383 }
384 }
a819a5f7 385 if(!sec->TECframe(nr,y,z)) try_again--;
46d29e70 386 }
387 }
388
389 return 1;
390}
391
392
393
394//_____________________________________________________________________________
bbf92647 395void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile)
46d29e70 396{
a819a5f7 397 // Opens a ROOT-file with TRD-clusters and reads the cluster-tree in
46d29e70 398
bbf92647 399 ReadClusters(fClusters, clusterfile);
400
401 // get geometry from the file with hits
402
403 TFile *fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(hitfile);
46d29e70 404 if (!fInputFile) {
405 printf("AliTRDtracker::Open -- ");
bbf92647 406 printf("Open the ALIROOT-file %s.\n",hitfile);
a819a5f7 407 fInputFile = new TFile(hitfile,"UPDATE");
46d29e70 408 }
409 else {
410 printf("AliTRDtracker::Open -- ");
bbf92647 411 printf("%s is already open.\n",hitfile);
46d29e70 412 }
413
414 // Get AliRun object from file or create it if not on file
46d29e70 415
bbf92647 416 gAlice = (AliRun*) fInputFile->Get("gAlice");
417 if (gAlice) {
418 printf("AliTRDtracker::GetEvent -- ");
419 printf("AliRun object found on file.\n");
420 }
421 else {
422 printf("AliTRDtracker::GetEvent -- ");
423 printf("Could not find AliRun object.\n");
424 }
46d29e70 425
a819a5f7 426 /*
46d29e70 427 // Import the Trees for the event nEvent in the file
428 Int_t nparticles = gAlice->GetEvent(fEvent);
429 cerr<<"nparticles = "<<nparticles<<endl;
430
431 if (nparticles <= 0) {
432 printf("AliTRDtracker::GetEvent -- ");
433 printf("No entries in the trees for event %d.\n",fEvent);
434 }
a819a5f7 435 */
46d29e70 436
437 AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
438 fGeom = TRD->GetGeometry();
439
46d29e70 440}
441
442
443//_____________________________________________________________________________
444void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec)
445{
446 // Fills clusters into TRD tracking_sectors
447 // Note that the numbering scheme for the TRD tracking_sectors
448 // differs from that of TRD sectors
449
bbf92647 450 for (Int_t i=0; i<AliTRDgeometry::Nsect(); i++) sec[i].SetUp();
46d29e70 451
452 // Sort clusters into AliTRDtimeBin's within AliTRDtrackSector's
453
a819a5f7 454 cerr<<"SetUpSectors: sorting clusters"<<endl;
46d29e70 455
456 Int_t ncl=fClusters->GetEntriesFast();
457 UInt_t index;
458 while (ncl--) {
a819a5f7 459 printf("\r %d left ",ncl);
46d29e70 460 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
461 Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin();
462 Int_t sector=fGeom->GetSector(detector);
463
a819a5f7 464 Int_t tracking_sector = AliTRDgeometry::kNsect - sector - 1;
46d29e70 465
466 Int_t tb=sec[sector].GetTimeBin(detector,local_time_bin);
467 index=ncl;
468 sec[tracking_sector][tb].InsertCluster(c,index);
a819a5f7 469
46d29e70 470 }
471 printf("\r\n");
472}
473
474
475//_____________________________________________________________________________
a819a5f7 476void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer,
477 AliTRDtrackingSector* fTrSec, Int_t turn,
478 TH1F *hs, TH1F *hd)
46d29e70 479{
480 // Creates track seeds using clusters in timeBins=i1,i2
481
bbf92647 482 Int_t i2 = inner, i1 = outer;
a819a5f7 483 Int_t ti[3], to[3];
484 Int_t nprim = 85600/2;
46d29e70 485
46d29e70 486
a819a5f7 487 TH1F *hsame = hs;
488 TH1F *hdiff = hd;
489 Bool_t match = false;
490 Int_t matched_index;
46d29e70 491
492 // find seeds
493
494 Double_t x[5], c[15];
dd56b762 495 Int_t max_sec=AliTRDgeometry::kNsect;
46d29e70 496
0e9c2ad5 497 Double_t alpha=AliTRDgeometry::GetAlpha();
498 Double_t shift=AliTRDgeometry::GetAlpha()/2.;
46d29e70 499 Double_t cs=cos(alpha), sn=sin(alpha);
a819a5f7 500 Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
46d29e70 501
502 Double_t x1 =fTrSec[0].GetX(i1);
503 Double_t xx2=fTrSec[0].GetX(i2);
504
a819a5f7 505
506 printf("\n");
507
508 if((turn != 1)&&(turn != 2)) {
509 printf("*** Error in MakeSeeds: unexpected turn = %d \n", turn);
510 return;
511 }
512
513
46d29e70 514 for (Int_t ns=0; ns<max_sec; ns++) {
515
a819a5f7 516 printf("\n MakeSeeds: sector %d \n", ns);
517
518 Int_t nl2=fTrSec[(ns-2+max_sec)%max_sec][i2];
46d29e70 519 Int_t nl=fTrSec[(ns-1+max_sec)%max_sec][i2];
520 Int_t nm=fTrSec[ns][i2];
521 Int_t nu=fTrSec[(ns+1)%max_sec][i2];
a819a5f7 522 Int_t nu2=fTrSec[(ns+2)%max_sec][i2];
46d29e70 523
524 AliTRDtimeBin& r1=fTrSec[ns][i1];
525
526 for (Int_t is=0; is < r1; is++) {
527 Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
a819a5f7 528 for(Int_t ii=0; ii<3; ii++) to[ii] = r1[is]->GetTrackIndex(ii);
46d29e70 529
a819a5f7 530 for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
531
46d29e70 532 const AliTRDcluster *cl;
a819a5f7 533 Double_t x2, y2, z2;
534 Double_t x3=0., y3=0.;
46d29e70 535
a819a5f7 536 if (js<nl2) {
537 if(turn != 2) continue;
538 AliTRDtimeBin& r2=fTrSec[(ns-2+max_sec)%max_sec][i2];
539 cl=r2[js];
540 y2=cl->GetY(); z2=cl->GetZ();
541 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
542
543 x2= xx2*cs2+y2*sn2;
544 y2=-xx2*sn2+y2*cs2;
545 }
546 else if (js<nl2+nl) {
547 if(turn != 1) continue;
46d29e70 548 AliTRDtimeBin& r2=fTrSec[(ns-1+max_sec)%max_sec][i2];
a819a5f7 549 cl=r2[js-nl2];
46d29e70 550 y2=cl->GetY(); z2=cl->GetZ();
a819a5f7 551 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
46d29e70 552
a819a5f7 553 x2= xx2*cs+y2*sn;
554 y2=-xx2*sn+y2*cs;
46d29e70 555
a819a5f7 556 }
557 else if (js<nl2+nl+nm) {
558 if(turn != 1) continue;
559 AliTRDtimeBin& r2=fTrSec[ns][i2];
560 cl=r2[js-nl2-nl];
561 x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
562 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
563 }
564 else if (js<nl2+nl+nm+nu) {
565 if(turn != 1) continue;
566 AliTRDtimeBin& r2=fTrSec[(ns+1)%max_sec][i2];
567 cl=r2[js-nl2-nl-nm];
568 y2=cl->GetY(); z2=cl->GetZ();
569 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
46d29e70 570
a819a5f7 571 x2=xx2*cs-y2*sn;
572 y2=xx2*sn+y2*cs;
46d29e70 573
a819a5f7 574 }
575 else {
576 if(turn != 2) continue;
577 AliTRDtimeBin& r2=fTrSec[(ns+2)%max_sec][i2];
578 cl=r2[js-nl2-nl-nm-nu];
579 y2=cl->GetY(); z2=cl->GetZ();
580 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
581
582 x2=xx2*cs2-y2*sn2;
583 y2=xx2*sn2+y2*cs2;
584 }
585
46d29e70 586
a819a5f7 587 match = false;
588 matched_index = -1;
589 for (Int_t ii=0; ii<3; ii++) {
590 // cerr<<"ti["<<ii<<"] = "<<ti[ii]<<"; to["<<ii<<"] = "<<to[ii]<<endl;
591 if(ti[ii] < 0) continue;
592 if(ti[ii] >= nprim) continue;
593 for (Int_t kk=0; kk<3; kk++) {
594 if(to[kk] < 0) continue;
595 if(to[kk] >= nprim) continue;
596 if(ti[ii] == to[kk]) {
597 //cerr<<"ti["<<ii<<"] = "<<ti[ii]<<" = "<<to[kk]<<" = to["<<kk<<"]"<<endl;
598 matched_index = ti[ii];
599 match = true;
600 }
601 }
602 }
46d29e70 603
a819a5f7 604 if(TMath::Abs(z1-z2) > fMaxSeedDeltaZ12) continue;
605
606 Double_t zz=z1 - z1/x1*(x1-x2);
607
bbf92647 608 if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue;
46d29e70 609
610 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
611 if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
612
613 x[0]=y1;
614 x[1]=z1;
615 x[2]=f1trd(x1,y1,x2,y2,x3,y3);
616
a819a5f7 617 if (TMath::Abs(x[2]) > fMaxSeedC) continue;
46d29e70 618
619 x[3]=f2trd(x1,y1,x2,y2,x3,y3);
620
621 if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue;
622
623 x[4]=f3trd(x1,y1,x2,y2,z1,z2);
624
bbf92647 625 if (TMath::Abs(x[4]) > fMaxSeedTan) continue;
46d29e70 626
627 Double_t a=asin(x[3]);
628 Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3]));
a819a5f7 629
bbf92647 630 if (TMath::Abs(zv)>fMaxSeedVertexZ) continue;
46d29e70 631
bbf92647 632 Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
633 Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2();
634 Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ;
46d29e70 635
636 Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
637 Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
638 Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
639 Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy;
640 Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy;
641 Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy;
642 Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy;
643 Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz;
644 Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy;
645 Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz;
646
647 c[0]=sy1;
648 c[1]=0.; c[2]=sz1;
649 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24;
650 c[6]=f30*sy1; c[7]=0.; c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24;
651 c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34;
652 c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22;
653 c[13]=f40*sy1*f30+f42*sy2*f32;
654 c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43;
655
656 UInt_t index=r1.GetIndex(is);
a819a5f7 657
658 AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
659
660 Int_t rc=FindProlongation(*track,fTrSec,ns,i2,matched_index,hsame,hdiff);
661
662 // if (match) hsame->Fill((Float_t) track->GetNclusters());
663 // else hdiff->Fill((Float_t) track->GetNclusters());
664 // delete track;
665 // continue;
46d29e70 666
bbf92647 667 if ((rc < 0) ||
668 (track->GetNclusters() < (i1-i2)*fMinClustersInSeed)) delete track;
46d29e70 669 else {
670 fSeeds->AddLast(track); fNseeds++;
a819a5f7 671 printf("\r found seed %d ", fNseeds);
46d29e70 672 }
673 }
674 }
675 }
676
46d29e70 677 fSeeds->Sort();
678}
679
a819a5f7 680//_____________________________________________________________________________
681void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename)
46d29e70 682{
a819a5f7 683 //
684 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
685 // from the file. The names of the cluster tree and branches
686 // should match the ones used in AliTRDclusterizer::WriteClusters()
687 //
46d29e70 688
a819a5f7 689 TDirectory *savedir=gDirectory;
690
691 TFile *file = TFile::Open(filename);
692 if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;}
693
694 TTree *ClusterTree = (TTree*)file->Get("ClusterTree");
46d29e70 695
a819a5f7 696 TObjArray *ClusterArray = new TObjArray(400);
697
698 ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray);
699
700 Int_t nEntries = (Int_t) ClusterTree->GetEntries();
701 printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
702
703 // Loop through all entries in the tree
704 Int_t nbytes;
705 AliTRDcluster *c = 0;
706
707 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
708
709 // Import the tree
710 nbytes += ClusterTree->GetEvent(iEntry);
711
712 // Get the number of points in the detector
713 Int_t nCluster = ClusterArray->GetEntriesFast();
714 printf("Read %d clusters from entry %d \n", nCluster, iEntry);
715
716 // Loop through all TRD digits
717 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
718 c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
719 AliTRDcluster *co = new AliTRDcluster(*c);
720 co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
721 array->AddLast(co);
722 delete ClusterArray->RemoveAt(iCluster);
723 }
724 }
725
726 file->Close();
727 delete ClusterArray;
728 savedir->cd();
729
730}
731
732//___________________________________________________________________
733void AliTRDtracker::FindTracks(AliTRDtrackingSector* fTrSec, TH1F *hs, TH1F *hd)
734{
735 //
736 // Finds tracks in TRD
737 //
738
739 TH1F *hsame = hs;
740 TH1F *hdiff = hd;
46d29e70 741
742 Int_t num_of_time_bins = fTrSec[0].GetNtimeBins();
743 Int_t nseed=fSeeds->GetEntriesFast();
744
745 Int_t nSeedClusters;
746 for (Int_t i=0; i<nseed; i++) {
747 cerr<<"FindTracks: seed "<<i+1<<" out of "<<nseed<<endl;
748
bbf92647 749 AliTRDtrack& t=*((AliTRDtrack*)fSeeds->UncheckedAt(i));
46d29e70 750
751 nSeedClusters = t.GetNclusters();
bbf92647 752 Double_t alpha=t.GetAlpha();
46d29e70 753
754 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
755 if (alpha < 0. ) alpha += 2.*TMath::Pi();
dd56b762 756 Int_t ns=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
46d29e70 757
a819a5f7 758 Int_t label = GetTrackLabel(t);
759
760 if (FindProlongation(t,fTrSec,ns,0,label,hsame,hdiff)) {
46d29e70 761 cerr<<"No of clusters in the track = "<<t.GetNclusters()<<endl;
bbf92647 762 if (t.GetNclusters() >= Int_t(fMinClustersInTrack*num_of_time_bins)) {
46d29e70 763 Int_t label = GetTrackLabel(t);
764 t.SetLabel(label);
a819a5f7 765 t.CookdEdx();
46d29e70 766 UseClusters(t);
bbf92647 767
768 AliTRDtrack *pt = new AliTRDtrack(t);
769 fTracks->AddLast(pt); fNtracks++;
770
46d29e70 771 cerr<<"found track "<<fNtracks<<endl;
772 }
773 }
bbf92647 774 delete fSeeds->RemoveAt(i);
a819a5f7 775 fNseeds--;
46d29e70 776 }
777}
778
779//__________________________________________________________________
bbf92647 780void AliTRDtracker::UseClusters(AliTRDtrack t) {
46d29e70 781 Int_t ncl=t.GetNclusters();
782 for (Int_t i=0; i<ncl; i++) {
783 Int_t index = t.GetClusterIndex(i);
784 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
785 c->Use();
786 }
787}
788
789//__________________________________________________________________
bbf92647 790Int_t AliTRDtracker::GetTrackLabel(AliTRDtrack t) {
46d29e70 791
792 Int_t label=123456789, index, i, j;
793 Int_t ncl=t.GetNclusters();
bbf92647 794 const Int_t range = AliTRDgeometry::kNplan * fGeom->GetTimeMax();
46d29e70 795 Bool_t label_added;
796
d1b06c24 797 // Int_t s[range][2];
798 Int_t **s = new Int_t* [range];
799 for (i=0; i<range; i++) {
800 s[i] = new Int_t[2];
801 }
46d29e70 802 for (i=0; i<range; i++) {
803 s[i][0]=-1;
804 s[i][1]=0;
805 }
806
807 Int_t t0,t1,t2;
808 for (i=0; i<ncl; i++) {
809 index=t.GetClusterIndex(i);
810 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
811 t0=c->GetTrackIndex(0);
812 t1=c->GetTrackIndex(1);
813 t2=c->GetTrackIndex(2);
814 }
815
816 for (i=0; i<ncl; i++) {
817 index=t.GetClusterIndex(i);
818 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
819 for (Int_t k=0; k<3; k++) {
820 label=c->GetTrackIndex(k);
57527628 821 label_added=kFALSE; j=0;
46d29e70 822 if (label >= 0) {
823 while ( (!label_added) && ( j < range ) ) {
824 if (s[j][0]==label || s[j][1]==0) {
825 s[j][0]=label;
826 s[j][1]=s[j][1]+1;
57527628 827 label_added=kTRUE;
46d29e70 828 }
829 j++;
830 }
831 }
832 }
833 }
834
46d29e70 835 Int_t max=0;
836 label = -123456789;
837
838 for (i=0; i<range; i++) {
839 if (s[i][1]>max) {
840 max=s[i][1]; label=s[i][0];
841 }
842 }
d1b06c24 843 delete []s;
bbf92647 844 if(max > ncl*fLabelFraction) return label;
46d29e70 845 else return -1;
846}
847
848//___________________________________________________________________
849
bbf92647 850Int_t AliTRDtracker::WriteTracks(const Char_t *filename) {
46d29e70 851
852 TDirectory *savedir=gDirectory;
853
a819a5f7 854 //TFile *out=TFile::Open(filename,"RECREATE");
855 TFile *out = (TFile*) gROOT->GetListOfFiles()->FindObject(filename);
856 if (!out) {
857 printf("AliTRDtracker::Open -- ");
858 printf("Open the ALIROOT-file %s.\n",filename);
859 out = new TFile(filename,"RECREATE");
860 }
861 else {
862 printf("AliTRDtracker::Open -- ");
863 printf("%s is already open.\n",filename);
864 }
46d29e70 865
866 TTree tracktree("TreeT","Tree with TRD tracks");
867
868 AliTRDtrack *iotrack=0;
869 tracktree.Branch("tracks","AliTRDtrack",&iotrack,32000,0);
870
871 Int_t ntracks=fTracks->GetEntriesFast();
872
873 for (Int_t i=0; i<ntracks; i++) {
874 AliTRDtrack *pt=(AliTRDtrack*)fTracks->UncheckedAt(i);
875 iotrack=pt;
876 tracktree.Fill();
877 cerr<<"WriteTracks: put track "<<i<<" in the tree"<<endl;
878 }
879
880 tracktree.Write();
881 out->Close();
882
883 savedir->cd();
884
885 cerr<<"WriteTracks: done"<<endl;
a819a5f7 886 return 1;
46d29e70 887}
888