When calculating a*a-b*b the form (a-b)*(a+b) is usually more numerically stable.
[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// //
c6f438c0 20// The standard TRD tracker //
21// Based on Kalman filltering approach //
029cd327 22// //
4551ea7c 23// Authors: //
24// M. Ivanov (Marian.Ivanov@cern.ch) //
4551ea7c 25// //
029cd327 26///////////////////////////////////////////////////////////////////////////////
27
46d29e70 28#include <TBranch.h>
a1e17193 29#include <TFile.h>
4551ea7c 30#include <TGraph.h>
a1e17193 31#include <TLinearFitter.h>
32#include <TObjArray.h>
33#include <TROOT.h>
34#include <TTree.h>
35#include <TTreeStream.h>
4551ea7c 36
af885e0f 37#include "AliESDEvent.h"
f3f19a6e 38#include "AliESDtrack.h"
4551ea7c 39#include "AliAlignObj.h"
40#include "AliRieman.h"
41#include "AliTrackPointArray.h"
46d29e70 42
46d29e70 43#include "AliTRDgeometry.h"
a5cadd36 44#include "AliTRDpadPlane.h"
0b2318ec 45#include "AliTRDgeometry.h"
46d29e70 46#include "AliTRDcluster.h"
47#include "AliTRDtrack.h"
75bd7f81 48#include "AliTRDseed.h"
3551db50 49#include "AliTRDcalibDB.h"
50#include "AliTRDCommonParam.h"
46d29e70 51#include "AliTRDtracker.h"
7b580082 52#include "AliTRDReconstructor.h"
d20df6fc 53#include "AliTRDrecoParam.h"
082a2283 54#include "AliTRDCalibraFillHisto.h"
55
7e448bcc 56ClassImp(AliTRDtracker)
46d29e70 57
b08a2e7e 58const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5; //
59const Float_t AliTRDtracker::fgkLabelFraction = 0.8; //
60const Double_t AliTRDtracker::fgkMaxChi2 = 12.0; //
0c349049 61const Double_t AliTRDtracker::fgkMaxSnp = 0.95; // Maximum local sine of the azimuthal angle
7e448bcc 62const Double_t AliTRDtracker::fgkMaxStep = 2.0; // Maximal step size in propagation
7ad19338 63
75bd7f81 64//_____________________________________________________________________________
3a039a31 65AliTRDtracker::AliTRDtracker(AliTRDReconstructor *rec)
4e009ce4 66 :AliTracker()
3a039a31 67 ,fReconstructor(rec)
a4cf6278 68 ,fGeom(0)
69 ,fNclusters(0)
70 ,fClusters(0)
71 ,fNseeds(0)
72 ,fSeeds(0)
73 ,fNtracks(0)
74 ,fTracks(0)
75 ,fTimeBinsPerPlane(0)
76 ,fAddTRDseeds(kFALSE)
77 ,fNoTilt(kFALSE)
4e009ce4 78 ,fHBackfit(0x0)
79 ,fHClSearch(0x0)
80 ,fHRefit(0x0)
81 ,fHX(0x0)
82 ,fHNCl(0x0)
83 ,fHNClTrack(0x0)
84 ,fHMinYPos(0x0)
85 ,fHMinYNeg(0x0)
86 ,fHMinZ(0x0)
87 ,fHMinD(0x0)
88 ,fHDeltaX(0x0)
89 ,fHXCl(0x0)
4e009ce4 90 ,fDebugStreamer(0)
89f05372 91{
75bd7f81 92 //
b7a0917f 93 // Default constructor
75bd7f81 94 //
b7a0917f 95
4551ea7c 96 for (Int_t i = 0; i < kTrackingSectors; i++) {
97 fTrSec[i] = 0;
98 }
75bd7f81 99
7e448bcc 100 InitLogHists();
4e009ce4 101
89f05372 102}
75bd7f81 103
104//_____________________________________________________________________________
6c94f330 105AliTRDtracker::AliTRDtracker(const AliTRDtracker &t)
106 :AliTracker(t)
3a039a31 107 ,fReconstructor(t.fReconstructor)
a4cf6278 108 ,fGeom(0)
109 ,fNclusters(0)
110 ,fClusters(0)
111 ,fNseeds(0)
112 ,fSeeds(0)
113 ,fNtracks(0)
114 ,fTracks(0)
115 ,fTimeBinsPerPlane(0)
116 ,fAddTRDseeds(kFALSE)
117 ,fNoTilt(kFALSE)
4e009ce4 118 ,fHBackfit(0x0)
119 ,fHClSearch(0x0)
120 ,fHRefit(0x0)
121 ,fHX(0x0)
122 ,fHNCl(0x0)
123 ,fHNClTrack(0x0)
124 ,fHMinYPos(0x0)
125 ,fHMinYNeg(0x0)
126 ,fHMinZ(0x0)
127 ,fHMinD(0x0)
128 ,fHDeltaX(0x0)
129 ,fHXCl(0x0)
ad670fba 130 ,fDebugStreamer(0)
6c94f330 131{
ad670fba 132 //
133 // Copy constructor
134 //
4e009ce4 135
ad670fba 136}
137
138//_____________________________________________________________________________
3a039a31 139AliTRDtracker::AliTRDtracker(const TFile */*geomfile*/, AliTRDReconstructor *rec)
ad670fba 140 :AliTracker()
3a039a31 141 ,fReconstructor(rec)
a4cf6278 142 ,fGeom(0)
143 ,fNclusters(0)
144 ,fClusters(new TObjArray(2000))
145 ,fNseeds(0)
146 ,fSeeds(new TObjArray(2000))
147 ,fNtracks(0)
148 ,fTracks(new TObjArray(1000))
149 ,fTimeBinsPerPlane(0)
150 ,fAddTRDseeds(kFALSE)
151 ,fNoTilt(kFALSE)
4e009ce4 152 ,fHBackfit(0x0)
153 ,fHClSearch(0x0)
154 ,fHRefit(0x0)
155 ,fHX(0x0)
156 ,fHNCl(0x0)
157 ,fHNClTrack(0x0)
158 ,fHMinYPos(0x0)
159 ,fHMinYNeg(0x0)
160 ,fHMinZ(0x0)
161 ,fHMinD(0x0)
162 ,fHDeltaX(0x0)
163 ,fHXCl(0x0)
ad670fba 164 ,fDebugStreamer(0)
46d29e70 165{
5443e65e 166 //
167 // Main constructor
168 //
b8dc2353 169
4551ea7c 170 TDirectory *savedir = gDirectory;
46d29e70 171
f162af62 172 fGeom = new AliTRDgeometry();
c630aafd 173
4551ea7c 174 for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
439c63c8 175 Int_t trS = geomS;
4551ea7c 176 fTrSec[trS] = new AliTRDtrackingSector(fGeom,geomS);
5443e65e 177 }
4551ea7c 178
f162af62 179 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0);
7ad19338 180 Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
4551ea7c 181 if (tiltAngle < 0.1) {
b8dc2353 182 fNoTilt = kTRUE;
183 }
184
fc88348b 185 if (!AliTRDcalibDB::Instance()) {
186 AliFatal("Could not get calibration object");
187 }
188 fTimeBinsPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
a819a5f7 189
4551ea7c 190 fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
46d29e70 191
9c9d2487 192 savedir->cd();
7e448bcc 193
194 InitLogHists();
75bd7f81 195
5443e65e 196}
46d29e70 197
75bd7f81 198//_____________________________________________________________________________
5443e65e 199AliTRDtracker::~AliTRDtracker()
46d29e70 200{
029cd327 201 //
202 // Destructor of AliTRDtracker
203 //
204
89f05372 205 if (fClusters) {
206 fClusters->Delete();
207 delete fClusters;
208 }
4551ea7c 209
89f05372 210 if (fTracks) {
211 fTracks->Delete();
212 delete fTracks;
213 }
4551ea7c 214
89f05372 215 if (fSeeds) {
216 fSeeds->Delete();
217 delete fSeeds;
218 }
4551ea7c 219
f162af62 220 if (fGeom) {
221 delete fGeom;
222 }
0a29d0f1 223
4551ea7c 224 for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
029cd327 225 delete fTrSec[geomS];
5443e65e 226 }
4551ea7c 227
7ad19338 228 if (fDebugStreamer) {
7ad19338 229 delete fDebugStreamer;
230 }
9c9d2487 231
75bd7f81 232}
59393e34 233
75bd7f81 234//_____________________________________________________________________________
235Int_t AliTRDtracker::LocalToGlobalID(Int_t lid)
236{
59393e34 237 //
75bd7f81 238 // Transform internal TRD ID to global detector ID
59393e34 239 //
75bd7f81 240
053767a4 241 Int_t isector = fGeom->GetSector(lid);
242 Int_t istack = fGeom->GetStack(lid);
243 Int_t ilayer = fGeom->GetLayer(lid);
4551ea7c 244
ae079791 245 AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
053767a4 246 switch (ilayer) {
59393e34 247 case 0:
ae079791 248 iLayer = AliGeomManager::kTRD1;
59393e34 249 break;
250 case 1:
ae079791 251 iLayer = AliGeomManager::kTRD2;
59393e34 252 break;
253 case 2:
ae079791 254 iLayer = AliGeomManager::kTRD3;
59393e34 255 break;
256 case 3:
ae079791 257 iLayer = AliGeomManager::kTRD4;
59393e34 258 break;
259 case 4:
ae079791 260 iLayer = AliGeomManager::kTRD5;
59393e34 261 break;
262 case 5:
ae079791 263 iLayer = AliGeomManager::kTRD6;
59393e34 264 break;
265 };
4551ea7c 266
053767a4 267 Int_t modId = isector * fGeom->Nstack() + istack;
ae079791 268 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,modId);
75bd7f81 269
59393e34 270 return volid;
75bd7f81 271
59393e34 272}
273
75bd7f81 274//_____________________________________________________________________________
275Int_t AliTRDtracker::GlobalToLocalID(Int_t gid)
276{
59393e34 277 //
75bd7f81 278 // Transform global detector ID to local detector ID
59393e34 279 //
75bd7f81 280
4551ea7c 281 Int_t modId = 0;
ae079791 282 AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(gid,modId);
4551ea7c 283
053767a4 284 Int_t isector = modId / fGeom->Nstack();
285 Int_t istack = modId % fGeom->Nstack();
286 Int_t iLayer = -1;
4551ea7c 287
59393e34 288 switch (layerId) {
ae079791 289 case AliGeomManager::kTRD1:
6c94f330 290 iLayer = 0;
59393e34 291 break;
ae079791 292 case AliGeomManager::kTRD2:
6c94f330 293 iLayer = 1;
59393e34 294 break;
ae079791 295 case AliGeomManager::kTRD3:
6c94f330 296 iLayer = 2;
59393e34 297 break;
ae079791 298 case AliGeomManager::kTRD4:
6c94f330 299 iLayer = 3;
59393e34 300 break;
ae079791 301 case AliGeomManager::kTRD5:
6c94f330 302 iLayer = 4;
59393e34 303 break;
ae079791 304 case AliGeomManager::kTRD6:
6c94f330 305 iLayer = 5;
59393e34 306 break;
307 default:
6c94f330 308 iLayer =-1;
59393e34 309 }
4551ea7c 310
311 if (iLayer < 0) {
312 return -1;
313 }
314
053767a4 315 Int_t lid = fGeom->GetDetector(iLayer,istack,isector);
75bd7f81 316
59393e34 317 return lid;
59393e34 318
75bd7f81 319}
59393e34 320
75bd7f81 321//_____________________________________________________________________________
75bd7f81 322Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track)
323{
9c9d2487 324 //
325 // Rotates the track when necessary
326 //
327
328 Double_t alpha = AliTRDgeometry::GetAlpha();
4551ea7c 329 Double_t y = track->GetY();
330 Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
9c9d2487 331
4551ea7c 332 if (y > ymax) {
4551ea7c 333 if (!track->Rotate( alpha)) {
334 return kFALSE;
335 }
336 }
337 else if (y < -ymax) {
4551ea7c 338 if (!track->Rotate(-alpha)) {
339 return kFALSE;
340 }
9c9d2487 341 }
342
343 return kTRUE;
9c9d2487 344
75bd7f81 345}
46e2d86c 346
75bd7f81 347//_____________________________________________________________________________
348AliTRDcluster *AliTRDtracker::GetCluster(AliTRDtrack *track, Int_t plane
349 , Int_t timebin, UInt_t &index)
350{
46e2d86c 351 //
75bd7f81 352 // Try to find cluster in the backup list
46e2d86c 353 //
75bd7f81 354
4551ea7c 355 AliTRDcluster *cl =0;
6c94f330 356 Int_t *indexes = track->GetBackupIndexes();
4551ea7c 357
358 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
359 if (indexes[i] == 0) {
360 break;
361 }
362 AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
363 if (!cli) {
364 break;
365 }
366 if (cli->GetLocalTimeBin() != timebin) {
367 continue;
368 }
053767a4 369 Int_t ilayer = fGeom->GetLayer(cli->GetDetector());
370 if (ilayer == plane) {
4551ea7c 371 cl = cli;
7ad19338 372 index = indexes[i];
46e2d86c 373 break;
374 }
375 }
75bd7f81 376
46e2d86c 377 return cl;
46e2d86c 378
75bd7f81 379}
3c625a9b 380
75bd7f81 381//_____________________________________________________________________________
4551ea7c 382Int_t AliTRDtracker::GetLastPlane(AliTRDtrack *track)
75bd7f81 383{
3c625a9b 384 //
75bd7f81 385 // Return last updated plane
386 //
387
4551ea7c 388 Int_t lastplane = 0;
389 Int_t *indexes = track->GetBackupIndexes();
390
391 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
392 AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
393 if (!cli) {
394 break;
395 }
053767a4 396 Int_t ilayer = fGeom->GetLayer(cli->GetDetector());
397 if (ilayer > lastplane) {
398 lastplane = ilayer;
3c625a9b 399 }
400 }
75bd7f81 401
3c625a9b 402 return lastplane;
75bd7f81 403
3c625a9b 404}
5443e65e 405
5443e65e 406//_____________________________________________________________________________
af885e0f 407Int_t AliTRDtracker::PropagateBack(AliESDEvent *event)
75bd7f81 408{
c630aafd 409 //
410 // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
411 // backpropagated by the TPC tracker. Each seed is first propagated
412 // to the TRD, and then its prolongation is searched in the TRD.
413 // If sufficiently long continuation of the track is found in the TRD
414 // the track is updated, otherwise it's stored as originaly defined
415 // by the TPC tracker.
416 //
417
e4f2f73d 418 Int_t found = 0; // number of tracks found
419 Float_t foundMin = 20.0;
4f1c04d3 420
e4f2f73d 421 Int_t nSeed = event->GetNumberOfTracks();
422 if(!nSeed){
423 // run stand alone tracking
3a039a31 424 if (fReconstructor->IsSeeding()) Clusters2Tracks(event);
e4f2f73d 425 return 0;
4f1c04d3 426 }
e4f2f73d 427
428 Float_t *quality = new Float_t[nSeed];
429 Int_t *index = new Int_t[nSeed];
430 for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
431 AliESDtrack *seed = event->GetTrack(iSeed);
432 Double_t covariance[15];
433 seed->GetExternalCovariance(covariance);
434 quality[iSeed] = covariance[0] + covariance[2];
7ad19338 435 }
e4f2f73d 436 // Sort tracks according to covariance of local Y and Z
437 TMath::Sort(nSeed,quality,index,kFALSE);
7e448bcc 438
e4f2f73d 439 // Backpropagate all seeds
440 for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
7e448bcc 441
e4f2f73d 442 // Get the seeds in sorted sequence
443 AliESDtrack *seed = event->GetTrack(index[iSeed]);
444 fHBackfit->Fill(0); // All seeds
445
446 // Check the seed status
447 ULong_t status = seed->GetStatus();
448 if ((status & AliESDtrack::kTPCout) == 0) {
449 fHBackfit->Fill(1); // TPC outer edge reached
450 continue;
451 }
452 if ((status & AliESDtrack::kTRDout) != 0) {
453 fHBackfit->Fill(2); // TRD outer edge reached (does this happen ?)
454 continue;
455 }
456
457 // Do the back prolongation
458 Int_t lbl = seed->GetLabel();
459 AliTRDtrack *track = new AliTRDtrack(*seed);
460 track->SetSeedLabel(lbl);
461 seed->UpdateTrackParams(track,AliESDtrack::kTRDbackup); // Make backup
462 fNseeds++;
463 Float_t p4 = track->GetC();
464 Int_t expectedClr = FollowBackProlongation(*track);
465 fHBackfit->Fill(3); // Back prolongation done
466 fHX->Fill(track->GetX());
467
468 if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
469 (track->Pt() > 0.8)) {
470
471 fHBackfit->Fill(4);
472
473 //
474 // Make backup for back propagation
475 //
476
477 Int_t foundClr = track->GetNumberOfClusters();
478 if (foundClr >= foundMin) {
479 track->CookdEdx();
480 track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
481 CookLabel(track,1 - fgkLabelFraction);
482 if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
483
484
485 // Sign only gold tracks
486 if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
487 if ((seed->GetKinkIndex(0) == 0) &&
488 (track->Pt() < 1.5)) {
489 UseClusters(track);
490 }
491 }
492 Bool_t isGold = kFALSE;
493
494 // Full gold track
495 if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
496 //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
497 if (track->GetBackupTrack()) {
498 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
499 }
500 isGold = kTRUE;
501 //fHBackfit->Fill()
502 }
503
504 // Almost gold track
505 if ((!isGold) && (track->GetNCross() == 0) &&
506 (track->GetChi2() / track->GetNumberOfClusters() < 7)) {
507 //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
508 if (track->GetBackupTrack()) {
509 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
510 }
511 isGold = kTRUE;
512 }
513
514 if ((!isGold) && (track->GetBackupTrack())) {
515 if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) && ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
516 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
517 isGold = kTRUE;
518 }
519 }
520
521 if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected()) > 0.4)) {
522 //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
523 }
524 }
525 }
526 /**/
527
528 /**/
529 // Debug part of tracking
530 TTreeSRedirector &cstream = *fDebugStreamer;
531 Int_t eventNrInFile = event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
3a039a31 532 if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0) {
e4f2f73d 533 if (track->GetBackupTrack()) {
534 cstream << "Tracks"
535 << "EventNrInFile=" << eventNrInFile
536 << "ESD.=" << seed
537 << "trd.=" << track
538 << "trdback.=" << track->GetBackupTrack()
539 << "\n";
540 }
541 else {
542 cstream << "Tracks"
543 << "EventNrInFile=" << eventNrInFile
544 << "ESD.=" << seed
545 << "trd.=" << track
546 << "trdback.=" << track
547 << "\n";
548 }
549 }
550 /**/
551
552 // Propagation to the TOF (I.Belikov)
553 if (track->GetStop() == kFALSE) {
554 fHBackfit->Fill(5);
555
556 Double_t xtof = 371.0;
557 Double_t xTOF0 = 370.0;
558
559 Double_t c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
560 if (TMath::Abs(c2) >= 0.99) {
561 fHBackfit->Fill(6);
562 delete track;
563 continue;
564 }
565
566 PropagateToX(*track,xTOF0,fgkMaxStep);
567
568 // Energy losses taken to the account - check one more time
569 c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
570 if (TMath::Abs(c2) >= 0.99) {
571 fHBackfit->Fill(7);
572 delete track;
573 continue;
574 }
575
576 //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
577 // fHBackfit->Fill(7);
578 //delete track;
579 // continue;
580 //}
581
582 Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
583 Double_t y;
584 track->GetYAt(xtof,GetBz(),y);
585 if (y > ymax) {
586 if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
587 fHBackfit->Fill(8);
588 delete track;
589 continue;
590 }
591 }
592 else if (y < -ymax) {
593 if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
594 fHBackfit->Fill(9);
595 delete track;
596 continue;
597 }
598 }
599
600 if (track->PropagateTo(xtof)) {
601 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
602 fHBackfit->Fill(10);
603
216b200f 604 seed->SetNumberOfTRDslices(AliTRDtrack::kNslice);
6984f7c1 605 for (Int_t i = 0; i < AliTRDtrack::kNplane; i++) {
606 for (Int_t j = 0; j < AliTRDtrack::kNslice; j++) {
607 seed->SetTRDslice(track->GetPIDsignals(i,j),i,j);
e4f2f73d 608 }
609 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
610 }
611 //seed->SetTRDtrack(new AliTRDtrack(*track));
612 if (track->GetNumberOfClusters() > foundMin) {
613 fHBackfit->Fill(11);
614 found++;
615 }
616 }
617 }
618 else {
619 fHBackfit->Fill(12);
620
621 if ((track->GetNumberOfClusters() > 15) &&
622 (track->GetNumberOfClusters() > 0.5*expectedClr)) {
623 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
624 fHBackfit->Fill(13);
625
626 //seed->SetStatus(AliESDtrack::kTRDStop);
216b200f 627
628 seed->SetNumberOfTRDslices(AliTRDtrack::kNslice);
6984f7c1 629 for (Int_t i = 0; i < AliTRDtrack::kNplane; i++) {
630 for (Int_t j = 0; j <AliTRDtrack::kNslice; j++) {
631 seed->SetTRDslice(track->GetPIDsignals(i,j),i,j);
e4f2f73d 632 }
633 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
634 }
635 //seed->SetTRDtrack(new AliTRDtrack(*track));
636 found++;
637 }
638 }
639
640 seed->SetTRDQuality(track->StatusForTOF());
641 seed->SetTRDBudget(track->GetBudget(0));
642
643 fHBackfit->Fill(14);
644 delete track;
4551ea7c 645 }
7e448bcc 646
e4f2f73d 647 AliInfo(Form("Number of seeds: %d",fNseeds));
648 AliInfo(Form("Number of back propagated TRD tracks: %d",found));
649
650 fSeeds->Clear();
651 fNseeds = 0;
652
653 delete [] index;
654 delete [] quality;
655
656 SaveLogHists();
7e448bcc 657
1e9bb598 658 return 0;
1e9bb598 659}
660
661//_____________________________________________________________________________
af885e0f 662Int_t AliTRDtracker::RefitInward(AliESDEvent *event)
1e9bb598 663{
664 //
665 // Refits tracks within the TRD. The ESD event is expected to contain seeds
666 // at the outer part of the TRD.
667 // The tracks are propagated to the innermost time bin
668 // of the TRD and the ESD event is updated
669 // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
670 //
671
0c349049 672 //Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
673 //Float_t foundMin = fgkMinClustersInTrack * timeBins;
4551ea7c 674 Int_t nseed = 0;
675 Int_t found = 0;
44dbae42 676 Int_t pidQ = 0;
4551ea7c 677 //Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
4f1c04d3 678 AliTRDtrack seed2;
6c94f330 679
bcb6fb78 680 // Calibration fill 2D
681 AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
682 if (!calibra) {
683 AliInfo("Could not get Calibra instance\n");
684 }
685
1e9bb598 686 Int_t n = event->GetNumberOfTracks();
4551ea7c 687 for (Int_t i = 0; i < n; i++) {
688
689 AliESDtrack *seed = event->GetTrack(i);
690 new (&seed2) AliTRDtrack(*seed);
7e448bcc 691 fHRefit->Fill(0);
692
4551ea7c 693 if (seed2.GetX() < 270.0) {
694 seed->UpdateTrackParams(&seed2,AliESDtrack::kTRDbackup); // Backup TPC track - only update
7e448bcc 695 fHRefit->Fill(1);
f4e9508c 696 continue;
697 }
698
4551ea7c 699 ULong_t status = seed->GetStatus();
700 if ((status & AliESDtrack::kTRDout) == 0) {
7e448bcc 701 fHRefit->Fill(2);
0dd7d129 702 continue;
703 }
4551ea7c 704 if ((status & AliESDtrack::kTRDin) != 0) {
7e448bcc 705 fHRefit->Fill(3);
0dd7d129 706 continue;
707 }
7e448bcc 708
6c94f330 709 nseed++;
7e448bcc 710 fHRefit->Fill(4);
6c94f330 711
4551ea7c 712 seed2.ResetCovariance(50.0);
6c94f330 713
4f1c04d3 714 AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
4551ea7c 715 Int_t *indexes2 = seed2.GetIndexes();
7bce990c 716 for (Int_t l = 0; l < AliTRDtrack::kNplane;++l) {
6984f7c1 717 for (Int_t j = 0; j < AliTRDtrack::kNslice;j++) {
7bce990c 718 pt->SetPIDsignals(seed2.GetPIDsignals(l,j),l,j);
6d45eaef 719 }
7bce990c 720 pt->SetPIDTimBin(seed2.GetPIDTimBin(l),l);
7ad19338 721 }
eab5961e 722
4551ea7c 723 Int_t *indexes3 = pt->GetBackupIndexes();
7bce990c 724 for (Int_t l = 0; l < 200;++l) {
725 if (indexes2[l] == 0) {
4551ea7c 726 break;
727 }
7bce990c 728 indexes3[l] = indexes2[l];
4551ea7c 729 }
730
43bfc8af 731 FollowProlongation(*pt);
0c349049 732 pt->CookdEdx();
733 pt->CookdEdxTimBin(seed->GetID());
7eccf2f3 734
735 //calculate PID methods
736 pt->SetPIDMethod(AliTRDtrack::kLQ);
0c349049 737 pt->CookPID(pidQ);
7eccf2f3 738 seed->SetTRDpid(pt->GetPID());
739 seed->SetTRDpidQuality(pidQ);
740
741 // update calibration
bcb6fb78 742 if(calibra->GetHisto2d()) calibra->UpdateHistograms(pt);
1e9bb598 743 found++;
33744e5d 744
4551ea7c 745 Double_t xTPC = 250.0;
43bfc8af 746 if (PropagateToX(*pt,xTPC,fgkMaxStep)) {
4551ea7c 747 seed->UpdateTrackParams(pt,AliESDtrack::kTRDrefit);
7e448bcc 748 fHRefit->Fill(5);
749
7bce990c 750 for (Int_t l = 0; l < AliTRDtrack::kNplane; ++l) {
6984f7c1 751 for (Int_t j = 0; j < AliTRDtrack::kNslice; j++) {
7bce990c 752 seed->SetTRDslice(pt->GetPIDsignals(l,j),l,j);
7eccf2f3 753 }
7bce990c 754 seed->SetTRDTimBin(pt->GetPIDTimBin(l),l);
7ad19338 755 }
7eccf2f3 756 } else {
4551ea7c 757 // If not prolongation to TPC - propagate without update
7e448bcc 758 fHRefit->Fill(5);
7bce990c 759 AliTRDtrack *seed2t = new AliTRDtrack(*seed);
760 seed2t->ResetCovariance(5.0);
761 AliTRDtrack *pt2 = new AliTRDtrack(*seed2t,seed2t->GetAlpha());
762 delete seed2t;
0c349049 763
59393e34 764 if (PropagateToX(*pt2,xTPC,fgkMaxStep)) {
0c349049 765 pt2->CookdEdx();
44dbae42 766 pt2->CookdEdxTimBin(seed->GetID());
7eccf2f3 767 seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit);
768 fHRefit->Fill(6);
7e448bcc 769
7bce990c 770 for (Int_t l = 0; l < AliTRDtrack::kNplane; ++l) {
6984f7c1 771 for (Int_t j = 0; j < AliTRDtrack::kNslice; j++) {
7bce990c 772 seed->SetTRDslice(pt2->GetPIDsignals(l,j),l,j);
7eccf2f3 773 }
7bce990c 774 seed->SetTRDTimBin(pt2->GetPIDTimBin(l),l);
7ad19338 775 }
eab5961e 776 }
306e3dc1 777
0c349049 778 // Add TRD track to ESDfriendTrack - maybe this tracks are
779 // not useful for post-processing - TODO make decision
3a039a31 780 if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0) {
306e3dc1 781 seed->AddCalibObject(new AliTRDtrack(*pt2/*, kTRUE*/));
782 }
7bed16a7 783 delete pt2;
0c349049 784
4551ea7c 785 }
786
306e3dc1 787 // Add TRD track to ESDfriendTrack
3a039a31 788 if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0) {
306e3dc1 789 seed->AddCalibObject(new AliTRDtrack(*pt/*, kTRUE*/));
790 }
1e9bb598 791 delete pt;
4551ea7c 792
eab5961e 793 }
1e9bb598 794
33744e5d 795 AliInfo(Form("Number of loaded seeds: %d",nseed));
796 AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
1e9bb598 797
7e448bcc 798 SaveLogHists();
0c349049 799
c630aafd 800 return 0;
0c349049 801
c630aafd 802}
803
75bd7f81 804//_____________________________________________________________________________
4551ea7c 805Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t)
8979685e 806{
0c349049 807 //
808 // Starting from current position on track=t this function tries
809 // to extrapolate the track up to timeBin=0 and to confirm prolongation
810 // if a close cluster is found. Returns the number of clusters
811 // expected to be found in sensitive layers
812 // GeoManager used to estimate mean density
813 //
4551ea7c 814
0c349049 815 Int_t sector;
816 Int_t lastplane = GetLastPlane(&t);
817 Double_t xx0 = 0.0;
818 Double_t xrho= 0.0;
819 Int_t expectedNumberOfClusters = 0;
4551ea7c 820
0c349049 821 for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
4551ea7c 822
0c349049 823 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
824 Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
4551ea7c 825
0c349049 826 // Propagate track close to the plane if neccessary
827 Double_t currentx = fTrSec[0]->GetLayer(rowlast)->GetX();
828 if (currentx < (-fgkMaxStep + t.GetX())) {
829 // Propagate closer to chamber - safety space fgkMaxStep
830 if (!PropagateToX(t,currentx+fgkMaxStep,fgkMaxStep)) {
831 break;
832 }
833 }
4551ea7c 834
0c349049 835 if (!AdjustSector(&t)) {
836 break;
837 }
4551ea7c 838
0c349049 839 // Get material budget
840 Double_t xyz0[3];
841 Double_t xyz1[3];
e4f2f73d 842 Double_t param[7];
0c349049 843 Double_t x;
844 Double_t y;
845 Double_t z;
4551ea7c 846
0c349049 847 // Starting global position
848 t.GetXYZ(xyz0);
849 // End global position
850 x = fTrSec[0]->GetLayer(row0)->GetX();
851 if (!t.GetProlongation(x,y,z)) {
852 break;
853 }
854 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
855 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
856 xyz1[2] = z;
857 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
858 xrho= param[0]*param[4];
859 xx0 = param[1]; // Get mean propagation parameters
860
861 // Flags for marking the track momentum and direction. The track is
862 // marked only if it has at least 1 cluster picked up in the current
863 // chamber.
864 Bool_t kUPDATED = kFALSE;
865 Bool_t kMARKED = kFALSE;
866
867 // Propagate and update
868 sector = t.GetSector();
869 //for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
870 for (Int_t itime = 0 ; itime < GetTimeBinsPerPlane(); itime++) {
871
872 // Mark track kinematics
873 if (itime > 10 && kUPDATED && !kMARKED) {
874 t.SetTrackSegmentDirMom(iplane);
875 kMARKED = kTRUE;
876 }
877
878 Int_t ilayer = GetGlobalTimeBin(0,iplane,itime);
879 expectedNumberOfClusters++;
880 t.SetNExpected(t.GetNExpected() + 1);
881 if (t.GetX() > 345.0) {
882 t.SetNExpectedLast(t.GetNExpectedLast() + 1);
883 }
884 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
885 AliTRDcluster *cl = 0;
886 UInt_t index = 0;
887 Double_t maxChi2 = fgkMaxChi2;
888 x = timeBin.GetX();
889
890 if (timeBin) {
891
892 AliTRDcluster *cl0 = timeBin[0];
893 if (!cl0) {
894 // No clusters in given time bin
895 continue;
4551ea7c 896 }
897
053767a4 898 Int_t layer = fGeom->GetLayer(cl0->GetDetector());
899 if (layer > lastplane) {
0c349049 900 continue;
901 }
902
903 Int_t timebin = cl0->GetLocalTimeBin();
053767a4 904 AliTRDcluster *cl2 = GetCluster(&t,layer,timebin,index);
0c349049 905
906 if (cl2) {
907 cl = cl2;
908 //Double_t h01 = GetTiltFactor(cl); //I.B's fix
909 //maxChi2=t.GetPredictedChi2(cl,h01);
910 }
911
912 if (cl) {
913
914 //if (cl->GetNPads()<5)
915 Double_t dxsample = timeBin.GetdX();
916 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
917 Double_t h01 = GetTiltFactor(cl);
918 Int_t det = cl->GetDetector();
053767a4 919 Int_t llayer = fGeom->GetLayer(det);
0c349049 920
921 if (t.GetX() > 345.0) {
922 t.SetNLast(t.GetNLast() + 1);
923 t.SetChi2Last(t.GetChi2Last() + maxChi2);
924 }
925
926 Double_t xcluster = cl->GetX();
927 t.PropagateTo(xcluster,xx0,xrho);
928 if (!AdjustSector(&t)) {
929 break; //I.B's fix
930 }
931
932 maxChi2 = t.GetPredictedChi2(cl,h01);
933 if (maxChi2 < 1e+10) {
053767a4 934 if (!t.UpdateMI(cl,maxChi2,index,h01,llayer)) {
0c349049 935 // ????
936 }
937 else {
938 //SetCluster(cl, GetNumberOfClusters()-1);
939 kUPDATED = kTRUE;
940 }
941 }
942
943 }
944
945 }
946
947 }
948
949 }
950
951 return expectedNumberOfClusters;
952
75bd7f81 953}
7b580082 954
75bd7f81 955//_____________________________________________________________________________
4551ea7c 956Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
8979685e 957{
75bd7f81 958 //
8979685e 959 // Starting from current radial position of track <t> this function
0c349049 960 // extrapolates the track up to the outer timebin and in the sensitive
8979685e 961 // layers confirms prolongation if a close cluster is found.
962 // Returns the number of clusters expected to be found in sensitive layers
0c349049 963 // Uses the geomanager for material description
75bd7f81 964 //
7e448bcc 965 // return number of assigned clusters ?
966 //
967
4551ea7c 968 Int_t sector;
0c349049 969
970 Double_t xx0 = 0.0;
971 Double_t xrho = 0.0;
972
973 Float_t ratio0 = 0.0;
974
4551ea7c 975 Int_t expectedNumberOfClusters = 0;
0c349049 976
8979685e 977 AliTRDtracklet tracklet;
33744e5d 978
0c349049 979 const Int_t kNclusters = 1000;
980 Int_t clusters[kNclusters];
981 for (Int_t i = 0; i < kNclusters; i++) {
982 clusters[i] = -1;
983 }
984
bcb6fb78 985// // Calibration fill 2D
986// AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
987// if (!calibra) {
988// AliInfo("Could not get Calibra instance\n");
989// }
990// if (calibra->GetMITracking()) {
991// calibra->ResetTrack();
992// }
77566f2a 993
0c349049 994 // Loop through the TRD planes
6984f7c1 995 for (Int_t iplane = 0; iplane < AliTRDtrack::kNplane; iplane++) {
4551ea7c 996
0c349049 997 Int_t hb = iplane * 10;
7e448bcc 998 fHClSearch->Fill(hb);
999
0c349049 1000 // Get the global time bin numbers for the first an last
1001 // local time bin of the given plane
4551ea7c 1002 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
1003 Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
0c349049 1004
1005 // Get the X coordinates of the propagation layer for the first time bin
4551ea7c 1006 Double_t currentx = fTrSec[0]->GetLayer(row0)->GetX();
1007 if (currentx < t.GetX()) {
7e448bcc 1008 fHClSearch->Fill(hb+1);
4551ea7c 1009 continue;
1010 }
1011
0c349049 1012 // Propagate closer to the current chamber if neccessary
4551ea7c 1013 if (currentx > (fgkMaxStep + t.GetX())) {
1014 if (!PropagateToX(t,currentx-fgkMaxStep,fgkMaxStep)) {
7e448bcc 1015 fHClSearch->Fill(hb+2);
4551ea7c 1016 break;
1017 }
1018 }
0c349049 1019
1020 // Rotate track to adjacent sector if neccessary
4551ea7c 1021 if (!AdjustSector(&t)) {
7e448bcc 1022 fHClSearch->Fill(hb+3);
4551ea7c 1023 break;
1024 }
0c349049 1025
1026 // Check whether azimuthal angle is getting too large
4551ea7c 1027 if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) {
7e448bcc 1028 fHClSearch->Fill(hb+4);
4551ea7c 1029 break;
8979685e 1030 }
4551ea7c 1031
4551ea7c 1032 Double_t xyz0[3];
1033 Double_t xyz1[3];
1034 Double_t param[7];
1035 Double_t x;
1036 Double_t y;
1037 Double_t z;
0c349049 1038 // Global start position (beginning of chamber)
4551ea7c 1039 t.GetXYZ(xyz0);
0c349049 1040 // X-position of the end of the chamber
7b580082 1041 x = fTrSec[0]->GetLayer(rowlast)->GetX();
0c349049 1042 // Get local Y and Z at the X-position of the end of the chamber
4551ea7c 1043 if (!t.GetProlongation(x,y,z)) {
7e448bcc 1044 fHClSearch->Fill(hb+5);
4551ea7c 1045 break;
1046 }
0c349049 1047 // Global end position (end of chamber)
4551ea7c 1048 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1049 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1050 xyz1[2] = z;
0c349049 1051
1052 // Calculate the mean material budget along the path inside the chamber
97ccb02b 1053 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
0c349049 1054 // The mean propagation parameters (density*length and radiation length)
97ccb02b 1055 xrho = param[0]*param[4];
0c349049 1056 xx0 = param[1];
4551ea7c 1057
0c349049 1058 // Find the clusters and tracklet along the path inside the chamber
1059 sector = t.GetSector();
4551ea7c 1060 Float_t ncl = FindClusters(sector,row0,rowlast,&t,clusters,tracklet);
7e448bcc 1061 fHNCl->Fill(tracklet.GetN());
1062
0c349049 1063 // Discard if in less than 1/3 of the available timebins
1064 // clusters are found
4551ea7c 1065 if (tracklet.GetN() < GetTimeBinsPerPlane()/3) {
7e448bcc 1066 fHClSearch->Fill(hb+6);
4551ea7c 1067 continue;
1068 }
1069
8979685e 1070 //
7b580082 1071 // Propagate and update track
8979685e 1072 //
4551ea7c 1073 for (Int_t itime = GetTimeBinsPerPlane()-1; itime >= 0; itime--) {
4551ea7c 1074 Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
8979685e 1075 expectedNumberOfClusters++;
27eaf44b 1076 t.SetNExpected(t.GetNExpected() + 1);
4551ea7c 1077 if (t.GetX() > 345.0) {
27eaf44b 1078 t.SetNExpectedLast(t.GetNExpectedLast() + 1);
4551ea7c 1079 }
1080 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
7e448bcc 1081 AliTRDcluster *cl = 0;
4551ea7c 1082 UInt_t index = 0;
1083 Double_t maxChi2 = fgkMaxChi2;
8979685e 1084 x = timeBin.GetX();
4551ea7c 1085
8979685e 1086 if (timeBin) {
e4f2f73d 1087 if (clusters[ilayer] > 0) {
1088 index = clusters[ilayer];
1089 cl = (AliTRDcluster *)GetCluster(index);
1090 //Double_t h01 = GetTiltFactor(cl); // I.B's fix
1091 //maxChi2=t.GetPredictedChi2(cl,h01); //
1092 }
8979685e 1093
1094 if (cl) {
e4f2f73d 1095 //if (cl->GetNPads() < 5)
1096 Double_t dxsample = timeBin.GetdX();
1097 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
1098 Double_t h01 = GetTiltFactor(cl);
1099 Int_t det = cl->GetDetector();
053767a4 1100 Int_t layer = fGeom->GetLayer(det);
e4f2f73d 1101 if (t.GetX() > 345.0) {
1102 t.SetNLast(t.GetNLast() + 1);
1103 t.SetChi2Last(t.GetChi2Last() + maxChi2);
1104 }
1105 Double_t xcluster = cl->GetX();
1106 t.PropagateTo(xcluster,xx0,xrho);
1107 maxChi2 = t.GetPredictedChi2(cl,h01);
1108
1109 if (maxChi2<1e+10)
053767a4 1110 if (!t.UpdateMI(cl,maxChi2,index,h01,layer)) {
e4f2f73d 1111 if (!t.Update(cl,maxChi2,index,h01)) {
1112 // ????
1113 }
1114 } // else SetCluster(cl, GetNumberOfClusters()-1); // A.Bercuci 25.07.07
43bfc8af 1115
4551ea7c 1116
bcb6fb78 1117// if (calibra->GetMITracking()) {
1118// calibra->UpdateHistograms(cl,&t);
1119// }
77566f2a 1120
e4f2f73d 1121 // Reset material budget if 2 consecutive gold
053767a4 1122 if (layer > 0) {
1123 if ((t.GetTracklets(layer).GetN() + t.GetTracklets(layer-1).GetN()) > 20) {
e4f2f73d 1124 t.SetBudget(2,0.0);
1125 }
1126 }
1127
1128 }
4551ea7c 1129
8979685e 1130 }
4551ea7c 1131
59393e34 1132 }
4551ea7c 1133
1134 ratio0 = ncl / Float_t(fTimeBinsPerPlane);
27eaf44b 1135 Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);
1136 if ((tracklet.GetChi2() < 18.0) &&
1137 (ratio0 > 0.8) &&
1138 (ratio1 > 0.6) &&
1139 (ratio0+ratio1 > 1.5) &&
1140 (t.GetNCross() == 0) &&
1141 (TMath::Abs(t.GetSnp()) < 0.85) &&
1142 (t.GetNumberOfClusters() > 20)){
7e448bcc 1143 //if (ratio0 > 0.8) {
4551ea7c 1144 t.MakeBackupTrack(); // Make backup of the track until is gold
59393e34 1145 }
7b580082 1146
8979685e 1147 }
46d29e70 1148
75bd7f81 1149 return expectedNumberOfClusters;
0c349049 1150
75bd7f81 1151}
1e9bb598 1152
75bd7f81 1153//_____________________________________________________________________________
4551ea7c 1154Int_t AliTRDtracker::PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep)
5443e65e 1155{
75bd7f81 1156 //
0c349049 1157 // Starting from current X-position of track <t> this function
5443e65e 1158 // extrapolates the track up to radial position <xToGo>.
1159 // Returns 1 if track reaches the plane, and 0 otherwise
75bd7f81 1160 //
1161
59393e34 1162 const Double_t kEpsilon = 0.00001;
0c349049 1163
1164 // Current track X-position
4551ea7c 1165 Double_t xpos = t.GetX();
4551ea7c 1166
0c349049 1167 // Direction: inward or outward
1168 Double_t dir = (xpos < xToGo) ? 1.0 : -1.0;
4551ea7c 1169
0c349049 1170 while (((xToGo - xpos) * dir) > kEpsilon) {
4551ea7c 1171
1172 Double_t xyz0[3];
1173 Double_t xyz1[3];
1174 Double_t param[7];
1175 Double_t x;
1176 Double_t y;
1177 Double_t z;
0c349049 1178
1179 // The next step size
1180 Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos),maxStep);
1181
1182 // Get the global position of the starting point
1183 t.GetXYZ(xyz0);
1184
1185 // X-position after next step
4551ea7c 1186 x = xpos + step;
1187
0c349049 1188 // Get local Y and Z at the X-position of the next step
4551ea7c 1189 if (!t.GetProlongation(x,y,z)) {
0c349049 1190 return 0; // No prolongation possible
4551ea7c 1191 }
1192
0c349049 1193 // The global position of the end point of this prolongation step
4551ea7c 1194 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1195 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1196 xyz1[2] = z;
1197
0c349049 1198 // Calculate the mean material budget between start and
1199 // end point of this prolongation step
1200 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1201
1202 // Propagate the track to the X-position after the next step
97ccb02b 1203 if (!t.PropagateTo(x,param[1],param[0]*param[4])) {
4551ea7c 1204 return 0;
1205 }
0c349049 1206
1207 // Rotate the track if necessary
59393e34 1208 AdjustSector(&t);
0c349049 1209
1210 // New track X-position
59393e34 1211 xpos = t.GetX();
4551ea7c 1212
46d29e70 1213 }
75bd7f81 1214
5443e65e 1215 return 1;
5443e65e 1216
59393e34 1217}
5443e65e 1218
46d29e70 1219//_____________________________________________________________________________
c630aafd 1220Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1221{
75bd7f81 1222 //
c630aafd 1223 // Fills clusters into TRD tracking_sectors
1224 // Note that the numbering scheme for the TRD tracking_sectors
1225 // differs from that of TRD sectors
75bd7f81 1226 //
1227
e4f2f73d 1228
1229 if (ReadClusters(fClusters, cTree)) {
0c349049 1230 AliError("Problem with reading the clusters !");
1231 return 1;
c630aafd 1232 }
0c349049 1233 Int_t ncl = fClusters->GetEntriesFast();
4551ea7c 1234 fNclusters = ncl;
1235 AliInfo(Form("Sorting %d clusters",ncl));
c630aafd 1236
1237 UInt_t index;
6c94f330 1238 while (ncl--) {
33744e5d 1239
4551ea7c 1240 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(ncl);
1241 Int_t detector = c->GetDetector();
1242 Int_t localTimeBin = c->GetLocalTimeBin();
1243 Int_t sector = fGeom->GetSector(detector);
053767a4 1244 Int_t layer = fGeom->GetLayer(detector);
439c63c8 1245 Int_t trackingSector = sector;
4551ea7c 1246
0c349049 1247 //if (c->GetQ() > 10) {
053767a4 1248 // Int_t stack = fGeom->GetStack(detector);
0c349049 1249 //}
c630aafd 1250
053767a4 1251 Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(layer,localTimeBin);
4551ea7c 1252 if (gtb < 0) {
1253 continue;
1254 }
053767a4 1255 Int_t trLayer = fTrSec[trackingSector]->GetLayerNumber(gtb);
c630aafd 1256
4551ea7c 1257 index = ncl;
33744e5d 1258
7e448bcc 1259 fHXCl->Fill(c->GetX());
33744e5d 1260
053767a4 1261 fTrSec[trackingSector]->GetLayer(trLayer)->SetX(c->GetX());
1262 fTrSec[trackingSector]->GetLayer(trLayer)->InsertCluster(c,index);
af26ce80 1263
4551ea7c 1264 }
75bd7f81 1265
c630aafd 1266 return 0;
75bd7f81 1267
c630aafd 1268}
1269
1270//_____________________________________________________________________________
b7a0917f 1271void AliTRDtracker::UnloadClusters()
5443e65e 1272{
0a29d0f1 1273 //
5443e65e 1274 // Clears the arrays of clusters and tracks. Resets sectors and timebins
0a29d0f1 1275 //
46d29e70 1276
4551ea7c 1277 Int_t i;
1278 Int_t nentr;
46d29e70 1279
5443e65e 1280 nentr = fClusters->GetEntriesFast();
4551ea7c 1281 for (i = 0; i < nentr; i++) {
1282 delete fClusters->RemoveAt(i);
1283 }
b7a0917f 1284 fNclusters = 0;
46d29e70 1285
5443e65e 1286 nentr = fSeeds->GetEntriesFast();
4551ea7c 1287 for (i = 0; i < nentr; i++) {
1288 delete fSeeds->RemoveAt(i);
1289 }
46d29e70 1290
5443e65e 1291 nentr = fTracks->GetEntriesFast();
4551ea7c 1292 for (i = 0; i < nentr; i++) {
1293 delete fTracks->RemoveAt(i);
1294 }
46d29e70 1295
053767a4 1296 Int_t nsec = AliTRDgeometry::kNsector;
5443e65e 1297 for (i = 0; i < nsec; i++) {
1298 for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1299 fTrSec[i]->GetLayer(pl)->Clear();
1300 }
1301 }
46d29e70 1302
5443e65e 1303}
a819a5f7 1304
75bd7f81 1305//_____________________________________________________________________________
e4f2f73d 1306Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *esd)
7ad19338 1307{
1308 //
1309 // Creates seeds using clusters between position inner plane and outer plane
1310 //
75bd7f81 1311
4551ea7c 1312 const Double_t kMaxTheta = 1.0;
1313 const Double_t kMaxPhi = 2.0;
1314
1315 const Double_t kRoad0y = 6.0; // Road for middle cluster
1316 const Double_t kRoad0z = 8.5; // Road for middle cluster
1317
1318 const Double_t kRoad1y = 2.0; // Road in y for seeded cluster
1319 const Double_t kRoad1z = 20.0; // Road in z for seeded cluster
1320
1321 const Double_t kRoad2y = 3.0; // Road in y for extrapolated cluster
1322 const Double_t kRoad2z = 20.0; // Road in z for extrapolated cluster
c6f438c0 1323 const Int_t kMaxSeed = 3000;
7ad19338 1324
053767a4 1325 Int_t maxSec = AliTRDgeometry::kNsector;
4551ea7c 1326
1327 // Linear fitters in planes
1328 TLinearFitter fitterTC(2,"hyp2"); // Fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
1329 TLinearFitter fitterT2(4,"hyp4"); // Fitting with tilting pads - kz not fixed
69b55c55 1330 fitterTC.StoreData(kTRUE);
1331 fitterT2.StoreData(kTRUE);
4551ea7c 1332 AliRieman rieman(1000); // Rieman fitter
1333 AliRieman rieman2(1000); // Rieman fitter
1334
1335 // Find the maximal and minimal layer for the planes
7ad19338 1336 Int_t layers[6][2];
4551ea7c 1337 AliTRDpropagationLayer *reflayers[6];
1338 for (Int_t i = 0; i < 6; i++) {
1339 layers[i][0] = 10000;
1340 layers[i][1] = 0;
1341 }
1342 for (Int_t ns = 0; ns < maxSec; ns++) {
1343 for (Int_t ilayer = 0; ilayer < fTrSec[ns]->GetNumberOfLayers(); ilayer++) {
053767a4 1344 AliTRDpropagationLayer &propLayer = *(fTrSec[ns]->GetLayer(ilayer));
1345 if (propLayer == 0) {
4551ea7c 1346 continue;
1347 }
053767a4 1348 Int_t det = propLayer[0]->GetDetector();
1349 Int_t layer = fGeom->GetLayer(det);
1350 if (ilayer < layers[layer][0]) {
1351 layers[layer][0] = ilayer;
4551ea7c 1352 }
053767a4 1353 if (ilayer > layers[layer][1]) {
1354 layers[layer][1] = ilayer;
4551ea7c 1355 }
7ad19338 1356 }
1357 }
4551ea7c 1358
f162af62 1359 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0);
6c94f330 1360 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
4551ea7c 1361 Double_t hL[6]; // Tilting angle
1362 Double_t xcl[6]; // X - position of reference cluster
1363 Double_t ycl[6]; // Y - position of reference cluster
1364 Double_t zcl[6]; // Z - position of reference cluster
1365
1366 AliTRDcluster *cl[6] = { 0, 0, 0, 0, 0, 0 }; // Seeding clusters
1367 Float_t padlength[6] = { 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 }; // Current pad-length
1368
1369 Double_t chi2R = 0.0;
1370 Double_t chi2Z = 0.0;
1371 Double_t chi2RF = 0.0;
1372 Double_t chi2ZF = 0.0;
1373
1374 Int_t nclusters; // Total number of clusters
1375 for (Int_t i = 0; i < 6; i++) {
1376 hL[i] = h01;
1377 if (i%2==1) {
1378 hL[i]*=-1.0;
1379 }
1380 }
1381
1382 // Registered seed
c6f438c0 1383 AliTRDseed *pseed = new AliTRDseed[kMaxSeed*6];
1384 AliTRDseed *seed[kMaxSeed];
4551ea7c 1385 for (Int_t iseed = 0; iseed < kMaxSeed; iseed++) {
1386 seed[iseed]= &pseed[iseed*6];
1387 }
69b55c55 1388 AliTRDseed *cseed = seed[0];
4551ea7c 1389
1390 Double_t seedquality[kMaxSeed];
1391 Double_t seedquality2[kMaxSeed];
1392 Double_t seedparams[kMaxSeed][7];
1393 Int_t seedlayer[kMaxSeed];
1394 Int_t registered = 0;
1395 Int_t sort[kMaxSeed];
1396
1397 //
1398 // Seeding part
1399 //
1400 for (Int_t ns = 0; ns < maxSec; ns++) { // Loop over sectors
1401 //for (Int_t ns = 0; ns < 5; ns++) { // Loop over sectors
1402
1403 registered = 0; // Reset registerd seed counter
1404 cseed = seed[registered];
1405 Float_t iter = 0.0;
1406
1407 for (Int_t sLayer = 2; sLayer >= 0; sLayer--) {
1408 //for (Int_t dseed = 5; dseed < 15; dseed += 3) {
1409
1410 iter += 1.0;
1411 Int_t dseed = 5 + Int_t(iter) * 3;
1412
69b55c55 1413 // Initialize seeding layers
4551ea7c 1414 for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
69b55c55 1415 reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
1416 xcl[ilayer] = reflayers[ilayer]->GetX();
4551ea7c 1417 }
1418
1419 Double_t xref = (xcl[sLayer+1] + xcl[sLayer+2]) * 0.5;
1420 AliTRDpropagationLayer &layer0 = *reflayers[sLayer+0];
1421 AliTRDpropagationLayer &layer1 = *reflayers[sLayer+1];
1422 AliTRDpropagationLayer &layer2 = *reflayers[sLayer+2];
1423 AliTRDpropagationLayer &layer3 = *reflayers[sLayer+3];
1424
1425 Int_t maxn3 = layer3;
1426 for (Int_t icl3 = 0; icl3 < maxn3; icl3++) {
33744e5d 1427
69b55c55 1428 AliTRDcluster *cl3 = layer3[icl3];
4551ea7c 1429 if (!cl3) {
1430 continue;
1431 }
1432 padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2() * 12.0);
69b55c55 1433 ycl[sLayer+3] = cl3->GetY();
1434 zcl[sLayer+3] = cl3->GetZ();
4551ea7c 1435 Float_t yymin0 = ycl[sLayer+3] - 1.0 - kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1436 Float_t yymax0 = ycl[sLayer+3] + 1.0 + kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1437 Int_t maxn0 = layer0;
1438
1439 for (Int_t icl0 = layer0.Find(yymin0); icl0 < maxn0; icl0++) {
1440
69b55c55 1441 AliTRDcluster *cl0 = layer0[icl0];
4551ea7c 1442 if (!cl0) {
1443 continue;
1444 }
1445 if (cl3->IsUsed() && cl0->IsUsed()) {
1446 continue;
1447 }
69b55c55 1448 ycl[sLayer+0] = cl0->GetY();
1449 zcl[sLayer+0] = cl0->GetZ();
4551ea7c 1450 if (ycl[sLayer+0] > yymax0) {
1451 break;
1452 }
1453 Double_t tanphi = (ycl[sLayer+3]-ycl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1454 if (TMath::Abs(tanphi) > kMaxPhi) {
1455 continue;
1456 }
1457 Double_t tantheta = (zcl[sLayer+3]-zcl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1458 if (TMath::Abs(tantheta) > kMaxTheta) {
1459 continue;
1460 }
1461 padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2() * 12.0);
1462
1463 // Expected position in 1 layer
1464 Double_t y1exp = ycl[sLayer+0] + (tanphi) * (xcl[sLayer+1]-xcl[sLayer+0]);
1465 Double_t z1exp = zcl[sLayer+0] + (tantheta) * (xcl[sLayer+1]-xcl[sLayer+0]);
1466 Float_t yymin1 = y1exp - kRoad0y - tanphi;
1467 Float_t yymax1 = y1exp + kRoad0y + tanphi;
1468 Int_t maxn1 = layer1;
1469
1470 for (Int_t icl1 = layer1.Find(yymin1); icl1 < maxn1; icl1++) {
1471
69b55c55 1472 AliTRDcluster *cl1 = layer1[icl1];
4551ea7c 1473 if (!cl1) {
1474 continue;
1475 }
69b55c55 1476 Int_t nusedCl = 0;
1477 if (cl3->IsUsed()) nusedCl++;
1478 if (cl0->IsUsed()) nusedCl++;
1479 if (cl1->IsUsed()) nusedCl++;
4551ea7c 1480 if (nusedCl > 1) {
1481 continue;
1482 }
69b55c55 1483 ycl[sLayer+1] = cl1->GetY();
1484 zcl[sLayer+1] = cl1->GetZ();
4551ea7c 1485 if (ycl[sLayer+1] > yymax1) {
1486 break;
1487 }
1488 if (TMath::Abs(ycl[sLayer+1]-y1exp) > kRoad0y+tanphi) {
1489 continue;
1490 }
1491 if (TMath::Abs(zcl[sLayer+1]-z1exp) > kRoad0z) {
1492 continue;
1493 }
1494 padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2() * 12.0);
1495
1496 Double_t y2exp = ycl[sLayer+0]+(tanphi) * (xcl[sLayer+2]-xcl[sLayer+0]) + (ycl[sLayer+1]-y1exp);
1497 Double_t z2exp = zcl[sLayer+0]+(tantheta) * (xcl[sLayer+2]-xcl[sLayer+0]);
1498 Int_t index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y,kRoad1z);
1499 if (index2 <= 0) {
1500 continue;
1501 }
1502 AliTRDcluster *cl2 = (AliTRDcluster *) GetCluster(index2);
1503 padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2() * 12.0);
1504 ycl[sLayer+2] = cl2->GetY();
1505 zcl[sLayer+2] = cl2->GetZ();
1506 if (TMath::Abs(cl2->GetZ()-z2exp) > kRoad0z) {
1507 continue;
1508 }
1509
69b55c55 1510 rieman.Reset();
1511 rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
1512 rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
1513 rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);
1514 rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
1515 rieman.Update();
4551ea7c 1516
1517 // Reset fitter
1518 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
69b55c55 1519 cseed[iLayer].Reset();
1520 }
4551ea7c 1521 chi2Z = 0.0;
1522 chi2R = 0.0;
1523
1524 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
ca1e1e5b 1525 cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1526 chi2Z += (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer])
1527 * (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer]);
1528 cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
1529 cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1530 chi2R += (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer])
1531 * (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer]);
1532 cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
69b55c55 1533 }
4551ea7c 1534 if (TMath::Sqrt(chi2R) > 1.0/iter) {
1535 continue;
69b55c55 1536 }
4551ea7c 1537 if (TMath::Sqrt(chi2Z) > 7.0/iter) {
1538 continue;
1539 }
1540
1541 Float_t minmax[2] = { -100.0, 100.0 };
1542 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1543 Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer] * 0.5
ca1e1e5b 1544 + 1.0 - cseed[sLayer+iLayer].GetZref(0);
4551ea7c 1545 if (max < minmax[1]) {
1546 minmax[1] = max;
1547 }
1548 Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer] * 0.5
ca1e1e5b 1549 - 1.0 - cseed[sLayer+iLayer].GetZref(0);
4551ea7c 1550 if (min > minmax[0]) {
1551 minmax[0] = min;
1552 }
1553 }
1554
69b55c55 1555 Bool_t isFake = kFALSE;
4551ea7c 1556 if (cl0->GetLabel(0) != cl3->GetLabel(0)) {
1557 isFake = kTRUE;
1558 }
1559 if (cl1->GetLabel(0) != cl3->GetLabel(0)) {
1560 isFake = kTRUE;
1561 }
1562 if (cl2->GetLabel(0) != cl3->GetLabel(0)) {
1563 isFake = kTRUE;
1564 }
1565
3a039a31 1566 if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0) {
4551ea7c 1567 if ((!isFake) || ((icl3%10) == 0)) { // Debugging print
1568 TTreeSRedirector &cstream = *fDebugStreamer;
1569 cstream << "Seeds0"
1570 << "isFake=" << isFake
1571 << "Cl0.=" << cl0
1572 << "Cl1.=" << cl1
1573 << "Cl2.=" << cl2
1574 << "Cl3.=" << cl3
1575 << "Xref=" << xref
1576 << "X0=" << xcl[sLayer+0]
1577 << "X1=" << xcl[sLayer+1]
1578 << "X2=" << xcl[sLayer+2]
1579 << "X3=" << xcl[sLayer+3]
1580 << "Y2exp=" << y2exp
1581 << "Z2exp=" << z2exp
1582 << "Chi2R=" << chi2R
1583 << "Chi2Z=" << chi2Z
1584 << "Seed0.=" << &cseed[sLayer+0]
1585 << "Seed1.=" << &cseed[sLayer+1]
1586 << "Seed2.=" << &cseed[sLayer+2]
1587 << "Seed3.=" << &cseed[sLayer+3]
1588 << "Zmin=" << minmax[0]
1589 << "Zmax=" << minmax[1]
1590 << "\n";
d337ef8d 1591 }
69b55c55 1592 }
33744e5d 1593
4551ea7c 1594 ////////////////////////////////////////////////////////////////////////////////////
33744e5d 1595 //
4551ea7c 1596 // Fit seeding part
33744e5d 1597 //
4551ea7c 1598 ////////////////////////////////////////////////////////////////////////////////////
33744e5d 1599
69b55c55 1600 cl[sLayer+0] = cl0;
1601 cl[sLayer+1] = cl1;
1602 cl[sLayer+2] = cl2;
1603 cl[sLayer+3] = cl3;
4551ea7c 1604 Bool_t isOK = kTRUE;
1605
1606 for (Int_t jLayer = 0; jLayer < 4; jLayer++) {
1607
ca1e1e5b 1608 cseed[sLayer+jLayer].SetTilt(hL[sLayer+jLayer]);
1609 cseed[sLayer+jLayer].SetPadLength(padlength[sLayer+jLayer]);
1610 cseed[sLayer+jLayer].SetX0(xcl[sLayer+jLayer]);
4551ea7c 1611
7bce990c 1612 for (Int_t jter = 0; jter < 2; jter++) {
4551ea7c 1613
69b55c55 1614 //
4551ea7c 1615 // In iteration 0 we try only one pad-row
1616 // If quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
69b55c55 1617 //
1618 AliTRDseed tseed = cseed[sLayer+jLayer];
4551ea7c 1619 Float_t roadz = padlength[sLayer+jLayer] * 0.5;
7bce990c 1620 if (jter > 0) {
4551ea7c 1621 roadz = padlength[sLayer+jLayer];
1622 }
1623
1624 Float_t quality = 10000.0;
1625
1626 for (Int_t iTime = 2; iTime < 20; iTime++) {
1627
1628 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
1629 Double_t dxlayer = layer.GetX() - xcl[sLayer+jLayer];
1630 Double_t zexp = cl[sLayer+jLayer]->GetZ();
1631
7bce990c 1632 if (jter > 0) {
4551ea7c 1633 // Try 2 pad-rows in second iteration
ca1e1e5b 1634 zexp = tseed.GetZref(0) + tseed.GetZref(1) * dxlayer;
4551ea7c 1635 if (zexp > cl[sLayer+jLayer]->GetZ()) {
1636 zexp = cl[sLayer+jLayer]->GetZ() + padlength[sLayer+jLayer]*0.5;
1637 }
1638 if (zexp < cl[sLayer+jLayer]->GetZ()) {
1639 zexp = cl[sLayer+jLayer]->GetZ() - padlength[sLayer+jLayer]*0.5;
1640 }
1641 }
1642
ca1e1e5b 1643 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
4551ea7c 1644 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
1645 if (index <= 0) {
1646 continue;
69b55c55 1647 }
7bce990c 1648 AliTRDcluster *clu = (AliTRDcluster *) GetCluster(index);
4551ea7c 1649
ca1e1e5b 1650 tseed.SetIndexes(iTime,index);
7bce990c 1651 tseed.SetClusters(iTime,clu); // Register cluster
1652 tseed.SetX(iTime,dxlayer); // Register cluster
1653 tseed.SetY(iTime,clu->GetY()); // Register cluster
1654 tseed.SetZ(iTime,clu->GetZ()); // Register cluster
4551ea7c 1655
1656 }
1657
69b55c55 1658 tseed.Update();
4551ea7c 1659
1660 // Count the number of clusters and distortions into quality
ca1e1e5b 1661 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1662 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
1663 + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
1664 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
7bce990c 1665 if ((jter == 0) && tseed.IsOK()) {
69b55c55 1666 cseed[sLayer+jLayer] = tseed;
4551ea7c 1667 quality = tquality;
1668 if (tquality < 5) {
1669 break;
1670 }
69b55c55 1671 }
4551ea7c 1672 if (tseed.IsOK() && (tquality < quality)) {
1673 cseed[sLayer+jLayer] = tseed;
1674 }
1675
7bce990c 1676 } // Loop: jter
4551ea7c 1677
1678 if (!cseed[sLayer+jLayer].IsOK()) {
69b55c55 1679 isOK = kFALSE;
1680 break;
4551ea7c 1681 }
1682
69b55c55 1683 cseed[sLayer+jLayer].CookLabels();
1684 cseed[sLayer+jLayer].UpdateUsed();
ca1e1e5b 1685 nusedCl += cseed[sLayer+jLayer].GetNUsed();
4551ea7c 1686 if (nusedCl > 25) {
69b55c55 1687 isOK = kFALSE;
1688 break;
4551ea7c 1689 }
1690
1691 } // Loop: jLayer
1692
1693 if (!isOK) {
1694 continue;
69b55c55 1695 }
4551ea7c 1696 nclusters = 0;
1697 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1698 if (cseed[sLayer+iLayer].IsOK()) {
ca1e1e5b 1699 nclusters += cseed[sLayer+iLayer].GetN2();
69b55c55 1700 }
1701 }
4551ea7c 1702
1703 // Iteration 0
69b55c55 1704 rieman.Reset();
4551ea7c 1705 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1706 rieman.AddPoint(xcl[sLayer+iLayer]
ca1e1e5b 1707 ,cseed[sLayer+iLayer].GetYfitR(0)
1708 ,cseed[sLayer+iLayer].GetZProb()
4551ea7c 1709 ,1
1710 ,10);
69b55c55 1711 }
1712 rieman.Update();
4551ea7c 1713
1714 chi2R = 0.0;
1715 chi2Z = 0.0;
1716
1717 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
ca1e1e5b 1718 cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1719 chi2R += (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0))
1720 * (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0));
1721 cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
1722 cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1723 chi2Z += (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz())
1724 * (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz());
1725 cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
69b55c55 1726 }
1727 Double_t curv = rieman.GetC();
4551ea7c 1728
69b55c55 1729 //
4551ea7c 1730 // Likelihoods
6c94f330 1731 //
ca1e1e5b 1732 Double_t sumda = TMath::Abs(cseed[sLayer+0].GetYfitR(1) - cseed[sLayer+0].GetYref(1))
1733 + TMath::Abs(cseed[sLayer+1].GetYfitR(1) - cseed[sLayer+1].GetYref(1))
1734 + TMath::Abs(cseed[sLayer+2].GetYfitR(1) - cseed[sLayer+2].GetYref(1))
1735 + TMath::Abs(cseed[sLayer+3].GetYfitR(1) - cseed[sLayer+3].GetYref(1));
4551ea7c 1736 Double_t likea = TMath::Exp(-sumda*10.6);
1737 Double_t likechi2 = 0.0000000001;
1738 if (chi2R < 0.5) {
1739 likechi2 += TMath::Exp(-TMath::Sqrt(chi2R) * 7.73);
1740 }
1741 Double_t likechi2z = TMath::Exp(-chi2Z * 0.088) / TMath::Exp(-chi2Z * 0.019);
1742 Double_t likeN = TMath::Exp(-(72 - nclusters) * 0.19);
1743 Double_t like = likea * likechi2 * likechi2z * likeN;
ca1e1e5b 1744 Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetYref(1) - 130.0*curv) * 1.9);
1745 Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetZref(1)
1746 - cseed[sLayer+0].GetZref(0) / xcl[sLayer+0]) * 5.9);
6c94f330 1747 Double_t likePrim = TMath::Max(likePrimY*likePrimZ,0.0005);
69b55c55 1748
4551ea7c 1749 seedquality[registered] = like;
1750 seedlayer[registered] = sLayer;
1751 if (TMath::Log(0.000000000000001 + like) < -15) {
1752 continue;
1753 }
69b55c55 1754 AliTRDseed seedb[6];
4551ea7c 1755 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
69b55c55 1756 seedb[iLayer] = cseed[iLayer];
1757 }
4551ea7c 1758
1759 ////////////////////////////////////////////////////////////////////////////////////
69b55c55 1760 //
4551ea7c 1761 // Full track fit part
33744e5d 1762 //
4551ea7c 1763 ////////////////////////////////////////////////////////////////////////////////////
1764
1765 Int_t nlayers = 0;
1766 Int_t nusedf = 0;
1767 Int_t findable = 0;
1768
69b55c55 1769 //
4551ea7c 1770 // Add new layers - avoid long extrapolation
69b55c55 1771 //
4551ea7c 1772 Int_t tLayer[2] = { 0, 0 };
1773 if (sLayer == 2) {
1774 tLayer[0] = 1;
1775 tLayer[1] = 0;
1776 }
1777 if (sLayer == 1) {
1778 tLayer[0] = 5;
1779 tLayer[1] = 0;
1780 }
1781 if (sLayer == 0) {
1782 tLayer[0] = 4;
1783 tLayer[1] = 5;
1784 }
1785
1786 for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
1787 Int_t jLayer = tLayer[iLayer]; // Set tracking layer
69b55c55 1788 cseed[jLayer].Reset();
ca1e1e5b 1789 cseed[jLayer].SetTilt(hL[jLayer]);
1790 cseed[jLayer].SetPadLength(padlength[jLayer]);
1791 cseed[jLayer].SetX0(xcl[jLayer]);
4551ea7c 1792 // Get pad length and rough cluster
ca1e1e5b 1793 Int_t indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].GetYref(0)
1794 ,cseed[jLayer].GetZref(0)
4551ea7c 1795 ,kRoad2y
1796 ,kRoad2z);
1797 if (indexdummy <= 0) {
1798 continue;
1799 }
1800 AliTRDcluster *cldummy = (AliTRDcluster *) GetCluster(indexdummy);
1801 padlength[jLayer] = TMath::Sqrt(cldummy->GetSigmaZ2() * 12.0);
69b55c55 1802 }
4551ea7c 1803 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
1804
1805 for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
1806
1807 Int_t jLayer = tLayer[iLayer]; // set tracking layer
1808 if ((jLayer == 0) && !(cseed[1].IsOK())) {
1809 continue; // break not allowed
1810 }
1811 if ((jLayer == 5) && !(cseed[4].IsOK())) {
1812 continue; // break not allowed
1813 }
ca1e1e5b 1814 Float_t zexp = cseed[jLayer].GetZref(0);
4551ea7c 1815 Double_t zroad = padlength[jLayer] * 0.5 + 1.0;
1816
7bce990c 1817 for (Int_t jter = 0; jter < 2; jter++) {
4551ea7c 1818
69b55c55 1819 AliTRDseed tseed = cseed[jLayer];
4551ea7c 1820 Float_t quality = 10000.0;
1821
1822 for (Int_t iTime = 2; iTime < 20; iTime++) {
1823 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
1824 Double_t dxlayer = layer.GetX()-xcl[jLayer];
ca1e1e5b 1825 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
4551ea7c 1826 Float_t yroad = kRoad1y;
1827 Int_t index = layer.FindNearestCluster(yexp,zexp,yroad,zroad);
1828 if (index <= 0) {
1829 continue;
1830 }
7bce990c 1831 AliTRDcluster *clu = (AliTRDcluster *) GetCluster(index);
ca1e1e5b 1832 tseed.SetIndexes(iTime,index);
7bce990c 1833 tseed.SetClusters(iTime,clu); // Register cluster
1834 tseed.SetX(iTime,dxlayer); // Register cluster
1835 tseed.SetY(iTime,clu->GetY()); // Register cluster
1836 tseed.SetZ(iTime,clu->GetZ()); // Register cluster
4551ea7c 1837 }
1838
69b55c55 1839 tseed.Update();
4551ea7c 1840 if (tseed.IsOK()) {
ca1e1e5b 1841 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1842 Float_t tquality = (18.0 - tseed.GetN2())/2.0 + TMath::Abs(dangle) / 0.1
1843 + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
1844 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 1845 if (tquality < quality) {
1846 cseed[jLayer] = tseed;
1847 quality = tquality;
69b55c55 1848 }
1849 }
4551ea7c 1850
1851 zroad *= 2.0;
1852
7bce990c 1853 } // Loop: jter
4551ea7c 1854
1855 if ( cseed[jLayer].IsOK()) {
69b55c55 1856 cseed[jLayer].CookLabels();
1857 cseed[jLayer].UpdateUsed();
ca1e1e5b 1858 nusedf += cseed[jLayer].GetNUsed();
4551ea7c 1859 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
69b55c55 1860 }
4551ea7c 1861
1862 } // Loop: iLayer
1863
1864 // Make copy
69b55c55 1865 AliTRDseed bseed[6];
4551ea7c 1866 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
69b55c55 1867 bseed[jLayer] = cseed[jLayer];
4551ea7c 1868 }
1869 Float_t lastquality = 10000.0;
1870 Float_t lastchi2 = 10000.0;
1871 Float_t chi2 = 1000.0;
1872
7bce990c 1873 for (Int_t jter = 0; jter < 4; jter++) {
4551ea7c 1874
1875 // Sort tracklets according "quality", try to "improve" 4 worst
1876 Float_t sumquality = 0.0;
69b55c55 1877 Float_t squality[6];
1878 Int_t sortindexes[6];
4551ea7c 1879
1880 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1881 if (bseed[jLayer].IsOK()) {
69b55c55 1882 AliTRDseed &tseed = bseed[jLayer];
ca1e1e5b 1883 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
1884 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1885 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
1886 + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
1887 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 1888 squality[jLayer] = tquality;
1889 }
1890 else {
1891 squality[jLayer] = -1.0;
ad670fba 1892 }
6c94f330 1893 sumquality +=squality[jLayer];
69b55c55 1894 }
1895
4551ea7c 1896 if ((sumquality >= lastquality) ||
1897 (chi2 > lastchi2)) {
1898 break;
1899 }
69b55c55 1900 lastquality = sumquality;
1901 lastchi2 = chi2;
7bce990c 1902 if (jter > 0) {
4551ea7c 1903 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
69b55c55 1904 cseed[jLayer] = bseed[jLayer];
1905 }
1906 }
1907 TMath::Sort(6,squality,sortindexes,kFALSE);
4551ea7c 1908
1909 for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
1910
6c94f330 1911 Int_t bLayer = sortindexes[jLayer];
1912 AliTRDseed tseed = bseed[bLayer];
4551ea7c 1913
1914 for (Int_t iTime = 2; iTime < 20; iTime++) {
1915
1916 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
1917 Double_t dxlayer = layer.GetX() - xcl[bLayer];
ca1e1e5b 1918 Double_t zexp = tseed.GetZref(0);
1919 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
4551ea7c 1920 Float_t roadz = padlength[bLayer] + 1;
ca1e1e5b 1921 if (TMath::Abs(tseed.GetZProb() - zexp) > 0.5*padlength[bLayer]) {
4551ea7c 1922 roadz = padlength[bLayer] * 0.5;
1923 }
ca1e1e5b 1924 if (tseed.GetZfit(1)*tseed.GetZref(1) < 0.0) {
4551ea7c 1925 roadz = padlength[bLayer] * 0.5;
1926 }
ca1e1e5b 1927 if (TMath::Abs(tseed.GetZProb() - zexp) < 0.1*padlength[bLayer]) {
1928 zexp = tseed.GetZProb();
4551ea7c 1929 roadz = padlength[bLayer] * 0.5;
1930 }
1931
ca1e1e5b 1932 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer - zcor;
4551ea7c 1933 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
1934 if (index <= 0) {
1935 continue;
69b55c55 1936 }
7bce990c 1937 AliTRDcluster *clu = (AliTRDcluster *) GetCluster(index);
4551ea7c 1938
ca1e1e5b 1939 tseed.SetIndexes(iTime,index);
7bce990c 1940 tseed.SetClusters(iTime,clu); // Register cluster
1941 tseed.SetX(iTime,dxlayer); // Register cluster
1942 tseed.SetY(iTime,clu->GetY()); // Register cluster
1943 tseed.SetZ(iTime,clu->GetZ()); // Register cluster
4551ea7c 1944
1945 }
1946
69b55c55 1947 tseed.Update();
c6f438c0 1948 if (tseed.IsOK()) {
ca1e1e5b 1949 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1950 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
1951 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0
1952 + TMath::Abs(dangle) / 0.1
1953 + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
1954 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 1955 if (tquality<squality[bLayer]) {
69b55c55 1956 bseed[bLayer] = tseed;
4551ea7c 1957 }
69b55c55 1958 }
4551ea7c 1959
1960 } // Loop: jLayer
1961
1962 chi2 = AliTRDseed::FitRiemanTilt(bseed,kTRUE);
1963
7bce990c 1964 } // Loop: jter
4551ea7c 1965
1966 nclusters = 0;
1967 nlayers = 0;
1968 findable = 0;
1969 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ca1e1e5b 1970 if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) {
69b55c55 1971 findable++;
4551ea7c 1972 }
1973 if (cseed[iLayer].IsOK()) {
ca1e1e5b 1974 nclusters += cseed[iLayer].GetN2();
69b55c55 1975 nlayers++;
1976 }
1977 }
4551ea7c 1978 if (nlayers < 3) {
1979 continue;
1980 }
69b55c55 1981 rieman.Reset();
4551ea7c 1982 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1983 if (cseed[iLayer].IsOK()) {
1984 rieman.AddPoint(xcl[iLayer]
ca1e1e5b 1985 ,cseed[iLayer].GetYfitR(0)
1986 ,cseed[iLayer].GetZProb()
4551ea7c 1987 ,1
1988 ,10);
1989 }
69b55c55 1990 }
1991 rieman.Update();
4551ea7c 1992
1993 chi2RF = 0.0;
1994 chi2ZF = 0.0;
1995 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1996 if (cseed[iLayer].IsOK()) {
ca1e1e5b 1997 cseed[iLayer].SetYref(0,rieman.GetYat(xcl[iLayer]));
1998 chi2RF += (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0))
1999 * (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0));
2000 cseed[iLayer].SetYref(1,rieman.GetDYat(xcl[iLayer]));
2001 cseed[iLayer].SetZref(0,rieman.GetZat(xcl[iLayer]));
2002 chi2ZF += (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz())
2003 * (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz());
2004 cseed[iLayer].SetZref(1,rieman.GetDZat(xcl[iLayer]));
69b55c55 2005 }
2006 }
4551ea7c 2007 chi2RF /= TMath::Max((nlayers - 3.0),1.0);
2008 chi2ZF /= TMath::Max((nlayers - 3.0),1.0);
2009 curv = rieman.GetC();
2010
2011 Double_t xref2 = (xcl[2] + xcl[3]) * 0.5; // Middle of the chamber
2012 Double_t dzmf = rieman.GetDZat(xref2);
2013 Double_t zmf = rieman.GetZat(xref2);
69b55c55 2014
69b55c55 2015 //
4551ea7c 2016 // Fit hyperplane
69b55c55 2017 //
4551ea7c 2018 Int_t npointsT = 0;
69b55c55 2019 fitterTC.ClearPoints();
2020 fitterT2.ClearPoints();
2021 rieman2.Reset();
4551ea7c 2022
2023 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2024
2025 if (!cseed[iLayer].IsOK()) {
2026 continue;
2027 }
2028
2029 for (Int_t itime = 0; itime < 25; itime++) {
2030
ca1e1e5b 2031 if (!cseed[iLayer].IsUsable(itime)) {
4551ea7c 2032 continue;
2033 }
2034 // X relative to the middle chamber
ca1e1e5b 2035 Double_t x = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;
2036 Double_t y = cseed[iLayer].GetY(itime);
2037 Double_t z = cseed[iLayer].GetZ(itime);
69b55c55 2038 // ExB correction to the correction
4551ea7c 2039 // Tilted rieman
69b55c55 2040 Double_t uvt[6];
4551ea7c 2041 // Global x
ca1e1e5b 2042 Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();
4551ea7c 2043 Double_t t = 1.0 / (x2*x2 + y*y);
2044 uvt[1] = t; // t
2045 uvt[0] = 2.0 * x2 * uvt[1]; // u
2046 uvt[2] = 2.0 * hL[iLayer] * uvt[1];
2047 uvt[3] = 2.0 * hL[iLayer] * x * uvt[1];
2048 uvt[4] = 2.0 * (y + hL[iLayer]*z) * uvt[1];
2049 Double_t error = 2.0 * 0.2 * uvt[1];
69b55c55 2050 fitterT2.AddPoint(uvt,uvt[4],error);
4551ea7c 2051
69b55c55 2052 //
4551ea7c 2053 // Constrained rieman
69b55c55 2054 //
ca1e1e5b 2055 z = cseed[iLayer].GetZ(itime);
4551ea7c 2056 uvt[0] = 2.0 * x2 * t; // u
2057 uvt[1] = 2.0 * hL[iLayer] * x2 * uvt[1];
2058 uvt[2] = 2.0 * (y + hL[iLayer] * (z - GetZ())) * t;
69b55c55 2059 fitterTC.AddPoint(uvt,uvt[2],error);
69b55c55 2060 rieman2.AddPoint(x2,y,z,1,10);
2061 npointsT++;
4551ea7c 2062
69b55c55 2063 }
4551ea7c 2064
2065 } // Loop: iLayer
2066
69b55c55 2067 rieman2.Update();
2068 fitterTC.Eval();
2069 fitterT2.Eval();
2070 Double_t rpolz0 = fitterT2.GetParameter(3);
2071 Double_t rpolz1 = fitterT2.GetParameter(4);
4551ea7c 2072
69b55c55 2073 //
4551ea7c 2074 // Linear fitter - not possible to make boundaries
2075 // Do not accept non possible z and dzdx combinations
69b55c55 2076 //
4551ea7c 2077 Bool_t acceptablez = kTRUE;
2078 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2079 if (cseed[iLayer].IsOK()) {
2080 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
ca1e1e5b 2081 if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > padlength[iLayer] * 0.5 + 1.0) {
69b55c55 2082 acceptablez = kFALSE;
4551ea7c 2083 }
69b55c55 2084 }
2085 }
4551ea7c 2086 if (!acceptablez) {
69b55c55 2087 fitterT2.FixParameter(3,zmf);
2088 fitterT2.FixParameter(4,dzmf);
2089 fitterT2.Eval();
2090 fitterT2.ReleaseParameter(3);
2091 fitterT2.ReleaseParameter(4);
2092 rpolz0 = fitterT2.GetParameter(3);
2093 rpolz1 = fitterT2.GetParameter(4);
2094 }
4551ea7c 2095
2096 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
2097 Double_t chi2TC = fitterTC.GetChisquare() / Float_t(npointsT);
69b55c55 2098 Double_t polz1c = fitterTC.GetParameter(2);
4551ea7c 2099 Double_t polz0c = polz1c * xref2;
6c94f330 2100 Double_t aC = fitterTC.GetParameter(0);
2101 Double_t bC = fitterTC.GetParameter(1);
4551ea7c 2102 Double_t cC = aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
6c94f330 2103 Double_t aR = fitterT2.GetParameter(0);
2104 Double_t bR = fitterT2.GetParameter(1);
2105 Double_t dR = fitterT2.GetParameter(2);
4551ea7c 2106 Double_t cR = 1.0 + bR*bR - dR*aR;
2107 Double_t dca = 0.0;
2108 if (cR > 0.0) {
2109 dca = -dR / (TMath::Sqrt(1.0 + bR*bR - dR*aR) + TMath::Sqrt(1.0 + bR*bR));
2110 cR = aR / TMath::Sqrt(cR);
69b55c55 2111 }
4551ea7c 2112
2113 Double_t chi2ZT2 = 0.0;
2114 Double_t chi2ZTC = 0.0;
2115 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2116 if (cseed[iLayer].IsOK()) {
2117 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
2118 Double_t zTC = polz0c + polz1c * (xcl[iLayer] - xref2);
ca1e1e5b 2119 chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
2120 chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
69b55c55 2121 }
2122 }
4551ea7c 2123 chi2ZT2 /= TMath::Max((nlayers - 3.0),1.0);
2124 chi2ZTC /= TMath::Max((nlayers - 3.0),1.0);
2125
2126 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
2127 Float_t sumdaf = 0.0;
2128 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2129 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2130 sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))
2131 / cseed[iLayer].GetSigmaY2());
4551ea7c 2132 }
2133 }
2134 sumdaf /= Float_t (nlayers - 2.0);
2135
69b55c55 2136 //
4551ea7c 2137 // Likelihoods for full track
69b55c55 2138 //
4551ea7c 2139 Double_t likezf = TMath::Exp(-chi2ZF * 0.14);
2140 Double_t likechi2C = TMath::Exp(-chi2TC * 0.677);
2141 Double_t likechi2TR = TMath::Exp(-chi2TR * 0.78);
2142 Double_t likeaf = TMath::Exp(-sumdaf * 3.23);
2143 seedquality2[registered] = likezf * likechi2TR * likeaf;
33744e5d 2144
4551ea7c 2145 // Still needed ????
69b55c55 2146// Bool_t isGold = kFALSE;
2147//
6c94f330 2148// if (nlayers == 6 && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE; // gold
2149// if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE; // gold
69b55c55 2150// if (isGold &&nusedf<10){
2151// for (Int_t jLayer=0;jLayer<6;jLayer++){
6c94f330 2152// if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
69b55c55 2153// seed[index][jLayer].UseClusters(); //sign gold
2154// }
2155// }
4551ea7c 2156
2157 Int_t index0 = 0;
2158 if (!cseed[0].IsOK()) {
69b55c55 2159 index0 = 1;
4551ea7c 2160 if (!cseed[1].IsOK()) {
2161 index0 = 2;
2162 }
69b55c55 2163 }
ca1e1e5b 2164 seedparams[registered][0] = cseed[index0].GetX0();
2165 seedparams[registered][1] = cseed[index0].GetYref(0);
2166 seedparams[registered][2] = cseed[index0].GetZref(0);
c6f438c0 2167 seedparams[registered][5] = cR;
ca1e1e5b 2168 seedparams[registered][3] = cseed[index0].GetX0() * cR - TMath::Sin(TMath::ATan(cseed[0].GetYref(1)));
2169 seedparams[registered][4] = cseed[index0].GetZref(1)
2170 / TMath::Sqrt(1.0 + cseed[index0].GetYref(1) * cseed[index0].GetYref(1));
69b55c55 2171 seedparams[registered][6] = ns;
4551ea7c 2172
2173 Int_t labels[12];
2174 Int_t outlab[24];
2175 Int_t nlab = 0;
2176 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2177 if (!cseed[iLayer].IsOK()) {
2178 continue;
2179 }
ca1e1e5b 2180 if (cseed[iLayer].GetLabels(0) >= 0) {
2181 labels[nlab] = cseed[iLayer].GetLabels(0);
69b55c55 2182 nlab++;
2183 }
ca1e1e5b 2184 if (cseed[iLayer].GetLabels(1) >= 0) {
2185 labels[nlab] = cseed[iLayer].GetLabels(1);
69b55c55 2186 nlab++;
4551ea7c 2187 }
69b55c55 2188 }
2189 Freq(nlab,labels,outlab,kFALSE);
4551ea7c 2190 Int_t label = outlab[0];
2191 Int_t frequency = outlab[1];
2192 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ca1e1e5b 2193 cseed[iLayer].SetFreq(frequency);
2194 cseed[iLayer].SetC(cR);
2195 cseed[iLayer].SetCC(cC);
2196 cseed[iLayer].SetChi2(chi2TR);
2197 cseed[iLayer].SetChi2Z(chi2ZF);
69b55c55 2198 }
33744e5d 2199
2200 // Debugging print
4551ea7c 2201 if (1 || (!isFake)) {
69b55c55 2202 Float_t zvertex = GetZ();
4551ea7c 2203 TTreeSRedirector &cstream = *fDebugStreamer;
3a039a31 2204 if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0) {
4551ea7c 2205 cstream << "Seeds1"
2206 << "isFake=" << isFake
2207 << "Vertex=" << zvertex
2208 << "Rieman2.=" << &rieman2
2209 << "Rieman.=" << &rieman
2210 << "Xref=" << xref
2211 << "X0=" << xcl[0]
2212 << "X1=" << xcl[1]
2213 << "X2=" << xcl[2]
2214 << "X3=" << xcl[3]
2215 << "X4=" << xcl[4]
2216 << "X5=" << xcl[5]
2217 << "Chi2R=" << chi2R
2218 << "Chi2Z=" << chi2Z
2219 << "Chi2RF=" << chi2RF // Chi2 of trackletes on full track
2220 << "Chi2ZF=" << chi2ZF // Chi2 z on tracklets on full track
2221 << "Chi2ZT2=" << chi2ZT2 // Chi2 z on tracklets on full track - rieman tilt
2222 << "Chi2ZTC=" << chi2ZTC // Chi2 z on tracklets on full track - rieman tilt const
2223 << "Chi2TR=" << chi2TR // Chi2 without vertex constrain
2224 << "Chi2TC=" << chi2TC // Chi2 with vertex constrain
2225 << "C=" << curv // Non constrained - no tilt correction
2226 << "DR=" << dR // DR parameter - tilt correction
2227 << "DCA=" << dca // DCA - tilt correction
2228 << "CR=" << cR // Non constrained curvature - tilt correction
2229 << "CC=" << cC // Constrained curvature
2230 << "Polz0=" << polz0c
2231 << "Polz1=" << polz1c
2232 << "RPolz0=" << rpolz0
2233 << "RPolz1=" << rpolz1
2234 << "Ncl=" << nclusters
2235 << "Nlayers=" << nlayers
2236 << "NUsedS=" << nusedCl
2237 << "NUsed=" << nusedf
2238 << "Findable=" << findable
2239 << "Like=" << like
2240 << "LikePrim=" << likePrim
2241 << "Likechi2C=" << likechi2C
2242 << "Likechi2TR=" << likechi2TR
2243 << "Likezf=" << likezf
2244 << "LikeF=" << seedquality2[registered]
2245 << "S0.=" << &cseed[0]
2246 << "S1.=" << &cseed[1]
2247 << "S2.=" << &cseed[2]
2248 << "S3.=" << &cseed[3]
2249 << "S4.=" << &cseed[4]
2250 << "S5.=" << &cseed[5]
2251 << "SB0.=" << &seedb[0]
2252 << "SB1.=" << &seedb[1]
2253 << "SB2.=" << &seedb[2]
2254 << "SB3.=" << &seedb[3]
2255 << "SB4.=" << &seedb[4]
2256 << "SB5.=" << &seedb[5]
2257 << "Label=" << label
2258 << "Freq=" << frequency
2259 << "sLayer=" << sLayer
2260 << "\n";
2261 }
69b55c55 2262 }
4551ea7c 2263
2264 if (registered<kMaxSeed - 1) {
69b55c55 2265 registered++;
2266 cseed = seed[registered];
2267 }
4551ea7c 2268
2269 } // End of loop over layer 1
2270
2271 } // End of loop over layer 0
2272
2273 } // End of loop over layer 3
2274
2275 } // End of loop over seeding time bins
2276
69b55c55 2277 //
4551ea7c 2278 // Choose best
69b55c55 2279 //
4551ea7c 2280
69b55c55 2281 TMath::Sort(registered,seedquality2,sort,kTRUE);
c6f438c0 2282 Bool_t signedseed[kMaxSeed];
4551ea7c 2283 for (Int_t i = 0; i < registered; i++) {
2284 signedseed[i] = kFALSE;
69b55c55 2285 }
4551ea7c 2286
7bce990c 2287 for (Int_t jter = 0; jter < 5; jter++) {
4551ea7c 2288
2289 for (Int_t iseed = 0; iseed < registered; iseed++) {
2290
69b55c55 2291 Int_t index = sort[iseed];
4551ea7c 2292 if (signedseed[index]) {
2293 continue;
2294 }
69b55c55 2295 Int_t labelsall[1000];
4551ea7c 2296 Int_t nlabelsall = 0;
2297 Int_t naccepted = 0;;
2298 Int_t sLayer = seedlayer[index];
2299 Int_t ncl = 0;
2300 Int_t nused = 0;
2301 Int_t nlayers = 0;
69b55c55 2302 Int_t findable = 0;
4551ea7c 2303
2304 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2305
ca1e1e5b 2306 if (TMath::Abs(seed[index][jLayer].GetYref(0) / xcl[jLayer]) < 0.15) {
69b55c55 2307 findable++;
4551ea7c 2308 }
2309 if (seed[index][jLayer].IsOK()) {
69b55c55 2310 seed[index][jLayer].UpdateUsed();
ca1e1e5b 2311 ncl +=seed[index][jLayer].GetN2();
2312 nused +=seed[index][jLayer].GetNUsed();
69b55c55 2313 nlayers++;
4551ea7c 2314 // Cooking label
2315 for (Int_t itime = 0; itime < 25; itime++) {
ca1e1e5b 2316 if (seed[index][jLayer].IsUsable(itime)) {
69b55c55 2317 naccepted++;
4551ea7c 2318 for (Int_t ilab = 0; ilab < 3; ilab++) {
ca1e1e5b 2319 Int_t tindex = seed[index][jLayer].GetClusters(itime)->GetLabel(ilab);
4551ea7c 2320 if (tindex >= 0) {
69b55c55 2321 labelsall[nlabelsall] = tindex;
2322 nlabelsall++;
2323 }
2324 }
2325 }
2326 }
2327 }
4551ea7c 2328
69b55c55 2329 }
4551ea7c 2330
2331 if (nused > 30) {
2332 continue;
69b55c55 2333 }
4551ea7c 2334
7bce990c 2335 if (jter == 0) {
4551ea7c 2336 if (nlayers < 6) {
2337 continue;
2338 }
2339 if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2340 continue; // Gold
2341 }
7ad19338 2342 }
4551ea7c 2343
7bce990c 2344 if (jter == 1) {
4551ea7c 2345 if (nlayers < findable) {
2346 continue;
2347 }
2348 if (TMath::Log(0.000000001+seedquality2[index]) < -4.0) {
2349 continue;
2350 }
69b55c55 2351 }
4551ea7c 2352
7bce990c 2353 if (jter == 2) {
4551ea7c 2354 if ((nlayers == findable) ||
2355 (nlayers == 6)) {
2356 continue;
2357 }
2358 if (TMath::Log(0.000000001+seedquality2[index]) < -6.0) {
2359 continue;
2360 }
69b55c55 2361 }
4551ea7c 2362
7bce990c 2363 if (jter == 3) {
4551ea7c 2364 if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2365 continue;
2366 }
69b55c55 2367 }
4551ea7c 2368
7bce990c 2369 if (jter == 4) {
4551ea7c 2370 if (TMath::Log(0.000000001+seedquality2[index]) - nused/(nlayers-3.0) < -15.0) {
2371 continue;
2372 }
2373 }
2374
69b55c55 2375 signedseed[index] = kTRUE;
4551ea7c 2376
2377 Int_t labels[1000];
2378 Int_t outlab[1000];
2379 Int_t nlab = 0;
2380 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2381 if (seed[index][iLayer].IsOK()) {
ca1e1e5b 2382 if (seed[index][iLayer].GetLabels(0) >= 0) {
2383 labels[nlab] = seed[index][iLayer].GetLabels(0);
69b55c55 2384 nlab++;
2385 }
ca1e1e5b 2386 if (seed[index][iLayer].GetLabels(1) >= 0) {
2387 labels[nlab] = seed[index][iLayer].GetLabels(1);
69b55c55 2388 nlab++;
4551ea7c 2389 }
2390 }
7ad19338 2391 }
69b55c55 2392 Freq(nlab,labels,outlab,kFALSE);
4551ea7c 2393 Int_t label = outlab[0];
2394 Int_t frequency = outlab[1];
69b55c55 2395 Freq(nlabelsall,labelsall,outlab,kFALSE);
4551ea7c 2396 Int_t label1 = outlab[0];
2397 Int_t label2 = outlab[2];
2398 Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
2399 Float_t ratio = Float_t(nused) / Float_t(ncl);
2400 if (ratio < 0.25) {
2401 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2402 if ((seed[index][jLayer].IsOK()) &&
ca1e1e5b 2403 (TMath::Abs(seed[index][jLayer].GetYfit(1) - seed[index][jLayer].GetYfit(1)) < 0.2)) {
4551ea7c 2404 seed[index][jLayer].UseClusters(); // Sign gold
2405 }
69b55c55 2406 }
7ad19338 2407 }
4551ea7c 2408
31fd97b2 2409 Int_t eventNrInFile = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
4551ea7c 2410 TTreeSRedirector &cstream = *fDebugStreamer;
2411
69b55c55 2412 //
4551ea7c 2413 // Register seed
69b55c55 2414 //
4551ea7c 2415 AliTRDtrack *track = RegisterSeed(seed[index],seedparams[index]);
2416 AliTRDtrack dummy;
2417 if (!track) {
2418 track = &dummy;
2419 }
2420 else {
69b55c55 2421 AliESDtrack esdtrack;
4551ea7c 2422 esdtrack.UpdateTrackParams(track,AliESDtrack::kTRDout);
69b55c55 2423 esdtrack.SetLabel(label);
2424 esd->AddTrack(&esdtrack);
3a039a31 2425 if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0) {
4551ea7c 2426 cstream << "Tracks"
31fd97b2 2427 << "EventNrInFile=" << eventNrInFile
4551ea7c 2428 << "ESD.=" << &esdtrack
2429 << "trd.=" << track
2430 << "trdback.=" << track
2431 << "\n";
2432 }
2433 }
2434
3a039a31 2435 if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0) {
4551ea7c 2436 cstream << "Seeds2"
7bce990c 2437 << "Iter=" << jter
4551ea7c 2438 << "Track.=" << track
2439 << "Like=" << seedquality[index]
2440 << "LikeF=" << seedquality2[index]
2441 << "S0.=" << &seed[index][0]
2442 << "S1.=" << &seed[index][1]
2443 << "S2.=" << &seed[index][2]
2444 << "S3.=" << &seed[index][3]
2445 << "S4.=" << &seed[index][4]
2446 << "S5.=" << &seed[index][5]
2447 << "Label=" << label
2448 << "Label1=" << label1
2449 << "Label2=" << label2
2450 << "FakeRatio=" << fakeratio
2451 << "Freq=" << frequency
2452 << "Ncl=" << ncl
2453 << "Nlayers=" << nlayers
2454 << "Findable=" << findable
2455 << "NUsed=" << nused
2456 << "sLayer=" << sLayer
31fd97b2 2457 << "EventNrInFile=" << eventNrInFile
4551ea7c 2458 << "\n";
7ad19338 2459 }
4551ea7c 2460
2461 } // Loop: iseed
2462
7bce990c 2463 } // Loop: jter
4551ea7c 2464
2465 } // End of loop over sectors
75bd7f81 2466
69b55c55 2467 delete [] pseed;
75bd7f81 2468
e4f2f73d 2469 return 0;
69b55c55 2470}
2471
a819a5f7 2472//_____________________________________________________________________________
4551ea7c 2473Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *clusterTree) const
46d29e70 2474{
a819a5f7 2475 //
0c349049 2476 // Reads AliTRDclusters from the file.
2477 // The names of the cluster tree and branches
a819a5f7 2478 // should match the ones used in AliTRDclusterizer::WriteClusters()
2479 //
75bd7f81 2480
4551ea7c 2481 Int_t nsize = Int_t(clusterTree->GetTotBytes() / (sizeof(AliTRDcluster)));
6c94f330 2482 TObjArray *clusterArray = new TObjArray(nsize+1000);
5443e65e 2483
4551ea7c 2484 TBranch *branch = clusterTree->GetBranch("TRDcluster");
c630aafd 2485 if (!branch) {
4551ea7c 2486 AliError("Can't get the branch !");
c630aafd 2487 return 1;
2488 }
029cd327 2489 branch->SetAddress(&clusterArray);
5443e65e 2490
a819a5f7 2491 // Loop through all entries in the tree
4551ea7c 2492 Int_t nEntries = (Int_t) clusterTree->GetEntries();
2493 Int_t nbytes = 0;
e4f2f73d 2494 AliTRDcluster *c = 0x0;
a819a5f7 2495 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2496
2497 // Import the tree
4551ea7c 2498 nbytes += clusterTree->GetEvent(iEntry);
5443e65e 2499
a819a5f7 2500 // Get the number of points in the detector
029cd327 2501 Int_t nCluster = clusterArray->GetEntriesFast();
a819a5f7 2502 // Loop through all TRD digits
2503 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
e4f2f73d 2504 if(!(c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster))) continue;
2505 array->AddLast(c);
2506 //printf("Add cluster 0x%x.\n", c);
4f1c04d3 2507 clusterArray->RemoveAt(iCluster);
a819a5f7 2508 }
4551ea7c 2509
a819a5f7 2510 }
029cd327 2511 delete clusterArray;
5443e65e 2512
c630aafd 2513 return 0;
75bd7f81 2514
a819a5f7 2515}
2516
75bd7f81 2517//_____________________________________________________________________________
4551ea7c 2518Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint &p) const
3551db50 2519{
2520 //
2521 // Get track space point with index i
2522 // Origin: C.Cheshkov
2523 //
2524
4551ea7c 2525 AliTRDcluster *cl = (AliTRDcluster *) fClusters->UncheckedAt(index);
053767a4 2526 Int_t idet = cl->GetDetector();
2527 Int_t isector = fGeom->GetSector(idet);
2528 Int_t istack = fGeom->GetStack(idet);
2529 Int_t ilayer = fGeom->GetLayer(idet);
3551db50 2530 Double_t local[3];
053767a4 2531 local[0] = GetX(isector,ilayer,cl->GetLocalTimeBin());
4551ea7c 2532 local[1] = cl->GetY();
2533 local[2] = cl->GetZ();
3551db50 2534 Double_t global[3];
2535 fGeom->RotateBack(idet,local,global);
2536 p.SetXYZ(global[0],global[1],global[2]);
053767a4 2537 AliGeomManager::ELayerID iGeoLayer = AliGeomManager::kTRD1;
2538 switch (ilayer) {
3551db50 2539 case 0:
053767a4 2540 iGeoLayer = AliGeomManager::kTRD1;
3551db50 2541 break;
2542 case 1:
053767a4 2543 iGeoLayer = AliGeomManager::kTRD2;
3551db50 2544 break;
2545 case 2:
053767a4 2546 iGeoLayer = AliGeomManager::kTRD3;
3551db50 2547 break;
2548 case 3:
053767a4 2549 iGeoLayer = AliGeomManager::kTRD4;
3551db50 2550 break;
2551 case 4:
053767a4 2552 iGeoLayer = AliGeomManager::kTRD5;
3551db50 2553 break;
2554 case 5:
053767a4 2555 iGeoLayer = AliGeomManager::kTRD6;
3551db50 2556 break;
2557 };
053767a4 2558 Int_t modId = isector * fGeom->Nstack() + istack;
2559 UShort_t volid = AliGeomManager::LayerToVolUID(iGeoLayer,modId);
3551db50 2560 p.SetVolumeID(volid);
2561
2562 return kTRUE;
2563
2564}
2565
75bd7f81 2566//_____________________________________________________________________________
4551ea7c 2567void AliTRDtracker::CookLabel(AliKalmanTrack *pt, Float_t wrong) const
029cd327 2568{
2569 //
2570 // This cooks a label. Mmmmh, smells good...
2571 //
46d29e70 2572
4551ea7c 2573 Int_t label = 123456789;
2574 Int_t index;
2575 Int_t i;
2576 Int_t j;
2577 Int_t ncl = pt->GetNumberOfClusters();
2578
2579 const Int_t kRange = fTrSec[0]->GetOuterTimeBin() + 1;
5443e65e 2580
029cd327 2581 Bool_t labelAdded;
46d29e70 2582
6c94f330 2583 Int_t **s = new Int_t* [kRange];
4551ea7c 2584 for (i = 0; i < kRange; i++) {
d1b06c24 2585 s[i] = new Int_t[2];
2586 }
4551ea7c 2587 for (i = 0; i < kRange; i++) {
2588 s[i][0] = -1;
2589 s[i][1] = 0;
46d29e70 2590 }
2591
4551ea7c 2592 Int_t t0;
2593 Int_t t1;
2594 Int_t t2;
2595
2596 for (i = 0; i < ncl; i++) {
2597 index = pt->GetClusterIndex(i);
2598 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
6c94f330 2599 t0=c->GetLabel(0);
2600 t1=c->GetLabel(1);
2601 t2=c->GetLabel(2);
46d29e70 2602 }
2603
4551ea7c 2604 for (i = 0; i < ncl; i++) {
2605 index = pt->GetClusterIndex(i);
2606 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2607 for (Int_t k = 0; k < 3; k++) {
2608 label = c->GetLabel(k);
2609 labelAdded = kFALSE;
2610 j = 0;
46d29e70 2611 if (label >= 0) {
4551ea7c 2612 while ((!labelAdded) && (j < kRange)) {
2613 if ((s[j][0] == label) ||
2614 (s[j][1] == 0)) {
2615 s[j][0] = label;
2616 s[j][1] = s[j][1] + 1;
2617 labelAdded = kTRUE;
a9814c09 2618 }
2619 j++;
2620 }
46d29e70 2621 }
2622 }
2623 }
2624
4551ea7c 2625 Int_t max = 0;
6c94f330 2626 label = -123456789;
46d29e70 2627
4551ea7c 2628 for (i = 0; i < kRange; i++) {
2629 if (s[i][1] > max) {
2630 max = s[i][1];
2631 label = s[i][0];
46d29e70 2632 }
2633 }
5443e65e 2634
4551ea7c 2635 for (i = 0; i < kRange; i++) {
6c94f330 2636 delete []s[i];
5443e65e 2637 }
2638
6c94f330 2639 delete []s;
5443e65e 2640
4551ea7c 2641 if ((1.0 - Float_t(max)/ncl) > wrong) {
2642 label = -label;
2643 }
5443e65e 2644
2645 pt->SetLabel(label);
2646
46d29e70 2647}
2648
75bd7f81 2649//_____________________________________________________________________________
4551ea7c 2650void AliTRDtracker::UseClusters(const AliKalmanTrack *t, Int_t from) const
029cd327 2651{
2652 //
2653 // Use clusters, but don't abuse them!
2654 //
75bd7f81 2655
4551ea7c 2656 const Float_t kmaxchi2 = 18;
2657 const Float_t kmincl = 10;
2658
2659 AliTRDtrack *track = (AliTRDtrack *) t;
2660
2661 Int_t ncl = t->GetNumberOfClusters();
2662 for (Int_t i = from; i < ncl; i++) {
2663 Int_t index = t->GetClusterIndex(i);
2664 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
053767a4 2665 Int_t ilayer = fGeom->GetLayer(c->GetDetector());
2666 if (track->GetTracklets(ilayer).GetChi2() > kmaxchi2) {
4551ea7c 2667 continue;
2668 }
053767a4 2669 if (track->GetTracklets(ilayer).GetN() < kmincl) {
4551ea7c 2670 continue;
2671 }
2672 if (!(c->IsUsed())) {
2673 c->Use();
2674 }
5443e65e 2675 }
5443e65e 2676
75bd7f81 2677}
5443e65e 2678
75bd7f81 2679//_____________________________________________________________________________
029cd327 2680Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
5443e65e 2681{
75bd7f81 2682 //
5443e65e 2683 // Parametrised "expected" error of the cluster reconstruction in Y
75bd7f81 2684 //
5443e65e 2685
2686 Double_t s = 0.08 * 0.08;
2687 return s;
75bd7f81 2688
5443e65e 2689}
2690
75bd7f81 2691//_____________________________________________________________________________
029cd327 2692Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
5443e65e 2693{
75bd7f81 2694 //
5443e65e 2695 // Parametrised "expected" error of the cluster reconstruction in Z
75bd7f81 2696 //
5443e65e 2697
4551ea7c 2698 Double_t s = 9.0 * 9.0 / 12.0;
5443e65e 2699 return s;
75bd7f81 2700
5443e65e 2701}
2702
75bd7f81 2703//_____________________________________________________________________________
029cd327 2704Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
5443e65e 2705{
2706 //
029cd327 2707 // Returns radial position which corresponds to time bin <localTB>
5443e65e 2708 // in tracking sector <sector> and plane <plane>
2709 //
2710
6c94f330 2711 Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
4551ea7c 2712 Int_t pl = fTrSec[sector]->GetLayerNumber(index);
2713
5443e65e 2714 return fTrSec[sector]->GetLayer(pl)->GetX();
2715
2716}
2717
5443e65e 2718
75bd7f81 2719//_____________________________________________________________________________
2720AliTRDtracker::AliTRDtrackingSector
ad670fba 2721 ::AliTRDtrackingSector(AliTRDgeometry *geo, Int_t gs)
2722 :fN(0)
2723 ,fGeom(geo)
2724 ,fGeomSector(gs)
0a29d0f1 2725{
2726 //
5443e65e 2727 // AliTRDtrackingSector Constructor
0a29d0f1 2728 //
75bd7f81 2729
7e448bcc 2730 AliTRDpadPlane *padPlane = 0;
2731 AliTRDpropagationLayer *ppl = 0;
a5cadd36 2732
ad670fba 2733 // Get holes description from geometry
053767a4 2734 Bool_t holes[AliTRDgeometry::kNstack];
2735 for (Int_t istack = 0; istack < AliTRDgeometry::kNstack; istack++) {
2736 holes[istack] = fGeom->IsHole(0,istack,gs);
3c625a9b 2737 }
3c625a9b 2738
ad670fba 2739 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
2740 fTimeBinIndex[i] = -1;
2741 }
5443e65e 2742
ad670fba 2743 Double_t x;
2744 Double_t dx;
2745 Double_t rho;
2746 Double_t radLength;
5443e65e 2747
ad670fba 2748 // Add layers for each of the planes
2749 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
a305677e 2750 //Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
5443e65e 2751
053767a4 2752 const Int_t kNstacks = AliTRDgeometry::Nstack();
a305677e 2753 Int_t tbIndex;
ad670fba 2754 Double_t ymax = 0;
2755 Double_t ymaxsensitive = 0;
053767a4 2756 Double_t *zc = new Double_t[kNstacks];
2757 Double_t *zmax = new Double_t[kNstacks];
2758 Double_t *zmaxsensitive = new Double_t[kNstacks];
3551db50 2759
053767a4 2760 for (Int_t layer = 0; layer < AliTRDgeometry::Nlayer(); layer++) {
ad670fba 2761
053767a4 2762 ymax = fGeom->GetChamberWidth(layer) / 2.0;
2763 padPlane = fGeom->GetPadPlane(layer,0);
ad670fba 2764 ymaxsensitive = (padPlane->GetColSize(1) * padPlane->GetNcols() - 4.0) / 2.0;
2765
053767a4 2766 for (Int_t st = 0; st < kNstacks; st++) {
2767 zmax[st] = fGeom->GetChamberLength(layer,st) / 2.0;
ad670fba 2768 Float_t pad = padPlane->GetRowSize(1);
053767a4 2769 Float_t row0 = fGeom->GetRow0(layer,st,0);
2770 Int_t nPads = fGeom->GetRowMax(layer,st,0);
2771 zmaxsensitive[st] = Float_t(nPads) * pad / 2.0;
2772 zc[st] = -(pad * nPads) / 2.0 + row0;
5443e65e 2773 }
2774
bcb6fb78 2775 AliTRDcalibDB *fCalibration = AliTRDcalibDB::Instance();
2776 dx = fCalibration->GetVdrift(0,0,0)
f162af62 2777 / AliTRDCommonParam::Instance()->GetSamplingFrequency();
ad670fba 2778 rho = 0.00295 * 0.85; //????
2779 radLength = 11.0;
5443e65e 2780
053767a4 2781 Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(layer);
a305677e 2782 //Double_t xbottom = x0 - dxDrift;
ad670fba 2783 //Double_t xtop = x0 + dxAmp;
bcb6fb78 2784
2785 //temporary !! (A.Bercuci)
053767a4 2786 Int_t T0 = (Int_t)fCalibration->GetT0Average(AliTRDgeometry::GetDetector(layer, 2, gs));
ad670fba 2787
59393e34 2788 Int_t nTimeBins = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
ad670fba 2789 for (Int_t iTime = 0; iTime < nTimeBins; iTime++) {
2790
2791 Double_t xlayer = iTime * dx - dxAmp;
2792 //if (xlayer<0) xlayer = dxAmp / 2.0;
59393e34 2793 x = x0 - xlayer;
ad670fba 2794
053767a4 2795 tbIndex = CookTimeBinIndex(layer,iTime);
2796 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex,layer);
3c625a9b 2797 ppl->SetYmax(ymax,ymaxsensitive);
ad670fba 2798 ppl->SetZ(zc,zmax,zmaxsensitive);
3c625a9b 2799 ppl->SetHoles(holes);
bcb6fb78 2800 if(iTime == T0) ppl->SetT0();
2801
2802 InsertLayer(ppl);
ad670fba 2803
5443e65e 2804 }
ad670fba 2805
5443e65e 2806 }
2807
5443e65e 2808 MapTimeBinLayers();
ad670fba 2809
029cd327 2810 delete [] zc;
2811 delete [] zmax;
4f1c04d3 2812 delete [] zmaxsensitive;
0a29d0f1 2813
46d29e70 2814}
2815
75bd7f81 2816//_____________________________________________________________________________
ad670fba 2817AliTRDtracker::AliTRDtrackingSector
2818 ::AliTRDtrackingSector(const AliTRDtrackingSector &/*t*/)
2819 :fN(0)
2820 ,fGeom(0)
2821 ,fGeomSector(0)
2822{
2823 //
2824 // Copy constructor
2825 //
2826
2827}
2828
2829//_____________________________________________________________________________
0c349049 2830AliTRDtracker::AliTRDtrackingSector
2831 ::~AliTRDtrackingSector()
2832{
2833 //
2834 // Destructor
2835 //
2836
2837 for (Int_t i = 0; i < fN; i++) {
2838 delete fLayers[i];
2839 }
2840
2841}
2842
2843//_____________________________________________________________________________
75bd7f81 2844Int_t AliTRDtracker::AliTRDtrackingSector
2845 ::CookTimeBinIndex(Int_t plane, Int_t localTB) const
5443e65e 2846{
2847 //
0c349049 2848 // Depending on the digitization parameters calculates global
2849 // (i.e. for a whole TRD stack of 6 planes) time bin index for
2850 // timebin <localTB> in plane <plane>
59393e34 2851 //
75bd7f81 2852
59393e34 2853 Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
4551ea7c 2854 Int_t gtb = (plane+1) * tbPerPlane - localTB - 1;
2855 if (localTB < 0) {
2856 return -1;
2857 }
2858 if (gtb < 0) {
2859 return -1;
2860 }
75bd7f81 2861
5443e65e 2862 return gtb;
5443e65e 2863
75bd7f81 2864}
5443e65e 2865
75bd7f81 2866//_____________________________________________________________________________
2867void AliTRDtracker::AliTRDtrackingSector
2868 ::MapTimeBinLayers()
5443e65e 2869{
2870 //
2871 // For all sensitive time bins sets corresponding layer index
2872 // in the array fTimeBins
2873 //
2874
2875 Int_t index;
2876
4551ea7c 2877 for (Int_t i = 0; i < fN; i++) {
2878
5443e65e 2879 index = fLayers[i]->GetTimeBinIndex();
6c94f330 2880
4551ea7c 2881 if (index < 0) {
5443e65e 2882 continue;
2883 }
4551ea7c 2884 if (index >= (Int_t) kMaxTimeBinIndex) {
2885 //AliWarning(Form("Index %d exceeds allowed maximum of %d!\n"
2886 // ,index,kMaxTimeBinIndex-1));
2887 continue;
2888 }
2889
5443e65e 2890 fTimeBinIndex[index] = i;
4551ea7c 2891
5443e65e 2892 }
5443e65e 2893
75bd7f81 2894}
5443e65e 2895
75bd7f81 2896//_____________________________________________________________________________
2897Int_t AliTRDtracker::AliTRDtrackingSector
2898 ::GetLayerNumber(Double_t x) const
5443e65e 2899{
2900 //
2901 // Returns the number of time bin which in radial position is closest to <x>
2902 //
2903
4551ea7c 2904 if (x >= fLayers[fN-1]->GetX()) {
2905 return fN - 1;
2906 }
2907 if (x <= fLayers[ 0]->GetX()) {
2908 return 0;
2909 }
2910
2911 Int_t b = 0;
2912 Int_t e = fN - 1;
2913 Int_t m = (b + e) / 2;
ad670fba 2914
4551ea7c 2915 for ( ; b < e; m = (b + e) / 2) {
2916 if (x > fLayers[m]->GetX()) {
2917 b = m + 1;
2918 }
2919 else {
2920 e = m;
2921 }
5443e65e 2922 }
75bd7f81 2923
4551ea7c 2924 if (TMath::Abs(x - fLayers[m]->GetX()) > TMath::Abs(x - fLayers[m+1]->GetX())) {
2925 return m + 1;
2926 }
2927 else {
2928 return m;
2929 }
5443e65e 2930
2931}
2932
75bd7f81 2933//_____________________________________________________________________________
2934Int_t AliTRDtracker::AliTRDtrackingSector
2935 ::GetInnerTimeBin() const
5443e65e 2936{
2937 //
2938 // Returns number of the innermost SENSITIVE propagation layer
2939 //
2940
2941 return GetLayerNumber(0);
5443e65e 2942
75bd7f81 2943}
5443e65e 2944
75bd7f81 2945//_____________________________________________________________________________
2946Int_t AliTRDtracker::AliTRDtrackingSector
2947 ::GetOuterTimeBin() const
5443e65e 2948{
2949 //
2950 // Returns number of the outermost SENSITIVE time bin
2951 //
2952
2953 return GetLayerNumber(GetNumberOfTimeBins() - 1);
5443e65e 2954
75bd7f81 2955}
5443e65e 2956
75bd7f81 2957//_____________________________________________________________________________
2958Int_t AliTRDtracker::AliTRDtrackingSector
2959 ::GetNumberOfTimeBins() const
5443e65e 2960{
2961 //
2962 // Returns number of SENSITIVE time bins
2963 //
2964
4551ea7c 2965 Int_t tb;
2966 Int_t layer;
2967
2968 for (tb = kMaxTimeBinIndex - 1; tb >= 0; tb--) {
5443e65e 2969 layer = GetLayerNumber(tb);
4551ea7c 2970 if (layer >= 0) {
2971 break;
2972 }
5443e65e 2973 }
75bd7f81 2974
4551ea7c 2975 return tb + 1;
5443e65e 2976
75bd7f81 2977}
5443e65e 2978
75bd7f81 2979//_____________________________________________________________________________
2980void AliTRDtracker::AliTRDtrackingSector
4551ea7c 2981 ::InsertLayer(AliTRDpropagationLayer *pl)
5443e65e 2982{
2983 //
2984 // Insert layer <pl> in fLayers array.
2985 // Layers are sorted according to X coordinate.
75bd7f81 2986 //
5443e65e 2987
4551ea7c 2988 if (fN == ((Int_t) kMaxLayersPerSector)) {
2989 //AliWarning("Too many layers !\n");
2990 return;
2991 }
2992
2993 if (fN == 0) {
2994 fLayers[fN++] = pl;
ad670fba 2995 return;
2996 }
2997
4551ea7c 2998 Int_t i = Find(pl->GetX());
2999
3000 memmove(fLayers+i+1,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
3001
3002 fLayers[i] = pl;
3003 fN++;
5443e65e 3004
3005}
3006
75bd7f81 3007//_____________________________________________________________________________
3008Int_t AliTRDtracker::AliTRDtrackingSector
3009 ::Find(Double_t x) const
5443e65e 3010{
3011 //
3012 // Returns index of the propagation layer nearest to X
3013 //
3014
4551ea7c 3015 if (x <= fLayers[0]->GetX()) {
3016 return 0;
3017 }
3018
3019 if (x > fLayers[fN-1]->GetX()) {
3020 return fN;
3021 }
3022
3023 Int_t b = 0;
3024 Int_t e = fN-1;
3025 Int_t m = (b + e) / 2;
3026
3027 for (; b < e; m = (b + e) / 2) {
3028 if (x > fLayers[m]->GetX()) {
3029 b = m + 1;
3030 }
3031 else {
3032 e = m;
3033 }
5443e65e 3034 }
7ad19338 3035
75bd7f81 3036 return m;
7ad19338 3037
75bd7f81 3038}
7ad19338 3039
75bd7f81 3040//_____________________________________________________________________________
4551ea7c 3041Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c)
75bd7f81 3042{
3043 //
3044 // Returns correction factor for tilted pads geometry
3045 //
5443e65e 3046
4551ea7c 3047 Int_t det = c->GetDetector();
053767a4 3048 Int_t layer = fGeom->GetLayer(det);
3049 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(layer,0);
4551ea7c 3050 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
b8dc2353 3051
4551ea7c 3052 if (fNoTilt) {
3053 h01 = 0;
3054 }
75bd7f81 3055