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