]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDtracker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update of tracking code for tilted pads
[u/mrichter/AliRoot.git] / TRD / AliTRDtracker.cxx
CommitLineData
46d29e70 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17$Log$
fd621f36 18Revision 1.19 2002/10/22 15:53:08 alibrary
19Introducing Riostream.h
20
a2cb5b3d 21Revision 1.18 2002/10/14 14:57:44 hristov
22Merging the VirtualMC branch to the main development branch (HEAD)
23
b9d0a01d 24Revision 1.14.6.2 2002/07/24 10:09:31 alibrary
25Updating VirtualMC
26
27Revision 1.17 2002/06/13 12:09:58 hristov
28Minor corrections
29
ca6c93d6 30Revision 1.16 2002/06/12 09:54:36 cblume
31Update of tracking code provided by Sergei
32
0a29d0f1 33Revision 1.14 2001/11/14 10:50:46 cblume
34Changes in digits IO. Add merging of summable digits
35
abaf1f1d 36Revision 1.13 2001/05/30 12:17:47 hristov
37Loop variables declared once
38
3ab6f951 39Revision 1.12 2001/05/28 17:07:58 hristov
40Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
41
71d9fa7b 42Revision 1.8 2000/12/20 13:00:44 cblume
43Modifications for the HP-compiler
44
d1b06c24 45Revision 1.7 2000/12/08 16:07:02 cblume
46Update of the tracking by Sergei
47
bbf92647 48Revision 1.6 2000/11/30 17:38:08 cblume
49Changes to get in line with new STEER and EVGEN
50
1819f4bb 51Revision 1.5 2000/11/14 14:40:27 cblume
52Correction for the Sun compiler (kTRUE and kFALSE)
53
57527628 54Revision 1.4 2000/11/10 14:57:52 cblume
55Changes in the geometry constants for the DEC compiler
56
dd56b762 57Revision 1.3 2000/10/15 23:40:01 cblume
58Remove AliTRDconst
59
0e9c2ad5 60Revision 1.2 2000/10/06 16:49:46 cblume
61Made Getters const
62
46d29e70 63Revision 1.1.2.2 2000/10/04 16:34:58 cblume
64Replace include files by forward declarations
65
66Revision 1.1.2.1 2000/09/22 14:47:52 cblume
67Add the tracking code
68
69*/
bbf92647 70
a2cb5b3d 71#include <Riostream.h>
46d29e70 72
73#include <TFile.h>
46d29e70 74#include <TBranch.h>
5443e65e 75#include <TTree.h>
76#include <TObjArray.h>
46d29e70 77
46d29e70 78#include "AliTRDgeometry.h"
5443e65e 79#include "AliTRDparameter.h"
80#include "AliTRDgeometryDetail.h"
46d29e70 81#include "AliTRDcluster.h"
82#include "AliTRDtrack.h"
5443e65e 83#include "../TPC/AliTPCtrack.h"
46d29e70 84
85#include "AliTRDtracker.h"
86
87ClassImp(AliTRDtracker)
88
5443e65e 89 const Float_t AliTRDtracker::fSeedDepth = 0.5;
90 const Float_t AliTRDtracker::fSeedStep = 0.10;
91 const Float_t AliTRDtracker::fSeedGap = 0.25;
92
93 const Float_t AliTRDtracker::fMaxSeedDeltaZ12 = 40.;
94 const Float_t AliTRDtracker::fMaxSeedDeltaZ = 25.;
95 const Float_t AliTRDtracker::fMaxSeedC = 0.0052;
96 const Float_t AliTRDtracker::fMaxSeedTan = 1.2;
97 const Float_t AliTRDtracker::fMaxSeedVertexZ = 150.;
a819a5f7 98
5443e65e 99 const Double_t AliTRDtracker::fSeedErrorSY = 0.2;
100 const Double_t AliTRDtracker::fSeedErrorSY3 = 2.5;
101 const Double_t AliTRDtracker::fSeedErrorSZ = 0.1;
bbf92647 102
5443e65e 103 const Float_t AliTRDtracker::fMinClustersInSeed = 0.7;
bbf92647 104
5443e65e 105 const Float_t AliTRDtracker::fMinClustersInTrack = 0.5;
106 const Float_t AliTRDtracker::fMinFractionOfFoundClusters = 0.8;
bbf92647 107
5443e65e 108 const Float_t AliTRDtracker::fSkipDepth = 0.05;
109 const Float_t AliTRDtracker::fLabelFraction = 0.8;
110 const Float_t AliTRDtracker::fWideRoad = 20.;
111
112 const Double_t AliTRDtracker::fMaxChi2 = 24.;
a819a5f7 113
bbf92647 114
46d29e70 115//____________________________________________________________________
5443e65e 116AliTRDtracker::AliTRDtracker(const TFile *geomfile)
46d29e70 117{
5443e65e 118 //
119 // Main constructor
120 //
46d29e70 121
5443e65e 122 Float_t fTzero = 0;
123 Int_t version = 2;
124
125 fAddTRDseeds = kFALSE;
bbf92647 126
5443e65e 127 fGeom = NULL;
128
129 TDirectory *savedir=gDirectory;
130 TFile *in=(TFile*)geomfile;
131 if (!in->IsOpen()) {
132 printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n");
133 printf(" DETAIL TRD geometry and DEFAULT TRD parameter will be used\n");
134 }
135 else {
136 in->cd();
137 in->ls();
138 fGeom = (AliTRDgeometry*) in->Get("TRDgeometry");
139 fPar = (AliTRDparameter*) in->Get("TRDparameter");
140 fGeom->Dump();
141 }
46d29e70 142
5443e65e 143 if(fGeom) {
144 // fTzero = geo->GetT0();
145 fTzero = 0.;
146 version = fGeom->IsVersion();
147 printf("Found geometry version %d on file \n", version);
148 }
149 else {
150 printf("AliTRDtracker::AliTRDtracker(): cann't find TRD geometry!\n");
151 printf(" DETAIL TRD geometry and DEFAULT TRD parameter will be used\n");
152 fGeom = new AliTRDgeometryDetail();
153 fPar = new AliTRDparameter();
154 }
155
156 savedir->cd();
46d29e70 157
5443e65e 158
159 // fGeom->SetT0(fTzero);
0a29d0f1 160
46d29e70 161 fEvent = 0;
bbf92647 162
46d29e70 163 fNclusters = 0;
164 fClusters = new TObjArray(2000);
165 fNseeds = 0;
5443e65e 166 fSeeds = new TObjArray(2000);
46d29e70 167 fNtracks = 0;
5443e65e 168 fTracks = new TObjArray(1000);
a819a5f7 169
5443e65e 170 for(Int_t geom_s = 0; geom_s < kTRACKING_SECTORS; geom_s++) {
171 Int_t tr_s = CookSectorIndex(geom_s);
172 fTrSec[tr_s] = new AliTRDtrackingSector(fGeom, geom_s, fPar);
173 }
a819a5f7 174
5443e65e 175 fSY2corr = 0.025;
176 fSZ2corr = 12.;
bbf92647 177
5443e65e 178 // calculate max gap on track
a819a5f7 179
5443e65e 180 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
181 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
a819a5f7 182
5443e65e 183 Double_t dx = (Double_t) fPar->GetTimeBinSize();
184 Int_t tbAmp = fPar->GetTimeBefore();
185 Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
186 Int_t tbDrift = fPar->GetTimeMax();
187 Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
a819a5f7 188
5443e65e 189 tbDrift = TMath::Min(tbDrift,maxDrift);
190 tbAmp = TMath::Min(tbAmp,maxAmp);
46d29e70 191
5443e65e 192 fTimeBinsPerPlane = tbAmp + tbDrift;
193 fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fSkipDepth);
46d29e70 194
5443e65e 195 fVocal = kFALSE;
0a29d0f1 196
5443e65e 197}
46d29e70 198
5443e65e 199//___________________________________________________________________
200AliTRDtracker::~AliTRDtracker()
46d29e70 201{
5443e65e 202 delete fClusters;
203 delete fTracks;
204 delete fSeeds;
205 delete fGeom;
206 delete fPar;
0a29d0f1 207
5443e65e 208 for(Int_t geom_s = 0; geom_s < kTRACKING_SECTORS; geom_s++) {
209 delete fTrSec[geom_s];
210 }
211}
46d29e70 212
213//_____________________________________________________________________
5443e65e 214inline Double_t f1trd(Double_t x1,Double_t y1,
215 Double_t x2,Double_t y2,
216 Double_t x3,Double_t y3)
46d29e70 217{
0a29d0f1 218 //
46d29e70 219 // Initial approximation of the track curvature
0a29d0f1 220 //
46d29e70 221 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
222 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
223 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
224 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
225 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
226
227 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
228
229 return -xr*yr/sqrt(xr*xr+yr*yr);
230}
231
232//_____________________________________________________________________
5443e65e 233inline Double_t f2trd(Double_t x1,Double_t y1,
234 Double_t x2,Double_t y2,
235 Double_t x3,Double_t y3)
46d29e70 236{
0a29d0f1 237 //
5443e65e 238 // Initial approximation of the track curvature times X coordinate
239 // of the center of curvature
0a29d0f1 240 //
46d29e70 241
242 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
243 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
244 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
245 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
246 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
247
248 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
249
250 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
251}
252
253//_____________________________________________________________________
5443e65e 254inline Double_t f3trd(Double_t x1,Double_t y1,
255 Double_t x2,Double_t y2,
256 Double_t z1,Double_t z2)
46d29e70 257{
0a29d0f1 258 //
46d29e70 259 // Initial approximation of the tangent of the track dip angle
0a29d0f1 260 //
46d29e70 261
262 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
263}
264
46d29e70 265//___________________________________________________________________
5443e65e 266Int_t AliTRDtracker::Clusters2Tracks(const TFile *inp, TFile *out)
46d29e70 267{
0a29d0f1 268 //
5443e65e 269 // Finds tracks within the TRD. File <inp> is expected to contain seeds
270 // at the outer part of the TRD. If <inp> is NULL, the seeds
271 // are found within the TRD if fAddTRDseeds is TRUE.
272 // The tracks are propagated to the innermost time bin
273 // of the TRD and stored in file <out>.
0a29d0f1 274 //
a819a5f7 275
5443e65e 276 LoadEvent();
277
278 TDirectory *savedir=gDirectory;
a819a5f7 279
5443e65e 280 char tname[100];
a819a5f7 281
5443e65e 282 if (!out->IsOpen()) {
283 cerr<<"AliTRDtracker::Clusters2Tracks(): output file is not open !\n";
284 return 1;
285 }
46d29e70 286
5443e65e 287 sprintf(tname,"seedTRDtoTPC_%d",fEvent);
288 TTree tpc_tree(tname,"Tree with seeds from TRD at outer TPC pad row");
289 AliTPCtrack *iotrack=0;
290 tpc_tree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
291
292 sprintf(tname,"TreeT%d_TRD",fEvent);
293 TTree trd_tree(tname,"TRD tracks at inner TRD time bin");
294 AliTRDtrack *iotrack_trd=0;
295 trd_tree.Branch("tracks","AliTRDtrack",&iotrack_trd,32000,0);
296
297 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
298 Float_t foundMin = fMinClustersInTrack * timeBins;
299
300 if (inp) {
301 TFile *in=(TFile*)inp;
302 if (!in->IsOpen()) {
303 cerr<<"AliTRDtracker::Clusters2Tracks(): file with seeds is not open !\n";
304 cerr<<" ... going for seeds finding inside the TRD\n";
305 }
306 else {
307 in->cd();
308 sprintf(tname,"TRDb_%d",fEvent);
309 TTree *seedTree=(TTree*)in->Get(tname);
310 if (!seedTree) {
311 cerr<<"AliTRDtracker::Clusters2Tracks(): ";
312 cerr<<"can't get a tree with track seeds !\n";
313 return 3;
314 }
315 AliTRDtrack *seed=new AliTRDtrack;
316 seedTree->SetBranchAddress("tracks",&seed);
317
318 Int_t n=(Int_t)seedTree->GetEntries();
319 for (Int_t i=0; i<n; i++) {
320 seedTree->GetEvent(i);
321 seed->ResetCovariance();
fd621f36 322 AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
323 fSeeds->AddLast(tr);
5443e65e 324 fNseeds++;
325 }
326 delete seed;
fd621f36 327 delete seedTree;
5443e65e 328 }
329 }
46d29e70 330
5443e65e 331 out->cd();
46d29e70 332
fd621f36 333
334
5443e65e 335 // find tracks from loaded seeds
a819a5f7 336
5443e65e 337 Int_t nseed=fSeeds->GetEntriesFast();
338 Int_t i, found = 0;
339 Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
340
341 for (i=0; i<nseed; i++) {
342 AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
343 FollowProlongation(t, innerTB);
344 if (t.GetNumberOfClusters() >= foundMin) {
345 UseClusters(&t);
346 CookLabel(pt, 1-fLabelFraction);
347 // t.CookdEdx();
348 }
349 iotrack_trd = pt;
350 trd_tree.Fill();
fd621f36 351 cout<<found++<<'\r';
5443e65e 352
353 if(PropagateToTPC(t)) {
354 AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
355 iotrack = tpc;
356 tpc_tree.Fill();
357 delete tpc;
358 }
359 delete fSeeds->RemoveAt(i);
360 fNseeds--;
361 }
a819a5f7 362
5443e65e 363 cout<<"Number of loaded seeds: "<<nseed<<endl;
364 cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
a819a5f7 365
5443e65e 366 // after tracks from loaded seeds are found and the corresponding
367 // clusters are used, look for additional seeds from TRD
46d29e70 368
5443e65e 369 if(fAddTRDseeds) {
370 // Find tracks for the seeds in the TRD
371 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
372
373 Int_t nSteps = (Int_t) (fSeedDepth / fSeedStep);
374 Int_t gap = (Int_t) (timeBins * fSeedGap);
375 Int_t step = (Int_t) (timeBins * fSeedStep);
376
377 // make a first turn with tight cut on initial curvature
378 for(Int_t turn = 1; turn <= 2; turn++) {
379 if(turn == 2) {
380 nSteps = (Int_t) (fSeedDepth / (3*fSeedStep));
381 step = (Int_t) (timeBins * (3*fSeedStep));
382 }
383 for(Int_t i=0; i<nSteps; i++) {
384 Int_t outer=timeBins-1-i*step;
385 Int_t inner=outer-gap;
46d29e70 386
5443e65e 387 nseed=fSeeds->GetEntriesFast();
5443e65e 388
389 MakeSeeds(inner, outer, turn);
390
391 nseed=fSeeds->GetEntriesFast();
fd621f36 392 printf("\n turn %d, step %d: number of seeds for TRD inward %d\n",
393 turn, i, nseed);
394
5443e65e 395 for (Int_t i=0; i<nseed; i++) {
396 AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
397 FollowProlongation(t,innerTB);
398 if (t.GetNumberOfClusters() >= foundMin) {
399 UseClusters(&t);
400 CookLabel(pt, 1-fLabelFraction);
401 t.CookdEdx();
fd621f36 402 cout<<found++<<'\r';
5443e65e 403 iotrack_trd = pt;
404 trd_tree.Fill();
405 if(PropagateToTPC(t)) {
406 AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
407 iotrack = tpc;
408 tpc_tree.Fill();
409 delete tpc;
410 }
411 }
412 delete fSeeds->RemoveAt(i);
413 fNseeds--;
414 }
46d29e70 415 }
5443e65e 416 }
417 }
418 tpc_tree.Write();
419 trd_tree.Write();
420
421 cout<<"Total number of found tracks: "<<found<<endl;
422
423 UnloadEvent();
424
425 savedir->cd();
426
427 return 0;
428}
429
430
431
432//_____________________________________________________________________________
433Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) {
434 //
435 // Reads seeds from file <inp>. The seeds are AliTPCtrack's found and
436 // backpropagated by the TPC tracker. Each seed is first propagated
437 // to the TRD, and then its prolongation is searched in the TRD.
438 // If sufficiently long continuation of the track is found in the TRD
439 // the track is updated, otherwise it's stored as originaly defined
440 // by the TPC tracker.
441 //
442
443
444 LoadEvent();
445
446 TDirectory *savedir=gDirectory;
447
448 TFile *in=(TFile*)inp;
449
450 if (!in->IsOpen()) {
451 cerr<<"AliTRDtracker::PropagateBack(): ";
452 cerr<<"file with back propagated TPC tracks is not open !\n";
453 return 1;
454 }
455
456 if (!out->IsOpen()) {
457 cerr<<"AliTRDtracker::PropagateBack(): ";
458 cerr<<"file for back propagated TRD tracks is not open !\n";
459 return 2;
460 }
461
462 in->cd();
463 char tname[100];
464 sprintf(tname,"seedsTPCtoTRD_%d",fEvent);
465 TTree *seedTree=(TTree*)in->Get(tname);
466 if (!seedTree) {
467 cerr<<"AliTRDtracker::PropagateBack(): ";
468 cerr<<"can't get a tree with seeds from TPC !\n";
469 cerr<<"check if your version of TPC tracker creates tree "<<tname<<"\n";
470 return 3;
471 }
46d29e70 472
5443e65e 473 AliTPCtrack *seed=new AliTPCtrack;
474 seedTree->SetBranchAddress("tracks",&seed);
475
476 Int_t n=(Int_t)seedTree->GetEntries();
477 for (Int_t i=0; i<n; i++) {
478 seedTree->GetEvent(i);
479 Int_t lbl = seed->GetLabel();
480 AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
481 tr->SetSeedLabel(lbl);
482 fSeeds->AddLast(tr);
483 fNseeds++;
484 }
a819a5f7 485
5443e65e 486 delete seed;
487 delete seedTree;
a819a5f7 488
5443e65e 489 out->cd();
a819a5f7 490
5443e65e 491 AliTPCtrack *otrack=0;
a819a5f7 492
5443e65e 493 sprintf(tname,"seedsTRDtoTOF1_%d",fEvent);
494 TTree tofTree1(tname,"Tracks back propagated through TPC and TRD");
495 tofTree1.Branch("tracks","AliTPCtrack",&otrack,32000,0);
a819a5f7 496
5443e65e 497 sprintf(tname,"seedsTRDtoTOF2_%d",fEvent);
498 TTree tofTree2(tname,"Tracks back propagated through TPC and TRD");
499 tofTree2.Branch("tracks","AliTPCtrack",&otrack,32000,0);
46d29e70 500
5443e65e 501 sprintf(tname,"seedsTRDtoPHOS_%d",fEvent);
502 TTree phosTree(tname,"Tracks back propagated through TPC and TRD");
503 phosTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
a819a5f7 504
5443e65e 505 sprintf(tname,"seedsTRDtoRICH_%d",fEvent);
506 TTree richTree(tname,"Tracks back propagated through TPC and TRD");
507 richTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
a819a5f7 508
5443e65e 509 sprintf(tname,"TRDb_%d",fEvent);
510 TTree trdTree(tname,"Back propagated TRD tracks at outer TRD time bin");
511 AliTRDtrack *otrack_trd=0;
512 trdTree.Branch("tracks","AliTRDtrack",&otrack_trd,32000,0);
513
514 Int_t found=0;
46d29e70 515
5443e65e 516 Int_t nseed=fSeeds->GetEntriesFast();
46d29e70 517
5443e65e 518 Float_t foundMin = fMinClustersInTrack * fTimeBinsPerPlane * fGeom->Nplan();
a819a5f7 519
5443e65e 520 Int_t outermost_tb = fTrSec[0]->GetOuterTimeBin();
a819a5f7 521
5443e65e 522 for (Int_t i=0; i<nseed; i++) {
a819a5f7 523
5443e65e 524 AliTRDtrack *ps=(AliTRDtrack*)fSeeds->UncheckedAt(i), &s=*ps;
525 Int_t expectedClr = FollowBackProlongation(s);
526 Int_t foundClr = s.GetNumberOfClusters();
527 Int_t last_tb = fTrSec[0]->GetLayerNumber(s.GetX());
a819a5f7 528
fd621f36 529 // printf("seed %d: found %d out of %d expected clusters, Min is %f\n",
530 // i, foundClr, expectedClr, foundMin);
a819a5f7 531
5443e65e 532 if (foundClr >= foundMin) {
533 s.CookdEdx();
534 CookLabel(ps, 1-fLabelFraction);
535 UseClusters(ps);
536 otrack_trd=ps;
537 trdTree.Fill();
538 cout<<found++<<'\r';
46d29e70 539 }
a819a5f7 540
5443e65e 541 if(((expectedClr < 10) && (last_tb == outermost_tb)) ||
542 ((expectedClr >= 10) &&
543 (((Float_t) foundClr) / ((Float_t) expectedClr) >=
544 fMinFractionOfFoundClusters) && (last_tb == outermost_tb))) {
545
546 Double_t x_tof = 375.5;
547
548 if(PropagateToOuterPlane(s,x_tof)) {
549 AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
550 otrack = pt;
551 tofTree1.Fill();
552 delete pt;
553
554 x_tof = 381.5;
555
556 if(PropagateToOuterPlane(s,x_tof)) {
557 AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
558 otrack = pt;
559 tofTree2.Fill();
560 delete pt;
561
562 Double_t x_phos = 460.;
563
564 if(PropagateToOuterPlane(s,x_phos)) {
565 AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
566 otrack = pt;
567 phosTree.Fill();
568 delete pt;
569
570 Double_t x_rich = 490+1.267;
571
572 if(PropagateToOuterPlane(s,x_rich)) {
573 AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
574 otrack = pt;
575 richTree.Fill();
576 delete pt;
577 }
578 }
579 }
580 }
46d29e70 581 }
582 }
5443e65e 583
584 tofTree1.Write();
585 tofTree2.Write();
586 phosTree.Write();
587 richTree.Write();
588 trdTree.Write();
46d29e70 589
5443e65e 590 savedir->cd();
591 cerr<<"Number of seeds: "<<nseed<<endl;
592 cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
46d29e70 593
5443e65e 594 UnloadEvent();
46d29e70 595
5443e65e 596 return 0;
46d29e70 597
5443e65e 598}
46d29e70 599
bbf92647 600
5443e65e 601//---------------------------------------------------------------------------
602Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
603{
604 // Starting from current position on track=t this function tries
605 // to extrapolate the track up to timeBin=0 and to confirm prolongation
606 // if a close cluster is found. Returns the number of clusters
607 // expected to be found in sensitive layers
bbf92647 608
5443e65e 609 Float_t wIndex, wTB, wChi2;
610 Float_t wYrt, wYclosest, wYcorrect, wYwindow;
611 Float_t wZrt, wZclosest, wZcorrect, wZwindow;
612 Float_t wPx, wPy, wPz, wC;
613 Double_t Px, Py, Pz;
614 Float_t wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
46d29e70 615
5443e65e 616 Int_t trackIndex = t.GetLabel();
46d29e70 617
5443e65e 618 Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
46d29e70 619
5443e65e 620 Int_t try_again=fMaxGap;
46d29e70 621
5443e65e 622 Double_t alpha=t.GetAlpha();
46d29e70 623
5443e65e 624 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
625 if (alpha < 0. ) alpha += 2.*TMath::Pi();
46d29e70 626
5443e65e 627 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
46d29e70 628
5443e65e 629 Double_t x0, rho, x, dx, y, ymax, z;
46d29e70 630
5443e65e 631 Int_t expectedNumberOfClusters = 0;
632 Bool_t lookForCluster;
46d29e70 633
5443e65e 634 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
46d29e70 635
5443e65e 636
637 for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) {
638
639 y = t.GetY(); z = t.GetZ();
640
641 // first propagate to the inner surface of the current time bin
642 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
643 x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
644 if(!t.PropagateTo(x,x0,rho,0.139)) break;
645 y = t.GetY();
646 ymax = x*TMath::Tan(0.5*alpha);
647 if (y > ymax) {
648 s = (s+1) % ns;
649 if (!t.Rotate(alpha)) break;
650 if(!t.PropagateTo(x,x0,rho,0.139)) break;
651 } else if (y <-ymax) {
652 s = (s-1+ns) % ns;
653 if (!t.Rotate(-alpha)) break;
654 if(!t.PropagateTo(x,x0,rho,0.139)) break;
655 }
46d29e70 656
5443e65e 657 y = t.GetY(); z = t.GetZ();
658
659 // now propagate to the middle plane of the next time bin
660 fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
661 x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
662 if(!t.PropagateTo(x,x0,rho,0.139)) break;
663 y = t.GetY();
664 ymax = x*TMath::Tan(0.5*alpha);
665 if (y > ymax) {
666 s = (s+1) % ns;
667 if (!t.Rotate(alpha)) break;
668 if(!t.PropagateTo(x,x0,rho,0.139)) break;
669 } else if (y <-ymax) {
670 s = (s-1+ns) % ns;
671 if (!t.Rotate(-alpha)) break;
672 if(!t.PropagateTo(x,x0,rho,0.139)) break;
673 }
46d29e70 674
46d29e70 675
5443e65e 676 if(lookForCluster) {
a819a5f7 677
46d29e70 678
5443e65e 679 expectedNumberOfClusters++;
46d29e70 680
5443e65e 681 wIndex = (Float_t) t.GetLabel();
682 wTB = nr;
46d29e70 683
46d29e70 684
5443e65e 685 AliTRDpropagationLayer& time_bin=*(fTrSec[s]->GetLayer(nr-1));
46d29e70 686
46d29e70 687
5443e65e 688 Double_t sy2=ExpectedSigmaY2(x,t.GetTgl(),t.GetPt());
46d29e70 689
46d29e70 690
5443e65e 691 Double_t sz2=ExpectedSigmaZ2(x,t.GetTgl());
46d29e70 692
5443e65e 693 Double_t road;
694 if((t.GetSigmaY2() + sy2) > 0) road=10.*sqrt(t.GetSigmaY2() + sy2);
695 else return expectedNumberOfClusters;
696
697 wYrt = (Float_t) y;
698 wZrt = (Float_t) z;
699 wYwindow = (Float_t) road;
700 t.GetPxPyPz(Px,Py,Pz);
701 wPx = (Float_t) Px;
702 wPy = (Float_t) Py;
703 wPz = (Float_t) Pz;
704 wC = (Float_t) t.GetC();
705 wSigmaC2 = (Float_t) t.GetSigmaC2();
706 wSigmaTgl2 = (Float_t) t.GetSigmaTgl2();
707 wSigmaY2 = (Float_t) t.GetSigmaY2();
708 wSigmaZ2 = (Float_t) t.GetSigmaZ2();
709 wChi2 = -1;
710
711 if (road>fWideRoad) {
712 if (t.GetNumberOfClusters()>4)
713 cerr<<t.GetNumberOfClusters()
714 <<"FindProlongation warning: Too broad road !\n";
715 return 0;
716 }
717
718 AliTRDcluster *cl=0;
719 UInt_t index=0;
720
721 Double_t max_chi2=fMaxChi2;
722
723 wYclosest = 12345678;
724 wYcorrect = 12345678;
725 wZclosest = 12345678;
726 wZcorrect = 12345678;
727 wZwindow = TMath::Sqrt(2.25 * 12 * sz2);
728
729 // Find the closest correct cluster for debugging purposes
730 if (time_bin) {
731 Float_t minDY = 1000000;
732 for (Int_t i=0; i<time_bin; i++) {
733 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
734 if((c->GetLabel(0) != trackIndex) &&
735 (c->GetLabel(1) != trackIndex) &&
736 (c->GetLabel(2) != trackIndex)) continue;
737 if(TMath::Abs(c->GetY() - y) > minDY) continue;
738 minDY = TMath::Abs(c->GetY() - y);
739 wYcorrect = c->GetY();
740 wZcorrect = c->GetZ();
fd621f36 741
742 Double_t h01 = GetTiltFactor(c);
743 wChi2 = t.GetPredictedChi2(c, h01);
5443e65e 744 }
745 }
46d29e70 746
5443e65e 747 // Now go for the real cluster search
a819a5f7 748
5443e65e 749 if (time_bin) {
a819a5f7 750
5443e65e 751 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
752 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
753 if (c->GetY() > y+road) break;
754 if (c->IsUsed() > 0) continue;
755 if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
fd621f36 756
757 Double_t h01 = GetTiltFactor(c);
758 Double_t chi2=t.GetPredictedChi2(c,h01);
5443e65e 759
760 if (chi2 > max_chi2) continue;
761 max_chi2=chi2;
762 cl=c;
763 index=time_bin.GetIndex(i);
764 }
765
766 if(!cl) {
767
768 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
769 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
770
771 if (c->GetY() > y+road) break;
772 if (c->IsUsed() > 0) continue;
773 if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
774
fd621f36 775 Double_t h01 = GetTiltFactor(c);
776 Double_t chi2=t.GetPredictedChi2(c, h01);
5443e65e 777
778 if (chi2 > max_chi2) continue;
779 max_chi2=chi2;
780 cl=c;
781 index=time_bin.GetIndex(i);
782 }
783 }
784
a819a5f7 785
5443e65e 786 if (cl) {
787 wYclosest = cl->GetY();
788 wZclosest = cl->GetZ();
fd621f36 789 Double_t h01 = GetTiltFactor(cl);
a819a5f7 790
5443e65e 791 t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
fd621f36 792 if(!t.Update(cl,max_chi2,index,h01)) {
5443e65e 793 if(!try_again--) return 0;
794 }
795 else try_again=fMaxGap;
796 }
797 else {
798 if (try_again==0) break;
799 try_again--;
800 }
46d29e70 801
5443e65e 802 /*
803 if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
804
805 printf(" %f", wIndex); //1
806 printf(" %f", wTB); //2
807 printf(" %f", wYrt); //3
808 printf(" %f", wYclosest); //4
809 printf(" %f", wYcorrect); //5
810 printf(" %f", wYwindow); //6
811 printf(" %f", wZrt); //7
812 printf(" %f", wZclosest); //8
813 printf(" %f", wZcorrect); //9
814 printf(" %f", wZwindow); //10
815 printf(" %f", wPx); //11
816 printf(" %f", wPy); //12
817 printf(" %f", wPz); //13
818 printf(" %f", wSigmaC2*1000000); //14
819 printf(" %f", wSigmaTgl2*1000); //15
820 printf(" %f", wSigmaY2); //16
821 // printf(" %f", wSigmaZ2); //17
822 printf(" %f", wChi2); //17
823 printf(" %f", wC); //18
824 printf("\n");
825 }
826 */
827 }
828 }
829 }
830 return expectedNumberOfClusters;
831
832
833}
a819a5f7 834
5443e65e 835//___________________________________________________________________
46d29e70 836
5443e65e 837Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
838{
839 // Starting from current radial position of track <t> this function
840 // extrapolates the track up to outer timebin and in the sensitive
841 // layers confirms prolongation if a close cluster is found.
842 // Returns the number of clusters expected to be found in sensitive layers
46d29e70 843
5443e65e 844 Float_t wIndex, wTB, wChi2;
845 Float_t wYrt, wYclosest, wYcorrect, wYwindow;
846 Float_t wZrt, wZclosest, wZcorrect, wZwindow;
847 Float_t wPx, wPy, wPz, wC;
848 Double_t Px, Py, Pz;
849 Float_t wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
46d29e70 850
5443e65e 851 Int_t trackIndex = t.GetLabel();
46d29e70 852
5443e65e 853 Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
a819a5f7 854
5443e65e 855 Int_t try_again=fMaxGap;
46d29e70 856
5443e65e 857 Double_t alpha=t.GetAlpha();
46d29e70 858
5443e65e 859 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
860 if (alpha < 0. ) alpha += 2.*TMath::Pi();
46d29e70 861
5443e65e 862 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
46d29e70 863
5443e65e 864 Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
46d29e70 865
5443e65e 866 Double_t x0, rho, x, dx, y, ymax, z;
a819a5f7 867
5443e65e 868 Bool_t lookForCluster;
a819a5f7 869
5443e65e 870 Int_t expectedNumberOfClusters = 0;
871 x = t.GetX();
46d29e70 872
5443e65e 873 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
46d29e70 874
46d29e70 875
5443e65e 876 for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB; nr++) {
46d29e70 877
5443e65e 878 y = t.GetY(); z = t.GetZ();
46d29e70 879
5443e65e 880 // first propagate to the outer surface of the current time bin
46d29e70 881
5443e65e 882 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
883 x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
46d29e70 884
5443e65e 885 if(!t.PropagateTo(x,x0,rho,0.139)) break;
46d29e70 886
5443e65e 887 y = t.GetY();
888 ymax = x*TMath::Tan(0.5*alpha);
a819a5f7 889
5443e65e 890 if (y > ymax) {
891 s = (s+1) % ns;
892 if (!t.Rotate(alpha)) break;
893 if(!t.PropagateTo(x,x0,rho,0.139)) break;
894 } else if (y <-ymax) {
895 s = (s-1+ns) % ns;
896 if (!t.Rotate(-alpha)) break;
897 if(!t.PropagateTo(x,x0,rho,0.139)) break;
898 }
899 y = t.GetY(); z = t.GetZ();
46d29e70 900
5443e65e 901 // now propagate to the middle plane of the next time bin
902 fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
46d29e70 903
5443e65e 904 x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
46d29e70 905
5443e65e 906 if(!t.PropagateTo(x,x0,rho,0.139)) break;
46d29e70 907
5443e65e 908 y = t.GetY();
a819a5f7 909
5443e65e 910 ymax = x*TMath::Tan(0.5*alpha);
a819a5f7 911
5443e65e 912 if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax);
913
914 if (y > ymax) {
915 s = (s+1) % ns;
916 if (!t.Rotate(alpha)) break;
917 if(!t.PropagateTo(x,x0,rho,0.139)) break;
918 } else if (y <-ymax) {
919 s = (s-1+ns) % ns;
920 if (!t.Rotate(-alpha)) break;
921 if(!t.PropagateTo(x,x0,rho,0.139)) break;
922 }
46d29e70 923
5443e65e 924 // printf("label %d, pl %d, lookForCluster %d \n",
925 // trackIndex, nr+1, lookForCluster);
926
927 if(lookForCluster) {
928 expectedNumberOfClusters++;
929
930 wIndex = (Float_t) t.GetLabel();
931 wTB = fTrSec[s]->GetLayer(nr+1)->GetTimeBinIndex();
932
933 AliTRDpropagationLayer& time_bin=*(fTrSec[s]->GetLayer(nr+1));
934 Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
935 Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
936 if((t.GetSigmaY2() + sy2) < 0) break;
937 Double_t road = 10.*sqrt(t.GetSigmaY2() + sy2);
938 Double_t y=t.GetY(), z=t.GetZ();
939
940 wYrt = (Float_t) y;
941 wZrt = (Float_t) z;
942 wYwindow = (Float_t) road;
943 t.GetPxPyPz(Px,Py,Pz);
944 wPx = (Float_t) Px;
945 wPy = (Float_t) Py;
946 wPz = (Float_t) Pz;
947 wC = (Float_t) t.GetC();
948 wSigmaC2 = (Float_t) t.GetSigmaC2();
949 wSigmaTgl2 = (Float_t) t.GetSigmaTgl2();
950 wSigmaY2 = (Float_t) t.GetSigmaY2();
951 wSigmaZ2 = (Float_t) t.GetSigmaZ2();
952 wChi2 = -1;
953
954 if (road>fWideRoad) {
955 if (t.GetNumberOfClusters()>4)
956 cerr<<t.GetNumberOfClusters()
957 <<"FindProlongation warning: Too broad road !\n";
958 return 0;
959 }
960
961 AliTRDcluster *cl=0;
962 UInt_t index=0;
963
964 Double_t max_chi2=fMaxChi2;
965
966 wYclosest = 12345678;
967 wYcorrect = 12345678;
968 wZclosest = 12345678;
969 wZcorrect = 12345678;
970 wZwindow = TMath::Sqrt(2.25 * 12 * sz2);
971
972 // Find the closest correct cluster for debugging purposes
973 if (time_bin) {
974 Float_t minDY = 1000000;
975 for (Int_t i=0; i<time_bin; i++) {
976 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
977 if((c->GetLabel(0) != trackIndex) &&
978 (c->GetLabel(1) != trackIndex) &&
979 (c->GetLabel(2) != trackIndex)) continue;
980 if(TMath::Abs(c->GetY() - y) > minDY) continue;
981 minDY = TMath::Abs(c->GetY() - y);
982 wYcorrect = c->GetY();
983 wZcorrect = c->GetZ();
fd621f36 984
985 Double_t h01 = GetTiltFactor(c);
986 wChi2 = t.GetPredictedChi2(c, h01);
46d29e70 987 }
5443e65e 988 }
46d29e70 989
5443e65e 990 // Now go for the real cluster search
46d29e70 991
5443e65e 992 if (time_bin) {
993
994 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
995 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
996 if (c->GetY() > y+road) break;
997 if (c->IsUsed() > 0) continue;
998 if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
fd621f36 999
1000 Double_t h01 = GetTiltFactor(c);
1001 Double_t chi2=t.GetPredictedChi2(c,h01);
5443e65e 1002
1003 if (chi2 > max_chi2) continue;
1004 max_chi2=chi2;
1005 cl=c;
1006 index=time_bin.GetIndex(i);
1007 }
1008
1009 if(!cl) {
1010
1011 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
1012 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
1013
1014 if (c->GetY() > y+road) break;
1015 if (c->IsUsed() > 0) continue;
1016 if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
1017
fd621f36 1018 Double_t h01 = GetTiltFactor(c);
1019 Double_t chi2=t.GetPredictedChi2(c,h01);
5443e65e 1020
1021 if (chi2 > max_chi2) continue;
1022 max_chi2=chi2;
1023 cl=c;
1024 index=time_bin.GetIndex(i);
1025 }
1026 }
1027
1028 if (cl) {
1029 wYclosest = cl->GetY();
1030 wZclosest = cl->GetZ();
1031
1032 t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
fd621f36 1033 Double_t h01 = GetTiltFactor(cl);
1034 if(!t.Update(cl,max_chi2,index,h01)) {
5443e65e 1035 if(!try_again--) return 0;
1036 }
1037 else try_again=fMaxGap;
1038 }
1039 else {
1040 if (try_again==0) break;
1041 try_again--;
1042 }
1043
1044 /*
1045 if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
1046
1047 printf(" %f", wIndex); //1
1048 printf(" %f", wTB); //2
1049 printf(" %f", wYrt); //3
1050 printf(" %f", wYclosest); //4
1051 printf(" %f", wYcorrect); //5
1052 printf(" %f", wYwindow); //6
1053 printf(" %f", wZrt); //7
1054 printf(" %f", wZclosest); //8
1055 printf(" %f", wZcorrect); //9
1056 printf(" %f", wZwindow); //10
1057 printf(" %f", wPx); //11
1058 printf(" %f", wPy); //12
1059 printf(" %f", wPz); //13
1060 printf(" %f", wSigmaC2*1000000); //14
1061 printf(" %f", wSigmaTgl2*1000); //15
1062 printf(" %f", wSigmaY2); //16
1063 // printf(" %f", wSigmaZ2); //17
1064 printf(" %f", wChi2); //17
1065 printf(" %f", wC); //18
1066 printf("\n");
1067 }
1068 */
1069 }
1070 }
1071 }
1072 return expectedNumberOfClusters;
1073}
1074
1075//___________________________________________________________________
1076
1077Int_t AliTRDtracker::PropagateToOuterPlane(AliTRDtrack& t, Double_t xToGo)
1078{
1079 // Starting from current radial position of track <t> this function
1080 // extrapolates the track up to radial position <xToGo>.
1081 // Returns 1 if track reaches the plane, and 0 otherwise
1082
1083 Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1084
1085 Double_t alpha=t.GetAlpha();
1086
1087 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1088 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1089
1090 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1091
1092 Bool_t lookForCluster;
1093 Double_t x0, rho, x, dx, y, ymax, z;
1094
1095 x = t.GetX();
1096
1097 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
1098
1099 Int_t plToGo = fTrSec[0]->GetLayerNumber(xToGo);
1100
1101 for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr<plToGo; nr++) {
1102
1103 y = t.GetY(); z = t.GetZ();
1104
1105 // first propagate to the outer surface of the current time bin
1106 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
1107 x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
1108 if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1109 y = t.GetY();
1110 ymax = x*TMath::Tan(0.5*alpha);
1111 if (y > ymax) {
1112 s = (s+1) % ns;
1113 if (!t.Rotate(alpha)) return 0;
1114 } else if (y <-ymax) {
1115 s = (s-1+ns) % ns;
1116 if (!t.Rotate(-alpha)) return 0;
1117 }
1118 if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1119
1120 y = t.GetY(); z = t.GetZ();
1121
1122 // now propagate to the middle plane of the next time bin
1123 fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
1124 x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
1125 if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1126 y = t.GetY();
1127 ymax = x*TMath::Tan(0.5*alpha);
1128 if (y > ymax) {
1129 s = (s+1) % ns;
1130 if (!t.Rotate(alpha)) return 0;
1131 } else if (y <-ymax) {
1132 s = (s-1+ns) % ns;
1133 if (!t.Rotate(-alpha)) return 0;
1134 }
1135 if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1136 }
1137 return 1;
1138}
1139
1140//___________________________________________________________________
1141
1142Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t)
1143{
1144 // Starting from current radial position of track <t> this function
1145 // extrapolates the track up to radial position of the outermost
1146 // padrow of the TPC.
1147 // Returns 1 if track reaches the TPC, and 0 otherwise
1148
1149 Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1150
1151 Double_t alpha=t.GetAlpha();
1152
1153 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1154 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1155
1156 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1157
1158 Bool_t lookForCluster;
1159 Double_t x0, rho, x, dx, y, ymax, z;
1160
1161 x = t.GetX();
1162
1163 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
1164
1165 Int_t plTPC = fTrSec[0]->GetLayerNumber(246.055);
1166
1167 for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr>plTPC; nr--) {
1168
1169 y = t.GetY(); z = t.GetZ();
1170
1171 // first propagate to the outer surface of the current time bin
1172 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
1173 x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
1174 if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1175 y = t.GetY();
1176 ymax = x*TMath::Tan(0.5*alpha);
1177 if (y > ymax) {
1178 s = (s+1) % ns;
1179 if (!t.Rotate(alpha)) return 0;
1180 } else if (y <-ymax) {
1181 s = (s-1+ns) % ns;
1182 if (!t.Rotate(-alpha)) return 0;
1183 }
1184 if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1185
1186 y = t.GetY(); z = t.GetZ();
1187
1188 // now propagate to the middle plane of the next time bin
1189 fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
1190 x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
1191 if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1192 y = t.GetY();
1193 ymax = x*TMath::Tan(0.5*alpha);
1194 if (y > ymax) {
1195 s = (s+1) % ns;
1196 if (!t.Rotate(alpha)) return 0;
1197 } else if (y <-ymax) {
1198 s = (s-1+ns) % ns;
1199 if (!t.Rotate(-alpha)) return 0;
1200 }
1201 if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1202 }
1203 return 1;
1204}
1205
1206
1207//_____________________________________________________________________________
1208void AliTRDtracker::LoadEvent()
1209{
1210 // Fills clusters into TRD tracking_sectors
1211 // Note that the numbering scheme for the TRD tracking_sectors
1212 // differs from that of TRD sectors
1213
1214 ReadClusters(fClusters);
1215 Int_t ncl=fClusters->GetEntriesFast();
1216 cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1217
1218 UInt_t index;
1219 while (ncl--) {
1220 printf("\r %d left ",ncl);
1221 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
1222 Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin();
1223 Int_t sector=fGeom->GetSector(detector);
1224 Int_t plane=fGeom->GetPlane(detector);
1225
1226 Int_t tracking_sector = CookSectorIndex(sector);
1227
1228 Int_t gtb = fTrSec[tracking_sector]->CookTimeBinIndex(plane,local_time_bin);
1229 Int_t layer = fTrSec[tracking_sector]->GetLayerNumber(gtb);
1230
1231 index=ncl;
1232 fTrSec[tracking_sector]->GetLayer(layer)->InsertCluster(c,index);
1233 }
1234 printf("\r\n");
1235
1236}
1237
1238//_____________________________________________________________________________
1239void AliTRDtracker::UnloadEvent()
1240{
1241 //
1242 // Clears the arrays of clusters and tracks. Resets sectors and timebins
1243 //
1244
1245 Int_t i, nentr;
1246
1247 nentr = fClusters->GetEntriesFast();
1248 for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i);
1249
1250 nentr = fSeeds->GetEntriesFast();
1251 for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i);
1252
1253 nentr = fTracks->GetEntriesFast();
1254 for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i);
1255
1256 Int_t nsec = AliTRDgeometry::kNsect;
1257
1258 for (i = 0; i < nsec; i++) {
1259 for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1260 fTrSec[i]->GetLayer(pl)->Clear();
1261 }
1262 }
1263
1264}
1265
1266//__________________________________________________________________________
1267void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn)
1268{
1269 // Creates track seeds using clusters in timeBins=i1,i2
1270
1271 if(turn > 2) {
1272 cerr<<"MakeSeeds: turn "<<turn<<" exceeds the limit of 2"<<endl;
1273 return;
1274 }
1275
1276 Double_t x[5], c[15];
1277 Int_t max_sec=AliTRDgeometry::kNsect;
1278
1279 Double_t alpha=AliTRDgeometry::GetAlpha();
1280 Double_t shift=AliTRDgeometry::GetAlpha()/2.;
1281 Double_t cs=cos(alpha), sn=sin(alpha);
1282 Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
1283
1284
1285 Int_t i2 = fTrSec[0]->GetLayerNumber(inner);
1286 Int_t i1 = fTrSec[0]->GetLayerNumber(outer);
1287
1288 Double_t x1 =fTrSec[0]->GetX(i1);
1289 Double_t xx2=fTrSec[0]->GetX(i2);
1290
1291 for (Int_t ns=0; ns<max_sec; ns++) {
1292
1293 Int_t nl2 = *(fTrSec[(ns-2+max_sec)%max_sec]->GetLayer(i2));
1294 Int_t nl=(*fTrSec[(ns-1+max_sec)%max_sec]->GetLayer(i2));
1295 Int_t nm=(*fTrSec[ns]->GetLayer(i2));
1296 Int_t nu=(*fTrSec[(ns+1)%max_sec]->GetLayer(i2));
1297 Int_t nu2=(*fTrSec[(ns+2)%max_sec]->GetLayer(i2));
1298
1299 AliTRDpropagationLayer& r1=*(fTrSec[ns]->GetLayer(i1));
1300
1301 for (Int_t is=0; is < r1; is++) {
1302 Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
1303
1304 for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
1305
1306 const AliTRDcluster *cl;
1307 Double_t x2, y2, z2;
1308 Double_t x3=0., y3=0.;
1309
1310 if (js<nl2) {
1311 if(turn != 2) continue;
1312 AliTRDpropagationLayer& r2=*(fTrSec[(ns-2+max_sec)%max_sec]->GetLayer(i2));
1313 cl=r2[js];
1314 y2=cl->GetY(); z2=cl->GetZ();
1315
1316 x2= xx2*cs2+y2*sn2;
1317 y2=-xx2*sn2+y2*cs2;
1318 }
1319 else if (js<nl2+nl) {
1320 if(turn != 1) continue;
1321 AliTRDpropagationLayer& r2=*(fTrSec[(ns-1+max_sec)%max_sec]->GetLayer(i2));
1322 cl=r2[js-nl2];
1323 y2=cl->GetY(); z2=cl->GetZ();
1324
1325 x2= xx2*cs+y2*sn;
1326 y2=-xx2*sn+y2*cs;
1327 }
1328 else if (js<nl2+nl+nm) {
1329 if(turn != 1) continue;
1330 AliTRDpropagationLayer& r2=*(fTrSec[ns]->GetLayer(i2));
1331 cl=r2[js-nl2-nl];
1332 x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
1333 }
1334 else if (js<nl2+nl+nm+nu) {
1335 if(turn != 1) continue;
1336 AliTRDpropagationLayer& r2=*(fTrSec[(ns+1)%max_sec]->GetLayer(i2));
1337 cl=r2[js-nl2-nl-nm];
1338 y2=cl->GetY(); z2=cl->GetZ();
1339
1340 x2=xx2*cs-y2*sn;
1341 y2=xx2*sn+y2*cs;
1342 }
1343 else {
1344 if(turn != 2) continue;
1345 AliTRDpropagationLayer& r2=*(fTrSec[(ns+2)%max_sec]->GetLayer(i2));
1346 cl=r2[js-nl2-nl-nm-nu];
1347 y2=cl->GetY(); z2=cl->GetZ();
1348
1349 x2=xx2*cs2-y2*sn2;
1350 y2=xx2*sn2+y2*cs2;
1351 }
1352
1353 if(TMath::Abs(z1-z2) > fMaxSeedDeltaZ12) continue;
1354
1355 Double_t zz=z1 - z1/x1*(x1-x2);
1356
1357 if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue;
1358
1359 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
1360 if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
1361
1362 x[0]=y1;
1363 x[1]=z1;
1364 x[2]=f1trd(x1,y1,x2,y2,x3,y3);
1365
1366 if (TMath::Abs(x[2]) > fMaxSeedC) continue;
1367
1368 x[3]=f2trd(x1,y1,x2,y2,x3,y3);
1369
1370 if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue;
1371
1372 x[4]=f3trd(x1,y1,x2,y2,z1,z2);
1373
1374 if (TMath::Abs(x[4]) > fMaxSeedTan) continue;
1375
1376 Double_t a=asin(x[3]);
1377 Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3]));
1378
1379 if (TMath::Abs(zv)>fMaxSeedVertexZ) continue;
1380
1381 Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
1382 Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2();
1383 Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ;
fd621f36 1384
1385 // Tilt changes
1386 Double_t h01 = GetTiltFactor(r1[is]);
1387 sy1=sy1+sz1*h01*h01;
1388 Double_t syz=sz1*h01;
1389 // end of tilt changes
5443e65e 1390
1391 Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1392 Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1393 Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1394 Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy;
1395 Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy;
1396 Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy;
1397 Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy;
1398 Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz;
1399 Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy;
1400 Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz;
1401
1402 c[0]=sy1;
fd621f36 1403 // c[1]=0.; c[2]=sz1;
1404 c[1]=syz; c[2]=sz1*100;
5443e65e 1405 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24;
1406 c[6]=f30*sy1; c[7]=0.; c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24;
1407 c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34;
1408 c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22;
1409 c[13]=f40*sy1*f30+f42*sy2*f32;
1410 c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43;
1411
1412 UInt_t index=r1.GetIndex(is);
1413
1414 AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
1415
1416 Int_t rc=FollowProlongation(*track, i2);
1417
1418 if ((rc < 1) ||
1419 (track->GetNumberOfClusters() <
1420 (outer-inner)*fMinClustersInSeed)) delete track;
1421 else {
1422 fSeeds->AddLast(track); fNseeds++;
1423 cerr<<"\r found seed "<<fNseeds;
1424 }
1425 }
1426 }
1427 }
1428}
1429
1430//_____________________________________________________________________________
1431void AliTRDtracker::ReadClusters(TObjArray *array, const TFile *inp)
1432{
1433 //
a819a5f7 1434 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
1435 // from the file. The names of the cluster tree and branches
1436 // should match the ones used in AliTRDclusterizer::WriteClusters()
1437 //
46d29e70 1438
a819a5f7 1439 TDirectory *savedir=gDirectory;
1440
5443e65e 1441 if (inp) {
1442 TFile *in=(TFile*)inp;
1443 if (!in->IsOpen()) {
1444 cerr<<"AliTRDtracker::ReadClusters(): input file is not open !\n";
1445 return;
1446 }
1447 else{
1448 in->cd();
1449 }
1450 }
a819a5f7 1451
abaf1f1d 1452 Char_t treeName[12];
1453 sprintf(treeName,"TreeR%d_TRD",fEvent);
5443e65e 1454 TTree *ClusterTree = (TTree*) gDirectory->Get(treeName);
1455
1456 TObjArray *ClusterArray = new TObjArray(400);
1457
1458 ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray);
1459
1460 Int_t nEntries = (Int_t) ClusterTree->GetEntries();
1461 printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
a819a5f7 1462
a819a5f7 1463 // Loop through all entries in the tree
1464 Int_t nbytes;
1465 AliTRDcluster *c = 0;
5443e65e 1466 printf("\n");
a819a5f7 1467
1468 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
1469
1470 // Import the tree
5443e65e 1471 nbytes += ClusterTree->GetEvent(iEntry);
1472
a819a5f7 1473 // Get the number of points in the detector
5443e65e 1474 Int_t nCluster = ClusterArray->GetEntriesFast();
1475 printf("\r Read %d clusters from entry %d", nCluster, iEntry);
1476
a819a5f7 1477 // Loop through all TRD digits
1478 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
5443e65e 1479 c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
a819a5f7 1480 AliTRDcluster *co = new AliTRDcluster(*c);
1481 co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
5443e65e 1482 Int_t ltb = co->GetLocalTimeBin();
1483 if(ltb != 0) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
1484
a819a5f7 1485 array->AddLast(co);
5443e65e 1486 delete ClusterArray->RemoveAt(iCluster);
a819a5f7 1487 }
1488 }
1489
5443e65e 1490 delete ClusterArray;
1491 savedir->cd();
1492
a819a5f7 1493}
1494
5443e65e 1495//______________________________________________________________________
1496void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename)
a819a5f7 1497{
1498 //
5443e65e 1499 // Reads AliTRDclusters from file <filename>. The names of the cluster
1500 // tree and branches should match the ones used in
1501 // AliTRDclusterizer::WriteClusters()
1502 // if <array> == 0, clusters are added into AliTRDtracker fCluster array
a819a5f7 1503 //
1504
5443e65e 1505 TDirectory *savedir=gDirectory;
46d29e70 1506
5443e65e 1507 TFile *file = TFile::Open(filename);
1508 if (!file->IsOpen()) {
1509 cerr<<"Can't open file with TRD clusters"<<endl;
1510 return;
1511 }
46d29e70 1512
5443e65e 1513 Char_t treeName[12];
1514 sprintf(treeName,"TreeR%d_TRD",fEvent);
1515 TTree *ClusterTree = (TTree*) gDirectory->Get(treeName);
46d29e70 1516
5443e65e 1517 if (!ClusterTree) {
1518 cerr<<"AliTRDtracker::ReadClusters(): ";
1519 cerr<<"can't get a tree with clusters !\n";
1520 return;
1521 }
46d29e70 1522
5443e65e 1523 TObjArray *ClusterArray = new TObjArray(400);
46d29e70 1524
5443e65e 1525 ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray);
46d29e70 1526
5443e65e 1527 Int_t nEntries = (Int_t) ClusterTree->GetEntries();
1528 cout<<"found "<<nEntries<<" in ClusterTree"<<endl;
a819a5f7 1529
5443e65e 1530 // Loop through all entries in the tree
1531 Int_t nbytes;
1532 AliTRDcluster *c = 0;
bbf92647 1533
5443e65e 1534 printf("\n");
bbf92647 1535
5443e65e 1536 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
46d29e70 1537
5443e65e 1538 // Import the tree
1539 nbytes += ClusterTree->GetEvent(iEntry);
0a29d0f1 1540
5443e65e 1541 // Get the number of points in the detector
1542 Int_t nCluster = ClusterArray->GetEntriesFast();
1543 printf("\r Read %d clusters from entry %d", nCluster, iEntry);
1544
1545 // Loop through all TRD digits
1546 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
1547 c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
1548 AliTRDcluster *co = new AliTRDcluster(*c);
1549 co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
1550 Int_t ltb = co->GetLocalTimeBin();
1551 if(ltb != 0) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
1552 array->AddLast(co);
1553 delete ClusterArray->RemoveAt(iCluster);
1554 }
46d29e70 1555 }
0a29d0f1 1556
5443e65e 1557 file->Close();
1558 delete ClusterArray;
1559 savedir->cd();
1560
1561}
1562
46d29e70 1563
1564//__________________________________________________________________
5443e65e 1565void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const {
46d29e70 1566
1567 Int_t label=123456789, index, i, j;
5443e65e 1568 Int_t ncl=pt->GetNumberOfClusters();
1569 const Int_t range = fTrSec[0]->GetOuterTimeBin()+1;
1570
1571 Bool_t label_added;
46d29e70 1572
5443e65e 1573 // Int_t s[range][2];
1574 Int_t **s = new Int_t* [range];
1575 for (i=0; i<range; i++) {
d1b06c24 1576 s[i] = new Int_t[2];
1577 }
5443e65e 1578 for (i=0; i<range; i++) {
46d29e70 1579 s[i][0]=-1;
1580 s[i][1]=0;
1581 }
1582
1583 Int_t t0,t1,t2;
1584 for (i=0; i<ncl; i++) {
5443e65e 1585 index=pt->GetClusterIndex(i);
46d29e70 1586 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
5443e65e 1587 t0=c->GetLabel(0);
1588 t1=c->GetLabel(1);
1589 t2=c->GetLabel(2);
46d29e70 1590 }
1591
1592 for (i=0; i<ncl; i++) {
5443e65e 1593 index=pt->GetClusterIndex(i);
46d29e70 1594 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
1595 for (Int_t k=0; k<3; k++) {
5443e65e 1596 label=c->GetLabel(k);
1597 label_added=kFALSE; j=0;
46d29e70 1598 if (label >= 0) {
5443e65e 1599 while ( (!label_added) && ( j < range ) ) {
46d29e70 1600 if (s[j][0]==label || s[j][1]==0) {
1601 s[j][0]=label;
1602 s[j][1]=s[j][1]+1;
5443e65e 1603 label_added=kTRUE;
46d29e70 1604 }
1605 j++;
1606 }
1607 }
1608 }
1609 }
1610
46d29e70 1611 Int_t max=0;
1612 label = -123456789;
1613
5443e65e 1614 for (i=0; i<range; i++) {
46d29e70 1615 if (s[i][1]>max) {
1616 max=s[i][1]; label=s[i][0];
1617 }
1618 }
5443e65e 1619
1620 for (i=0; i<range; i++) {
1621 delete []s[i];
1622 }
1623
d1b06c24 1624 delete []s;
5443e65e 1625
1626 if ((1.- Float_t(max)/ncl) > wrong) label=-label;
1627
1628 pt->SetLabel(label);
1629
46d29e70 1630}
1631
5443e65e 1632
1633//__________________________________________________________________
ca6c93d6 1634void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const {
5443e65e 1635 Int_t ncl=t->GetNumberOfClusters();
1636 for (Int_t i=from; i<ncl; i++) {
1637 Int_t index = t->GetClusterIndex(i);
1638 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
1639 c->Use();
1640 }
1641}
1642
1643
1644//_____________________________________________________________________
1645Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt)
1646{
1647 // Parametrised "expected" error of the cluster reconstruction in Y
1648
1649 Double_t s = 0.08 * 0.08;
1650 return s;
1651}
1652
1653//_____________________________________________________________________
1654Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl)
0a29d0f1 1655{
5443e65e 1656 // Parametrised "expected" error of the cluster reconstruction in Z
1657
1658 Double_t s = 6 * 6 /12.;
1659 return s;
1660}
1661
1662
1663//_____________________________________________________________________
1664Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t local_tb) const
1665{
1666 //
1667 // Returns radial position which corresponds to time bin <local_tb>
1668 // in tracking sector <sector> and plane <plane>
1669 //
1670
1671 Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, local_tb);
1672 Int_t pl = fTrSec[sector]->GetLayerNumber(index);
1673 return fTrSec[sector]->GetLayer(pl)->GetX();
1674
1675}
1676
1677
1678//_______________________________________________________
1679AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x,
1680 Double_t dx, Double_t rho, Double_t x0, Int_t tb_index)
1681{
0a29d0f1 1682 //
5443e65e 1683 // AliTRDpropagationLayer constructor
0a29d0f1 1684 //
46d29e70 1685
5443e65e 1686 fN = 0; fX = x; fdX = dx; fRho = rho; fX0 = x0;
1687 fClusters = NULL; fIndex = NULL; fTimeBinIndex = tb_index;
1688
46d29e70 1689
5443e65e 1690 for(Int_t i=0; i < (Int_t) kZONES; i++) {
1691 fZc[i]=0; fZmax[i] = 0;
a819a5f7 1692 }
5443e65e 1693
1694 fYmax = 0;
1695
1696 if(fTimeBinIndex >= 0) {
ca6c93d6 1697 fClusters = new AliTRDcluster*[kMAX_CLUSTER_PER_TIME_BIN];
1698 fIndex = new UInt_t[kMAX_CLUSTER_PER_TIME_BIN];
a819a5f7 1699 }
46d29e70 1700
5443e65e 1701 fHole = kFALSE;
1702 fHoleZc = 0;
1703 fHoleZmax = 0;
1704 fHoleYc = 0;
1705 fHoleYmax = 0;
1706 fHoleRho = 0;
1707 fHoleX0 = 0;
1708
1709}
1710
1711//_______________________________________________________
1712void AliTRDtracker::AliTRDpropagationLayer::SetHole(
1713 Double_t Zmax, Double_t Ymax, Double_t rho,
1714 Double_t x0, Double_t Yc, Double_t Zc)
1715{
1716 //
1717 // Sets hole in the layer
1718 //
1719
1720 fHole = kTRUE;
1721 fHoleZc = Zc;
1722 fHoleZmax = Zmax;
1723 fHoleYc = Yc;
1724 fHoleYmax = Ymax;
1725 fHoleRho = rho;
1726 fHoleX0 = x0;
1727}
1728
46d29e70 1729
5443e65e 1730//_______________________________________________________
1731AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs, AliTRDparameter* par)
1732{
1733 //
1734 // AliTRDtrackingSector Constructor
1735 //
1736
1737 fGeom = geo;
1738 fPar = par;
1739 fGeomSector = gs;
1740 fTzeroShift = 0.13;
1741 fN = 0;
1742
1743 for(UInt_t i=0; i < kMAX_TIME_BIN_INDEX; i++) fTimeBinIndex[i] = -1;
1744
1745
1746 AliTRDpropagationLayer* ppl;
1747
1748 Double_t x, xin, xout, dx, rho, x0;
1749 Int_t steps;
46d29e70 1750
5443e65e 1751 // set time bins in the gas of the TPC
46d29e70 1752
5443e65e 1753 xin = 246.055; xout = 254.055; steps = 20; dx = (xout-xin)/steps;
1754 rho = 0.9e-3; x0 = 28.94;
1755
1756 for(Int_t i=0; i<steps; i++) {
1757 x = xin + i*dx + dx/2;
1758 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
1759 InsertLayer(ppl);
46d29e70 1760 }
1761
5443e65e 1762 // set time bins in the outer field cage vessel
46d29e70 1763
5443e65e 1764 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar
1765 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1766 InsertLayer(ppl);
46d29e70 1767
5443e65e 1768 dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; x0 = 44.86; // prepreg
1769 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1770 InsertLayer(ppl);
1771
1772 dx = 2.; xin = xout; xout = xin + dx; rho = 1.45*0.02; x0 = 41.28; // Nomex
1773 steps = 5; dx = (xout - xin)/steps;
1774 for(Int_t i=0; i<steps; i++) {
1775 x = xin + i*dx + dx/2;
1776 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
1777 InsertLayer(ppl);
1778 }
1779
1780 dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; x0 = 44.86; // prepreg
1781 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1782 InsertLayer(ppl);
1783
1784 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar
1785 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1786 InsertLayer(ppl);
0a29d0f1 1787
5443e65e 1788
1789 // set time bins in CO2
1790
1791 xin = xout; xout = 275.0;
1792 steps = 50; dx = (xout - xin)/steps;
1793 rho = 1.977e-3; x0 = 36.2;
1794
1795 for(Int_t i=0; i<steps; i++) {
1796 x = xin + i*dx + dx/2;
1797 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
1798 InsertLayer(ppl);
1799 }
1800
1801 // set time bins in the outer containment vessel
1802
1803 dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; x0 = 24.01; // Al
1804 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1805 InsertLayer(ppl);
1806
1807 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar
1808 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1809 InsertLayer(ppl);
1810
1811 dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; x0 = 44.86; // prepreg
1812 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1813 InsertLayer(ppl);
1814
1815 dx = 3.; xin = xout; xout = xin + dx; rho = 1.45*0.02; x0 = 41.28; // Nomex
1816 steps = 10; dx = (xout - xin)/steps;
1817 for(Int_t i=0; i<steps; i++) {
1818 x = xin + i*dx + dx/2;
1819 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
1820 InsertLayer(ppl);
1821 }
1822
1823 dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; x0 = 44.86; // prepreg
1824 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1825 InsertLayer(ppl);
1826
1827 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar
1828 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1829 InsertLayer(ppl);
1830
1831 dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; x0 = 24.01; // Al
1832 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1833 InsertLayer(ppl);
1834
1835 Double_t xtrd = (Double_t) fGeom->Rmin();
1836
1837 // add layers between TPC and TRD (Air temporarily)
1838 xin = xout; xout = xtrd;
1839 steps = 50; dx = (xout - xin)/steps;
1840 rho = 1.2e-3; x0 = 36.66;
1841
1842 for(Int_t i=0; i<steps; i++) {
1843 x = xin + i*dx + dx/2;
1844 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
1845 InsertLayer(ppl);
1846 }
1847
1848
1849 Double_t alpha=AliTRDgeometry::GetAlpha();
1850
1851 // add layers for each of the planes
1852
1853 Double_t dxRo = (Double_t) fGeom->CroHght(); // Rohacell
1854 Double_t dxSpace = (Double_t) fGeom->Cspace(); // Spacing between planes
1855 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
1856 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
1857 Double_t dxRad = (Double_t) fGeom->CraHght(); // Radiator
1858 Double_t dxTEC = dxRad + dxDrift + dxAmp + dxRo;
1859 Double_t dxPlane = dxTEC + dxSpace;
1860
1861 Int_t tbBefore = (Int_t) (dxAmp/fPar->GetTimeBinSize());
1862 if(tbBefore > fPar->GetTimeBefore()) tbBefore = fPar->GetTimeBefore();
1863
1864 Int_t tb, tb_index;
1865 const Int_t nChambers = AliTRDgeometry::Ncham();
ca6c93d6 1866 Double_t Ymax = 0, holeYmax = 0;
1867 Double_t * Zc = new Double_t[nChambers];
1868 Double_t * Zmax = new Double_t[nChambers];
5443e65e 1869 Double_t holeZmax = 1000.; // the whole sector is missing
1870
1871
1872 for(Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
1873
1874 // Radiator
1875 xin = xtrd + plane * dxPlane; xout = xin + dxRad;
1876 steps = 12; dx = (xout - xin)/steps; rho = 0.074; x0 = 40.6;
1877 for(Int_t i=0; i<steps; i++) {
1878 x = xin + i*dx + dx/2;
1879 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
1880 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1881 holeYmax = x*TMath::Tan(0.5*alpha);
1882 ppl->SetHole(holeYmax, holeZmax);
1883 }
1884 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1885 holeYmax = x*TMath::Tan(0.5*alpha);
1886 ppl->SetHole(holeYmax, holeZmax);
1887 }
1888 InsertLayer(ppl);
1889 }
1890
1891 Ymax = fGeom->GetChamberWidth(plane)/2;
1892 for(Int_t ch = 0; ch < nChambers; ch++) {
1893 Zmax[ch] = fGeom->GetChamberLength(plane,ch)/2;
1894 Float_t pad = fPar->GetRowPadSize(plane,ch,0);
1895 Float_t row0 = fPar->GetRow0(plane,ch,0);
1896 Int_t nPads = fPar->GetRowMax(plane,ch,0);
1897 Zc[ch] = (pad * nPads)/2 + row0 - pad/2;
1898 }
1899
1900 dx = fPar->GetTimeBinSize();
1901 rho = 0.00295 * 0.85; x0 = 11.0;
1902
1903 Double_t x0 = (Double_t) fPar->GetTime0(plane);
1904 Double_t xbottom = x0 - dxDrift;
1905 Double_t xtop = x0 + dxAmp;
1906
1907 // Amplification region
1908
1909 steps = (Int_t) (dxAmp/dx);
1910
1911 for(tb = 0; tb < steps; tb++) {
1912 x = x0 + tb * dx + dx/2;
1913 tb_index = CookTimeBinIndex(plane, -tb-1);
1914 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index);
1915 ppl->SetYmax(Ymax);
1916 for(Int_t ch = 0; ch < nChambers; ch++) {
1917 ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
1918 }
1919 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1920 holeYmax = x*TMath::Tan(0.5*alpha);
1921 ppl->SetHole(holeYmax, holeZmax);
1922 }
1923 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1924 holeYmax = x*TMath::Tan(0.5*alpha);
1925 ppl->SetHole(holeYmax, holeZmax);
1926 }
1927 InsertLayer(ppl);
1928 }
1929 tb_index = CookTimeBinIndex(plane, -steps);
1930 x = (x + dx/2 + xtop)/2;
1931 dx = 2*(xtop-x);
1932 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index);
1933 ppl->SetYmax(Ymax);
1934 for(Int_t ch = 0; ch < nChambers; ch++) {
1935 ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
1936 }
1937 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1938 holeYmax = x*TMath::Tan(0.5*alpha);
1939 ppl->SetHole(holeYmax, holeZmax);
1940 }
1941 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1942 holeYmax = x*TMath::Tan(0.5*alpha);
1943 ppl->SetHole(holeYmax, holeZmax);
1944 }
1945 InsertLayer(ppl);
1946
1947 // Drift region
1948 dx = fPar->GetTimeBinSize();
1949 steps = (Int_t) (dxDrift/dx);
1950
1951 for(tb = 0; tb < steps; tb++) {
1952 x = x0 - tb * dx - dx/2;
1953 tb_index = CookTimeBinIndex(plane, tb);
1954
1955 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index);
1956 ppl->SetYmax(Ymax);
1957 for(Int_t ch = 0; ch < nChambers; ch++) {
1958 ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
1959 }
1960 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1961 holeYmax = x*TMath::Tan(0.5*alpha);
1962 ppl->SetHole(holeYmax, holeZmax);
1963 }
1964 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1965 holeYmax = x*TMath::Tan(0.5*alpha);
1966 ppl->SetHole(holeYmax, holeZmax);
1967 }
1968 InsertLayer(ppl);
1969 }
1970 tb_index = CookTimeBinIndex(plane, steps);
1971 x = (x - dx/2 + xbottom)/2;
1972 dx = 2*(x-xbottom);
1973 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index);
1974 ppl->SetYmax(Ymax);
1975 for(Int_t ch = 0; ch < nChambers; ch++) {
1976 ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
1977 }
1978 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1979 holeYmax = x*TMath::Tan(0.5*alpha);
1980 ppl->SetHole(holeYmax, holeZmax);
1981 }
1982 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1983 holeYmax = x*TMath::Tan(0.5*alpha);
1984 ppl->SetHole(holeYmax, holeZmax);
1985 }
1986 InsertLayer(ppl);
1987
1988 // Pad Plane
1989 xin = xtop; dx = 0.025; xout = xin + dx; rho = 1.7; x0 = 33.0;
1990 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1991 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1992 holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha);
1993 ppl->SetHole(holeYmax, holeZmax);
1994 }
1995 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1996 holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha);
1997 ppl->SetHole(holeYmax, holeZmax);
1998 }
1999 InsertLayer(ppl);
2000
2001 // Rohacell
2002 xin = xout; xout = xtrd + (plane + 1) * dxPlane - dxSpace;
2003 steps = 5; dx = (xout - xin)/steps; rho = 0.074; x0 = 40.6;
2004 for(Int_t i=0; i<steps; i++) {
2005 x = xin + i*dx + dx/2;
2006 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
2007 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2008 holeYmax = x*TMath::Tan(0.5*alpha);
2009 ppl->SetHole(holeYmax, holeZmax);
2010 }
2011 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2012 holeYmax = x*TMath::Tan(0.5*alpha);
2013 ppl->SetHole(holeYmax, holeZmax);
2014 }
2015 InsertLayer(ppl);
2016 }
2017
2018 // Space between the chambers, air
2019 xin = xout; xout = xtrd + (plane + 1) * dxPlane;
2020 steps = 5; dx = (xout - xin)/steps; rho = 1.29e-3; x0 = 36.66;
2021 for(Int_t i=0; i<steps; i++) {
2022 x = xin + i*dx + dx/2;
2023 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
2024 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2025 holeYmax = x*TMath::Tan(0.5*alpha);
2026 ppl->SetHole(holeYmax, holeZmax);
2027 }
2028 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2029 holeYmax = x*TMath::Tan(0.5*alpha);
2030 ppl->SetHole(holeYmax, holeZmax);
2031 }
2032 InsertLayer(ppl);
2033 }
2034 }
2035
2036 // Space between the TRD and RICH
2037 Double_t xRICH = 500.;
2038 xin = xout; xout = xRICH;
2039 steps = 200; dx = (xout - xin)/steps; rho = 1.29e-3; x0 = 36.66;
2040 for(Int_t i=0; i<steps; i++) {
2041 x = xin + i*dx + dx/2;
2042 ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
2043 InsertLayer(ppl);
2044 }
2045
2046 MapTimeBinLayers();
ca6c93d6 2047 delete [] Zc;
2048 delete [] Zmax;
5443e65e 2049
2050}
2051
2052//______________________________________________________
2053
2054Int_t AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t local_tb) const
2055{
2056 //
2057 // depending on the digitization parameters calculates "global"
2058 // time bin index for timebin <local_tb> in plane <plane>
2059 //
2060
2061 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
2062 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
2063 Double_t dx = (Double_t) fPar->GetTimeBinSize();
2064
2065 Int_t tbAmp = fPar->GetTimeBefore();
2066 Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
2067 Int_t tbDrift = fPar->GetTimeMax();
2068 Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
2069
2070 Int_t tb_per_plane = TMath::Min(tbAmp,maxAmp) + TMath::Min(tbDrift,maxDrift);
2071
2072 Int_t gtb = (plane+1) * tb_per_plane - local_tb - 1;
2073
2074 if((local_tb < 0) &&
2075 (TMath::Abs(local_tb) > TMath::Min(tbAmp,maxAmp))) return -1;
2076 if(local_tb >= TMath::Min(tbDrift,maxDrift)) return -1;
2077
2078 return gtb;
2079
2080
2081}
2082
2083//______________________________________________________
2084
2085void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers()
2086{
2087 //
2088 // For all sensitive time bins sets corresponding layer index
2089 // in the array fTimeBins
2090 //
2091
2092 Int_t index;
2093
2094 for(Int_t i = 0; i < fN; i++) {
2095 index = fLayers[i]->GetTimeBinIndex();
2096
2097 // printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX());
2098
2099 if(index < 0) continue;
2100 if(index >= (Int_t) kMAX_TIME_BIN_INDEX) {
2101 printf("*** AliTRDtracker::MapTimeBinLayers: \n");
2102 printf(" index %d exceeds allowed maximum of %d!\n",
2103 index, kMAX_TIME_BIN_INDEX-1);
2104 continue;
2105 }
2106 fTimeBinIndex[index] = i;
2107 }
2108
2109 Double_t x1, dx1, x2, dx2, gap;
2110
2111 for(Int_t i = 0; i < fN-1; i++) {
2112 x1 = fLayers[i]->GetX();
2113 dx1 = fLayers[i]->GetdX();
2114 x2 = fLayers[i+1]->GetX();
2115 dx2 = fLayers[i+1]->GetdX();
2116 gap = (x2 - dx2/2) - (x1 + dx1/2);
2117 if(gap < -0.01) {
2118 printf("*** warning: layers %d and %d are overlayed:\n",i,i+1);
2119 printf(" %f + %f + %f > %f\n", x1, dx1/2, dx2/2, x2);
2120 }
2121 if(gap > 0.01) {
2122 printf("*** warning: layers %d and %d have a large gap:\n",i,i+1);
2123 printf(" (%f - %f) - (%f + %f) = %f\n",
2124 x2, dx2/2, x1, dx1, gap);
2125 }
2126 }
2127}
2128
2129
2130//______________________________________________________
2131
2132
2133Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const
2134{
2135 //
2136 // Returns the number of time bin which in radial position is closest to <x>
2137 //
2138
2139 if(x >= fLayers[fN-1]->GetX()) return fN-1;
2140 if(x <= fLayers[0]->GetX()) return 0;
2141
2142 Int_t b=0, e=fN-1, m=(b+e)/2;
2143 for (; b<e; m=(b+e)/2) {
2144 if (x > fLayers[m]->GetX()) b=m+1;
2145 else e=m;
2146 }
2147 if(TMath::Abs(x - fLayers[m]->GetX()) >
2148 TMath::Abs(x - fLayers[m+1]->GetX())) return m+1;
2149 else return m;
2150
2151}
2152
2153//______________________________________________________
2154
2155Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const
2156{
2157 //
2158 // Returns number of the innermost SENSITIVE propagation layer
2159 //
2160
2161 return GetLayerNumber(0);
2162}
2163
2164//______________________________________________________
2165
2166Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const
2167{
2168 //
2169 // Returns number of the outermost SENSITIVE time bin
2170 //
2171
2172 return GetLayerNumber(GetNumberOfTimeBins() - 1);
46d29e70 2173}
2174
5443e65e 2175//______________________________________________________
2176
2177Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const
2178{
2179 //
2180 // Returns number of SENSITIVE time bins
2181 //
2182
2183 Int_t tb, layer;
2184 for(tb = kMAX_TIME_BIN_INDEX-1; tb >=0; tb--) {
2185 layer = GetLayerNumber(tb);
2186 if(layer>=0) break;
2187 }
2188 return tb+1;
2189}
2190
2191//______________________________________________________
2192
2193void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl)
2194{
2195 //
2196 // Insert layer <pl> in fLayers array.
2197 // Layers are sorted according to X coordinate.
2198
2199 if ( fN == ((Int_t) kMAX_LAYERS_PER_SECTOR)) {
2200 printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
2201 return;
2202 }
2203 if (fN==0) {fLayers[fN++] = pl; return;}
2204 Int_t i=Find(pl->GetX());
2205
2206 memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
2207 fLayers[i]=pl; fN++;
2208
2209}
2210
2211//______________________________________________________
2212
2213Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const
2214{
2215 //
2216 // Returns index of the propagation layer nearest to X
2217 //
2218
2219 if (x <= fLayers[0]->GetX()) return 0;
2220 if (x > fLayers[fN-1]->GetX()) return fN;
2221 Int_t b=0, e=fN-1, m=(b+e)/2;
2222 for (; b<e; m=(b+e)/2) {
2223 if (x > fLayers[m]->GetX()) b=m+1;
2224 else e=m;
2225 }
2226 return m;
2227}
2228
2229//______________________________________________________
2230
2231void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
2232 Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &x0,
2233 Bool_t &lookForCluster) const
2234{
2235 //
2236 // Returns radial step <dx>, density <rho>, rad. length <x0>,
2237 // and sensitivity <lookForCluster> in point <y,z>
2238 //
2239
2240 dx = fdX;
2241 rho = fRho;
2242 x0 = fX0;
2243 lookForCluster = kFALSE;
2244
2245 // check dead regions
2246 if(fTimeBinIndex >= 0) {
2247 for(Int_t ch = 0; ch < (Int_t) kZONES; ch++) {
2248 if(TMath::Abs(z - fZc[ch]) < fZmax[ch])
2249 lookForCluster = kTRUE;
2250 }
2251 if(TMath::Abs(y) > fYmax) lookForCluster = kFALSE;
2252 if(!lookForCluster) {
2253 // rho = 1.7; x0 = 33.0; // G10
2254 }
2255 }
2256
2257 // check hole
2258 if(fHole && (TMath::Abs(y - fHoleYc) < fHoleYmax) &&
2259 (TMath::Abs(z - fHoleZc) < fHoleZmax)) {
2260 lookForCluster = kFALSE;
2261 rho = fHoleRho;
2262 x0 = fHoleX0;
2263 }
2264
2265 return;
2266}
2267
2268//______________________________________________________
2269
2270void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c,
2271 UInt_t index) {
2272
2273// Insert cluster in cluster array.
2274// Clusters are sorted according to Y coordinate.
2275
2276 if(fTimeBinIndex < 0) {
2277 printf("*** attempt to insert cluster into non-sensitive time bin!\n");
2278 return;
2279 }
2280
2281 if (fN== (Int_t) kMAX_CLUSTER_PER_TIME_BIN) {
2282 printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n");
2283 return;
2284 }
2285 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
2286 Int_t i=Find(c->GetY());
2287 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
2288 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
2289 fIndex[i]=index; fClusters[i]=c; fN++;
2290}
2291
2292//______________________________________________________
2293
2294Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Double_t y) const {
2295
2296// Returns index of the cluster nearest in Y
2297
2298 if (y <= fClusters[0]->GetY()) return 0;
2299 if (y > fClusters[fN-1]->GetY()) return fN;
2300 Int_t b=0, e=fN-1, m=(b+e)/2;
2301 for (; b<e; m=(b+e)/2) {
2302 if (y > fClusters[m]->GetY()) b=m+1;
2303 else e=m;
2304 }
2305 return m;
2306}
2307
fd621f36 2308//---------------------------------------------------------
5443e65e 2309
fd621f36 2310Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
2311//
2312// Returns correction factor for tilted pads geometry
2313//
5443e65e 2314
fd621f36 2315 Double_t h01 = sin(TMath::Pi() / 180.0 * fPar->GetTiltingAngle());
2316 Int_t det = c->GetDetector();
2317 Int_t plane = fGeom->GetPlane(det);
2318 if((plane == 0) || (plane == 2) || (plane == 4)) h01=-h01;
2319
2320 return h01;
2321}
5443e65e 2322
2323
2324
2325
2326
2327