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