Replacing fabs by TMath::Abs
[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
c630aafd 16/*
17$Log$
029cd327 18Revision 1.30 2003/09/18 09:06:07 cblume
19Geometry update, Removal of compiler warnings
20
73ae7b59 21Revision 1.29 2003/07/22 15:56:14 hristov
22Implementing ESD functionality in the NewIO (Yu.Belikov)
23
c630aafd 24Revision 1.27.2.2 2003/07/14 09:19:33 hristov
25TOF included in the combined PID (Yu.Belikov)
26
27Revision 1.27.2.1 2003/07/11 10:53:01 hristov
28Inward refit for TPC and TRD in the ESD schema (T.Kuhr)
29
30Revision 1.27 2003/05/27 17:46:13 hristov
31TRD PID included in the ESD schema (T.Kuhr)
32
33Revision 1.26 2003/04/10 10:36:54 hristov
34Code for unified TPC/TRD tracking (S.Radomski)
35
36Revision 1.25 2003/03/19 17:14:11 hristov
37Load/UnloadClusters added to the base class and the derived classes changed correspondingly. Possibility to give 2 input files for ITS and TPC tracks in PropagateBack. TRD tracker uses fEventN from the base class (T.Kuhr)
38
39Revision 1.24 2003/02/19 09:02:28 hristov
40Track time measurement (S.Radomski)
41
42Revision 1.23 2003/02/10 14:06:10 cblume
43Add tracking without tilted pads as option
44
45Revision 1.22 2003/01/30 15:19:58 cblume
46New set of parameters
47
48Revision 1.21 2003/01/27 16:34:49 cblume
49Update of tracking by Sergei and Chuncheng
50
51Revision 1.20 2002/11/07 15:52:09 cblume
52Update of tracking code for tilted pads
53
54Revision 1.19 2002/10/22 15:53:08 alibrary
55Introducing Riostream.h
56
57Revision 1.18 2002/10/14 14:57:44 hristov
58Merging the VirtualMC branch to the main development branch (HEAD)
59
60Revision 1.14.6.2 2002/07/24 10:09:31 alibrary
61Updating VirtualMC
62
63Revision 1.17 2002/06/13 12:09:58 hristov
64Minor corrections
65
66Revision 1.16 2002/06/12 09:54:36 cblume
67Update of tracking code provided by Sergei
68
69Revision 1.14 2001/11/14 10:50:46 cblume
70Changes in digits IO. Add merging of summable digits
71
72Revision 1.13 2001/05/30 12:17:47 hristov
73Loop variables declared once
74
75Revision 1.12 2001/05/28 17:07:58 hristov
76Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
77
78Revision 1.8 2000/12/20 13:00:44 cblume
79Modifications for the HP-compiler
80
81Revision 1.7 2000/12/08 16:07:02 cblume
82Update of the tracking by Sergei
83
84Revision 1.6 2000/11/30 17:38:08 cblume
85Changes to get in line with new STEER and EVGEN
86
87Revision 1.5 2000/11/14 14:40:27 cblume
88Correction for the Sun compiler (kTRUE and kFALSE)
89
90Revision 1.4 2000/11/10 14:57:52 cblume
91Changes in the geometry constants for the DEC compiler
92
93Revision 1.3 2000/10/15 23:40:01 cblume
94Remove AliTRDconst
95
96Revision 1.2 2000/10/06 16:49:46 cblume
97Made Getters const
98
99Revision 1.1.2.2 2000/10/04 16:34:58 cblume
100Replace include files by forward declarations
101
102Revision 1.1.2.1 2000/09/22 14:47:52 cblume
103Add the tracking code
104
105*/
bbf92647 106
029cd327 107///////////////////////////////////////////////////////////////////////////////
108// //
109// The standard TRD tracker //
110// //
111///////////////////////////////////////////////////////////////////////////////
112
a2cb5b3d 113#include <Riostream.h>
46d29e70 114
115#include <TFile.h>
46d29e70 116#include <TBranch.h>
5443e65e 117#include <TTree.h>
c630aafd 118#include <TObjArray.h>
46d29e70 119
46d29e70 120#include "AliTRDgeometry.h"
5443e65e 121#include "AliTRDparameter.h"
122#include "AliTRDgeometryDetail.h"
46d29e70 123#include "AliTRDcluster.h"
124#include "AliTRDtrack.h"
c630aafd 125#include "AliTRDPartID.h"
5443e65e 126#include "../TPC/AliTPCtrack.h"
46d29e70 127
128#include "AliTRDtracker.h"
129
130ClassImp(AliTRDtracker)
131
029cd327 132 const Float_t AliTRDtracker::fgkSeedDepth = 0.5;
133 const Float_t AliTRDtracker::fgkSeedStep = 0.10;
134 const Float_t AliTRDtracker::fgkSeedGap = 0.25;
5443e65e 135
029cd327 136 const Float_t AliTRDtracker::fgkMaxSeedDeltaZ12 = 40.;
137 const Float_t AliTRDtracker::fgkMaxSeedDeltaZ = 25.;
138 const Float_t AliTRDtracker::fgkMaxSeedC = 0.0052;
139 const Float_t AliTRDtracker::fgkMaxSeedTan = 1.2;
140 const Float_t AliTRDtracker::fgkMaxSeedVertexZ = 150.;
5443e65e 141
029cd327 142 const Double_t AliTRDtracker::fgkSeedErrorSY = 0.2;
143 const Double_t AliTRDtracker::fgkSeedErrorSY3 = 2.5;
144 const Double_t AliTRDtracker::fgkSeedErrorSZ = 0.1;
a819a5f7 145
029cd327 146 const Float_t AliTRDtracker::fgkMinClustersInSeed = 0.7;
bbf92647 147
029cd327 148 const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5;
149 const Float_t AliTRDtracker::fgkMinFractionOfFoundClusters = 0.8;
bbf92647 150
029cd327 151 const Float_t AliTRDtracker::fgkSkipDepth = 0.3;
152 const Float_t AliTRDtracker::fgkLabelFraction = 0.8;
153 const Float_t AliTRDtracker::fgkWideRoad = 20.;
bbf92647 154
029cd327 155 const Double_t AliTRDtracker::fgkMaxChi2 = 12.;
a819a5f7 156
029cd327 157const Int_t AliTRDtracker::fgkFirstPlane = 5;
158const Int_t AliTRDtracker::fgkLastPlane = 17;
9c9d2487 159
bbf92647 160
46d29e70 161//____________________________________________________________________
c630aafd 162AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
46d29e70 163{
5443e65e 164 //
165 // Main constructor
166 //
46d29e70 167
c630aafd 168 //Float_t fTzero = 0;
b8dc2353 169
5443e65e 170 fAddTRDseeds = kFALSE;
5443e65e 171 fGeom = NULL;
b8dc2353 172 fNoTilt = kFALSE;
5443e65e 173
174 TDirectory *savedir=gDirectory;
175 TFile *in=(TFile*)geomfile;
176 if (!in->IsOpen()) {
177 printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n");
178 printf(" DETAIL TRD geometry and DEFAULT TRD parameter will be used\n");
179 }
180 else {
181 in->cd();
c630aafd 182// in->ls();
5443e65e 183 fGeom = (AliTRDgeometry*) in->Get("TRDgeometry");
184 fPar = (AliTRDparameter*) in->Get("TRDparameter");
c630aafd 185// fGeom->Dump();
5443e65e 186 }
46d29e70 187
5443e65e 188 if(fGeom) {
189 // fTzero = geo->GetT0();
b8dc2353 190 printf("Found geometry version %d on file \n", fGeom->IsVersion());
5443e65e 191 }
192 else {
c630aafd 193 printf("AliTRDtracker::AliTRDtracker(): can't find TRD geometry!\n");
194 printf("The DETAIL TRD geometry will be used\n");
195 fGeom = new AliTRDgeometryDetail();
196 }
197
198 if (!fPar) {
199 printf("AliTRDtracker::AliTRDtracker(): can't find TRD parameter!\n");
200 printf("The DEFAULT TRD parameter will be used\n");
5443e65e 201 fPar = new AliTRDparameter();
202 }
c630aafd 203 fPar->ReInit();
5443e65e 204
205 savedir->cd();
46d29e70 206
5443e65e 207
208 // fGeom->SetT0(fTzero);
0a29d0f1 209
46d29e70 210 fNclusters = 0;
211 fClusters = new TObjArray(2000);
212 fNseeds = 0;
5443e65e 213 fSeeds = new TObjArray(2000);
46d29e70 214 fNtracks = 0;
5443e65e 215 fTracks = new TObjArray(1000);
46d29e70 216
029cd327 217 for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
218 Int_t trS = CookSectorIndex(geomS);
219 fTrSec[trS] = new AliTRDtrackingSector(fGeom, geomS, fPar);
5443e65e 220 }
a819a5f7 221
029cd327 222 Float_t tiltAngle = TMath::Abs(fPar->GetTiltingAngle());
223 if(tiltAngle < 0.1) {
b8dc2353 224 fNoTilt = kTRUE;
225 }
226
227 fSY2corr = 0.2;
228 fSZ2corr = 120.;
229
029cd327 230 if(fNoTilt && (tiltAngle > 0.1)) fSY2corr = fSY2corr + tiltAngle * 0.05;
b8dc2353 231
a819a5f7 232
5443e65e 233 // calculate max gap on track
a819a5f7 234
5443e65e 235 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
236 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
bbf92647 237
5443e65e 238 Double_t dx = (Double_t) fPar->GetTimeBinSize();
239 Int_t tbAmp = fPar->GetTimeBefore();
240 Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
b3a5a838 241 if(kTRUE) maxAmp = 0; // intentional until we change the parameter class
5443e65e 242 Int_t tbDrift = fPar->GetTimeMax();
243 Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
a819a5f7 244
5443e65e 245 tbDrift = TMath::Min(tbDrift,maxDrift);
246 tbAmp = TMath::Min(tbAmp,maxAmp);
a819a5f7 247
5443e65e 248 fTimeBinsPerPlane = tbAmp + tbDrift;
029cd327 249 fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fgkSkipDepth);
a819a5f7 250
5443e65e 251 fVocal = kFALSE;
46d29e70 252
9c9d2487 253
254 // Barrel Tracks [SR, 03.04.2003]
255
256 fBarrelFile = 0;
257 fBarrelTree = 0;
258 fBarrelArray = 0;
259 fBarrelTrack = 0;
260
261 savedir->cd();
5443e65e 262}
46d29e70 263
5443e65e 264//___________________________________________________________________
265AliTRDtracker::~AliTRDtracker()
46d29e70 266{
029cd327 267 //
268 // Destructor of AliTRDtracker
269 //
270
5443e65e 271 delete fClusters;
272 delete fTracks;
273 delete fSeeds;
274 delete fGeom;
275 delete fPar;
0a29d0f1 276
029cd327 277 for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
278 delete fTrSec[geomS];
5443e65e 279 }
280}
46d29e70 281
282//_____________________________________________________________________
9c9d2487 283
284void AliTRDtracker::SetBarrelTree(const char *mode) {
285 //
286 //
287 //
288
289 if (!IsStoringBarrel()) return;
290
291 TDirectory *sav = gDirectory;
292 if (!fBarrelFile) fBarrelFile = new TFile("AliBarrelTracks.root", "UPDATE");
293
294 char buff[40];
295 sprintf(buff, "BarrelTRD_%d_%s", GetEventNumber(), mode);
296
297 fBarrelFile->cd();
298 fBarrelTree = new TTree(buff, "Barrel TPC tracks");
299
029cd327 300 Int_t nRefs = fgkLastPlane - fgkFirstPlane + 1;
9c9d2487 301
302 if (!fBarrelArray) fBarrelArray = new TClonesArray("AliBarrelTrack", nRefs);
303 for(Int_t i=0; i<nRefs; i++) new((*fBarrelArray)[i]) AliBarrelTrack();
304
305 fBarrelTree->Branch("tracks", &fBarrelArray);
306 sav->cd();
307}
308
309//_____________________________________________________________________
310
311void AliTRDtracker::StoreBarrelTrack(AliTRDtrack *ps, Int_t refPlane, Int_t isIn) {
312 //
313 //
314 //
315
316 if (!IsStoringBarrel()) return;
317
318 static Int_t nClusters;
319 static Int_t nWrong;
320 static Double_t chi2;
321 static Int_t index;
322 static Bool_t wasLast = kTRUE;
323
324 Int_t newClusters, newWrong;
325 Double_t newChi2;
326
327 if (wasLast) {
328
329 fBarrelArray->Clear();
330 nClusters = nWrong = 0;
331 chi2 = 0.0;
332 index = 0;
333 wasLast = kFALSE;
334 }
335
336 fBarrelTrack = (AliBarrelTrack*)(*fBarrelArray)[index++];
337 ps->GetBarrelTrack(fBarrelTrack);
338
339 newClusters = ps->GetNumberOfClusters() - nClusters;
340 newWrong = ps->GetNWrong() - nWrong;
341 newChi2 = ps->GetChi2() - chi2;
342
343 nClusters = ps->GetNumberOfClusters();
344 nWrong = ps->GetNWrong();
345 chi2 = ps->GetChi2();
346
029cd327 347 if (refPlane != fgkLastPlane) {
9c9d2487 348 fBarrelTrack->SetNClusters(newClusters, newChi2);
349 fBarrelTrack->SetNWrongClusters(newWrong);
350 } else {
351 wasLast = kTRUE;
352 }
353
354 fBarrelTrack->SetRefPlane(refPlane, isIn);
355}
356
357//_____________________________________________________________________
358
359Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track) {
360 //
361 // Rotates the track when necessary
362 //
363
364 Double_t alpha = AliTRDgeometry::GetAlpha();
365 Double_t y = track->GetY();
366 Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
367
c630aafd 368 //Int_t ns = AliTRDgeometry::kNsect;
9c9d2487 369 //Int_t s=Int_t(track->GetAlpha()/alpha)%ns;
370
371 if (y > ymax) {
372 //s = (s+1) % ns;
373 if (!track->Rotate(alpha)) return kFALSE;
374 } else if (y <-ymax) {
375 //s = (s-1+ns) % ns;
376 if (!track->Rotate(-alpha)) return kFALSE;
377 }
378
379 return kTRUE;
380}
381
382//_____________________________________________________________________
5443e65e 383inline Double_t f1trd(Double_t x1,Double_t y1,
a9814c09 384 Double_t x2,Double_t y2,
385 Double_t x3,Double_t y3)
46d29e70 386{
0a29d0f1 387 //
46d29e70 388 // Initial approximation of the track curvature
0a29d0f1 389 //
46d29e70 390 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
391 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
392 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
393 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
394 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
395
396 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
397
398 return -xr*yr/sqrt(xr*xr+yr*yr);
399}
400
401//_____________________________________________________________________
5443e65e 402inline Double_t f2trd(Double_t x1,Double_t y1,
a9814c09 403 Double_t x2,Double_t y2,
404 Double_t x3,Double_t y3)
46d29e70 405{
0a29d0f1 406 //
5443e65e 407 // Initial approximation of the track curvature times X coordinate
408 // of the center of curvature
0a29d0f1 409 //
46d29e70 410
411 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
412 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
413 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
414 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
415 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
416
417 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
418
419 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
420}
421
422//_____________________________________________________________________
5443e65e 423inline Double_t f3trd(Double_t x1,Double_t y1,
a9814c09 424 Double_t x2,Double_t y2,
425 Double_t z1,Double_t z2)
46d29e70 426{
0a29d0f1 427 //
46d29e70 428 // Initial approximation of the tangent of the track dip angle
0a29d0f1 429 //
46d29e70 430
431 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
432}
433
46d29e70 434//___________________________________________________________________
5443e65e 435Int_t AliTRDtracker::Clusters2Tracks(const TFile *inp, TFile *out)
46d29e70 436{
0a29d0f1 437 //
5443e65e 438 // Finds tracks within the TRD. File <inp> is expected to contain seeds
439 // at the outer part of the TRD. If <inp> is NULL, the seeds
440 // are found within the TRD if fAddTRDseeds is TRUE.
441 // The tracks are propagated to the innermost time bin
442 // of the TRD and stored in file <out>.
0a29d0f1 443 //
a819a5f7 444
5443e65e 445 LoadEvent();
446
447 TDirectory *savedir=gDirectory;
a819a5f7 448
5443e65e 449 char tname[100];
a819a5f7 450
5443e65e 451 if (!out->IsOpen()) {
452 cerr<<"AliTRDtracker::Clusters2Tracks(): output file is not open !\n";
453 return 1;
454 }
46d29e70 455
e24ea474 456 sprintf(tname,"seedTRDtoTPC_%d",GetEventNumber());
029cd327 457 TTree tpcTree(tname,"Tree with seeds from TRD at outer TPC pad row");
5443e65e 458 AliTPCtrack *iotrack=0;
029cd327 459 tpcTree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
5443e65e 460
e24ea474 461 sprintf(tname,"TreeT%d_TRD",GetEventNumber());
029cd327 462 TTree trdTree(tname,"TRD tracks at inner TRD time bin");
463 AliTRDtrack *iotrackTRD=0;
464 trdTree.Branch("tracks","AliTRDtrack",&iotrackTRD,32000,0);
5443e65e 465
466 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
029cd327 467 Float_t foundMin = fgkMinClustersInTrack * timeBins;
5443e65e 468
469 if (inp) {
470 TFile *in=(TFile*)inp;
471 if (!in->IsOpen()) {
c630aafd 472 cerr<<
473 "AliTRDtracker::Clusters2Tracks(): file with seeds is not open !\n";
5443e65e 474 cerr<<" ... going for seeds finding inside the TRD\n";
475 }
476 else {
477 in->cd();
e24ea474 478 sprintf(tname,"TRDb_%d",GetEventNumber());
5443e65e 479 TTree *seedTree=(TTree*)in->Get(tname);
480 if (!seedTree) {
a9814c09 481 cerr<<"AliTRDtracker::Clusters2Tracks(): ";
482 cerr<<"can't get a tree with track seeds !\n";
483 return 3;
5443e65e 484 }
485 AliTRDtrack *seed=new AliTRDtrack;
486 seedTree->SetBranchAddress("tracks",&seed);
487
488 Int_t n=(Int_t)seedTree->GetEntries();
489 for (Int_t i=0; i<n; i++) {
a9814c09 490 seedTree->GetEvent(i);
491 seed->ResetCovariance();
492 AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
493 fSeeds->AddLast(tr);
494 fNseeds++;
5443e65e 495 }
496 delete seed;
fd621f36 497 delete seedTree;
5443e65e 498 }
499 }
46d29e70 500
5443e65e 501 out->cd();
46d29e70 502
fd621f36 503
5443e65e 504 // find tracks from loaded seeds
a819a5f7 505
5443e65e 506 Int_t nseed=fSeeds->GetEntriesFast();
507 Int_t i, found = 0;
508 Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
509
510 for (i=0; i<nseed; i++) {
511 AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
512 FollowProlongation(t, innerTB);
513 if (t.GetNumberOfClusters() >= foundMin) {
514 UseClusters(&t);
029cd327 515 CookLabel(pt, 1-fgkLabelFraction);
5443e65e 516 // t.CookdEdx();
517 }
029cd327 518 iotrackTRD = pt;
519 trdTree.Fill();
e24ea474 520 found++;
521// cout<<found<<'\r';
5443e65e 522
523 if(PropagateToTPC(t)) {
524 AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
525 iotrack = tpc;
029cd327 526 tpcTree.Fill();
5443e65e 527 delete tpc;
528 }
529 delete fSeeds->RemoveAt(i);
530 fNseeds--;
531 }
a819a5f7 532
5443e65e 533 cout<<"Number of loaded seeds: "<<nseed<<endl;
534 cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
a819a5f7 535
5443e65e 536 // after tracks from loaded seeds are found and the corresponding
537 // clusters are used, look for additional seeds from TRD
538
539 if(fAddTRDseeds) {
540 // Find tracks for the seeds in the TRD
541 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
542
029cd327 543 Int_t nSteps = (Int_t) (fgkSeedDepth / fgkSeedStep);
544 Int_t gap = (Int_t) (timeBins * fgkSeedGap);
545 Int_t step = (Int_t) (timeBins * fgkSeedStep);
5443e65e 546
547 // make a first turn with tight cut on initial curvature
548 for(Int_t turn = 1; turn <= 2; turn++) {
549 if(turn == 2) {
029cd327 550 nSteps = (Int_t) (fgkSeedDepth / (3*fgkSeedStep));
551 step = (Int_t) (timeBins * (3*fgkSeedStep));
5443e65e 552 }
553 for(Int_t i=0; i<nSteps; i++) {
a9814c09 554 Int_t outer=timeBins-1-i*step;
555 Int_t inner=outer-gap;
5443e65e 556
a9814c09 557 nseed=fSeeds->GetEntriesFast();
5443e65e 558
a9814c09 559 MakeSeeds(inner, outer, turn);
5443e65e 560
a9814c09 561 nseed=fSeeds->GetEntriesFast();
562 printf("\n turn %d, step %d: number of seeds for TRD inward %d\n",
563 turn, i, nseed);
564
565 for (Int_t i=0; i<nseed; i++) {
566 AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
567 FollowProlongation(t,innerTB);
568 if (t.GetNumberOfClusters() >= foundMin) {
569 UseClusters(&t);
029cd327 570 CookLabel(pt, 1-fgkLabelFraction);
a9814c09 571 t.CookdEdx();
c630aafd 572 found++;
e24ea474 573// cout<<found<<'\r';
029cd327 574 iotrackTRD = pt;
575 trdTree.Fill();
a9814c09 576 if(PropagateToTPC(t)) {
577 AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
578 iotrack = tpc;
029cd327 579 tpcTree.Fill();
a9814c09 580 delete tpc;
581 }
582 }
583 delete fSeeds->RemoveAt(i);
584 fNseeds--;
585 }
5443e65e 586 }
a819a5f7 587 }
5443e65e 588 }
029cd327 589 tpcTree.Write();
590 trdTree.Write();
5443e65e 591
592 cout<<"Total number of found tracks: "<<found<<endl;
593
594 UnloadEvent();
595
596 savedir->cd();
597
598 return 0;
599}
c630aafd 600
601//___________________________________________________________________
602Int_t AliTRDtracker::Clusters2Tracks(AliESD* event)
603{
604 //
605 // Finds tracks within the TRD. The ESD event is expected to contain seeds
606 // at the outer part of the TRD. The seeds
607 // are found within the TRD if fAddTRDseeds is TRUE.
608 // The tracks are propagated to the innermost time bin
609 // of the TRD and the ESD event is updated
610 //
611
612 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
029cd327 613 Float_t foundMin = fgkMinClustersInTrack * timeBins;
c630aafd 614 Int_t nseed = 0;
615 Int_t found = 0;
616 Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
617
618 Int_t n = event->GetNumberOfTracks();
619 for (Int_t i=0; i<n; i++) {
620 AliESDtrack* seed=event->GetTrack(i);
621 ULong_t status=seed->GetStatus();
622 if ( (status & AliESDtrack::kTRDout ) == 0 ) continue;
623 if ( (status & AliESDtrack::kTRDin) != 0 ) continue;
624 nseed++;
625
626 AliTRDtrack* seed2 = new AliTRDtrack(*seed);
627 seed2->ResetCovariance();
628 AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
629 AliTRDtrack &t=*pt;
630 FollowProlongation(t, innerTB);
631 if (t.GetNumberOfClusters() >= foundMin) {
632 UseClusters(&t);
029cd327 633 CookLabel(pt, 1-fgkLabelFraction);
c630aafd 634 // t.CookdEdx();
635 }
636 found++;
637// cout<<found<<'\r';
638
639 if(PropagateToTPC(t)) {
640 seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
641 }
642 delete seed2;
643 delete pt;
644 }
645
646 cout<<"Number of loaded seeds: "<<nseed<<endl;
647 cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
648
649 // after tracks from loaded seeds are found and the corresponding
650 // clusters are used, look for additional seeds from TRD
651
652 if(fAddTRDseeds) {
653 // Find tracks for the seeds in the TRD
654 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
655
029cd327 656 Int_t nSteps = (Int_t) (fgkSeedDepth / fgkSeedStep);
657 Int_t gap = (Int_t) (timeBins * fgkSeedGap);
658 Int_t step = (Int_t) (timeBins * fgkSeedStep);
c630aafd 659
660 // make a first turn with tight cut on initial curvature
661 for(Int_t turn = 1; turn <= 2; turn++) {
662 if(turn == 2) {
029cd327 663 nSteps = (Int_t) (fgkSeedDepth / (3*fgkSeedStep));
664 step = (Int_t) (timeBins * (3*fgkSeedStep));
c630aafd 665 }
666 for(Int_t i=0; i<nSteps; i++) {
667 Int_t outer=timeBins-1-i*step;
668 Int_t inner=outer-gap;
669
670 nseed=fSeeds->GetEntriesFast();
671
672 MakeSeeds(inner, outer, turn);
673
674 nseed=fSeeds->GetEntriesFast();
675 printf("\n turn %d, step %d: number of seeds for TRD inward %d\n",
676 turn, i, nseed);
677
678 for (Int_t i=0; i<nseed; i++) {
679 AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
680 FollowProlongation(t,innerTB);
681 if (t.GetNumberOfClusters() >= foundMin) {
682 UseClusters(&t);
029cd327 683 CookLabel(pt, 1-fgkLabelFraction);
c630aafd 684 t.CookdEdx();
685 found++;
686// cout<<found<<'\r';
687 if(PropagateToTPC(t)) {
688 AliESDtrack track;
689 track.UpdateTrackParams(pt,AliESDtrack::kTRDin);
690 event->AddTrack(&track);
691 }
692 }
693 delete fSeeds->RemoveAt(i);
694 fNseeds--;
695 }
696 }
697 }
698 }
699
700 cout<<"Total number of found tracks: "<<found<<endl;
701
702 return 0;
703}
5443e65e 704
c630aafd 705
5443e65e 706
707//_____________________________________________________________________________
708Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) {
709 //
710 // Reads seeds from file <inp>. The seeds are AliTPCtrack's found and
711 // backpropagated by the TPC tracker. Each seed is first propagated
712 // to the TRD, and then its prolongation is searched in the TRD.
713 // If sufficiently long continuation of the track is found in the TRD
714 // the track is updated, otherwise it's stored as originaly defined
715 // by the TPC tracker.
716 //
717
5443e65e 718 LoadEvent();
719
720 TDirectory *savedir=gDirectory;
721
722 TFile *in=(TFile*)inp;
723
724 if (!in->IsOpen()) {
725 cerr<<"AliTRDtracker::PropagateBack(): ";
726 cerr<<"file with back propagated TPC tracks is not open !\n";
727 return 1;
728 }
729
730 if (!out->IsOpen()) {
731 cerr<<"AliTRDtracker::PropagateBack(): ";
732 cerr<<"file for back propagated TRD tracks is not open !\n";
733 return 2;
734 }
735
736 in->cd();
737 char tname[100];
e24ea474 738 sprintf(tname,"seedsTPCtoTRD_%d",GetEventNumber());
5443e65e 739 TTree *seedTree=(TTree*)in->Get(tname);
740 if (!seedTree) {
741 cerr<<"AliTRDtracker::PropagateBack(): ";
742 cerr<<"can't get a tree with seeds from TPC !\n";
743 cerr<<"check if your version of TPC tracker creates tree "<<tname<<"\n";
744 return 3;
745 }
746
747 AliTPCtrack *seed=new AliTPCtrack;
748 seedTree->SetBranchAddress("tracks",&seed);
749
750 Int_t n=(Int_t)seedTree->GetEntries();
751 for (Int_t i=0; i<n; i++) {
752 seedTree->GetEvent(i);
753 Int_t lbl = seed->GetLabel();
754 AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
755 tr->SetSeedLabel(lbl);
756 fSeeds->AddLast(tr);
757 fNseeds++;
758 }
759
760 delete seed;
761 delete seedTree;
762
763 out->cd();
764
765 AliTPCtrack *otrack=0;
766
e24ea474 767 sprintf(tname,"seedsTRDtoTOF1_%d",GetEventNumber());
5443e65e 768 TTree tofTree1(tname,"Tracks back propagated through TPC and TRD");
769 tofTree1.Branch("tracks","AliTPCtrack",&otrack,32000,0);
770
e24ea474 771 sprintf(tname,"seedsTRDtoTOF2_%d",GetEventNumber());
5443e65e 772 TTree tofTree2(tname,"Tracks back propagated through TPC and TRD");
773 tofTree2.Branch("tracks","AliTPCtrack",&otrack,32000,0);
774
e24ea474 775 sprintf(tname,"seedsTRDtoPHOS_%d",GetEventNumber());
5443e65e 776 TTree phosTree(tname,"Tracks back propagated through TPC and TRD");
777 phosTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
778
e24ea474 779 sprintf(tname,"seedsTRDtoRICH_%d",GetEventNumber());
5443e65e 780 TTree richTree(tname,"Tracks back propagated through TPC and TRD");
781 richTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
782
e24ea474 783 sprintf(tname,"TRDb_%d",GetEventNumber());
5443e65e 784 TTree trdTree(tname,"Back propagated TRD tracks at outer TRD time bin");
029cd327 785 AliTRDtrack *otrackTRD=0;
786 trdTree.Branch("tracks","AliTRDtrack",&otrackTRD,32000,0);
5443e65e 787
9c9d2487 788 if (IsStoringBarrel()) SetBarrelTree("back");
789 out->cd();
5443e65e 790
9c9d2487 791 Int_t found=0;
5443e65e 792 Int_t nseed=fSeeds->GetEntriesFast();
793
029cd327 794 // Float_t foundMin = fgkMinClustersInTrack * fTimeBinsPerPlane * fGeom->Nplan();
9c9d2487 795 Float_t foundMin = 40;
5443e65e 796
029cd327 797 Int_t outermostTB = fTrSec[0]->GetOuterTimeBin();
5443e65e 798
799 for (Int_t i=0; i<nseed; i++) {
800
801 AliTRDtrack *ps=(AliTRDtrack*)fSeeds->UncheckedAt(i), &s=*ps;
802 Int_t expectedClr = FollowBackProlongation(s);
9c9d2487 803
804 if (IsStoringBarrel()) {
029cd327 805 StoreBarrelTrack(ps, fgkLastPlane, kTrackBack);
9c9d2487 806 fBarrelTree->Fill();
807 }
808
5443e65e 809 Int_t foundClr = s.GetNumberOfClusters();
029cd327 810 Int_t lastTB = fTrSec[0]->GetLayerNumber(s.GetX());
5443e65e 811
fd621f36 812 // printf("seed %d: found %d out of %d expected clusters, Min is %f\n",
a9814c09 813 // i, foundClr, expectedClr, foundMin);
5443e65e 814
815 if (foundClr >= foundMin) {
a9814c09 816 if(foundClr >= 2) {
817 s.CookdEdx();
029cd327 818 CookLabel(ps, 1-fgkLabelFraction);
a9814c09 819 UseClusters(ps);
820 }
0d5b5c27 821
822 // Propagate to outer reference plane [SR, GSI, 18.02.2003]
823 ps->PropagateTo(364.8);
029cd327 824 otrackTRD=ps;
5443e65e 825 trdTree.Fill();
e24ea474 826 found++;
827// cout<<found<<'\r';
46d29e70 828 }
829
029cd327 830 if(((expectedClr < 10) && (lastTB == outermostTB)) ||
5443e65e 831 ((expectedClr >= 10) &&
a9814c09 832 (((Float_t) foundClr) / ((Float_t) expectedClr) >=
029cd327 833 fgkMinFractionOfFoundClusters) && (lastTB == outermostTB))) {
46d29e70 834
029cd327 835 Double_t xTOF = 375.5;
5443e65e 836
029cd327 837 if(PropagateToOuterPlane(s,xTOF)) {
a9814c09 838 AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
839 otrack = pt;
840 tofTree1.Fill();
841 delete pt;
5443e65e 842
029cd327 843 xTOF = 381.5;
5443e65e 844
029cd327 845 if(PropagateToOuterPlane(s,xTOF)) {
a9814c09 846 AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
847 otrack = pt;
848 tofTree2.Fill();
849 delete pt;
850
029cd327 851 Double_t xPHOS = 460.;
a9814c09 852
029cd327 853 if(PropagateToOuterPlane(s,xPHOS)) {
a9814c09 854 AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
855 otrack = pt;
856 phosTree.Fill();
857 delete pt;
858
029cd327 859 Double_t xRICH = 490+1.267;
a9814c09 860
029cd327 861 if(PropagateToOuterPlane(s,xRICH)) {
a9814c09 862 AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
863 otrack = pt;
864 richTree.Fill();
865 delete pt;
866 }
867 }
868 }
5443e65e 869 }
870 }
871 }
872
9c9d2487 873
874 out->cd();
5443e65e 875 tofTree1.Write();
876 tofTree2.Write();
877 phosTree.Write();
878 richTree.Write();
879 trdTree.Write();
880
9c9d2487 881 if (IsStoringBarrel()) { // [SR, 03.04.2003]
882 fBarrelFile->cd();
883 fBarrelTree->Write();
884 fBarrelFile->Flush();
885 }
886
5443e65e 887 savedir->cd();
888 cerr<<"Number of seeds: "<<nseed<<endl;
889 cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
46d29e70 890
5443e65e 891 UnloadEvent();
a819a5f7 892
5443e65e 893 return 0;
894
895}
896
c630aafd 897//_____________________________________________________________________________
898Int_t AliTRDtracker::PropagateBack(AliESD* event) {
899 //
900 // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
901 // backpropagated by the TPC tracker. Each seed is first propagated
902 // to the TRD, and then its prolongation is searched in the TRD.
903 // If sufficiently long continuation of the track is found in the TRD
904 // the track is updated, otherwise it's stored as originaly defined
905 // by the TPC tracker.
906 //
907
908 Int_t found=0;
909 Float_t foundMin = 40;
910
911 Int_t n = event->GetNumberOfTracks();
912 for (Int_t i=0; i<n; i++) {
913 AliESDtrack* seed=event->GetTrack(i);
914 ULong_t status=seed->GetStatus();
915 if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
916 if ( (status & AliESDtrack::kTRDout) != 0 ) continue;
917
918 Int_t lbl = seed->GetLabel();
919 AliTRDtrack *track = new AliTRDtrack(*seed);
920 track->SetSeedLabel(lbl);
921 fNseeds++;
922
923 /*Int_t expectedClr = */FollowBackProlongation(*track);
924
925 Int_t foundClr = track->GetNumberOfClusters();
926 if (foundClr >= foundMin) {
927 if(foundClr >= 2) {
928 track->CookdEdx();
029cd327 929// CookLabel(track, 1-fgkLabelFraction);
c630aafd 930 UseClusters(track);
931 }
932
933 // Propagate to outer reference plane [SR, GSI, 18.02.2003]
934// track->PropagateTo(364.8); why?
935
936 //seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
937 //found++;
938 }
939
940 if (track->PropagateTo(376.)) { //Propagation to the TOF (I.Belikov)
941 seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
942 found++;
943 }
944
945 }
946
947 cerr<<"Number of seeds: "<<fNseeds<<endl;
948 cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
949
950 return 0;
951
952}
953
5443e65e 954
955//---------------------------------------------------------------------------
956Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
957{
958 // Starting from current position on track=t this function tries
959 // to extrapolate the track up to timeBin=0 and to confirm prolongation
960 // if a close cluster is found. Returns the number of clusters
961 // expected to be found in sensitive layers
962
963 Float_t wIndex, wTB, wChi2;
964 Float_t wYrt, wYclosest, wYcorrect, wYwindow;
965 Float_t wZrt, wZclosest, wZcorrect, wZwindow;
966 Float_t wPx, wPy, wPz, wC;
029cd327 967 Double_t px, py, pz;
5443e65e 968 Float_t wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
969
970 Int_t trackIndex = t.GetLabel();
971
972 Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
973
029cd327 974 Int_t tryAgain=fMaxGap;
5443e65e 975
976 Double_t alpha=t.GetAlpha();
c630aafd 977 alpha = TVector2::Phi_0_2pi(alpha);
5443e65e 978
979 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
029cd327 980 Double_t radLength, rho, x, dx, y, ymax, z;
5443e65e 981
982 Int_t expectedNumberOfClusters = 0;
983 Bool_t lookForCluster;
984
985 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
986
987
988 for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) {
989
990 y = t.GetY(); z = t.GetZ();
991
992 // first propagate to the inner surface of the current time bin
029cd327 993 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
5443e65e 994 x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
029cd327 995 if(!t.PropagateTo(x,radLength,rho)) break;
5443e65e 996 y = t.GetY();
997 ymax = x*TMath::Tan(0.5*alpha);
998 if (y > ymax) {
999 s = (s+1) % ns;
1000 if (!t.Rotate(alpha)) break;
029cd327 1001 if(!t.PropagateTo(x,radLength,rho)) break;
5443e65e 1002 } else if (y <-ymax) {
1003 s = (s-1+ns) % ns;
1004 if (!t.Rotate(-alpha)) break;
029cd327 1005 if(!t.PropagateTo(x,radLength,rho)) break;
5443e65e 1006 }
1007
1008 y = t.GetY(); z = t.GetZ();
1009
1010 // now propagate to the middle plane of the next time bin
029cd327 1011 fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
5443e65e 1012 x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
029cd327 1013 if(!t.PropagateTo(x,radLength,rho)) break;
5443e65e 1014 y = t.GetY();
1015 ymax = x*TMath::Tan(0.5*alpha);
1016 if (y > ymax) {
1017 s = (s+1) % ns;
1018 if (!t.Rotate(alpha)) break;
029cd327 1019 if(!t.PropagateTo(x,radLength,rho)) break;
5443e65e 1020 } else if (y <-ymax) {
1021 s = (s-1+ns) % ns;
1022 if (!t.Rotate(-alpha)) break;
029cd327 1023 if(!t.PropagateTo(x,radLength,rho)) break;
5443e65e 1024 }
a819a5f7 1025
a819a5f7 1026
5443e65e 1027 if(lookForCluster) {
a819a5f7 1028
5443e65e 1029 expectedNumberOfClusters++;
5443e65e 1030 wIndex = (Float_t) t.GetLabel();
1031 wTB = nr;
a819a5f7 1032
029cd327 1033 AliTRDpropagationLayer& timeBin=*(fTrSec[s]->GetLayer(nr-1));
46d29e70 1034
5443e65e 1035 Double_t sy2=ExpectedSigmaY2(x,t.GetTgl(),t.GetPt());
5443e65e 1036 Double_t sz2=ExpectedSigmaZ2(x,t.GetTgl());
1037
1038 Double_t road;
1039 if((t.GetSigmaY2() + sy2) > 0) road=10.*sqrt(t.GetSigmaY2() + sy2);
1040 else return expectedNumberOfClusters;
1041
1042 wYrt = (Float_t) y;
1043 wZrt = (Float_t) z;
1044 wYwindow = (Float_t) road;
029cd327 1045 t.GetPxPyPz(px,py,pz);
1046 wPx = (Float_t) px;
1047 wPy = (Float_t) py;
1048 wPz = (Float_t) pz;
5443e65e 1049 wC = (Float_t) t.GetC();
1050 wSigmaC2 = (Float_t) t.GetSigmaC2();
1051 wSigmaTgl2 = (Float_t) t.GetSigmaTgl2();
1052 wSigmaY2 = (Float_t) t.GetSigmaY2();
1053 wSigmaZ2 = (Float_t) t.GetSigmaZ2();
1054 wChi2 = -1;
a819a5f7 1055
029cd327 1056 if (road>fgkWideRoad) {
a9814c09 1057 if (t.GetNumberOfClusters()>4)
1058 cerr<<t.GetNumberOfClusters()
1059 <<"FindProlongation warning: Too broad road !\n";
1060 return 0;
5443e65e 1061 }
1062
1063 AliTRDcluster *cl=0;
1064 UInt_t index=0;
1065
029cd327 1066 Double_t maxChi2=fgkMaxChi2;
5443e65e 1067
1068 wYclosest = 12345678;
1069 wYcorrect = 12345678;
1070 wZclosest = 12345678;
1071 wZcorrect = 12345678;
1072 wZwindow = TMath::Sqrt(2.25 * 12 * sz2);
1073
1074 // Find the closest correct cluster for debugging purposes
029cd327 1075 if (timeBin) {
a9814c09 1076 Float_t minDY = 1000000;
029cd327 1077 for (Int_t i=0; i<timeBin; i++) {
1078 AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
a9814c09 1079 if((c->GetLabel(0) != trackIndex) &&
1080 (c->GetLabel(1) != trackIndex) &&
1081 (c->GetLabel(2) != trackIndex)) continue;
1082 if(TMath::Abs(c->GetY() - y) > minDY) continue;
1083 minDY = TMath::Abs(c->GetY() - y);
1084 wYcorrect = c->GetY();
1085 wZcorrect = c->GetZ();
1086
1087 Double_t h01 = GetTiltFactor(c);
1088 wChi2 = t.GetPredictedChi2(c, h01);
1089 }
5443e65e 1090 }
a819a5f7 1091
5443e65e 1092 // Now go for the real cluster search
a819a5f7 1093
029cd327 1094 if (timeBin) {
5443e65e 1095
029cd327 1096 for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
1097 AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
a9814c09 1098 if (c->GetY() > y+road) break;
1099 if (c->IsUsed() > 0) continue;
1100 if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
1101
1102 Double_t h01 = GetTiltFactor(c);
1103 Double_t chi2=t.GetPredictedChi2(c,h01);
1104
029cd327 1105 if (chi2 > maxChi2) continue;
1106 maxChi2=chi2;
a9814c09 1107 cl=c;
029cd327 1108 index=timeBin.GetIndex(i);
a9814c09 1109 }
1110
1111 if(!cl) {
1112
029cd327 1113 for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
1114 AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
a9814c09 1115
1116 if (c->GetY() > y+road) break;
1117 if (c->IsUsed() > 0) continue;
1118 if((c->GetZ()-z)*(c->GetZ()-z) > 12 * sz2) continue;
1119
1120 Double_t h01 = GetTiltFactor(c);
1121 Double_t chi2=t.GetPredictedChi2(c, h01);
1122
029cd327 1123 if (chi2 > maxChi2) continue;
1124 maxChi2=chi2;
a9814c09 1125 cl=c;
029cd327 1126 index=timeBin.GetIndex(i);
a9814c09 1127 }
1128 }
1129
1130
1131 if (cl) {
1132 wYclosest = cl->GetY();
1133 wZclosest = cl->GetZ();
1134 Double_t h01 = GetTiltFactor(cl);
1135
1136 t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
029cd327 1137 if(!t.Update(cl,maxChi2,index,h01)) {
1138 if(!tryAgain--) return 0;
a9814c09 1139 }
029cd327 1140 else tryAgain=fMaxGap;
a9814c09 1141 }
1142 else {
029cd327 1143 if (tryAgain==0) break;
1144 tryAgain--;
a9814c09 1145 }
1146
1147 /*
1148 if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
1149
1150 printf(" %f", wIndex); //1
1151 printf(" %f", wTB); //2
1152 printf(" %f", wYrt); //3
1153 printf(" %f", wYclosest); //4
1154 printf(" %f", wYcorrect); //5
1155 printf(" %f", wYwindow); //6
1156 printf(" %f", wZrt); //7
1157 printf(" %f", wZclosest); //8
1158 printf(" %f", wZcorrect); //9
1159 printf(" %f", wZwindow); //10
1160 printf(" %f", wPx); //11
1161 printf(" %f", wPy); //12
1162 printf(" %f", wPz); //13
1163 printf(" %f", wSigmaC2*1000000); //14
1164 printf(" %f", wSigmaTgl2*1000); //15
1165 printf(" %f", wSigmaY2); //16
1166 // printf(" %f", wSigmaZ2); //17
1167 printf(" %f", wChi2); //17
1168 printf(" %f", wC); //18
1169 printf("\n");
1170 }
1171 */
a819a5f7 1172 }
5443e65e 1173 }
1174 }
1175 return expectedNumberOfClusters;
1176
1177
1178}
1179
1180//___________________________________________________________________
1181
1182Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
1183{
1184 // Starting from current radial position of track <t> this function
1185 // extrapolates the track up to outer timebin and in the sensitive
1186 // layers confirms prolongation if a close cluster is found.
1187 // Returns the number of clusters expected to be found in sensitive layers
1188
1189 Float_t wIndex, wTB, wChi2;
1190 Float_t wYrt, wYclosest, wYcorrect, wYwindow;
1191 Float_t wZrt, wZclosest, wZcorrect, wZwindow;
1192 Float_t wPx, wPy, wPz, wC;
029cd327 1193 Double_t px, py, pz;
5443e65e 1194 Float_t wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
1195
1196 Int_t trackIndex = t.GetLabel();
029cd327 1197 Int_t tryAgain=fMaxGap;
5443e65e 1198
1199 Double_t alpha=t.GetAlpha();
9c9d2487 1200 TVector2::Phi_0_2pi(alpha);
5443e65e 1201
9c9d2487 1202 Int_t s;
5443e65e 1203
1204 Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
029cd327 1205 Double_t radLength, rho, x, dx, y, ymax = 0, z;
5443e65e 1206 Bool_t lookForCluster;
1207
1208 Int_t expectedNumberOfClusters = 0;
1209 x = t.GetX();
1210
1211 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
1212
029cd327 1213 Int_t nRefPlane = fgkFirstPlane;
9c9d2487 1214 Bool_t isNewLayer = kFALSE;
5443e65e 1215
9c9d2487 1216 Double_t chi2;
1217 Double_t minDY;
5443e65e 1218
9c9d2487 1219 for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB+1; nr++) {
1220
1221 y = t.GetY();
1222 z = t.GetZ();
5443e65e 1223
1224 // first propagate to the outer surface of the current time bin
1225
9c9d2487 1226 s = t.GetSector();
029cd327 1227 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
9c9d2487 1228 x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2;
1229 y = t.GetY();
1230 z = t.GetZ();
5443e65e 1231
029cd327 1232 if(!t.PropagateTo(x,radLength,rho)) break;
9c9d2487 1233 if (!AdjustSector(&t)) break;
1234 s = t.GetSector();
029cd327 1235 if (!t.PropagateTo(x,radLength,rho)) break;
5443e65e 1236
1237 y = t.GetY();
9c9d2487 1238 z = t.GetZ();
5443e65e 1239
9c9d2487 1240 // Barrel Tracks [SR, 04.04.2003]
5443e65e 1241
9c9d2487 1242 s = t.GetSector();
1243 if (fTrSec[s]->GetLayer(nr)->IsSensitive() !=
1244 fTrSec[s]->GetLayer(nr+1)->IsSensitive() ) {
5443e65e 1245
c630aafd 1246// if (IsStoringBarrel()) StoreBarrelTrack(&t, nRefPlane++, kTrackBack);
9c9d2487 1247 }
5443e65e 1248
9c9d2487 1249 if (fTrSec[s]->GetLayer(nr-1)->IsSensitive() &&
1250 ! fTrSec[s]->GetLayer(nr)->IsSensitive()) {
1251 isNewLayer = kTRUE;
1252 } else {isNewLayer = kFALSE;}
5443e65e 1253
1254 y = t.GetY();
9c9d2487 1255 z = t.GetZ();
5443e65e 1256
9c9d2487 1257 // now propagate to the middle plane of the next time bin
029cd327 1258 fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
5443e65e 1259
9c9d2487 1260 x = fTrSec[s]->GetLayer(nr+1)->GetX();
029cd327 1261 if(!t.PropagateTo(x,radLength,rho)) break;
9c9d2487 1262 if (!AdjustSector(&t)) break;
1263 s = t.GetSector();
029cd327 1264 if(!t.PropagateTo(x,radLength,rho)) break;
5443e65e 1265
9c9d2487 1266 y = t.GetY();
1267 z = t.GetZ();
1268
1269 if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax);
5443e65e 1270 // printf("label %d, pl %d, lookForCluster %d \n",
a9814c09 1271 // trackIndex, nr+1, lookForCluster);
5443e65e 1272
1273 if(lookForCluster) {
1274 expectedNumberOfClusters++;
1275
1276 wIndex = (Float_t) t.GetLabel();
1277 wTB = fTrSec[s]->GetLayer(nr+1)->GetTimeBinIndex();
1278
029cd327 1279 AliTRDpropagationLayer& timeBin=*(fTrSec[s]->GetLayer(nr+1));
5443e65e 1280 Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
1281 Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
1282 if((t.GetSigmaY2() + sy2) < 0) break;
1283 Double_t road = 10.*sqrt(t.GetSigmaY2() + sy2);
1284 Double_t y=t.GetY(), z=t.GetZ();
1285
1286 wYrt = (Float_t) y;
1287 wZrt = (Float_t) z;
1288 wYwindow = (Float_t) road;
029cd327 1289 t.GetPxPyPz(px,py,pz);
1290 wPx = (Float_t) px;
1291 wPy = (Float_t) py;
1292 wPz = (Float_t) pz;
5443e65e 1293 wC = (Float_t) t.GetC();
1294 wSigmaC2 = (Float_t) t.GetSigmaC2();
1295 wSigmaTgl2 = (Float_t) t.GetSigmaTgl2();
1296 wSigmaY2 = (Float_t) t.GetSigmaY2();
1297 wSigmaZ2 = (Float_t) t.GetSigmaZ2();
1298 wChi2 = -1;
a819a5f7 1299
029cd327 1300 if (road>fgkWideRoad) {
a9814c09 1301 if (t.GetNumberOfClusters()>4)
1302 cerr<<t.GetNumberOfClusters()
1303 <<"FindProlongation warning: Too broad road !\n";
1304 return 0;
5443e65e 1305 }
1306
1307 AliTRDcluster *cl=0;
1308 UInt_t index=0;
1309
029cd327 1310 Double_t maxChi2=fgkMaxChi2;
5443e65e 1311
9c9d2487 1312 if (isNewLayer) {
1313 road = 3 * road;
1314 //sz2 = 3 * sz2;
029cd327 1315 maxChi2 = 10 * fgkMaxChi2;
9c9d2487 1316 }
1317
029cd327 1318 if (nRefPlane == fgkFirstPlane) maxChi2 = 20 * fgkMaxChi2;
1319 if (nRefPlane == fgkFirstPlane+2) maxChi2 = 15 * fgkMaxChi2;
1320 if (t.GetNRotate() > 0) maxChi2 = 3 * maxChi2;
9c9d2487 1321
1322
5443e65e 1323 wYclosest = 12345678;
1324 wYcorrect = 12345678;
1325 wZclosest = 12345678;
1326 wZcorrect = 12345678;
1327 wZwindow = TMath::Sqrt(2.25 * 12 * sz2);
1328
1329 // Find the closest correct cluster for debugging purposes
029cd327 1330 if (timeBin) {
9c9d2487 1331 minDY = 1000000;
029cd327 1332 for (Int_t i=0; i<timeBin; i++) {
1333 AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
a9814c09 1334 if((c->GetLabel(0) != trackIndex) &&
1335 (c->GetLabel(1) != trackIndex) &&
1336 (c->GetLabel(2) != trackIndex)) continue;
1337 if(TMath::Abs(c->GetY() - y) > minDY) continue;
9c9d2487 1338 //minDY = TMath::Abs(c->GetY() - y);
1339 minDY = c->GetY() - y;
a9814c09 1340 wYcorrect = c->GetY();
1341 wZcorrect = c->GetZ();
1342
1343 Double_t h01 = GetTiltFactor(c);
1344 wChi2 = t.GetPredictedChi2(c, h01);
1345 }
5443e65e 1346 }
1347
1348 // Now go for the real cluster search
1349
029cd327 1350 if (timeBin) {
a819a5f7 1351
029cd327 1352 for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
1353 AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
a9814c09 1354 if (c->GetY() > y+road) break;
9c9d2487 1355 if (c->IsUsed() > 0) continue;
a9814c09 1356 if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
1357
1358 Double_t h01 = GetTiltFactor(c);
9c9d2487 1359 chi2=t.GetPredictedChi2(c,h01);
a9814c09 1360
029cd327 1361 if (chi2 > maxChi2) continue;
1362 maxChi2=chi2;
a9814c09 1363 cl=c;
029cd327 1364 index=timeBin.GetIndex(i);
9c9d2487 1365
1366 //check is correct
1367 if((c->GetLabel(0) != trackIndex) &&
1368 (c->GetLabel(1) != trackIndex) &&
1369 (c->GetLabel(2) != trackIndex)) t.AddNWrong();
a9814c09 1370 }
5443e65e 1371
a9814c09 1372 if(!cl) {
1373
029cd327 1374 for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
1375 AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
a9814c09 1376
1377 if (c->GetY() > y+road) break;
9c9d2487 1378 if (c->IsUsed() > 0) continue;
a9814c09 1379 if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
1380
1381 Double_t h01 = GetTiltFactor(c);
9c9d2487 1382 chi2=t.GetPredictedChi2(c,h01);
a9814c09 1383
029cd327 1384 if (chi2 > maxChi2) continue;
1385 maxChi2=chi2;
a9814c09 1386 cl=c;
029cd327 1387 index=timeBin.GetIndex(i);
a9814c09 1388 }
1389 }
1390
1391 if (cl) {
1392 wYclosest = cl->GetY();
1393 wZclosest = cl->GetZ();
1394
1395 t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
1396 Double_t h01 = GetTiltFactor(cl);
029cd327 1397 if(!t.Update(cl,maxChi2,index,h01)) {
1398 if(!tryAgain--) return 0;
a9814c09 1399 }
029cd327 1400 else tryAgain=fMaxGap;
a9814c09 1401 }
1402 else {
029cd327 1403 if (tryAgain==0) break;
1404 tryAgain--;
9c9d2487 1405
1406 //if (minDY < 1000000 && isNewLayer)
1407 //cout << "\t" << nRefPlane << "\t" << "\t" << t.GetNRotate() << "\t" <<
029cd327 1408 // road << "\t" << minDY << "\t" << chi2 << "\t" << wChi2 << "\t" << maxChi2 << endl;
9c9d2487 1409
a9814c09 1410 }
1411
9c9d2487 1412 isNewLayer = kFALSE;
1413
a9814c09 1414 /*
1415 if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
1416
1417 printf(" %f", wIndex); //1
1418 printf(" %f", wTB); //2
1419 printf(" %f", wYrt); //3
1420 printf(" %f", wYclosest); //4
1421 printf(" %f", wYcorrect); //5
1422 printf(" %f", wYwindow); //6
1423 printf(" %f", wZrt); //7
1424 printf(" %f", wZclosest); //8
1425 printf(" %f", wZcorrect); //9
1426 printf(" %f", wZwindow); //10
1427 printf(" %f", wPx); //11
1428 printf(" %f", wPy); //12
1429 printf(" %f", wPz); //13
1430 printf(" %f", wSigmaC2*1000000); //14
1431 printf(" %f", wSigmaTgl2*1000); //15
1432 printf(" %f", wSigmaY2); //16
1433 // printf(" %f", wSigmaZ2); //17
1434 printf(" %f", wChi2); //17
1435 printf(" %f", wC); //18
1436 printf("\n");
1437 }
1438 */
46d29e70 1439 }
5443e65e 1440 }
46d29e70 1441 }
5443e65e 1442 return expectedNumberOfClusters;
9c9d2487 1443
1444
5443e65e 1445}
46d29e70 1446
5443e65e 1447//___________________________________________________________________
46d29e70 1448
5443e65e 1449Int_t AliTRDtracker::PropagateToOuterPlane(AliTRDtrack& t, Double_t xToGo)
1450{
1451 // Starting from current radial position of track <t> this function
1452 // extrapolates the track up to radial position <xToGo>.
1453 // Returns 1 if track reaches the plane, and 0 otherwise
46d29e70 1454
5443e65e 1455 Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
46d29e70 1456
5443e65e 1457 Double_t alpha=t.GetAlpha();
46d29e70 1458
5443e65e 1459 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1460 if (alpha < 0. ) alpha += 2.*TMath::Pi();
bbf92647 1461
5443e65e 1462 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
bbf92647 1463
5443e65e 1464 Bool_t lookForCluster;
029cd327 1465 Double_t radLength, rho, x, dx, y, ymax, z;
46d29e70 1466
5443e65e 1467 x = t.GetX();
46d29e70 1468
5443e65e 1469 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
1470
1471 Int_t plToGo = fTrSec[0]->GetLayerNumber(xToGo);
46d29e70 1472
5443e65e 1473 for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr<plToGo; nr++) {
46d29e70 1474
5443e65e 1475 y = t.GetY(); z = t.GetZ();
1476
1477 // first propagate to the outer surface of the current time bin
029cd327 1478 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
5443e65e 1479 x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
029cd327 1480 if(!t.PropagateTo(x,radLength,rho)) return 0;
5443e65e 1481 y = t.GetY();
1482 ymax = x*TMath::Tan(0.5*alpha);
1483 if (y > ymax) {
1484 s = (s+1) % ns;
1485 if (!t.Rotate(alpha)) return 0;
1486 } else if (y <-ymax) {
1487 s = (s-1+ns) % ns;
1488 if (!t.Rotate(-alpha)) return 0;
1489 }
029cd327 1490 if(!t.PropagateTo(x,radLength,rho)) return 0;
5443e65e 1491
1492 y = t.GetY(); z = t.GetZ();
1493
1494 // now propagate to the middle plane of the next time bin
029cd327 1495 fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
5443e65e 1496 x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
029cd327 1497 if(!t.PropagateTo(x,radLength,rho)) return 0;
5443e65e 1498 y = t.GetY();
1499 ymax = x*TMath::Tan(0.5*alpha);
1500 if (y > ymax) {
1501 s = (s+1) % ns;
1502 if (!t.Rotate(alpha)) return 0;
1503 } else if (y <-ymax) {
1504 s = (s-1+ns) % ns;
1505 if (!t.Rotate(-alpha)) return 0;
1506 }
029cd327 1507 if(!t.PropagateTo(x,radLength,rho)) return 0;
46d29e70 1508 }
5443e65e 1509 return 1;
1510}
46d29e70 1511
5443e65e 1512//___________________________________________________________________
46d29e70 1513
5443e65e 1514Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t)
1515{
1516 // Starting from current radial position of track <t> this function
1517 // extrapolates the track up to radial position of the outermost
1518 // padrow of the TPC.
1519 // Returns 1 if track reaches the TPC, and 0 otherwise
1520
c630aafd 1521 //Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
5443e65e 1522
1523 Double_t alpha=t.GetAlpha();
c630aafd 1524 alpha = TVector2::Phi_0_2pi(alpha);
5443e65e 1525
1526 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1527
1528 Bool_t lookForCluster;
029cd327 1529 Double_t radLength, rho, x, dx, y, /*ymax,*/ z;
5443e65e 1530
1531 x = t.GetX();
1532
1533 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
5443e65e 1534 Int_t plTPC = fTrSec[0]->GetLayerNumber(246.055);
1535
1536 for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr>plTPC; nr--) {
1537
9c9d2487 1538 y = t.GetY();
1539 z = t.GetZ();
5443e65e 1540
1541 // first propagate to the outer surface of the current time bin
029cd327 1542 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
9c9d2487 1543 x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2;
1544
029cd327 1545 if(!t.PropagateTo(x,radLength,rho)) return 0;
9c9d2487 1546 AdjustSector(&t);
029cd327 1547 if(!t.PropagateTo(x,radLength,rho)) return 0;
5443e65e 1548
9c9d2487 1549 y = t.GetY();
1550 z = t.GetZ();
5443e65e 1551
1552 // now propagate to the middle plane of the next time bin
029cd327 1553 fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
9c9d2487 1554 x = fTrSec[s]->GetLayer(nr-1)->GetX();
1555
029cd327 1556 if(!t.PropagateTo(x,radLength,rho)) return 0;
9c9d2487 1557 AdjustSector(&t);
029cd327 1558 if(!t.PropagateTo(x,radLength,rho)) return 0;
5443e65e 1559 }
1560 return 1;
1561}
46d29e70 1562
5443e65e 1563void AliTRDtracker::LoadEvent()
46d29e70 1564{
1565 // Fills clusters into TRD tracking_sectors
1566 // Note that the numbering scheme for the TRD tracking_sectors
1567 // differs from that of TRD sectors
46d29e70 1568
5443e65e 1569 ReadClusters(fClusters);
46d29e70 1570 Int_t ncl=fClusters->GetEntriesFast();
5443e65e 1571 cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1572
46d29e70 1573 UInt_t index;
1574 while (ncl--) {
e24ea474 1575// printf("\r %d left ",ncl);
46d29e70 1576 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
029cd327 1577 Int_t detector=c->GetDetector();
1578 Int_t localTimeBin=c->GetLocalTimeBin();
46d29e70 1579 Int_t sector=fGeom->GetSector(detector);
5443e65e 1580 Int_t plane=fGeom->GetPlane(detector);
46d29e70 1581
029cd327 1582 Int_t trackingSector = CookSectorIndex(sector);
46d29e70 1583
029cd327 1584 Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
b3a5a838 1585 if(gtb < 0) continue;
029cd327 1586 Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
a819a5f7 1587
5443e65e 1588 index=ncl;
029cd327 1589 fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
46d29e70 1590 }
1591 printf("\r\n");
46d29e70 1592
5443e65e 1593}
46d29e70 1594
1595//_____________________________________________________________________________
c630aafd 1596Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1597{
1598 // Fills clusters into TRD tracking_sectors
1599 // Note that the numbering scheme for the TRD tracking_sectors
1600 // differs from that of TRD sectors
1601
1602 if (ReadClusters(fClusters,cTree)) {
1603 Error("LoadClusters","Problem with reading the clusters !");
1604 return 1;
1605 }
1606 Int_t ncl=fClusters->GetEntriesFast();
1607 cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1608
1609 UInt_t index;
1610 while (ncl--) {
1611// printf("\r %d left ",ncl);
1612 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
029cd327 1613 Int_t detector=c->GetDetector();
1614 Int_t localTimeBin=c->GetLocalTimeBin();
c630aafd 1615 Int_t sector=fGeom->GetSector(detector);
1616 Int_t plane=fGeom->GetPlane(detector);
1617
029cd327 1618 Int_t trackingSector = CookSectorIndex(sector);
c630aafd 1619
029cd327 1620 Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
c630aafd 1621 if(gtb < 0) continue;
029cd327 1622 Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
c630aafd 1623
1624 index=ncl;
029cd327 1625 fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
c630aafd 1626 }
1627 printf("\r\n");
1628
1629 return 0;
1630}
1631
1632//_____________________________________________________________________________
5443e65e 1633void AliTRDtracker::UnloadEvent()
1634{
0a29d0f1 1635 //
5443e65e 1636 // Clears the arrays of clusters and tracks. Resets sectors and timebins
0a29d0f1 1637 //
46d29e70 1638
5443e65e 1639 Int_t i, nentr;
46d29e70 1640
5443e65e 1641 nentr = fClusters->GetEntriesFast();
1642 for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i);
46d29e70 1643
5443e65e 1644 nentr = fSeeds->GetEntriesFast();
1645 for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i);
46d29e70 1646
5443e65e 1647 nentr = fTracks->GetEntriesFast();
1648 for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i);
46d29e70 1649
5443e65e 1650 Int_t nsec = AliTRDgeometry::kNsect;
46d29e70 1651
5443e65e 1652 for (i = 0; i < nsec; i++) {
1653 for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1654 fTrSec[i]->GetLayer(pl)->Clear();
1655 }
1656 }
46d29e70 1657
5443e65e 1658}
a819a5f7 1659
5443e65e 1660//__________________________________________________________________________
1661void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn)
1662{
1663 // Creates track seeds using clusters in timeBins=i1,i2
a819a5f7 1664
5443e65e 1665 if(turn > 2) {
1666 cerr<<"MakeSeeds: turn "<<turn<<" exceeds the limit of 2"<<endl;
a819a5f7 1667 return;
1668 }
1669
5443e65e 1670 Double_t x[5], c[15];
029cd327 1671 Int_t maxSec=AliTRDgeometry::kNsect;
5443e65e 1672
1673 Double_t alpha=AliTRDgeometry::GetAlpha();
1674 Double_t shift=AliTRDgeometry::GetAlpha()/2.;
1675 Double_t cs=cos(alpha), sn=sin(alpha);
1676 Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
1677
1678
1679 Int_t i2 = fTrSec[0]->GetLayerNumber(inner);
1680 Int_t i1 = fTrSec[0]->GetLayerNumber(outer);
1681
1682 Double_t x1 =fTrSec[0]->GetX(i1);
1683 Double_t xx2=fTrSec[0]->GetX(i2);
1684
029cd327 1685 for (Int_t ns=0; ns<maxSec; ns++) {
5443e65e 1686
029cd327 1687 Int_t nl2 = *(fTrSec[(ns-2+maxSec)%maxSec]->GetLayer(i2));
1688 Int_t nl=(*fTrSec[(ns-1+maxSec)%maxSec]->GetLayer(i2));
5443e65e 1689 Int_t nm=(*fTrSec[ns]->GetLayer(i2));
029cd327 1690 Int_t nu=(*fTrSec[(ns+1)%maxSec]->GetLayer(i2));
1691 Int_t nu2=(*fTrSec[(ns+2)%maxSec]->GetLayer(i2));
5443e65e 1692
1693 AliTRDpropagationLayer& r1=*(fTrSec[ns]->GetLayer(i1));
1694
46d29e70 1695 for (Int_t is=0; is < r1; is++) {
1696 Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
5443e65e 1697
a819a5f7 1698 for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
a9814c09 1699
1700 const AliTRDcluster *cl;
1701 Double_t x2, y2, z2;
1702 Double_t x3=0., y3=0.;
1703
1704 if (js<nl2) {
1705 if(turn != 2) continue;
029cd327 1706 AliTRDpropagationLayer& r2=*(fTrSec[(ns-2+maxSec)%maxSec]->GetLayer(i2));
a9814c09 1707 cl=r2[js];
1708 y2=cl->GetY(); z2=cl->GetZ();
1709
1710 x2= xx2*cs2+y2*sn2;
1711 y2=-xx2*sn2+y2*cs2;
1712 }
1713 else if (js<nl2+nl) {
1714 if(turn != 1) continue;
029cd327 1715 AliTRDpropagationLayer& r2=*(fTrSec[(ns-1+maxSec)%maxSec]->GetLayer(i2));
a9814c09 1716 cl=r2[js-nl2];
1717 y2=cl->GetY(); z2=cl->GetZ();
1718
1719 x2= xx2*cs+y2*sn;
1720 y2=-xx2*sn+y2*cs;
1721 }
1722 else if (js<nl2+nl+nm) {
1723 if(turn != 1) continue;
1724 AliTRDpropagationLayer& r2=*(fTrSec[ns]->GetLayer(i2));
1725 cl=r2[js-nl2-nl];
1726 x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
1727 }
1728 else if (js<nl2+nl+nm+nu) {
1729 if(turn != 1) continue;
029cd327 1730 AliTRDpropagationLayer& r2=*(fTrSec[(ns+1)%maxSec]->GetLayer(i2));
a9814c09 1731 cl=r2[js-nl2-nl-nm];
1732 y2=cl->GetY(); z2=cl->GetZ();
1733
1734 x2=xx2*cs-y2*sn;
1735 y2=xx2*sn+y2*cs;
1736 }
1737 else {
1738 if(turn != 2) continue;
029cd327 1739 AliTRDpropagationLayer& r2=*(fTrSec[(ns+2)%maxSec]->GetLayer(i2));
a9814c09 1740 cl=r2[js-nl2-nl-nm-nu];
1741 y2=cl->GetY(); z2=cl->GetZ();
1742
1743 x2=xx2*cs2-y2*sn2;
1744 y2=xx2*sn2+y2*cs2;
1745 }
1746
029cd327 1747 if(TMath::Abs(z1-z2) > fgkMaxSeedDeltaZ12) continue;
a9814c09 1748
1749 Double_t zz=z1 - z1/x1*(x1-x2);
1750
029cd327 1751 if (TMath::Abs(zz-z2)>fgkMaxSeedDeltaZ) continue;
a9814c09 1752
1753 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
1754 if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
1755
1756 x[0]=y1;
1757 x[1]=z1;
1758 x[4]=f1trd(x1,y1,x2,y2,x3,y3);
1759
029cd327 1760 if (TMath::Abs(x[4]) > fgkMaxSeedC) continue;
a9814c09 1761
1762 x[2]=f2trd(x1,y1,x2,y2,x3,y3);
1763
1764 if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1765
1766 x[3]=f3trd(x1,y1,x2,y2,z1,z2);
1767
029cd327 1768 if (TMath::Abs(x[3]) > fgkMaxSeedTan) continue;
a9814c09 1769
1770 Double_t a=asin(x[2]);
1771 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
1772
029cd327 1773 if (TMath::Abs(zv)>fgkMaxSeedVertexZ) continue;
a9814c09 1774
1775 Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
1776 Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2();
029cd327 1777 Double_t sy3=fgkSeedErrorSY3, sy=fgkSeedErrorSY, sz=fgkSeedErrorSZ;
a9814c09 1778
1779 // Tilt changes
1780 Double_t h01 = GetTiltFactor(r1[is]);
029cd327 1781 Double_t xuFactor = 100.;
b8dc2353 1782 if(fNoTilt) {
1783 h01 = 0;
029cd327 1784 xuFactor = 1;
b8dc2353 1785 }
1786
fd621f36 1787 sy1=sy1+sz1*h01*h01;
b3a5a838 1788 Double_t syz=sz1*(-h01);
a9814c09 1789 // end of tilt changes
1790
b3a5a838 1791 Double_t f40=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1792 Double_t f42=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1793 Double_t f43=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1794 Double_t f20=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1795 Double_t f22=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1796 Double_t f23=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1797 Double_t f30=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1798 Double_t f31=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1799 Double_t f32=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1800 Double_t f34=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1801
a9814c09 1802
b3a5a838 1803 c[0]=sy1;
a9814c09 1804 // c[1]=0.; c[2]=sz1;
029cd327 1805 c[1]=syz; c[2]=sz1*xuFactor;
b3a5a838 1806 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1807 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
1808 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1809 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1810 c[13]=f30*sy1*f40+f32*sy2*f42;
1811 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
a9814c09 1812
1813 UInt_t index=r1.GetIndex(is);
1814
1815 AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
1816
1817 Int_t rc=FollowProlongation(*track, i2);
1818
1819 if ((rc < 1) ||
1820 (track->GetNumberOfClusters() <
029cd327 1821 (outer-inner)*fgkMinClustersInSeed)) delete track;
a9814c09 1822 else {
1823 fSeeds->AddLast(track); fNseeds++;
e24ea474 1824// cerr<<"\r found seed "<<fNseeds;
a9814c09 1825 }
46d29e70 1826 }
1827 }
1828 }
5443e65e 1829}
46d29e70 1830
a819a5f7 1831//_____________________________________________________________________________
c630aafd 1832Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *ClusterTree)
46d29e70 1833{
a819a5f7 1834 //
1835 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
1836 // from the file. The names of the cluster tree and branches
1837 // should match the ones used in AliTRDclusterizer::WriteClusters()
1838 //
029cd327 1839 TObjArray *clusterArray = new TObjArray(400);
5443e65e 1840
c630aafd 1841 TBranch *branch=ClusterTree->GetBranch("TRDcluster");
1842 if (!branch) {
1843 Error("ReadClusters","Can't get the branch !");
1844 return 1;
1845 }
029cd327 1846 branch->SetAddress(&clusterArray);
5443e65e 1847
1848 Int_t nEntries = (Int_t) ClusterTree->GetEntries();
1849 printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
a819a5f7 1850
a819a5f7 1851 // Loop through all entries in the tree
1852 Int_t nbytes;
1853 AliTRDcluster *c = 0;
5443e65e 1854 printf("\n");
a819a5f7 1855
1856 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
1857
1858 // Import the tree
5443e65e 1859 nbytes += ClusterTree->GetEvent(iEntry);
1860
a819a5f7 1861 // Get the number of points in the detector
029cd327 1862 Int_t nCluster = clusterArray->GetEntriesFast();
e24ea474 1863// printf("\r Read %d clusters from entry %d", nCluster, iEntry);
5443e65e 1864
a819a5f7 1865 // Loop through all TRD digits
1866 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
029cd327 1867 c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
a819a5f7 1868 AliTRDcluster *co = new AliTRDcluster(*c);
1869 co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
5443e65e 1870 Int_t ltb = co->GetLocalTimeBin();
b8dc2353 1871 if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
1872 else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
a819a5f7 1873 array->AddLast(co);
029cd327 1874 delete clusterArray->RemoveAt(iCluster);
a819a5f7 1875 }
1876 }
1877
029cd327 1878 delete clusterArray;
5443e65e 1879
c630aafd 1880 return 0;
a819a5f7 1881}
1882
5443e65e 1883//______________________________________________________________________
1884void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename)
a819a5f7 1885{
1886 //
5443e65e 1887 // Reads AliTRDclusters from file <filename>. The names of the cluster
1888 // tree and branches should match the ones used in
1889 // AliTRDclusterizer::WriteClusters()
1890 // if <array> == 0, clusters are added into AliTRDtracker fCluster array
a819a5f7 1891 //
1892
5443e65e 1893 TDirectory *savedir=gDirectory;
46d29e70 1894
5443e65e 1895 TFile *file = TFile::Open(filename);
1896 if (!file->IsOpen()) {
1897 cerr<<"Can't open file with TRD clusters"<<endl;
1898 return;
1899 }
46d29e70 1900
5443e65e 1901 Char_t treeName[12];
e24ea474 1902 sprintf(treeName,"TreeR%d_TRD",GetEventNumber());
029cd327 1903 TTree *clusterTree = (TTree*) gDirectory->Get(treeName);
46d29e70 1904
029cd327 1905 if (!clusterTree) {
5443e65e 1906 cerr<<"AliTRDtracker::ReadClusters(): ";
1907 cerr<<"can't get a tree with clusters !\n";
1908 return;
1909 }
46d29e70 1910
029cd327 1911 TObjArray *clusterArray = new TObjArray(400);
46d29e70 1912
029cd327 1913 clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
46d29e70 1914
029cd327 1915 Int_t nEntries = (Int_t) clusterTree->GetEntries();
1916 cout<<"found "<<nEntries<<" in clusterTree"<<endl;
a819a5f7 1917
5443e65e 1918 // Loop through all entries in the tree
1919 Int_t nbytes;
1920 AliTRDcluster *c = 0;
bbf92647 1921
5443e65e 1922 printf("\n");
bbf92647 1923
5443e65e 1924 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
46d29e70 1925
5443e65e 1926 // Import the tree
029cd327 1927 nbytes += clusterTree->GetEvent(iEntry);
0a29d0f1 1928
5443e65e 1929 // Get the number of points in the detector
029cd327 1930 Int_t nCluster = clusterArray->GetEntriesFast();
b3a5a838 1931 printf("\n Read %d clusters from entry %d", nCluster, iEntry);
5443e65e 1932
1933 // Loop through all TRD digits
1934 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
029cd327 1935 c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
5443e65e 1936 AliTRDcluster *co = new AliTRDcluster(*c);
1937 co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
1938 Int_t ltb = co->GetLocalTimeBin();
b8dc2353 1939 if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
1940 else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
5443e65e 1941 array->AddLast(co);
029cd327 1942 delete clusterArray->RemoveAt(iCluster);
5443e65e 1943 }
46d29e70 1944 }
0a29d0f1 1945
5443e65e 1946 file->Close();
029cd327 1947 delete clusterArray;
5443e65e 1948 savedir->cd();
1949
1950}
1951
c630aafd 1952void AliTRDtracker::ReadClusters(TObjArray *array, const TFile *inp)
1953{
1954 //
1955 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
1956 // from the file. The names of the cluster tree and branches
1957 // should match the ones used in AliTRDclusterizer::WriteClusters()
1958 //
1959
1960 TDirectory *savedir=gDirectory;
1961
1962 if (inp) {
1963 TFile *in=(TFile*)inp;
1964 if (!in->IsOpen()) {
1965 cerr<<"AliTRDtracker::ReadClusters(): input file is not open !\n";
1966 return;
1967 }
1968 else{
1969 in->cd();
1970 }
1971 }
1972
1973 Char_t treeName[12];
1974 sprintf(treeName,"TreeR%d_TRD",GetEventNumber());
029cd327 1975 TTree *clusterTree = (TTree*) gDirectory->Get(treeName);
c630aafd 1976
029cd327 1977 TObjArray *clusterArray = new TObjArray(400);
c630aafd 1978
029cd327 1979 clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
c630aafd 1980
029cd327 1981 Int_t nEntries = (Int_t) clusterTree->GetEntries();
1982 printf("found %d entries in %s.\n",nEntries,clusterTree->GetName());
c630aafd 1983
1984 // Loop through all entries in the tree
1985 Int_t nbytes;
1986 AliTRDcluster *c = 0;
1987 printf("\n");
1988
1989 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
1990
1991 // Import the tree
029cd327 1992 nbytes += clusterTree->GetEvent(iEntry);
c630aafd 1993
1994 // Get the number of points in the detector
029cd327 1995 Int_t nCluster = clusterArray->GetEntriesFast();
c630aafd 1996// printf("\r Read %d clusters from entry %d", nCluster, iEntry);
1997
1998 // Loop through all TRD digits
1999 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
029cd327 2000 c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
c630aafd 2001 AliTRDcluster *co = new AliTRDcluster(*c);
2002 co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
2003 Int_t ltb = co->GetLocalTimeBin();
2004 if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
2005 else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
2006 array->AddLast(co);
029cd327 2007 delete clusterArray->RemoveAt(iCluster);
c630aafd 2008 }
2009 }
2010
029cd327 2011 delete clusterArray;
c630aafd 2012 savedir->cd();
2013
2014}
46d29e70 2015
2016//__________________________________________________________________
029cd327 2017void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const
2018{
2019 //
2020 // This cooks a label. Mmmmh, smells good...
2021 //
46d29e70 2022
2023 Int_t label=123456789, index, i, j;
5443e65e 2024 Int_t ncl=pt->GetNumberOfClusters();
029cd327 2025 const Int_t kRange = fTrSec[0]->GetOuterTimeBin()+1;
5443e65e 2026
029cd327 2027 Bool_t labelAdded;
46d29e70 2028
029cd327 2029 // Int_t s[kRange][2];
2030 Int_t **s = new Int_t* [kRange];
2031 for (i=0; i<kRange; i++) {
d1b06c24 2032 s[i] = new Int_t[2];
2033 }
029cd327 2034 for (i=0; i<kRange; i++) {
46d29e70 2035 s[i][0]=-1;
2036 s[i][1]=0;
2037 }
2038
2039 Int_t t0,t1,t2;
2040 for (i=0; i<ncl; i++) {
5443e65e 2041 index=pt->GetClusterIndex(i);
46d29e70 2042 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
5443e65e 2043 t0=c->GetLabel(0);
2044 t1=c->GetLabel(1);
2045 t2=c->GetLabel(2);
46d29e70 2046 }
2047
2048 for (i=0; i<ncl; i++) {
5443e65e 2049 index=pt->GetClusterIndex(i);
46d29e70 2050 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2051 for (Int_t k=0; k<3; k++) {
5443e65e 2052 label=c->GetLabel(k);
029cd327 2053 labelAdded=kFALSE; j=0;
46d29e70 2054 if (label >= 0) {
029cd327 2055 while ( (!labelAdded) && ( j < kRange ) ) {
a9814c09 2056 if (s[j][0]==label || s[j][1]==0) {
2057 s[j][0]=label;
2058 s[j][1]=s[j][1]+1;
029cd327 2059 labelAdded=kTRUE;
a9814c09 2060 }
2061 j++;
2062 }
46d29e70 2063 }
2064 }
2065 }
2066
46d29e70 2067 Int_t max=0;
2068 label = -123456789;
2069
029cd327 2070 for (i=0; i<kRange; i++) {
46d29e70 2071 if (s[i][1]>max) {
2072 max=s[i][1]; label=s[i][0];
2073 }
2074 }
5443e65e 2075
029cd327 2076 for (i=0; i<kRange; i++) {
5443e65e 2077 delete []s[i];
2078 }
2079
d1b06c24 2080 delete []s;
5443e65e 2081
2082 if ((1.- Float_t(max)/ncl) > wrong) label=-label;
2083
2084 pt->SetLabel(label);
2085
46d29e70 2086}
2087
c630aafd 2088
5443e65e 2089//__________________________________________________________________
029cd327 2090void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const
2091{
2092 //
2093 // Use clusters, but don't abuse them!
2094 //
2095
5443e65e 2096 Int_t ncl=t->GetNumberOfClusters();
2097 for (Int_t i=from; i<ncl; i++) {
2098 Int_t index = t->GetClusterIndex(i);
2099 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2100 c->Use();
2101 }
2102}
2103
2104
2105//_____________________________________________________________________
029cd327 2106Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
5443e65e 2107{
2108 // Parametrised "expected" error of the cluster reconstruction in Y
2109
2110 Double_t s = 0.08 * 0.08;
2111 return s;
2112}
2113
2114//_____________________________________________________________________
029cd327 2115Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
5443e65e 2116{
2117 // Parametrised "expected" error of the cluster reconstruction in Z
2118
a9814c09 2119 Double_t s = 9 * 9 /12.;
5443e65e 2120 return s;
2121}
2122
5443e65e 2123//_____________________________________________________________________
029cd327 2124Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
5443e65e 2125{
2126 //
029cd327 2127 // Returns radial position which corresponds to time bin <localTB>
5443e65e 2128 // in tracking sector <sector> and plane <plane>
2129 //
2130
029cd327 2131 Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
5443e65e 2132 Int_t pl = fTrSec[sector]->GetLayerNumber(index);
2133 return fTrSec[sector]->GetLayer(pl)->GetX();
2134
2135}
2136
c630aafd 2137
5443e65e 2138//_______________________________________________________
2139AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x,
029cd327 2140 Double_t dx, Double_t rho, Double_t radLength, Int_t tbIndex)
5443e65e 2141{
2142 //
2143 // AliTRDpropagationLayer constructor
2144 //
2145
029cd327 2146 fN = 0; fX = x; fdX = dx; fRho = rho; fX0 = radLength;
2147 fClusters = NULL; fIndex = NULL; fTimeBinIndex = tbIndex;
5443e65e 2148
2149
029cd327 2150 for(Int_t i=0; i < (Int_t) kZones; i++) {
5443e65e 2151 fZc[i]=0; fZmax[i] = 0;
2152 }
2153
2154 fYmax = 0;
2155
2156 if(fTimeBinIndex >= 0) {
029cd327 2157 fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
2158 fIndex = new UInt_t[kMaxClusterPerTimeBin];
5443e65e 2159 }
2160
2161 fHole = kFALSE;
2162 fHoleZc = 0;
2163 fHoleZmax = 0;
2164 fHoleYc = 0;
2165 fHoleYmax = 0;
2166 fHoleRho = 0;
2167 fHoleX0 = 0;
2168
2169}
2170
2171//_______________________________________________________
2172void AliTRDtracker::AliTRDpropagationLayer::SetHole(
a9814c09 2173 Double_t Zmax, Double_t Ymax, Double_t rho,
029cd327 2174 Double_t radLength, Double_t Yc, Double_t Zc)
5443e65e 2175{
2176 //
2177 // Sets hole in the layer
2178 //
2179
2180 fHole = kTRUE;
2181 fHoleZc = Zc;
2182 fHoleZmax = Zmax;
2183 fHoleYc = Yc;
2184 fHoleYmax = Ymax;
2185 fHoleRho = rho;
029cd327 2186 fHoleX0 = radLength;
5443e65e 2187}
2188
2189
2190//_______________________________________________________
2191AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs, AliTRDparameter* par)
0a29d0f1 2192{
2193 //
5443e65e 2194 // AliTRDtrackingSector Constructor
0a29d0f1 2195 //
46d29e70 2196
5443e65e 2197 fGeom = geo;
2198 fPar = par;
2199 fGeomSector = gs;
2200 fTzeroShift = 0.13;
2201 fN = 0;
46d29e70 2202
029cd327 2203 for(UInt_t i=0; i < kMaxTimeBinIndex; i++) fTimeBinIndex[i] = -1;
5443e65e 2204
2205
2206 AliTRDpropagationLayer* ppl;
2207
029cd327 2208 Double_t x, xin, xout, dx, rho, radLength;
5443e65e 2209 Int_t steps;
2210
2211 // set time bins in the gas of the TPC
2212
2213 xin = 246.055; xout = 254.055; steps = 20; dx = (xout-xin)/steps;
029cd327 2214 rho = 0.9e-3; radLength = 28.94;
5443e65e 2215
2216 for(Int_t i=0; i<steps; i++) {
2217 x = xin + i*dx + dx/2;
029cd327 2218 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2219 InsertLayer(ppl);
a819a5f7 2220 }
5443e65e 2221
2222 // set time bins in the outer field cage vessel
2223
029cd327 2224 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2225 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2226 InsertLayer(ppl);
2227
029cd327 2228 dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2229 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2230 InsertLayer(ppl);
2231
029cd327 2232 dx = 2.; xin = xout; xout = xin + dx; rho = 1.45*0.02; radLength = 41.28; // Nomex
5443e65e 2233 steps = 5; dx = (xout - xin)/steps;
2234 for(Int_t i=0; i<steps; i++) {
2235 x = xin + i*dx + dx/2;
029cd327 2236 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2237 InsertLayer(ppl);
a819a5f7 2238 }
46d29e70 2239
029cd327 2240 dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2241 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2242 InsertLayer(ppl);
2243
029cd327 2244 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2245 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2246 InsertLayer(ppl);
46d29e70 2247
46d29e70 2248
5443e65e 2249 // set time bins in CO2
46d29e70 2250
5443e65e 2251 xin = xout; xout = 275.0;
2252 steps = 50; dx = (xout - xin)/steps;
029cd327 2253 rho = 1.977e-3; radLength = 36.2;
5443e65e 2254
2255 for(Int_t i=0; i<steps; i++) {
2256 x = xin + i*dx + dx/2;
029cd327 2257 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2258 InsertLayer(ppl);
46d29e70 2259 }
2260
5443e65e 2261 // set time bins in the outer containment vessel
46d29e70 2262
029cd327 2263 dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; radLength = 24.01; // Al
2264 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2265 InsertLayer(ppl);
46d29e70 2266
029cd327 2267 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2268 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2269 InsertLayer(ppl);
2270
029cd327 2271 dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2272 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2273 InsertLayer(ppl);
2274
029cd327 2275 dx = 3.; xin = xout; xout = xin + dx; rho = 1.45*0.02; radLength = 41.28; // Nomex
5443e65e 2276 steps = 10; dx = (xout - xin)/steps;
2277 for(Int_t i=0; i<steps; i++) {
2278 x = xin + i*dx + dx/2;
029cd327 2279 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2280 InsertLayer(ppl);
2281 }
2282
029cd327 2283 dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2284 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2285 InsertLayer(ppl);
2286
029cd327 2287 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2288 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2289 InsertLayer(ppl);
2290
029cd327 2291 dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; radLength = 24.01; // Al
2292 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2293 InsertLayer(ppl);
2294
2295 Double_t xtrd = (Double_t) fGeom->Rmin();
2296
2297 // add layers between TPC and TRD (Air temporarily)
2298 xin = xout; xout = xtrd;
2299 steps = 50; dx = (xout - xin)/steps;
029cd327 2300 rho = 1.2e-3; radLength = 36.66;
5443e65e 2301
2302 for(Int_t i=0; i<steps; i++) {
2303 x = xin + i*dx + dx/2;
029cd327 2304 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2305 InsertLayer(ppl);
2306 }
2307
2308
2309 Double_t alpha=AliTRDgeometry::GetAlpha();
2310
2311 // add layers for each of the planes
2312
2313 Double_t dxRo = (Double_t) fGeom->CroHght(); // Rohacell
2314 Double_t dxSpace = (Double_t) fGeom->Cspace(); // Spacing between planes
2315 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
2316 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
2317 Double_t dxRad = (Double_t) fGeom->CraHght(); // Radiator
2318 Double_t dxTEC = dxRad + dxDrift + dxAmp + dxRo;
2319 Double_t dxPlane = dxTEC + dxSpace;
2320
029cd327 2321 Int_t tb, tbIndex;
2322 const Int_t kNchambers = AliTRDgeometry::Ncham();
2323 Double_t ymax = 0, holeYmax = 0;
2324 Double_t *zc = new Double_t[kNchambers];
2325 Double_t *zmax = new Double_t[kNchambers];
5443e65e 2326 Double_t holeZmax = 1000.; // the whole sector is missing
2327
5443e65e 2328 for(Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
2329
2330 // Radiator
2331 xin = xtrd + plane * dxPlane; xout = xin + dxRad;
029cd327 2332 steps = 12; dx = (xout - xin)/steps; rho = 0.074; radLength = 40.6;
5443e65e 2333 for(Int_t i=0; i<steps; i++) {
2334 x = xin + i*dx + dx/2;
029cd327 2335 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2336 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
a9814c09 2337 holeYmax = x*TMath::Tan(0.5*alpha);
2338 ppl->SetHole(holeYmax, holeZmax);
5443e65e 2339 }
2340 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
a9814c09 2341 holeYmax = x*TMath::Tan(0.5*alpha);
2342 ppl->SetHole(holeYmax, holeZmax);
5443e65e 2343 }
2344 InsertLayer(ppl);
2345 }
2346
029cd327 2347 ymax = fGeom->GetChamberWidth(plane)/2;
2348 for(Int_t ch = 0; ch < kNchambers; ch++) {
2349 zmax[ch] = fGeom->GetChamberLength(plane,ch)/2;
5443e65e 2350 Float_t pad = fPar->GetRowPadSize(plane,ch,0);
2351 Float_t row0 = fPar->GetRow0(plane,ch,0);
2352 Int_t nPads = fPar->GetRowMax(plane,ch,0);
029cd327 2353 zc[ch] = (pad * nPads)/2 + row0 - pad/2;
5443e65e 2354 }
2355
2356 dx = fPar->GetTimeBinSize();
029cd327 2357 rho = 0.00295 * 0.85; radLength = 11.0;
5443e65e 2358
2359 Double_t x0 = (Double_t) fPar->GetTime0(plane);
2360 Double_t xbottom = x0 - dxDrift;
2361 Double_t xtop = x0 + dxAmp;
2362
2363 // Amplification region
2364
2365 steps = (Int_t) (dxAmp/dx);
2366
2367 for(tb = 0; tb < steps; tb++) {
2368 x = x0 + tb * dx + dx/2;
029cd327 2369 tbIndex = CookTimeBinIndex(plane, -tb-1);
2370 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
2371 ppl->SetYmax(ymax);
2372 for(Int_t ch = 0; ch < kNchambers; ch++) {
2373 ppl->SetZmax(ch, zc[ch], zmax[ch]);
5443e65e 2374 }
2375 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
a9814c09 2376 holeYmax = x*TMath::Tan(0.5*alpha);
2377 ppl->SetHole(holeYmax, holeZmax);
5443e65e 2378 }
2379 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
a9814c09 2380 holeYmax = x*TMath::Tan(0.5*alpha);
2381 ppl->SetHole(holeYmax, holeZmax);
5443e65e 2382 }
2383 InsertLayer(ppl);
2384 }
029cd327 2385 tbIndex = CookTimeBinIndex(plane, -steps);
5443e65e 2386 x = (x + dx/2 + xtop)/2;
2387 dx = 2*(xtop-x);
029cd327 2388 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
2389 ppl->SetYmax(ymax);
2390 for(Int_t ch = 0; ch < kNchambers; ch++) {
2391 ppl->SetZmax(ch, zc[ch], zmax[ch]);
5443e65e 2392 }
2393 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2394 holeYmax = x*TMath::Tan(0.5*alpha);
2395 ppl->SetHole(holeYmax, holeZmax);
2396 }
2397 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2398 holeYmax = x*TMath::Tan(0.5*alpha);
2399 ppl->SetHole(holeYmax, holeZmax);
2400 }
2401 InsertLayer(ppl);
2402
2403 // Drift region
2404 dx = fPar->GetTimeBinSize();
2405 steps = (Int_t) (dxDrift/dx);
2406
2407 for(tb = 0; tb < steps; tb++) {
2408 x = x0 - tb * dx - dx/2;
029cd327 2409 tbIndex = CookTimeBinIndex(plane, tb);
5443e65e 2410
029cd327 2411 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
2412 ppl->SetYmax(ymax);
2413 for(Int_t ch = 0; ch < kNchambers; ch++) {
2414 ppl->SetZmax(ch, zc[ch], zmax[ch]);
5443e65e 2415 }
2416 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
a9814c09 2417 holeYmax = x*TMath::Tan(0.5*alpha);
2418 ppl->SetHole(holeYmax, holeZmax);
5443e65e 2419 }
2420 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
a9814c09 2421 holeYmax = x*TMath::Tan(0.5*alpha);
2422 ppl->SetHole(holeYmax, holeZmax);
5443e65e 2423 }
2424 InsertLayer(ppl);
2425 }
029cd327 2426 tbIndex = CookTimeBinIndex(plane, steps);
5443e65e 2427 x = (x - dx/2 + xbottom)/2;
2428 dx = 2*(x-xbottom);
029cd327 2429 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
2430 ppl->SetYmax(ymax);
2431 for(Int_t ch = 0; ch < kNchambers; ch++) {
2432 ppl->SetZmax(ch, zc[ch], zmax[ch]);
5443e65e 2433 }
2434 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2435 holeYmax = x*TMath::Tan(0.5*alpha);
2436 ppl->SetHole(holeYmax, holeZmax);
2437 }
2438 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2439 holeYmax = x*TMath::Tan(0.5*alpha);
2440 ppl->SetHole(holeYmax, holeZmax);
2441 }
2442 InsertLayer(ppl);
2443
2444 // Pad Plane
029cd327 2445 xin = xtop; dx = 0.025; xout = xin + dx; rho = 1.7; radLength = 33.0;
2446 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2447 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2448 holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha);
2449 ppl->SetHole(holeYmax, holeZmax);
2450 }
2451 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2452 holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha);
2453 ppl->SetHole(holeYmax, holeZmax);
2454 }
2455 InsertLayer(ppl);
2456
2457 // Rohacell
2458 xin = xout; xout = xtrd + (plane + 1) * dxPlane - dxSpace;
029cd327 2459 steps = 5; dx = (xout - xin)/steps; rho = 0.074; radLength = 40.6;
5443e65e 2460 for(Int_t i=0; i<steps; i++) {
2461 x = xin + i*dx + dx/2;
029cd327 2462 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2463 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
a9814c09 2464 holeYmax = x*TMath::Tan(0.5*alpha);
2465 ppl->SetHole(holeYmax, holeZmax);
5443e65e 2466 }
2467 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
a9814c09 2468 holeYmax = x*TMath::Tan(0.5*alpha);
2469 ppl->SetHole(holeYmax, holeZmax);
5443e65e 2470 }
2471 InsertLayer(ppl);
2472 }
2473
2474 // Space between the chambers, air
2475 xin = xout; xout = xtrd + (plane + 1) * dxPlane;
029cd327 2476 steps = 5; dx = (xout - xin)/steps; rho = 1.29e-3; radLength = 36.66;
5443e65e 2477 for(Int_t i=0; i<steps; i++) {
2478 x = xin + i*dx + dx/2;
029cd327 2479 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2480 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
a9814c09 2481 holeYmax = x*TMath::Tan(0.5*alpha);
2482 ppl->SetHole(holeYmax, holeZmax);
5443e65e 2483 }
2484 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
a9814c09 2485 holeYmax = x*TMath::Tan(0.5*alpha);
2486 ppl->SetHole(holeYmax, holeZmax);
5443e65e 2487 }
2488 InsertLayer(ppl);
2489 }
2490 }
2491
2492 // Space between the TRD and RICH
2493 Double_t xRICH = 500.;
2494 xin = xout; xout = xRICH;
029cd327 2495 steps = 200; dx = (xout - xin)/steps; rho = 1.29e-3; radLength = 36.66;
5443e65e 2496 for(Int_t i=0; i<steps; i++) {
2497 x = xin + i*dx + dx/2;
029cd327 2498 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2499 InsertLayer(ppl);
2500 }
2501
2502 MapTimeBinLayers();
029cd327 2503 delete [] zc;
2504 delete [] zmax;
0a29d0f1 2505
46d29e70 2506}
2507
5443e65e 2508//______________________________________________________
2509
029cd327 2510Int_t AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t localTB) const
5443e65e 2511{
2512 //
2513 // depending on the digitization parameters calculates "global"
029cd327 2514 // time bin index for timebin <localTB> in plane <plane>
5443e65e 2515 //
2516
2517 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
2518 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
2519 Double_t dx = (Double_t) fPar->GetTimeBinSize();
2520
2521 Int_t tbAmp = fPar->GetTimeBefore();
2522 Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
b3a5a838 2523 if(kTRUE) maxAmp = 0; // intentional until we change parameter class
5443e65e 2524 Int_t tbDrift = fPar->GetTimeMax();
2525 Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
2526
029cd327 2527 Int_t tbPerPlane = TMath::Min(tbAmp,maxAmp) + TMath::Min(tbDrift,maxDrift);
5443e65e 2528
029cd327 2529 Int_t gtb = (plane+1) * tbPerPlane - localTB - 1 - TMath::Min(tbAmp,maxAmp);
5443e65e 2530
029cd327 2531 if((localTB < 0) &&
2532 (TMath::Abs(localTB) > TMath::Min(tbAmp,maxAmp))) return -1;
2533 if(localTB >= TMath::Min(tbDrift,maxDrift)) return -1;
5443e65e 2534
2535 return gtb;
2536
2537
2538}
2539
2540//______________________________________________________
2541
2542void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers()
2543{
2544 //
2545 // For all sensitive time bins sets corresponding layer index
2546 // in the array fTimeBins
2547 //
2548
2549 Int_t index;
2550
2551 for(Int_t i = 0; i < fN; i++) {
2552 index = fLayers[i]->GetTimeBinIndex();
2553
2554 // printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX());
2555
2556 if(index < 0) continue;
029cd327 2557 if(index >= (Int_t) kMaxTimeBinIndex) {
5443e65e 2558 printf("*** AliTRDtracker::MapTimeBinLayers: \n");
2559 printf(" index %d exceeds allowed maximum of %d!\n",
029cd327 2560 index, kMaxTimeBinIndex-1);
5443e65e 2561 continue;
2562 }
2563 fTimeBinIndex[index] = i;
2564 }
2565
2566 Double_t x1, dx1, x2, dx2, gap;
2567
2568 for(Int_t i = 0; i < fN-1; i++) {
2569 x1 = fLayers[i]->GetX();
2570 dx1 = fLayers[i]->GetdX();
2571 x2 = fLayers[i+1]->GetX();
2572 dx2 = fLayers[i+1]->GetdX();
2573 gap = (x2 - dx2/2) - (x1 + dx1/2);
2574 if(gap < -0.01) {
2575 printf("*** warning: layers %d and %d are overlayed:\n",i,i+1);
2576 printf(" %f + %f + %f > %f\n", x1, dx1/2, dx2/2, x2);
2577 }
2578 if(gap > 0.01) {
2579 printf("*** warning: layers %d and %d have a large gap:\n",i,i+1);
2580 printf(" (%f - %f) - (%f + %f) = %f\n",
a9814c09 2581 x2, dx2/2, x1, dx1, gap);
5443e65e 2582 }
2583 }
2584}
2585
2586
2587//______________________________________________________
2588
2589
2590Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const
2591{
2592 //
2593 // Returns the number of time bin which in radial position is closest to <x>
2594 //
2595
2596 if(x >= fLayers[fN-1]->GetX()) return fN-1;
2597 if(x <= fLayers[0]->GetX()) return 0;
2598
2599 Int_t b=0, e=fN-1, m=(b+e)/2;
2600 for (; b<e; m=(b+e)/2) {
2601 if (x > fLayers[m]->GetX()) b=m+1;
2602 else e=m;
2603 }
2604 if(TMath::Abs(x - fLayers[m]->GetX()) >
2605 TMath::Abs(x - fLayers[m+1]->GetX())) return m+1;
2606 else return m;
2607
2608}
2609
2610//______________________________________________________
2611
2612Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const
2613{
2614 //
2615 // Returns number of the innermost SENSITIVE propagation layer
2616 //
2617
2618 return GetLayerNumber(0);
2619}
2620
2621//______________________________________________________
2622
2623Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const
2624{
2625 //
2626 // Returns number of the outermost SENSITIVE time bin
2627 //
2628
2629 return GetLayerNumber(GetNumberOfTimeBins() - 1);
2630}
2631
2632//______________________________________________________
2633
2634Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const
2635{
2636 //
2637 // Returns number of SENSITIVE time bins
2638 //
2639
2640 Int_t tb, layer;
029cd327 2641 for(tb = kMaxTimeBinIndex-1; tb >=0; tb--) {
5443e65e 2642 layer = GetLayerNumber(tb);
2643 if(layer>=0) break;
2644 }
2645 return tb+1;
2646}
2647
2648//______________________________________________________
2649
2650void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl)
2651{
2652 //
2653 // Insert layer <pl> in fLayers array.
2654 // Layers are sorted according to X coordinate.
2655
029cd327 2656 if ( fN == ((Int_t) kMaxLayersPerSector)) {
5443e65e 2657 printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
2658 return;
2659 }
2660 if (fN==0) {fLayers[fN++] = pl; return;}
2661 Int_t i=Find(pl->GetX());
2662
2663 memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
2664 fLayers[i]=pl; fN++;
2665
2666}
2667
2668//______________________________________________________
2669
2670Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const
2671{
2672 //
2673 // Returns index of the propagation layer nearest to X
2674 //
2675
2676 if (x <= fLayers[0]->GetX()) return 0;
2677 if (x > fLayers[fN-1]->GetX()) return fN;
2678 Int_t b=0, e=fN-1, m=(b+e)/2;
2679 for (; b<e; m=(b+e)/2) {
2680 if (x > fLayers[m]->GetX()) b=m+1;
2681 else e=m;
2682 }
2683 return m;
2684}
2685
2686//______________________________________________________
2687
2688void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
029cd327 2689 Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &radLength,
a9814c09 2690 Bool_t &lookForCluster) const
5443e65e 2691{
2692 //
029cd327 2693 // Returns radial step <dx>, density <rho>, rad. length <radLength>,
5443e65e 2694 // and sensitivity <lookForCluster> in point <y,z>
2695 //
2696
2697 dx = fdX;
2698 rho = fRho;
029cd327 2699 radLength = fX0;
5443e65e 2700 lookForCluster = kFALSE;
2701
2702 // check dead regions
2703 if(fTimeBinIndex >= 0) {
029cd327 2704 for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
5443e65e 2705 if(TMath::Abs(z - fZc[ch]) < fZmax[ch])
a9814c09 2706 lookForCluster = kTRUE;
029cd327 2707 // else { rho = 1.7; radLength = 33.0; } // G10
5443e65e 2708 }
2709 if(TMath::Abs(y) > fYmax) lookForCluster = kFALSE;
2710 if(!lookForCluster) {
029cd327 2711 // rho = 1.7; radLength = 33.0; // G10
5443e65e 2712 }
2713 }
2714
2715 // check hole
2716 if(fHole && (TMath::Abs(y - fHoleYc) < fHoleYmax) &&
2717 (TMath::Abs(z - fHoleZc) < fHoleZmax)) {
2718 lookForCluster = kFALSE;
2719 rho = fHoleRho;
029cd327 2720 radLength = fHoleX0;
5443e65e 2721 }
2722
2723 return;
2724}
2725
2726//______________________________________________________
2727
2728void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c,
a9814c09 2729 UInt_t index) {
5443e65e 2730
2731// Insert cluster in cluster array.
2732// Clusters are sorted according to Y coordinate.
2733
2734 if(fTimeBinIndex < 0) {
2735 printf("*** attempt to insert cluster into non-sensitive time bin!\n");
2736 return;
2737 }
2738
029cd327 2739 if (fN== (Int_t) kMaxClusterPerTimeBin) {
5443e65e 2740 printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n");
2741 return;
2742 }
2743 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
2744 Int_t i=Find(c->GetY());
2745 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
2746 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
2747 fIndex[i]=index; fClusters[i]=c; fN++;
2748}
2749
2750//______________________________________________________
2751
2752Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Double_t y) const {
2753
2754// Returns index of the cluster nearest in Y
2755
2756 if (y <= fClusters[0]->GetY()) return 0;
2757 if (y > fClusters[fN-1]->GetY()) return fN;
2758 Int_t b=0, e=fN-1, m=(b+e)/2;
2759 for (; b<e; m=(b+e)/2) {
2760 if (y > fClusters[m]->GetY()) b=m+1;
2761 else e=m;
2762 }
2763 return m;
2764}
2765
fd621f36 2766//---------------------------------------------------------
5443e65e 2767
fd621f36 2768Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
2769//
2770// Returns correction factor for tilted pads geometry
2771//
5443e65e 2772
fd621f36 2773 Double_t h01 = sin(TMath::Pi() / 180.0 * fPar->GetTiltingAngle());
2774 Int_t det = c->GetDetector();
2775 Int_t plane = fGeom->GetPlane(det);
b3a5a838 2776
2777 if((plane == 1) || (plane == 3) || (plane == 5)) h01=-h01;
b8dc2353 2778
2779 if(fNoTilt) h01 = 0;
fd621f36 2780
2781 return h01;
2782}
5443e65e 2783
c630aafd 2784
2785
2786
2787
2788