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