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