Latest version
[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$
0a29d0f1 18Revision 1.14 2001/11/14 10:50:46 cblume
19Changes in digits IO. Add merging of summable digits
20
abaf1f1d 21Revision 1.13 2001/05/30 12:17:47 hristov
22Loop variables declared once
23
3ab6f951 24Revision 1.12 2001/05/28 17:07:58 hristov
25Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
26
71d9fa7b 27Revision 1.8 2000/12/20 13:00:44 cblume
28Modifications for the HP-compiler
29
d1b06c24 30Revision 1.7 2000/12/08 16:07:02 cblume
31Update of the tracking by Sergei
32
bbf92647 33Revision 1.6 2000/11/30 17:38:08 cblume
34Changes to get in line with new STEER and EVGEN
35
1819f4bb 36Revision 1.5 2000/11/14 14:40:27 cblume
37Correction for the Sun compiler (kTRUE and kFALSE)
38
57527628 39Revision 1.4 2000/11/10 14:57:52 cblume
40Changes in the geometry constants for the DEC compiler
41
dd56b762 42Revision 1.3 2000/10/15 23:40:01 cblume
43Remove AliTRDconst
44
0e9c2ad5 45Revision 1.2 2000/10/06 16:49:46 cblume
46Made Getters const
47
46d29e70 48Revision 1.1.2.2 2000/10/04 16:34:58 cblume
49Replace include files by forward declarations
50
51Revision 1.1.2.1 2000/09/22 14:47:52 cblume
52Add the tracking code
53
54*/
bbf92647 55
0a29d0f1 56////////////////////////////////////////////////////////////////////////////////
57// //
58// The TRD tracker //
59// //
60////////////////////////////////////////////////////////////////////////////////
61
1819f4bb 62#include <iostream.h>
46d29e70 63
64#include <TFile.h>
65#include <TROOT.h>
66#include <TBranch.h>
67#include <TTree.h>
68
69#include "AliRun.h"
46d29e70 70#include "AliTRD.h"
46d29e70 71#include "AliTRDgeometry.h"
72#include "AliTRDrecPoint.h"
73#include "AliTRDcluster.h"
74#include "AliTRDtrack.h"
75#include "AliTRDtrackingSector.h"
76#include "AliTRDtimeBin.h"
77
78#include "AliTRDtracker.h"
79
80ClassImp(AliTRDtracker)
81
0a29d0f1 82 const Float_t AliTRDtracker::fgkSeedDepth = 0.5;
83 const Float_t AliTRDtracker::fgkSeedStep = 0.05;
84 const Float_t AliTRDtracker::fgkSeedGap = 0.25;
a819a5f7 85
0a29d0f1 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.;
bbf92647 91
0a29d0f1 92 const Double_t AliTRDtracker::fgkSeedErrorSY = 0.2;
93 const Double_t AliTRDtracker::fgkSeedErrorSY3 = 2.5;
94 const Double_t AliTRDtracker::fgkSeedErrorSZ = 0.1;
bbf92647 95
0a29d0f1 96 const Float_t AliTRDtracker::fgkMinClustersInSeed = 0.7;
bbf92647 97
0a29d0f1 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.;
a819a5f7 102
0a29d0f1 103 const Double_t AliTRDtracker::fgkMaxChi2 = 24.;
bbf92647 104
46d29e70 105//____________________________________________________________________
106AliTRDtracker::AliTRDtracker()
107{
108 //
109 // Default constructor
110 //
111
46d29e70 112 fEvent = 0;
46d29e70 113 fGeom = NULL;
bbf92647 114
46d29e70 115 fNclusters = 0;
116 fClusters = NULL;
117 fNseeds = 0;
118 fSeeds = NULL;
119 fNtracks = 0;
120 fTracks = NULL;
121
a819a5f7 122 fSY2corr = 0.025;
46d29e70 123}
124
125//____________________________________________________________________
126AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title)
127 :TNamed(name, title)
128{
0a29d0f1 129 //
130 // TRD tracker contructor
131 //
132
46d29e70 133 fEvent = 0;
46d29e70 134 fGeom = NULL;
bbf92647 135
46d29e70 136 fNclusters = 0;
137 fClusters = new TObjArray(2000);
138 fNseeds = 0;
139 fSeeds = new TObjArray(20000);
140 fNtracks = 0;
141 fTracks = new TObjArray(10000);
142
a819a5f7 143 fSY2corr = 0.025;
46d29e70 144}
145
46d29e70 146//___________________________________________________________________
147AliTRDtracker::~AliTRDtracker()
148{
0a29d0f1 149 //
150 // Destructor
151 //
152
46d29e70 153 delete fClusters;
154 delete fTracks;
155 delete fSeeds;
156 delete fGeom;
0a29d0f1 157
46d29e70 158}
159
bbf92647 160//___________________________________________________________________
a819a5f7 161void AliTRDtracker::Clusters2Tracks(TH1F *hs, TH1F *hd)
bbf92647 162{
0a29d0f1 163 //
164 // Do the trackfinding
165 //
166
3ab6f951 167 Int_t i;
168
bbf92647 169 Int_t inner, outer;
170 Int_t fTotalNofTimeBins = fGeom->GetTimeMax() * AliTRDgeometry::Nplan();
0a29d0f1 171 Int_t nSteps = (Int_t) (fgkSeedDepth / fgkSeedStep);
172 Int_t gap = (Int_t) (fTotalNofTimeBins * fgkSeedGap);
173 Int_t step = (Int_t) (fTotalNofTimeBins * fgkSeedStep);
a819a5f7 174
175
176 // nSteps = 1;
177
178 if (!fClusters) return;
179
180 AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect];
181 SetUpSectors(fTrSec);
182
183 // make a first turn looking for seed ends in the same (n,n)
184 // and in the adjacent sectors (n,n+1)
bbf92647 185
3ab6f951 186 for(i=0; i<nSteps; i++) {
bbf92647 187 printf("step %d out of %d \n", i+1, nSteps);
a819a5f7 188 outer=fTotalNofTimeBins-1-i*step; inner=outer-gap;
189 MakeSeeds(inner,outer, fTrSec, 1, hs, hd);
190 FindTracks(fTrSec, hs, hd);
bbf92647 191 }
a819a5f7 192
193 // make a second turn looking for seed ends in next-to-adjacent
194 // sectors (n,n+2)
195
3ab6f951 196 for(i=0; i<nSteps; i++) {
a819a5f7 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);
201 }
202
bbf92647 203}
46d29e70 204
205//_____________________________________________________________________
0a29d0f1 206Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) const
46d29e70 207{
0a29d0f1 208 //
46d29e70 209 // Parametrised "expected" error of the cluster reconstruction in Y
0a29d0f1 210 //
46d29e70 211
a819a5f7 212 Double_t s = 0.08 * 0.08;
46d29e70 213 return s;
0a29d0f1 214
46d29e70 215}
216
217//_____________________________________________________________________
0a29d0f1 218Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl) const
46d29e70 219{
0a29d0f1 220 //
46d29e70 221 // Parametrised "expected" error of the cluster reconstruction in Z
0a29d0f1 222 //
46d29e70 223
a819a5f7 224 Double_t s = 6 * 6 /12.;
46d29e70 225 return s;
0a29d0f1 226
46d29e70 227}
228
229//_____________________________________________________________________
bdbd0f7a 230Double_t f1trd(Double_t x1,Double_t y1,
a819a5f7 231 Double_t x2,Double_t y2,
232 Double_t x3,Double_t y3)
46d29e70 233{
0a29d0f1 234 //
46d29e70 235 // Initial approximation of the track curvature
236 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
0a29d0f1 237 //
46d29e70 238
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));
244
245 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
246
247 return -xr*yr/sqrt(xr*xr+yr*yr);
0a29d0f1 248
46d29e70 249}
250
251//_____________________________________________________________________
bdbd0f7a 252Double_t f2trd(Double_t x1,Double_t y1,
a819a5f7 253 Double_t x2,Double_t y2,
254 Double_t x3,Double_t y3)
46d29e70 255{
0a29d0f1 256 //
46d29e70 257 // Initial approximation of the track curvature times center of curvature
258 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
0a29d0f1 259 //
46d29e70 260
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));
266
267 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
268
269 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
0a29d0f1 270
46d29e70 271}
272
273//_____________________________________________________________________
bdbd0f7a 274Double_t f3trd(Double_t x1,Double_t y1,
a819a5f7 275 Double_t x2,Double_t y2,
276 Double_t z1,Double_t z2)
46d29e70 277{
0a29d0f1 278 //
46d29e70 279 // Initial approximation of the tangent of the track dip angle
280 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
0a29d0f1 281 //
46d29e70 282
283 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
0a29d0f1 284
46d29e70 285}
286
287
288//___________________________________________________________________
289
bbf92647 290Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec,
0a29d0f1 291 Int_t s, Int_t rf, Int_t matchedIndex,
a819a5f7 292 TH1F *hs, TH1F *hd)
46d29e70 293{
0a29d0f1 294 //
46d29e70 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.
0a29d0f1 299 //
a819a5f7 300
301 // TH1F *hsame = hs;
302 // TH1F *hdiff = hd;
303
304 // Bool_t good_match;
305
0a29d0f1 306 const Int_t kTimeBinsToSkip=Int_t(fgkSkipDepth*sec->GetNtimeBins());
307 Int_t tryAgain=kTimeBinsToSkip;
46d29e70 308
0e9c2ad5 309 Double_t alpha=AliTRDgeometry::GetAlpha();
46d29e70 310
311 Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5);
312
a819a5f7 313 Double_t x0, rho;
314
46d29e70 315 for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) {
a819a5f7 316
317 Double_t x=sec->GetX(nr);
318 Double_t ymax=x*TMath::Tan(0.5*alpha);
319
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
323 }
324 if(TMath::Abs(x - t.GetX()) > 3) {
325 rho = 0.0559; x0 = 55.6; // radiator
326 }
327 if (!t.PropagateTo(x,x0,rho,0.139)) {
46d29e70 328 cerr<<"Can't propagate to x = "<<x<<endl;
329 return 0;
330 }
331
0a29d0f1 332 AliTRDtimeBin& timeBin=sec[s][nr];
bbf92647 333 Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
334 Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
a819a5f7 335 Double_t road=25.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
46d29e70 336
0a29d0f1 337 if (road>fgkWideRoad) {
46d29e70 338 if (t.GetNclusters() > 4) {
339 cerr<<t.GetNclusters()<<" FindProlongation: Road is too wide !\n";
340 }
341 return 0;
342 }
343
a819a5f7 344 AliTRDcluster *cl=0;
345 UInt_t index=0;
346 // Int_t ncl = 0;
347
0a29d0f1 348 Double_t maxChi2=fgkMaxChi2;
a819a5f7 349
0a29d0f1 350 if (timeBin) {
a819a5f7 351
0a29d0f1 352 for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
353 AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
a819a5f7 354
355 // good_match = kFALSE;
0a29d0f1 356 // if((c->GetTrackIndex(0) == matchedIndex) ||
357 // (c->GetTrackIndex(1) == matchedIndex) ||
358 // (c->GetTrackIndex(2) == matchedIndex)) good_match = kTRUE;
a819a5f7 359 // if(hsame) hsame->Fill(TMath::Abs(c->GetY()-y)/road);
360 // if(hdiff) hdiff->Fill(road);
361
46d29e70 362 if (c->GetY() > y+road) break;
363 if (c->IsUsed() > 0) continue;
364
a819a5f7 365 // if(good_match) hsame->Fill(TMath::Abs(c->GetZ()-z));
366 // else hdiff->Fill(TMath::Abs(c->GetZ()-z));
367
368 // if(!good_match) continue;
369
370 if((c->GetZ()-z)*(c->GetZ()-z) > 3 * 12 * sz2) continue;
46d29e70 371
372 Double_t chi2=t.GetPredictedChi2(c);
373
0a29d0f1 374 // if((c->GetTrackIndex(0) == matchedIndex) ||
375 // (c->GetTrackIndex(1) == matchedIndex) ||
376 // (c->GetTrackIndex(2) == matchedIndex))
a819a5f7 377 // hdiff->Fill(chi2);
378
379 // ncl++;
380
0a29d0f1 381 if (chi2 > maxChi2) continue;
382 maxChi2=chi2;
46d29e70 383 cl=c;
0a29d0f1 384 index=timeBin.GetIndex(i);
46d29e70 385 }
a819a5f7 386
387 if(!cl) {
388
0a29d0f1 389 for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
390 AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
a819a5f7 391
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;
395
396 Double_t chi2=t.GetPredictedChi2(c);
397
398 // ncl++;
399
0a29d0f1 400 if (chi2 > maxChi2) continue;
401 maxChi2=chi2;
a819a5f7 402 cl=c;
0a29d0f1 403 index=timeBin.GetIndex(i);
a819a5f7 404 }
405 }
406
46d29e70 407 }
a819a5f7 408
46d29e70 409 if (cl) {
410
0a29d0f1 411 t.Update(cl,maxChi2,index);
46d29e70 412
0a29d0f1 413 tryAgain=kTimeBinsToSkip;
46d29e70 414 } else {
0a29d0f1 415 if (tryAgain==0) break;
46d29e70 416 if (y > ymax) {
a819a5f7 417 // cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl;
46d29e70 418 s = (s+1) % ns;
419 if (!t.Rotate(alpha)) {
420 cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
421 return 0;
422 }
423 } else if (y <-ymax) {
a819a5f7 424 // cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl;
46d29e70 425 s = (s-1+ns) % ns;
426 if (!t.Rotate(-alpha)) {
427 cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
428 return 0;
429 }
430 }
0a29d0f1 431 if(!sec->TECframe(nr,y,z)) tryAgain--;
46d29e70 432 }
433 }
434
435 return 1;
436}
437
438
439
440//_____________________________________________________________________________
bbf92647 441void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile)
46d29e70 442{
0a29d0f1 443 //
a819a5f7 444 // Opens a ROOT-file with TRD-clusters and reads the cluster-tree in
0a29d0f1 445 //
46d29e70 446
bbf92647 447 ReadClusters(fClusters, clusterfile);
448
449 // get geometry from the file with hits
450
451 TFile *fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(hitfile);
46d29e70 452 if (!fInputFile) {
453 printf("AliTRDtracker::Open -- ");
bbf92647 454 printf("Open the ALIROOT-file %s.\n",hitfile);
a819a5f7 455 fInputFile = new TFile(hitfile,"UPDATE");
46d29e70 456 }
457 else {
458 printf("AliTRDtracker::Open -- ");
bbf92647 459 printf("%s is already open.\n",hitfile);
46d29e70 460 }
461
462 // Get AliRun object from file or create it if not on file
46d29e70 463
bbf92647 464 gAlice = (AliRun*) fInputFile->Get("gAlice");
465 if (gAlice) {
466 printf("AliTRDtracker::GetEvent -- ");
467 printf("AliRun object found on file.\n");
468 }
469 else {
470 printf("AliTRDtracker::GetEvent -- ");
471 printf("Could not find AliRun object.\n");
472 }
46d29e70 473
a819a5f7 474 /*
46d29e70 475 // Import the Trees for the event nEvent in the file
476 Int_t nparticles = gAlice->GetEvent(fEvent);
477 cerr<<"nparticles = "<<nparticles<<endl;
478
479 if (nparticles <= 0) {
480 printf("AliTRDtracker::GetEvent -- ");
481 printf("No entries in the trees for event %d.\n",fEvent);
482 }
a819a5f7 483 */
46d29e70 484
0a29d0f1 485 AliTRD *trd = (AliTRD*) gAlice->GetDetector("TRD");
486 fGeom = trd->GetGeometry();
46d29e70 487
46d29e70 488}
489
490
491//_____________________________________________________________________________
492void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec)
493{
0a29d0f1 494 //
46d29e70 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
0a29d0f1 498 //
46d29e70 499
bbf92647 500 for (Int_t i=0; i<AliTRDgeometry::Nsect(); i++) sec[i].SetUp();
46d29e70 501
502 // Sort clusters into AliTRDtimeBin's within AliTRDtrackSector's
503
a819a5f7 504 cerr<<"SetUpSectors: sorting clusters"<<endl;
46d29e70 505
506 Int_t ncl=fClusters->GetEntriesFast();
507 UInt_t index;
508 while (ncl--) {
a819a5f7 509 printf("\r %d left ",ncl);
46d29e70 510 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
0a29d0f1 511 Int_t detector=c->GetDetector(), localTimeBin=c->GetLocalTimeBin();
46d29e70 512 Int_t sector=fGeom->GetSector(detector);
513
0a29d0f1 514 Int_t trackingSector = AliTRDgeometry::kNsect - sector - 1;
46d29e70 515
0a29d0f1 516 Int_t tb=sec[sector].GetTimeBin(detector,localTimeBin);
46d29e70 517 index=ncl;
0a29d0f1 518 sec[trackingSector][tb].InsertCluster(c,index);
a819a5f7 519
46d29e70 520 }
521 printf("\r\n");
522}
523
524
525//_____________________________________________________________________________
a819a5f7 526void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer,
527 AliTRDtrackingSector* fTrSec, Int_t turn,
528 TH1F *hs, TH1F *hd)
46d29e70 529{
0a29d0f1 530 //
46d29e70 531 // Creates track seeds using clusters in timeBins=i1,i2
0a29d0f1 532 //
46d29e70 533
bbf92647 534 Int_t i2 = inner, i1 = outer;
a819a5f7 535 Int_t ti[3], to[3];
536 Int_t nprim = 85600/2;
46d29e70 537
46d29e70 538
a819a5f7 539 TH1F *hsame = hs;
540 TH1F *hdiff = hd;
541 Bool_t match = false;
0a29d0f1 542 Int_t matchedIndex;
46d29e70 543
544 // find seeds
545
546 Double_t x[5], c[15];
0a29d0f1 547 Int_t maxSec=AliTRDgeometry::kNsect;
46d29e70 548
0e9c2ad5 549 Double_t alpha=AliTRDgeometry::GetAlpha();
550 Double_t shift=AliTRDgeometry::GetAlpha()/2.;
46d29e70 551 Double_t cs=cos(alpha), sn=sin(alpha);
a819a5f7 552 Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
46d29e70 553
554 Double_t x1 =fTrSec[0].GetX(i1);
555 Double_t xx2=fTrSec[0].GetX(i2);
556
a819a5f7 557
558 printf("\n");
559
560 if((turn != 1)&&(turn != 2)) {
561 printf("*** Error in MakeSeeds: unexpected turn = %d \n", turn);
562 return;
563 }
564
565
0a29d0f1 566 for (Int_t ns=0; ns<maxSec; ns++) {
46d29e70 567
a819a5f7 568 printf("\n MakeSeeds: sector %d \n", ns);
569
0a29d0f1 570 Int_t nl2=fTrSec[(ns-2+maxSec)%maxSec][i2];
571 Int_t nl=fTrSec[(ns-1+maxSec)%maxSec][i2];
46d29e70 572 Int_t nm=fTrSec[ns][i2];
0a29d0f1 573 Int_t nu=fTrSec[(ns+1)%maxSec][i2];
574 Int_t nu2=fTrSec[(ns+2)%maxSec][i2];
46d29e70 575
576 AliTRDtimeBin& r1=fTrSec[ns][i1];
577
578 for (Int_t is=0; is < r1; is++) {
579 Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
a819a5f7 580 for(Int_t ii=0; ii<3; ii++) to[ii] = r1[is]->GetTrackIndex(ii);
46d29e70 581
a819a5f7 582 for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
583
46d29e70 584 const AliTRDcluster *cl;
a819a5f7 585 Double_t x2, y2, z2;
586 Double_t x3=0., y3=0.;
46d29e70 587
a819a5f7 588 if (js<nl2) {
589 if(turn != 2) continue;
0a29d0f1 590 AliTRDtimeBin& r2=fTrSec[(ns-2+maxSec)%maxSec][i2];
a819a5f7 591 cl=r2[js];
592 y2=cl->GetY(); z2=cl->GetZ();
593 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
594
595 x2= xx2*cs2+y2*sn2;
596 y2=-xx2*sn2+y2*cs2;
597 }
598 else if (js<nl2+nl) {
599 if(turn != 1) continue;
0a29d0f1 600 AliTRDtimeBin& r2=fTrSec[(ns-1+maxSec)%maxSec][i2];
a819a5f7 601 cl=r2[js-nl2];
46d29e70 602 y2=cl->GetY(); z2=cl->GetZ();
a819a5f7 603 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
46d29e70 604
a819a5f7 605 x2= xx2*cs+y2*sn;
606 y2=-xx2*sn+y2*cs;
46d29e70 607
a819a5f7 608 }
609 else if (js<nl2+nl+nm) {
610 if(turn != 1) continue;
611 AliTRDtimeBin& r2=fTrSec[ns][i2];
612 cl=r2[js-nl2-nl];
613 x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
614 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
615 }
616 else if (js<nl2+nl+nm+nu) {
617 if(turn != 1) continue;
0a29d0f1 618 AliTRDtimeBin& r2=fTrSec[(ns+1)%maxSec][i2];
a819a5f7 619 cl=r2[js-nl2-nl-nm];
620 y2=cl->GetY(); z2=cl->GetZ();
621 for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
46d29e70 622
a819a5f7 623 x2=xx2*cs-y2*sn;
624 y2=xx2*sn+y2*cs;
46d29e70 625
a819a5f7 626 }
627 else {
628 if(turn != 2) continue;
0a29d0f1 629 AliTRDtimeBin& r2=fTrSec[(ns+2)%maxSec][i2];
a819a5f7 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);
633
634 x2=xx2*cs2-y2*sn2;
635 y2=xx2*sn2+y2*cs2;
636 }
637
46d29e70 638
a819a5f7 639 match = false;
0a29d0f1 640 matchedIndex = -1;
a819a5f7 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;
0a29d0f1 650 matchedIndex = ti[ii];
a819a5f7 651 match = true;
652 }
653 }
654 }
46d29e70 655
0a29d0f1 656 if(TMath::Abs(z1-z2) > fgkMaxSeedDeltaZ12) continue;
a819a5f7 657
658 Double_t zz=z1 - z1/x1*(x1-x2);
659
0a29d0f1 660 if (TMath::Abs(zz-z2)>fgkMaxSeedDeltaZ) continue;
46d29e70 661
662 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
663 if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
664
665 x[0]=y1;
666 x[1]=z1;
667 x[2]=f1trd(x1,y1,x2,y2,x3,y3);
668
0a29d0f1 669 if (TMath::Abs(x[2]) > fgkMaxSeedC) continue;
46d29e70 670
671 x[3]=f2trd(x1,y1,x2,y2,x3,y3);
672
673 if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue;
674
675 x[4]=f3trd(x1,y1,x2,y2,z1,z2);
676
0a29d0f1 677 if (TMath::Abs(x[4]) > fgkMaxSeedTan) continue;
46d29e70 678
679 Double_t a=asin(x[3]);
680 Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3]));
a819a5f7 681
0a29d0f1 682 if (TMath::Abs(zv)>fgkMaxSeedVertexZ) continue;
46d29e70 683
bbf92647 684 Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
685 Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2();
0a29d0f1 686 Double_t sy3=fgkSeedErrorSY3, sy=fgkSeedErrorSY, sz=fgkSeedErrorSZ;
46d29e70 687
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;
698
699 c[0]=sy1;
700 c[1]=0.; c[2]=sz1;
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;
707
708 UInt_t index=r1.GetIndex(is);
a819a5f7 709
710 AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
711
0a29d0f1 712 Int_t rc=FindProlongation(*track,fTrSec,ns,i2,matchedIndex,hsame,hdiff);
a819a5f7 713
714 // if (match) hsame->Fill((Float_t) track->GetNclusters());
715 // else hdiff->Fill((Float_t) track->GetNclusters());
716 // delete track;
717 // continue;
46d29e70 718
bbf92647 719 if ((rc < 0) ||
0a29d0f1 720 (track->GetNclusters() < (i1-i2)*fgkMinClustersInSeed)) delete track;
46d29e70 721 else {
722 fSeeds->AddLast(track); fNseeds++;
a819a5f7 723 printf("\r found seed %d ", fNseeds);
46d29e70 724 }
725 }
726 }
727 }
728
46d29e70 729 fSeeds->Sort();
730}
731
a819a5f7 732//_____________________________________________________________________________
733void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename)
46d29e70 734{
a819a5f7 735 //
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()
739 //
46d29e70 740
a819a5f7 741 TDirectory *savedir=gDirectory;
742
743 TFile *file = TFile::Open(filename);
744 if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;}
745
abaf1f1d 746 Char_t treeName[12];
747 sprintf(treeName,"TreeR%d_TRD",fEvent);
0a29d0f1 748 TTree *clusterTree = (TTree*) file->Get(treeName);
46d29e70 749
0a29d0f1 750 TObjArray *clusterArray = new TObjArray(400);
a819a5f7 751
0a29d0f1 752 clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
a819a5f7 753
0a29d0f1 754 Int_t nEntries = (Int_t) clusterTree->GetEntries();
755 printf("found %d entries in %s.\n",nEntries,clusterTree->GetName());
a819a5f7 756
757 // Loop through all entries in the tree
758 Int_t nbytes;
759 AliTRDcluster *c = 0;
760
761 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
762
763 // Import the tree
0a29d0f1 764 nbytes += clusterTree->GetEvent(iEntry);
a819a5f7 765
766 // Get the number of points in the detector
0a29d0f1 767 Int_t nCluster = clusterArray->GetEntriesFast();
a819a5f7 768 printf("Read %d clusters from entry %d \n", nCluster, iEntry);
769
770 // Loop through all TRD digits
771 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
0a29d0f1 772 c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
a819a5f7 773 AliTRDcluster *co = new AliTRDcluster(*c);
774 co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
775 array->AddLast(co);
0a29d0f1 776 delete clusterArray->RemoveAt(iCluster);
a819a5f7 777 }
778 }
779
780 file->Close();
0a29d0f1 781 delete clusterArray;
a819a5f7 782 savedir->cd();
783
784}
785
786//___________________________________________________________________
787void AliTRDtracker::FindTracks(AliTRDtrackingSector* fTrSec, TH1F *hs, TH1F *hd)
788{
789 //
790 // Finds tracks in TRD
791 //
792
793 TH1F *hsame = hs;
794 TH1F *hdiff = hd;
46d29e70 795
0a29d0f1 796 Int_t numOfTimeBins = fTrSec[0].GetNtimeBins();
46d29e70 797 Int_t nseed=fSeeds->GetEntriesFast();
798
799 Int_t nSeedClusters;
800 for (Int_t i=0; i<nseed; i++) {
801 cerr<<"FindTracks: seed "<<i+1<<" out of "<<nseed<<endl;
802
bbf92647 803 AliTRDtrack& t=*((AliTRDtrack*)fSeeds->UncheckedAt(i));
46d29e70 804
805 nSeedClusters = t.GetNclusters();
bbf92647 806 Double_t alpha=t.GetAlpha();
46d29e70 807
808 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
809 if (alpha < 0. ) alpha += 2.*TMath::Pi();
dd56b762 810 Int_t ns=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
46d29e70 811
a819a5f7 812 Int_t label = GetTrackLabel(t);
813
814 if (FindProlongation(t,fTrSec,ns,0,label,hsame,hdiff)) {
46d29e70 815 cerr<<"No of clusters in the track = "<<t.GetNclusters()<<endl;
0a29d0f1 816 if (t.GetNclusters() >= Int_t(fgkMinClustersInTrack*numOfTimeBins)) {
46d29e70 817 Int_t label = GetTrackLabel(t);
818 t.SetLabel(label);
a819a5f7 819 t.CookdEdx();
46d29e70 820 UseClusters(t);
bbf92647 821
822 AliTRDtrack *pt = new AliTRDtrack(t);
823 fTracks->AddLast(pt); fNtracks++;
824
46d29e70 825 cerr<<"found track "<<fNtracks<<endl;
826 }
827 }
bbf92647 828 delete fSeeds->RemoveAt(i);
a819a5f7 829 fNseeds--;
46d29e70 830 }
831}
832
833//__________________________________________________________________
0a29d0f1 834void AliTRDtracker::UseClusters(AliTRDtrack t)
835{
836 //
837 // Mark used cluster
838 //
839
46d29e70 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);
844 c->Use();
845 }
0a29d0f1 846
46d29e70 847}
848
849//__________________________________________________________________
0a29d0f1 850Int_t AliTRDtracker::GetTrackLabel(AliTRDtrack t)
851{
852 //
853 // Get MC label
854 //
46d29e70 855
856 Int_t label=123456789, index, i, j;
857 Int_t ncl=t.GetNclusters();
0a29d0f1 858 const Int_t kRange = AliTRDgeometry::kNplan * fGeom->GetTimeMax();
859 Bool_t labelAdded;
46d29e70 860
0a29d0f1 861 // Int_t s[kRange][2];
862 Int_t **s = new Int_t* [kRange];
863 for (i=0; i<kRange; i++) {
d1b06c24 864 s[i] = new Int_t[2];
865 }
0a29d0f1 866 for (i=0; i<kRange; i++) {
46d29e70 867 s[i][0]=-1;
868 s[i][1]=0;
869 }
870
871 Int_t t0,t1,t2;
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);
878 }
879
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);
0a29d0f1 885 labelAdded=kFALSE; j=0;
46d29e70 886 if (label >= 0) {
0a29d0f1 887 while ( (!labelAdded) && ( j < kRange ) ) {
46d29e70 888 if (s[j][0]==label || s[j][1]==0) {
889 s[j][0]=label;
890 s[j][1]=s[j][1]+1;
0a29d0f1 891 labelAdded=kTRUE;
46d29e70 892 }
893 j++;
894 }
895 }
896 }
897 }
898
46d29e70 899 Int_t max=0;
900 label = -123456789;
901
0a29d0f1 902 for (i=0; i<kRange; i++) {
46d29e70 903 if (s[i][1]>max) {
904 max=s[i][1]; label=s[i][0];
905 }
906 }
d1b06c24 907 delete []s;
0a29d0f1 908 if(max > ncl*fgkLabelFraction) return label;
46d29e70 909 else return -1;
910}
911
912//___________________________________________________________________
0a29d0f1 913Int_t AliTRDtracker::WriteTracks(const Char_t *filename)
914{
915 //
916 // Write the tracks to the output file
917 //
46d29e70 918
919 TDirectory *savedir=gDirectory;
920
a819a5f7 921 //TFile *out=TFile::Open(filename,"RECREATE");
922 TFile *out = (TFile*) gROOT->GetListOfFiles()->FindObject(filename);
923 if (!out) {
924 printf("AliTRDtracker::Open -- ");
925 printf("Open the ALIROOT-file %s.\n",filename);
926 out = new TFile(filename,"RECREATE");
927 }
928 else {
929 printf("AliTRDtracker::Open -- ");
930 printf("%s is already open.\n",filename);
931 }
46d29e70 932
abaf1f1d 933 Char_t treeName[12];
934 sprintf(treeName,"TreeT%d_TRD",fEvent);
935 TTree tracktree(treeName,"Tree with TRD tracks");
46d29e70 936
937 AliTRDtrack *iotrack=0;
938 tracktree.Branch("tracks","AliTRDtrack",&iotrack,32000,0);
939
940 Int_t ntracks=fTracks->GetEntriesFast();
941
942 for (Int_t i=0; i<ntracks; i++) {
943 AliTRDtrack *pt=(AliTRDtrack*)fTracks->UncheckedAt(i);
944 iotrack=pt;
945 tracktree.Fill();
946 cerr<<"WriteTracks: put track "<<i<<" in the tree"<<endl;
947 }
948
949 tracktree.Write();
950 out->Close();
951
952 savedir->cd();
953
954 cerr<<"WriteTracks: done"<<endl;
a819a5f7 955 return 1;
0a29d0f1 956
46d29e70 957}
958