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