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