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