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