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