]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDtracker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coding convention + var inizialization corrected
[u/mrichter/AliRoot.git] / TRD / AliTRDtracker.cxx
CommitLineData
46d29e70 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
0fa7dfa7 16/* $Id$ */
bbf92647 17
029cd327 18///////////////////////////////////////////////////////////////////////////////
19// //
c6f438c0 20// The standard TRD tracker //
21// Based on Kalman filltering approach //
029cd327 22// //
4551ea7c 23// Authors: //
24// M. Ivanov (Marian.Ivanov@cern.ch) //
25// Y. Belikov (Jouri.Belikov@cern.ch) //
26// //
029cd327 27///////////////////////////////////////////////////////////////////////////////
28
7e448bcc 29
a2cb5b3d 30#include <Riostream.h>
46d29e70 31#include <TBranch.h>
a1e17193 32#include <TFile.h>
4551ea7c 33#include <TGraph.h>
7e448bcc 34#include <TH1D.h>
35#include <TH2D.h>
a1e17193 36#include <TLinearFitter.h>
37#include <TObjArray.h>
38#include <TROOT.h>
39#include <TTree.h>
40#include <TTreeStream.h>
4551ea7c 41
af885e0f 42#include "AliESDEvent.h"
f3f19a6e 43#include "AliESDtrack.h"
4551ea7c 44#include "AliAlignObj.h"
45#include "AliRieman.h"
46#include "AliTrackPointArray.h"
46d29e70 47
46d29e70 48#include "AliTRDgeometry.h"
a5cadd36 49#include "AliTRDpadPlane.h"
0b2318ec 50#include "AliTRDgeometry.h"
46d29e70 51#include "AliTRDcluster.h"
52#include "AliTRDtrack.h"
75bd7f81 53#include "AliTRDseed.h"
3551db50 54#include "AliTRDcalibDB.h"
55#include "AliTRDCommonParam.h"
46d29e70 56#include "AliTRDtracker.h"
7b580082 57#include "AliTRDReconstructor.h"
082a2283 58#include "AliTRDCalibraFillHisto.h"
59
7e448bcc 60ClassImp(AliTRDtracker)
46d29e70 61
b08a2e7e 62const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5; //
63const Float_t AliTRDtracker::fgkLabelFraction = 0.8; //
64const Double_t AliTRDtracker::fgkMaxChi2 = 12.0; //
7e448bcc 65const Double_t AliTRDtracker::fgkMaxSnp = 0.95; // Corresponds to tan = 3
66const Double_t AliTRDtracker::fgkMaxStep = 2.0; // Maximal step size in propagation
7ad19338 67
75bd7f81 68//_____________________________________________________________________________
4551ea7c 69AliTRDtracker::AliTRDtracker()
4e009ce4 70 :AliTracker()
71 ,fHBackfit(0x0)
72 ,fHClSearch(0x0)
73 ,fHRefit(0x0)
74 ,fHX(0x0)
75 ,fHNCl(0x0)
76 ,fHNClTrack(0x0)
77 ,fHMinYPos(0x0)
78 ,fHMinYNeg(0x0)
79 ,fHMinZ(0x0)
80 ,fHMinD(0x0)
81 ,fHDeltaX(0x0)
82 ,fHXCl(0x0)
83 ,fGeom(0)
84 ,fNclusters(0)
85 ,fClusters(0)
86 ,fNseeds(0)
87 ,fSeeds(0)
88 ,fNtracks(0)
89 ,fTracks(0)
90 ,fTimeBinsPerPlane(0)
91 ,fAddTRDseeds(kFALSE)
92 ,fNoTilt(kFALSE)
93 ,fDebugStreamer(0)
89f05372 94{
75bd7f81 95 //
b7a0917f 96 // Default constructor
75bd7f81 97 //
b7a0917f 98
4551ea7c 99 for (Int_t i = 0; i < kTrackingSectors; i++) {
100 fTrSec[i] = 0;
101 }
102 for (Int_t j = 0; j < 5; j++) {
103 for (Int_t k = 0; k < 18; k++) {
104 fHoles[j][k] = kFALSE;
105 }
106 }
75bd7f81 107
7e448bcc 108 InitLogHists();
4e009ce4 109
89f05372 110}
75bd7f81 111
112//_____________________________________________________________________________
6c94f330 113AliTRDtracker::AliTRDtracker(const AliTRDtracker &t)
114 :AliTracker(t)
4e009ce4 115 ,fHBackfit(0x0)
116 ,fHClSearch(0x0)
117 ,fHRefit(0x0)
118 ,fHX(0x0)
119 ,fHNCl(0x0)
120 ,fHNClTrack(0x0)
121 ,fHMinYPos(0x0)
122 ,fHMinYNeg(0x0)
123 ,fHMinZ(0x0)
124 ,fHMinD(0x0)
125 ,fHDeltaX(0x0)
126 ,fHXCl(0x0)
ad670fba 127 ,fGeom(0)
128 ,fNclusters(0)
129 ,fClusters(0)
130 ,fNseeds(0)
131 ,fSeeds(0)
132 ,fNtracks(0)
133 ,fTracks(0)
134 ,fTimeBinsPerPlane(0)
135 ,fAddTRDseeds(kFALSE)
136 ,fNoTilt(kFALSE)
137 ,fDebugStreamer(0)
6c94f330 138{
ad670fba 139 //
140 // Copy constructor
141 //
4e009ce4 142
ad670fba 143}
144
145//_____________________________________________________________________________
eccb46a3 146AliTRDtracker::AliTRDtracker(const TFile */*geomfile*/)
ad670fba 147 :AliTracker()
4e009ce4 148 ,fHBackfit(0x0)
149 ,fHClSearch(0x0)
150 ,fHRefit(0x0)
151 ,fHX(0x0)
152 ,fHNCl(0x0)
153 ,fHNClTrack(0x0)
154 ,fHMinYPos(0x0)
155 ,fHMinYNeg(0x0)
156 ,fHMinZ(0x0)
157 ,fHMinD(0x0)
158 ,fHDeltaX(0x0)
159 ,fHXCl(0x0)
ad670fba 160 ,fGeom(0)
161 ,fNclusters(0)
162 ,fClusters(new TObjArray(2000))
163 ,fNseeds(0)
164 ,fSeeds(new TObjArray(2000))
165 ,fNtracks(0)
166 ,fTracks(new TObjArray(1000))
167 ,fTimeBinsPerPlane(0)
168 ,fAddTRDseeds(kFALSE)
169 ,fNoTilt(kFALSE)
170 ,fDebugStreamer(0)
46d29e70 171{
5443e65e 172 //
173 // Main constructor
174 //
b8dc2353 175
4551ea7c 176 TDirectory *savedir = gDirectory;
46d29e70 177
f162af62 178 fGeom = new AliTRDgeometry();
c6f438c0 179 fGeom->ReadGeoMatrices();
c630aafd 180
4551ea7c 181 for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
439c63c8 182 Int_t trS = geomS;
4551ea7c 183 fTrSec[trS] = new AliTRDtrackingSector(fGeom,geomS);
184 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
eccb46a3 185 // Could also go ...
4551ea7c 186 fHoles[icham][trS] = fGeom->IsHole(0,icham,geomS);
3c625a9b 187 }
5443e65e 188 }
4551ea7c 189
f162af62 190 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0);
7ad19338 191 Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
4551ea7c 192 if (tiltAngle < 0.1) {
b8dc2353 193 fNoTilt = kTRUE;
194 }
195
fc88348b 196 if (!AliTRDcalibDB::Instance()) {
197 AliFatal("Could not get calibration object");
198 }
199 fTimeBinsPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
46d29e70 200
4551ea7c 201 fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
0a29d0f1 202
9c9d2487 203 savedir->cd();
7e448bcc 204
205 InitLogHists();
75bd7f81 206
5443e65e 207}
46d29e70 208
75bd7f81 209//_____________________________________________________________________________
5443e65e 210AliTRDtracker::~AliTRDtracker()
46d29e70 211{
029cd327 212 //
213 // Destructor of AliTRDtracker
214 //
215
89f05372 216 if (fClusters) {
217 fClusters->Delete();
218 delete fClusters;
219 }
4551ea7c 220
89f05372 221 if (fTracks) {
222 fTracks->Delete();
223 delete fTracks;
224 }
4551ea7c 225
89f05372 226 if (fSeeds) {
227 fSeeds->Delete();
228 delete fSeeds;
229 }
4551ea7c 230
f162af62 231 if (fGeom) {
232 delete fGeom;
233 }
0a29d0f1 234
4551ea7c 235 for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
029cd327 236 delete fTrSec[geomS];
5443e65e 237 }
4551ea7c 238
7ad19338 239 if (fDebugStreamer) {
7ad19338 240 delete fDebugStreamer;
241 }
9c9d2487 242
75bd7f81 243}
59393e34 244
75bd7f81 245//_____________________________________________________________________________
246Int_t AliTRDtracker::LocalToGlobalID(Int_t lid)
247{
59393e34 248 //
75bd7f81 249 // Transform internal TRD ID to global detector ID
59393e34 250 //
75bd7f81 251
4551ea7c 252 Int_t isector = fGeom->GetSector(lid);
253 Int_t ichamber = fGeom->GetChamber(lid);
254 Int_t iplan = fGeom->GetPlane(lid);
255
ae079791 256 AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
59393e34 257 switch (iplan) {
258 case 0:
ae079791 259 iLayer = AliGeomManager::kTRD1;
59393e34 260 break;
261 case 1:
ae079791 262 iLayer = AliGeomManager::kTRD2;
59393e34 263 break;
264 case 2:
ae079791 265 iLayer = AliGeomManager::kTRD3;
59393e34 266 break;
267 case 3:
ae079791 268 iLayer = AliGeomManager::kTRD4;
59393e34 269 break;
270 case 4:
ae079791 271 iLayer = AliGeomManager::kTRD5;
59393e34 272 break;
273 case 5:
ae079791 274 iLayer = AliGeomManager::kTRD6;
59393e34 275 break;
276 };
4551ea7c 277
278 Int_t modId = isector * fGeom->Ncham() + ichamber;
ae079791 279 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,modId);
75bd7f81 280
59393e34 281 return volid;
75bd7f81 282
59393e34 283}
284
75bd7f81 285//_____________________________________________________________________________
286Int_t AliTRDtracker::GlobalToLocalID(Int_t gid)
287{
59393e34 288 //
75bd7f81 289 // Transform global detector ID to local detector ID
59393e34 290 //
75bd7f81 291
4551ea7c 292 Int_t modId = 0;
ae079791 293 AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(gid,modId);
4551ea7c 294
295 Int_t isector = modId / fGeom->Ncham();
296 Int_t ichamber = modId % fGeom->Ncham();
297 Int_t iLayer = -1;
298
59393e34 299 switch (layerId) {
ae079791 300 case AliGeomManager::kTRD1:
6c94f330 301 iLayer = 0;
59393e34 302 break;
ae079791 303 case AliGeomManager::kTRD2:
6c94f330 304 iLayer = 1;
59393e34 305 break;
ae079791 306 case AliGeomManager::kTRD3:
6c94f330 307 iLayer = 2;
59393e34 308 break;
ae079791 309 case AliGeomManager::kTRD4:
6c94f330 310 iLayer = 3;
59393e34 311 break;
ae079791 312 case AliGeomManager::kTRD5:
6c94f330 313 iLayer = 4;
59393e34 314 break;
ae079791 315 case AliGeomManager::kTRD6:
6c94f330 316 iLayer = 5;
59393e34 317 break;
318 default:
6c94f330 319 iLayer =-1;
59393e34 320 }
4551ea7c 321
322 if (iLayer < 0) {
323 return -1;
324 }
325
59393e34 326 Int_t lid = fGeom->GetDetector(iLayer,ichamber,isector);
75bd7f81 327
59393e34 328 return lid;
59393e34 329
75bd7f81 330}
59393e34 331
75bd7f81 332//_____________________________________________________________________________
333Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track)
334{
9c9d2487 335 //
336 // Rotates the track when necessary
337 //
338
339 Double_t alpha = AliTRDgeometry::GetAlpha();
4551ea7c 340 Double_t y = track->GetY();
341 Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
9c9d2487 342
4551ea7c 343 // Is this still needed ????
c630aafd 344 //Int_t ns = AliTRDgeometry::kNsect;
9c9d2487 345 //Int_t s=Int_t(track->GetAlpha()/alpha)%ns;
346
4551ea7c 347 if (y > ymax) {
9c9d2487 348 //s = (s+1) % ns;
4551ea7c 349 if (!track->Rotate( alpha)) {
350 return kFALSE;
351 }
352 }
353 else if (y < -ymax) {
9c9d2487 354 //s = (s-1+ns) % ns;
4551ea7c 355 if (!track->Rotate(-alpha)) {
356 return kFALSE;
357 }
9c9d2487 358 }
359
360 return kTRUE;
9c9d2487 361
75bd7f81 362}
46e2d86c 363
75bd7f81 364//_____________________________________________________________________________
365AliTRDcluster *AliTRDtracker::GetCluster(AliTRDtrack *track, Int_t plane
366 , Int_t timebin, UInt_t &index)
367{
46e2d86c 368 //
75bd7f81 369 // Try to find cluster in the backup list
46e2d86c 370 //
75bd7f81 371
4551ea7c 372 AliTRDcluster *cl =0;
6c94f330 373 Int_t *indexes = track->GetBackupIndexes();
4551ea7c 374
375 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
376 if (indexes[i] == 0) {
377 break;
378 }
379 AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
380 if (!cli) {
381 break;
382 }
383 if (cli->GetLocalTimeBin() != timebin) {
384 continue;
385 }
46e2d86c 386 Int_t iplane = fGeom->GetPlane(cli->GetDetector());
4551ea7c 387 if (iplane == plane) {
388 cl = cli;
7ad19338 389 index = indexes[i];
46e2d86c 390 break;
391 }
392 }
75bd7f81 393
46e2d86c 394 return cl;
46e2d86c 395
75bd7f81 396}
3c625a9b 397
75bd7f81 398//_____________________________________________________________________________
4551ea7c 399Int_t AliTRDtracker::GetLastPlane(AliTRDtrack *track)
75bd7f81 400{
3c625a9b 401 //
75bd7f81 402 // Return last updated plane
403 //
404
4551ea7c 405 Int_t lastplane = 0;
406 Int_t *indexes = track->GetBackupIndexes();
407
408 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
409 AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
410 if (!cli) {
411 break;
412 }
3c625a9b 413 Int_t iplane = fGeom->GetPlane(cli->GetDetector());
4551ea7c 414 if (iplane > lastplane) {
3c625a9b 415 lastplane = iplane;
416 }
417 }
75bd7f81 418
3c625a9b 419 return lastplane;
75bd7f81 420
3c625a9b 421}
75bd7f81 422
423//_____________________________________________________________________________
af885e0f 424Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *event)
c630aafd 425{
426 //
427 // Finds tracks within the TRD. The ESD event is expected to contain seeds
428 // at the outer part of the TRD. The seeds
429 // are found within the TRD if fAddTRDseeds is TRUE.
430 // The tracks are propagated to the innermost time bin
431 // of the TRD and the ESD event is updated
432 //
433
4551ea7c 434 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
029cd327 435 Float_t foundMin = fgkMinClustersInTrack * timeBins;
4551ea7c 436 Int_t nseed = 0;
437 Int_t found = 0;
438 //Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
c630aafd 439
440 Int_t n = event->GetNumberOfTracks();
4551ea7c 441 for (Int_t i = 0; i < n; i++) {
442
443 AliESDtrack *seed = event->GetTrack(i);
444 ULong_t status = seed->GetStatus();
445 if ((status & AliESDtrack::kTRDout) == 0) {
446 continue;
447 }
448 if ((status & AliESDtrack::kTRDin) != 0) {
449 continue;
450 }
c630aafd 451 nseed++;
7ad19338 452
4551ea7c 453 AliTRDtrack *seed2 = new AliTRDtrack(*seed);
46e2d86c 454 //seed2->ResetCovariance();
4551ea7c 455 AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
456 AliTRDtrack &t = *pt;
7b580082 457 FollowProlongation(t);
c630aafd 458 if (t.GetNumberOfClusters() >= foundMin) {
459 UseClusters(&t);
4551ea7c 460 CookLabel(pt,1 - fgkLabelFraction);
461 //t.CookdEdx();
c630aafd 462 }
463 found++;
c630aafd 464
4551ea7c 465 Double_t xTPC = 250.0;
59393e34 466 if (PropagateToX(t,xTPC,fgkMaxStep)) {
c630aafd 467 seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
468 }
469 delete seed2;
470 delete pt;
4551ea7c 471
472 }
c630aafd 473
33744e5d 474 AliInfo(Form("Number of loaded seeds: %d",nseed));
475 AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
476 AliInfo(Form("Total number of found tracks: %d",found));
c630aafd 477
478 return 0;
75bd7f81 479
c630aafd 480}
5443e65e 481
c630aafd 482//_____________________________________________________________________________
af885e0f 483Int_t AliTRDtracker::PropagateBack(AliESDEvent *event)
75bd7f81 484{
c630aafd 485 //
486 // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
487 // backpropagated by the TPC tracker. Each seed is first propagated
488 // to the TRD, and then its prolongation is searched in the TRD.
489 // If sufficiently long continuation of the track is found in the TRD
490 // the track is updated, otherwise it's stored as originaly defined
491 // by the TPC tracker.
492 //
493
7e448bcc 494 Int_t found = 0; // number of tracks found
4551ea7c 495 Float_t foundMin = 20.0;
496 Int_t n = event->GetNumberOfTracks();
497
498 // Sort tracks
499 Float_t *quality = new Float_t[n];
500 Int_t *index = new Int_t[n];
501 for (Int_t i = 0; i < n; i++) {
502 AliESDtrack *seed = event->GetTrack(i);
4f1c04d3 503 Double_t covariance[15];
504 seed->GetExternalCovariance(covariance);
505 quality[i] = covariance[0]+covariance[2];
7e448bcc 506 //quality[i] = covariance[0];
4f1c04d3 507 }
508 TMath::Sort(n,quality,index,kFALSE);
4f1c04d3 509
4551ea7c 510 for (Int_t i = 0; i < n; i++) {
511
512 //AliESDtrack *seed = event->GetTrack(i);
513 AliESDtrack *seed = event->GetTrack(index[i]);
7e448bcc 514 fHBackfit->Fill(0);
c630aafd 515
4551ea7c 516 ULong_t status = seed->GetStatus();
517 if ((status & AliESDtrack::kTPCout) == 0) {
7e448bcc 518 fHBackfit->Fill(1);
4551ea7c 519 continue;
520 }
7e448bcc 521
4551ea7c 522 if ((status & AliESDtrack::kTRDout) != 0) {
7e448bcc 523 fHBackfit->Fill(2);
4551ea7c 524 continue;
525 }
526
527 Int_t lbl = seed->GetLabel();
c630aafd 528 AliTRDtrack *track = new AliTRDtrack(*seed);
529 track->SetSeedLabel(lbl);
4551ea7c 530 seed->UpdateTrackParams(track,AliESDtrack::kTRDbackup); // Make backup
c630aafd 531 fNseeds++;
4551ea7c 532 Float_t p4 = track->GetC();
f6625211 533 Int_t expectedClr = FollowBackProlongation(*track);
7e448bcc 534
535 fHBackfit->Fill(3);
536 fHX->Fill(track->GetX());
537
538
539 // store the last measurement
540 /*
541 fHNClTrack->Fill(track->GetNumberOfClusters());
542 if (track->GetNumberOfClusters() >= foundMin) {
543
544 fHBackfit->Fill(4);
545 track->CookdEdx();
546 CookdEdxTimBin(*track);
547 CookLabel(track,1 - fgkLabelFraction);
548 if (track->GetBackupTrack()) {
549 //fHBackfit->Fill(5);
550 UseClusters(track->GetBackupTrack());
551 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
552 }
553 }
554 */
555
556 /**/
557 // inter-tracks competition ???
4551ea7c 558 if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
6c23ffed 559 (track->Pt() > 0.8)) {
7e448bcc 560
561 fHBackfit->Fill(4);
562
4f1c04d3 563 //
4551ea7c 564 // Make backup for back propagation
4f1c04d3 565 //
7e448bcc 566
4f1c04d3 567 Int_t foundClr = track->GetNumberOfClusters();
568 if (foundClr >= foundMin) {
569 track->CookdEdx();
43bfc8af 570 track->CookdEdxTimBin(); // A.Bercuci 25.07.07
4551ea7c 571 CookLabel(track,1 - fgkLabelFraction);
572 if (track->GetBackupTrack()) {
573 UseClusters(track->GetBackupTrack());
574 }
575
576 // Sign only gold tracks
577 if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
7e448bcc 578 if ((seed->GetKinkIndex(0) == 0) &&
6c23ffed 579 (track->Pt() < 1.5)) {
4551ea7c 580 UseClusters(track);
581 }
4f1c04d3 582 }
583 Bool_t isGold = kFALSE;
584
4551ea7c 585 // Full gold track
586 if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
587 //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
588 if (track->GetBackupTrack()) {
589 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
7e448bcc 590
4551ea7c 591 }
4f1c04d3 592 isGold = kTRUE;
7e448bcc 593 //fHBackfit->Fill()
4f1c04d3 594 }
4551ea7c 595
596 // Almost gold track
597 if ((!isGold) &&
598 (track->GetNCross() == 0) &&
599 (track->GetChi2() / track->GetNumberOfClusters() < 7)) {
600 //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
601 if (track->GetBackupTrack()) {
602 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
603 }
f4e9508c 604 isGold = kTRUE;
605 }
4551ea7c 606
607 if ((!isGold) &&
608 (track->GetBackupTrack())) {
7e448bcc 609 if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) &&
610 ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
4551ea7c 611 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
4f1c04d3 612 isGold = kTRUE;
613 }
614 }
4551ea7c 615
7e448bcc 616 if ((track->StatusForTOF() > 0) &&
617 (track->GetNCross() == 0) &&
27eaf44b 618 (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected()) > 0.4)) {
6c94f330 619 //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
7ad19338 620 }
4551ea7c 621
16d9fbba 622 }
c630aafd 623 }
7e448bcc 624 /**/
4551ea7c 625
7e448bcc 626 /**/
8979685e 627 // Debug part of tracking
4551ea7c 628 TTreeSRedirector &cstream = *fDebugStreamer;
31fd97b2 629 Int_t eventNrInFile = event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
4551ea7c 630 if (AliTRDReconstructor::StreamLevel() > 0) {
631 if (track->GetBackupTrack()) {
632 cstream << "Tracks"
31fd97b2 633 << "EventNrInFile=" << eventNrInFile
4551ea7c 634 << "ESD.=" << seed
635 << "trd.=" << track
636 << "trdback.=" << track->GetBackupTrack()
637 << "\n";
638 }
639 else {
640 cstream << "Tracks"
31fd97b2 641 << "EventNrInFile=" << eventNrInFile
4551ea7c 642 << "ESD.=" << seed
643 << "trd.=" << track
644 << "trdback.=" << track
645 << "\n";
d337ef8d 646 }
8979685e 647 }
7e448bcc 648 /**/
4551ea7c 649
650 // Propagation to the TOF (I.Belikov)
651 if (track->GetStop() == kFALSE) {
7e448bcc 652 fHBackfit->Fill(5);
653
4551ea7c 654 Double_t xtof = 371.0;
7e448bcc 655 Double_t xTOF0 = 370.0;
656
4551ea7c 657 Double_t c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
658 if (TMath::Abs(c2) >= 0.99) {
7e448bcc 659
660 fHBackfit->Fill(6);
c5a8e3df 661 delete track;
662 continue;
663 }
7e448bcc 664
59393e34 665 PropagateToX(*track,xTOF0,fgkMaxStep);
4551ea7c 666
667 // Energy losses taken to the account - check one more time
668 c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
669 if (TMath::Abs(c2) >= 0.99) {
7e448bcc 670
671 fHBackfit->Fill(7);
4f1c04d3 672 delete track;
673 continue;
674 }
7e448bcc 675
676 //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
677 // fHBackfit->Fill(7);
678 //delete track;
679 // continue;
680 //}
4f1c04d3 681
4551ea7c 682 Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
683 Double_t y;
684 track->GetYAt(xtof,GetBz(),y);
7e448bcc 685 if (y > ymax) {
4551ea7c 686 if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
7e448bcc 687 fHBackfit->Fill(8);
7ac6fa52 688 delete track;
7bed16a7 689 continue;
7ac6fa52 690 }
4551ea7c 691 }
692 else if (y < -ymax) {
7ac6fa52 693 if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
7e448bcc 694 fHBackfit->Fill(9);
7ac6fa52 695 delete track;
7bed16a7 696 continue;
7ac6fa52 697 }
3c625a9b 698 }
7e448bcc 699
3c625a9b 700 if (track->PropagateTo(xtof)) {
4551ea7c 701 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
7e448bcc 702 fHBackfit->Fill(10);
703
4551ea7c 704 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
705 for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
6d45eaef 706 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
707 }
708 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
eab5961e 709 }
4551ea7c 710 //seed->SetTRDtrack(new AliTRDtrack(*track));
711 if (track->GetNumberOfClusters() > foundMin) {
7e448bcc 712 fHBackfit->Fill(11);
4551ea7c 713 found++;
714 }
3c625a9b 715 }
4551ea7c 716
717 }
718 else {
7e448bcc 719
720 fHBackfit->Fill(12);
721
4551ea7c 722 if ((track->GetNumberOfClusters() > 15) &&
723 (track->GetNumberOfClusters() > 0.5*expectedClr)) {
7e448bcc 724
4551ea7c 725 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
7e448bcc 726 fHBackfit->Fill(13);
727
16d9fbba 728 //seed->SetStatus(AliESDtrack::kTRDStop);
4551ea7c 729 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
730 for (Int_t j = 0; j <AliESDtrack::kNSlice; j++) {
6d45eaef 731 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
732 }
733 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
eab5961e 734 }
7ad19338 735 //seed->SetTRDtrack(new AliTRDtrack(*track));
3c625a9b 736 found++;
737 }
4551ea7c 738
1e9bb598 739 }
4551ea7c 740
7ad19338 741 seed->SetTRDQuality(track->StatusForTOF());
27eaf44b 742 seed->SetTRDBudget(track->GetBudget(0));
8979685e 743
7e448bcc 744 fHBackfit->Fill(14);
d9b8978b 745 delete track;
c630aafd 746 }
ad670fba 747
33744e5d 748 AliInfo(Form("Number of seeds: %d",fNseeds));
749 AliInfo(Form("Number of back propagated TRD tracks: %d",found));
6c94f330 750
33744e5d 751 // New seeding
752 if (AliTRDReconstructor::SeedingOn()) {
4551ea7c 753 MakeSeedsMI(3,5,event);
33744e5d 754 }
755
756 fSeeds->Clear();
4551ea7c 757 fNseeds = 0;
7ad19338 758
4f1c04d3 759 delete [] index;
760 delete [] quality;
761
7e448bcc 762 SaveLogHists();
763
1e9bb598 764 return 0;
1e9bb598 765}
766
767//_____________________________________________________________________________
af885e0f 768Int_t AliTRDtracker::RefitInward(AliESDEvent *event)
1e9bb598 769{
770 //
771 // Refits tracks within the TRD. The ESD event is expected to contain seeds
772 // at the outer part of the TRD.
773 // The tracks are propagated to the innermost time bin
774 // of the TRD and the ESD event is updated
775 // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
776 //
777
4551ea7c 778 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
1e9bb598 779 Float_t foundMin = fgkMinClustersInTrack * timeBins;
4551ea7c 780 Int_t nseed = 0;
781 Int_t found = 0;
782 //Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
4f1c04d3 783 AliTRDtrack seed2;
6c94f330 784
1e9bb598 785 Int_t n = event->GetNumberOfTracks();
4551ea7c 786 for (Int_t i = 0; i < n; i++) {
787
788 AliESDtrack *seed = event->GetTrack(i);
789 new (&seed2) AliTRDtrack(*seed);
7e448bcc 790 fHRefit->Fill(0);
791
4551ea7c 792 if (seed2.GetX() < 270.0) {
793 seed->UpdateTrackParams(&seed2,AliESDtrack::kTRDbackup); // Backup TPC track - only update
7e448bcc 794 fHRefit->Fill(1);
f4e9508c 795 continue;
796 }
797
4551ea7c 798 ULong_t status = seed->GetStatus();
799 if ((status & AliESDtrack::kTRDout) == 0) {
7e448bcc 800 fHRefit->Fill(2);
0dd7d129 801 continue;
802 }
4551ea7c 803 if ((status & AliESDtrack::kTRDin) != 0) {
7e448bcc 804 fHRefit->Fill(3);
0dd7d129 805 continue;
806 }
7e448bcc 807
6c94f330 808 nseed++;
7e448bcc 809 fHRefit->Fill(4);
6c94f330 810
4551ea7c 811 seed2.ResetCovariance(50.0);
6c94f330 812
4f1c04d3 813 AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
4551ea7c 814 Int_t *indexes2 = seed2.GetIndexes();
815 for (Int_t i = 0; i < AliESDtrack::kNPlane;i++) {
816 for (Int_t j = 0; j < AliESDtrack::kNSlice;j++) {
6d45eaef 817 pt->SetPIDsignals(seed2.GetPIDsignals(i,j),i,j);
818 }
7ad19338 819 pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i);
820 }
eab5961e 821
4551ea7c 822 Int_t *indexes3 = pt->GetBackupIndexes();
823 for (Int_t i = 0; i < 200;i++) {
824 if (indexes2[i] == 0) {
825 break;
826 }
46e2d86c 827 indexes3[i] = indexes2[i];
4551ea7c 828 }
829
46e2d86c 830 //AliTRDtrack *pt = seed2;
306e3dc1 831 //AliTRDtrack &t = *pt;
43bfc8af 832 FollowProlongation(*pt);
833 if (pt->GetNumberOfClusters() >= foundMin) {
4551ea7c 834 //UseClusters(&t);
6c94f330 835 //CookLabel(pt, 1-fgkLabelFraction);
43bfc8af 836 pt->CookdEdx();
837 pt->CookdEdxTimBin();
306e3dc1 838 pt->CookPID(seed);
839 //pt->Calibrate(); // slot for calibration
1e9bb598 840 }
841 found++;
33744e5d 842
4551ea7c 843 Double_t xTPC = 250.0;
43bfc8af 844 if (PropagateToX(*pt,xTPC,fgkMaxStep)) {
7e448bcc 845
4551ea7c 846 seed->UpdateTrackParams(pt,AliESDtrack::kTRDrefit);
7e448bcc 847 fHRefit->Fill(5);
848
4551ea7c 849 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
850 for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
6d45eaef 851 seed->SetTRDsignals(pt->GetPIDsignals(i,j),i,j);
852 }
7ad19338 853 seed->SetTRDTimBin(pt->GetPIDTimBin(i),i);
854 }
4551ea7c 855 }
856 else {
857 // If not prolongation to TPC - propagate without update
7e448bcc 858 fHRefit->Fill(5);
4551ea7c 859 AliTRDtrack *seed2 = new AliTRDtrack(*seed);
860 seed2->ResetCovariance(5.0);
861 AliTRDtrack *pt2 = new AliTRDtrack(*seed2,seed2->GetAlpha());
7bed16a7 862 delete seed2;
59393e34 863 if (PropagateToX(*pt2,xTPC,fgkMaxStep)) {
4551ea7c 864 pt2->CookdEdx( );
306e3dc1 865 pt2->CookdEdxTimBin();
866 seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit);
867 fHRefit->Fill(6);
7e448bcc 868
4551ea7c 869 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
870 for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
6d45eaef 871 seed->SetTRDsignals(pt2->GetPIDsignals(i,j),i,j);
43bfc8af 872 }
7ad19338 873 seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
874 }
eab5961e 875 }
306e3dc1 876
877 // Add TRD track to ESDfriendTrack - maybe this tracks are not useful for post-processing - TODO make decission
878 if (AliTRDReconstructor::StreamLevel() > 0) {
879 seed->AddCalibObject(new AliTRDtrack(*pt2/*, kTRUE*/));
880 }
7bed16a7 881 delete pt2;
4551ea7c 882 }
883
306e3dc1 884 // Add TRD track to ESDfriendTrack
885 if (AliTRDReconstructor::StreamLevel() > 0) {
886 seed->AddCalibObject(new AliTRDtrack(*pt/*, kTRUE*/));
887 }
1e9bb598 888 delete pt;
4551ea7c 889
eab5961e 890 }
1e9bb598 891
33744e5d 892 AliInfo(Form("Number of loaded seeds: %d",nseed));
893 AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
1e9bb598 894
7e448bcc 895 SaveLogHists();
c630aafd 896 return 0;
c630aafd 897}
898
75bd7f81 899//_____________________________________________________________________________
4551ea7c 900Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t)
8979685e 901{
43bfc8af 902 //
903 // Starting from current position on track=t this function tries
904 // to extrapolate the track up to timeBin=0 and to confirm prolongation
905 // if a close cluster is found. Returns the number of clusters
906 // expected to be found in sensitive layers
907 // GeoManager used to estimate mean density
908 //
4551ea7c 909
43bfc8af 910 Int_t sector;
911 Int_t lastplane = GetLastPlane(&t);
97ccb02b 912 Double_t xx0 = 0.0;
913 Double_t xrho= 0.0;
43bfc8af 914 Int_t expectedNumberOfClusters = 0;
4551ea7c 915
43bfc8af 916 for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
4551ea7c 917
43bfc8af 918 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
919 Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
4551ea7c 920
43bfc8af 921 //
922 // Propagate track close to the plane if neccessary
923 //
924 Double_t currentx = fTrSec[0]->GetLayer(rowlast)->GetX();
925 if (currentx < (-fgkMaxStep + t.GetX())) {
926 // Propagate closer to chamber - safety space fgkMaxStep
927 if (!PropagateToX(t,currentx+fgkMaxStep,fgkMaxStep)) {
928 break;
929 }
930 }
4551ea7c 931
43bfc8af 932 if (!AdjustSector(&t)) {
933 break;
934 }
4551ea7c 935
43bfc8af 936 //
937 // Get material budget
938 //
939 Double_t xyz0[3];
940 Double_t xyz1[3];
941 Double_t param[7];
942 Double_t x;
943 Double_t y;
944 Double_t z;
945
946 // Starting global position
947 t.GetXYZ(xyz0);
948 // End global position
949 x = fTrSec[0]->GetLayer(row0)->GetX();
950 if (!t.GetProlongation(x,y,z)) {
951 break;
952 }
953 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
954 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
955 xyz1[2] = z;
97ccb02b 956 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
957 xrho= param[0]*param[4];
958 xx0 = param[1]; // Get mean propagation parameters
43bfc8af 959 // A.Bercuci 25.07.07
960 // Flags for marking the track momentum and direction. The track is
961 // marked only if it has at least 1 cluster picked up in the current
962 // chamber.
963 Bool_t kUPDATED = kFALSE, kMARKED = kFALSE;
4551ea7c 964
43bfc8af 965 //
966 // Propagate and update
967 //
968 sector = t.GetSector();
969 //for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
970 for (Int_t itime = 0 ; itime < GetTimeBinsPerPlane(); itime++) {
971 // A.Bercuci 25.07.07
972 // Mark track kinematics
973 if(itime > 10 && kUPDATED && !kMARKED){
974 t.SetTrackSegmentDirMom(iplane);
975 kMARKED = kTRUE;
976 }
977
978 Int_t ilayer = GetGlobalTimeBin(0,iplane,itime);
979 expectedNumberOfClusters++;
980 t.SetNExpected(t.GetNExpected() + 1);
981 if (t.GetX() > 345.0) {
982 t.SetNExpectedLast(t.GetNExpectedLast() + 1);
983 }
984 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
985 AliTRDcluster *cl = 0;
986 UInt_t index = 0;
987 Double_t maxChi2 = fgkMaxChi2;
988 x = timeBin.GetX();
989
990 if (timeBin) {
991
992 AliTRDcluster *cl0 = timeBin[0];
993 if (!cl0) {
994 // No clusters in given time bin
995 continue;
996 }
997
998 Int_t plane = fGeom->GetPlane(cl0->GetDetector());
999 if (plane > lastplane) {
1000 continue;
1001 }
1002
1003 Int_t timebin = cl0->GetLocalTimeBin();
1004 AliTRDcluster *cl2 = GetCluster(&t,plane,timebin,index);
1005
1006 if (cl2) {
1007 cl = cl2;
1008 //Double_t h01 = GetTiltFactor(cl); //I.B's fix
1009 //maxChi2=t.GetPredictedChi2(cl,h01);
1010 }
1011 if (cl) {
1012 //if (cl->GetNPads()<5)
1013 Double_t dxsample = timeBin.GetdX();
1014 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
1015 Double_t h01 = GetTiltFactor(cl);
1016 Int_t det = cl->GetDetector();
1017 Int_t plane = fGeom->GetPlane(det);
1018 if (t.GetX() > 345.0) {
1019 t.SetNLast(t.GetNLast() + 1);
1020 t.SetChi2Last(t.GetChi2Last() + maxChi2);
1021 }
1022
1023 Double_t xcluster = cl->GetX();
97ccb02b 1024 t.PropagateTo(xcluster,xx0,xrho);
43bfc8af 1025
1026 if (!AdjustSector(&t)) {
1027 break; //I.B's fix
1028 }
1029 maxChi2 = t.GetPredictedChi2(cl,h01);
1030
1031 if (maxChi2<1e+10)
1032 if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1033 // ????
1034 } else {
1035 // A.Bercuci 25.07.07
1036 //SetCluster(cl, GetNumberOfClusters()-1);
1037 kUPDATED = kTRUE;
1038 }
1039 }
1040 }
1041 }
4551ea7c 1042 }
1043
43bfc8af 1044 return expectedNumberOfClusters;
75bd7f81 1045}
7b580082 1046
75bd7f81 1047//_____________________________________________________________________________
4551ea7c 1048Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
8979685e 1049{
75bd7f81 1050 //
8979685e 1051 // Starting from current radial position of track <t> this function
1052 // extrapolates the track up to outer timebin and in the sensitive
1053 // layers confirms prolongation if a close cluster is found.
1054 // Returns the number of clusters expected to be found in sensitive layers
1055 // Use GEO manager for material Description
75bd7f81 1056 //
7e448bcc 1057 // return number of assigned clusters ?
1058 //
1059
59393e34 1060
4551ea7c 1061 Int_t sector;
1062 Int_t clusters[1000];
97ccb02b 1063 Double_t xx0 = 0.0;
1064 Double_t xrho = 0.0;
4551ea7c 1065 Int_t expectedNumberOfClusters = 0;
1066 Float_t ratio0 = 0.0;
8979685e 1067 AliTRDtracklet tracklet;
33744e5d 1068
77566f2a 1069 // Calibration fill 2D
082a2283 1070 AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
77566f2a 1071 if (!calibra) {
1072 AliInfo("Could not get Calibra instance\n");
1073 }
8ec526a4 1074 if (calibra->GetMITracking()) {
1075 calibra->ResetTrack();
77566f2a 1076 }
1077
33744e5d 1078 for (Int_t i = 0; i < 1000; i++) {
1079 clusters[i] = -1;
1080 }
1081
4551ea7c 1082 for (Int_t iplane = 0; iplane < AliESDtrack::kNPlane; iplane++) {
1083
7e448bcc 1084 int hb = iplane * 10;
1085 fHClSearch->Fill(hb);
1086
4551ea7c 1087 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
1088 Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
1089 Double_t currentx = fTrSec[0]->GetLayer(row0)->GetX();
1090 if (currentx < t.GetX()) {
7e448bcc 1091 fHClSearch->Fill(hb+1);
4551ea7c 1092 continue;
1093 }
1094
6c94f330 1095 //
4551ea7c 1096 // Propagate closer to chamber if neccessary
6c94f330 1097 //
4551ea7c 1098 if (currentx > (fgkMaxStep + t.GetX())) {
1099 if (!PropagateToX(t,currentx-fgkMaxStep,fgkMaxStep)) {
7e448bcc 1100 fHClSearch->Fill(hb+2);
4551ea7c 1101 break;
1102 }
1103 }
1104 if (!AdjustSector(&t)) {
7e448bcc 1105 fHClSearch->Fill(hb+3);
4551ea7c 1106 break;
8979685e 1107 }
4551ea7c 1108 if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) {
7e448bcc 1109 fHClSearch->Fill(hb+4);
4551ea7c 1110 break;
1111 }
1112
8979685e 1113 //
4551ea7c 1114 // Get material budget inside of chamber
59393e34 1115 //
4551ea7c 1116 Double_t xyz0[3];
1117 Double_t xyz1[3];
1118 Double_t param[7];
1119 Double_t x;
1120 Double_t y;
1121 Double_t z;
1122 // Starting global position
1123 t.GetXYZ(xyz0);
1124 // End global position
7b580082 1125 x = fTrSec[0]->GetLayer(rowlast)->GetX();
4551ea7c 1126 if (!t.GetProlongation(x,y,z)) {
7e448bcc 1127 fHClSearch->Fill(hb+5);
4551ea7c 1128 break;
1129 }
1130 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1131 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1132 xyz1[2] = z;
97ccb02b 1133 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1134 xrho = param[0]*param[4];
1135 xx0 = param[1]; // Get mean propagation parameters
4551ea7c 1136
8979685e 1137 //
7b580082 1138 // Find clusters
8979685e 1139 //
6c94f330 1140 sector = t.GetSector();
4551ea7c 1141 Float_t ncl = FindClusters(sector,row0,rowlast,&t,clusters,tracklet);
7e448bcc 1142 fHNCl->Fill(tracklet.GetN());
1143
4551ea7c 1144 if (tracklet.GetN() < GetTimeBinsPerPlane()/3) {
7e448bcc 1145 fHClSearch->Fill(hb+6);
4551ea7c 1146 continue;
1147 }
1148
8979685e 1149 //
7b580082 1150 // Propagate and update track
8979685e 1151 //
4551ea7c 1152 for (Int_t itime = GetTimeBinsPerPlane()-1; itime >= 0; itime--) {
1153
1154 Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
8979685e 1155 expectedNumberOfClusters++;
27eaf44b 1156 t.SetNExpected(t.GetNExpected() + 1);
4551ea7c 1157 if (t.GetX() > 345.0) {
27eaf44b 1158 t.SetNExpectedLast(t.GetNExpectedLast() + 1);
4551ea7c 1159 }
1160 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
7e448bcc 1161 AliTRDcluster *cl = 0;
4551ea7c 1162 UInt_t index = 0;
1163 Double_t maxChi2 = fgkMaxChi2;
8979685e 1164 x = timeBin.GetX();
4551ea7c 1165
8979685e 1166 if (timeBin) {
4551ea7c 1167
1168 if (clusters[ilayer] > 0) {
8979685e 1169 index = clusters[ilayer];
4551ea7c 1170 cl = (AliTRDcluster *)GetCluster(index);
33744e5d 1171 //Double_t h01 = GetTiltFactor(cl); // I.B's fix
1172 //maxChi2=t.GetPredictedChi2(cl,h01); //
8979685e 1173 }
1174
1175 if (cl) {
4551ea7c 1176
1177 //if (cl->GetNPads() < 5)
59393e34 1178 Double_t dxsample = timeBin.GetdX();
1179 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
4551ea7c 1180 Double_t h01 = GetTiltFactor(cl);
1181 Int_t det = cl->GetDetector();
1182 Int_t plane = fGeom->GetPlane(det);
1183 if (t.GetX() > 345.0) {
27eaf44b 1184 t.SetNLast(t.GetNLast() + 1);
1185 t.SetChi2Last(t.GetChi2Last() + maxChi2);
8979685e 1186 }
59393e34 1187 Double_t xcluster = cl->GetX();
97ccb02b 1188 t.PropagateTo(xcluster,xx0,xrho);
4551ea7c 1189 maxChi2 = t.GetPredictedChi2(cl,h01);
6c94f330 1190
c2239ca8 1191 if (maxChi2<1e+10)
1192 if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1193 if (!t.Update(cl,maxChi2,index,h01)) {
1194 // ????
1195 }
43bfc8af 1196 } // else SetCluster(cl, GetNumberOfClusters()-1); // A.Bercuci 25.07.07
1197
4551ea7c 1198
8ec526a4 1199 if (calibra->GetMITracking()) {
77566f2a 1200 calibra->UpdateHistograms(cl,&t);
1201 }
1202
4551ea7c 1203 // Reset material budget if 2 consecutive gold
1204 if (plane > 0) {
27eaf44b 1205 if ((t.GetTracklets(plane).GetN() + t.GetTracklets(plane-1).GetN()) > 20) {
1206 t.SetBudget(2,0.0);
4551ea7c 1207 }
1208 }
1209
1210 }
1211
8979685e 1212 }
4551ea7c 1213
59393e34 1214 }
4551ea7c 1215
1216 ratio0 = ncl / Float_t(fTimeBinsPerPlane);
27eaf44b 1217 Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);
1218 if ((tracklet.GetChi2() < 18.0) &&
1219 (ratio0 > 0.8) &&
1220 (ratio1 > 0.6) &&
1221 (ratio0+ratio1 > 1.5) &&
1222 (t.GetNCross() == 0) &&
1223 (TMath::Abs(t.GetSnp()) < 0.85) &&
1224 (t.GetNumberOfClusters() > 20)){
7e448bcc 1225 //if (ratio0 > 0.8) {
4551ea7c 1226 t.MakeBackupTrack(); // Make backup of the track until is gold
59393e34 1227 }
7b580082 1228
8979685e 1229 }
5443e65e 1230
75bd7f81 1231 return expectedNumberOfClusters;
75bd7f81 1232}
1e9bb598 1233
75bd7f81 1234//_____________________________________________________________________________
4551ea7c 1235Int_t AliTRDtracker::PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep)
5443e65e 1236{
75bd7f81 1237 //
5443e65e 1238 // Starting from current radial position of track <t> this function
1239 // extrapolates the track up to radial position <xToGo>.
1240 // Returns 1 if track reaches the plane, and 0 otherwise
75bd7f81 1241 //
1242
59393e34 1243 const Double_t kEpsilon = 0.00001;
4551ea7c 1244 //Double_t tanmax = TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
1245 Double_t xpos = t.GetX();
1246 Double_t dir = (xpos<xToGo) ? 1.0 : -1.0;
1247
1248 while (((xToGo-xpos)*dir) > kEpsilon) {
1249
1250 Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos),maxStep);
1251
1252 Double_t xyz0[3];
1253 Double_t xyz1[3];
1254 Double_t param[7];
1255 Double_t x;
1256 Double_t y;
1257 Double_t z;
1258 // Starting global position
1259 t.GetXYZ(xyz0);
1260 x = xpos + step;
1261
1262 if (!t.GetProlongation(x,y,z)) {
1263 return 0; // No prolongation
1264 }
1265
1266 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1267 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1268 xyz1[2] = z;
1269
97ccb02b 1270 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1271 if (!t.PropagateTo(x,param[1],param[0]*param[4])) {
4551ea7c 1272 return 0;
1273 }
59393e34 1274 AdjustSector(&t);
1275 xpos = t.GetX();
4551ea7c 1276
5443e65e 1277 }
75bd7f81 1278
5443e65e 1279 return 1;
5443e65e 1280
59393e34 1281}
5443e65e 1282
c630aafd 1283//_____________________________________________________________________________
1284Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1285{
75bd7f81 1286 //
c630aafd 1287 // Fills clusters into TRD tracking_sectors
1288 // Note that the numbering scheme for the TRD tracking_sectors
1289 // differs from that of TRD sectors
75bd7f81 1290 //
1291
c630aafd 1292 if (ReadClusters(fClusters,cTree)) {
4551ea7c 1293 AliError("Problem with reading the clusters !");
c630aafd 1294 return 1;
1295 }
4551ea7c 1296 Int_t ncl = fClusters->GetEntriesFast();
1297 fNclusters = ncl;
1298 AliInfo(Form("Sorting %d clusters",ncl));
c630aafd 1299
1300 UInt_t index;
4551ea7c 1301 for (Int_t ichamber = 0; ichamber < 5; ichamber++) {
1302 for (Int_t isector = 0; isector < 18; isector++) {
1303 fHoles[ichamber][isector] = kTRUE;
3c625a9b 1304 }
4551ea7c 1305 }
ad670fba 1306
6c94f330 1307 while (ncl--) {
33744e5d 1308
4551ea7c 1309 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(ncl);
1310 Int_t detector = c->GetDetector();
1311 Int_t localTimeBin = c->GetLocalTimeBin();
1312 Int_t sector = fGeom->GetSector(detector);
1313 Int_t plane = fGeom->GetPlane(detector);
439c63c8 1314 Int_t trackingSector = sector;
4551ea7c 1315
7e448bcc 1316 if (c->GetQ() > 10) {
3c625a9b 1317 Int_t chamber = fGeom->GetChamber(detector);
4551ea7c 1318 fHoles[chamber][trackingSector] = kFALSE;
3c625a9b 1319 }
c630aafd 1320
af26ce80 1321 Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
4551ea7c 1322 if (gtb < 0) {
1323 continue;
1324 }
029cd327 1325 Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
c630aafd 1326
4551ea7c 1327 index = ncl;
33744e5d 1328
7e448bcc 1329 fHXCl->Fill(c->GetX());
33744e5d 1330
7e448bcc 1331 fTrSec[trackingSector]->GetLayer(layer)->SetX(c->GetX());
1332 fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
af26ce80 1333
4551ea7c 1334 }
75bd7f81 1335
c630aafd 1336 return 0;
75bd7f81 1337
c630aafd 1338}
1339
5443e65e 1340//_____________________________________________________________________________
b7a0917f 1341void AliTRDtracker::UnloadClusters()
5443e65e 1342{
1343 //
1344 // Clears the arrays of clusters and tracks. Resets sectors and timebins
1345 //
1346
4551ea7c 1347 Int_t i;
1348 Int_t nentr;
5443e65e 1349
1350 nentr = fClusters->GetEntriesFast();
4551ea7c 1351 for (i = 0; i < nentr; i++) {
1352 delete fClusters->RemoveAt(i);
1353 }
b7a0917f 1354 fNclusters = 0;
5443e65e 1355
1356 nentr = fSeeds->GetEntriesFast();
4551ea7c 1357 for (i = 0; i < nentr; i++) {
1358 delete fSeeds->RemoveAt(i);
1359 }
5443e65e 1360
1361 nentr = fTracks->GetEntriesFast();
4551ea7c 1362 for (i = 0; i < nentr; i++) {
1363 delete fTracks->RemoveAt(i);
1364 }
5443e65e 1365
1366 Int_t nsec = AliTRDgeometry::kNsect;
5443e65e 1367 for (i = 0; i < nsec; i++) {
1368 for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1369 fTrSec[i]->GetLayer(pl)->Clear();
1370 }
1371 }
1372
1373}
1374
75bd7f81 1375//_____________________________________________________________________________
af885e0f 1376void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESDEvent *esd)
7ad19338 1377{
1378 //
1379 // Creates seeds using clusters between position inner plane and outer plane
1380 //
75bd7f81 1381
4551ea7c 1382 const Double_t kMaxTheta = 1.0;
1383 const Double_t kMaxPhi = 2.0;
1384
1385 const Double_t kRoad0y = 6.0; // Road for middle cluster
1386 const Double_t kRoad0z = 8.5; // Road for middle cluster
1387
1388 const Double_t kRoad1y = 2.0; // Road in y for seeded cluster
1389 const Double_t kRoad1z = 20.0; // Road in z for seeded cluster
1390
1391 const Double_t kRoad2y = 3.0; // Road in y for extrapolated cluster
1392 const Double_t kRoad2z = 20.0; // Road in z for extrapolated cluster
c6f438c0 1393 const Int_t kMaxSeed = 3000;
7ad19338 1394
4551ea7c 1395 Int_t maxSec = AliTRDgeometry::kNsect;
1396
1397 // Linear fitters in planes
1398 TLinearFitter fitterTC(2,"hyp2"); // Fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
1399 TLinearFitter fitterT2(4,"hyp4"); // Fitting with tilting pads - kz not fixed
69b55c55 1400 fitterTC.StoreData(kTRUE);
1401 fitterT2.StoreData(kTRUE);
4551ea7c 1402 AliRieman rieman(1000); // Rieman fitter
1403 AliRieman rieman2(1000); // Rieman fitter
1404
1405 // Find the maximal and minimal layer for the planes
7ad19338 1406 Int_t layers[6][2];
4551ea7c 1407 AliTRDpropagationLayer *reflayers[6];
1408 for (Int_t i = 0; i < 6; i++) {
1409 layers[i][0] = 10000;
1410 layers[i][1] = 0;
1411 }
1412 for (Int_t ns = 0; ns < maxSec; ns++) {
1413 for (Int_t ilayer = 0; ilayer < fTrSec[ns]->GetNumberOfLayers(); ilayer++) {
1414 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(ilayer));
1415 if (layer == 0) {
1416 continue;
1417 }
7ad19338 1418 Int_t det = layer[0]->GetDetector();
1419 Int_t plane = fGeom->GetPlane(det);
4551ea7c 1420 if (ilayer < layers[plane][0]) {
1421 layers[plane][0] = ilayer;
1422 }
1423 if (ilayer > layers[plane][1]) {
1424 layers[plane][1] = ilayer;
1425 }
7ad19338 1426 }
1427 }
4551ea7c 1428
f162af62 1429 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0);
6c94f330 1430 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
4551ea7c 1431 Double_t hL[6]; // Tilting angle
1432 Double_t xcl[6]; // X - position of reference cluster
1433 Double_t ycl[6]; // Y - position of reference cluster
1434 Double_t zcl[6]; // Z - position of reference cluster
1435
1436 AliTRDcluster *cl[6] = { 0, 0, 0, 0, 0, 0 }; // Seeding clusters
1437 Float_t padlength[6] = { 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 }; // Current pad-length
1438
1439 Double_t chi2R = 0.0;
1440 Double_t chi2Z = 0.0;
1441 Double_t chi2RF = 0.0;
1442 Double_t chi2ZF = 0.0;
1443
1444 Int_t nclusters; // Total number of clusters
1445 for (Int_t i = 0; i < 6; i++) {
1446 hL[i] = h01;
1447 if (i%2==1) {
1448 hL[i]*=-1.0;
1449 }
1450 }
1451
1452 // Registered seed
c6f438c0 1453 AliTRDseed *pseed = new AliTRDseed[kMaxSeed*6];
1454 AliTRDseed *seed[kMaxSeed];
4551ea7c 1455 for (Int_t iseed = 0; iseed < kMaxSeed; iseed++) {
1456 seed[iseed]= &pseed[iseed*6];
1457 }
69b55c55 1458 AliTRDseed *cseed = seed[0];
4551ea7c 1459
1460 Double_t seedquality[kMaxSeed];
1461 Double_t seedquality2[kMaxSeed];
1462 Double_t seedparams[kMaxSeed][7];
1463 Int_t seedlayer[kMaxSeed];
1464 Int_t registered = 0;
1465 Int_t sort[kMaxSeed];
1466
1467 //
1468 // Seeding part
1469 //
1470 for (Int_t ns = 0; ns < maxSec; ns++) { // Loop over sectors
1471 //for (Int_t ns = 0; ns < 5; ns++) { // Loop over sectors
1472
1473 registered = 0; // Reset registerd seed counter
1474 cseed = seed[registered];
1475 Float_t iter = 0.0;
1476
1477 for (Int_t sLayer = 2; sLayer >= 0; sLayer--) {
1478 //for (Int_t dseed = 5; dseed < 15; dseed += 3) {
1479
1480 iter += 1.0;
1481 Int_t dseed = 5 + Int_t(iter) * 3;
1482
69b55c55 1483 // Initialize seeding layers
4551ea7c 1484 for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
69b55c55 1485 reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
1486 xcl[ilayer] = reflayers[ilayer]->GetX();
4551ea7c 1487 }
1488
1489 Double_t xref = (xcl[sLayer+1] + xcl[sLayer+2]) * 0.5;
1490 AliTRDpropagationLayer &layer0 = *reflayers[sLayer+0];
1491 AliTRDpropagationLayer &layer1 = *reflayers[sLayer+1];
1492 AliTRDpropagationLayer &layer2 = *reflayers[sLayer+2];
1493 AliTRDpropagationLayer &layer3 = *reflayers[sLayer+3];
1494
1495 Int_t maxn3 = layer3;
1496 for (Int_t icl3 = 0; icl3 < maxn3; icl3++) {
33744e5d 1497
69b55c55 1498 AliTRDcluster *cl3 = layer3[icl3];
4551ea7c 1499 if (!cl3) {
1500 continue;
1501 }
1502 padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2() * 12.0);
69b55c55 1503 ycl[sLayer+3] = cl3->GetY();
1504 zcl[sLayer+3] = cl3->GetZ();
4551ea7c 1505 Float_t yymin0 = ycl[sLayer+3] - 1.0 - kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1506 Float_t yymax0 = ycl[sLayer+3] + 1.0 + kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1507 Int_t maxn0 = layer0;
1508
1509 for (Int_t icl0 = layer0.Find(yymin0); icl0 < maxn0; icl0++) {
1510
69b55c55 1511 AliTRDcluster *cl0 = layer0[icl0];
4551ea7c 1512 if (!cl0) {
1513 continue;
1514 }
1515 if (cl3->IsUsed() && cl0->IsUsed()) {
1516 continue;
1517 }
69b55c55 1518 ycl[sLayer+0] = cl0->GetY();
1519 zcl[sLayer+0] = cl0->GetZ();
4551ea7c 1520 if (ycl[sLayer+0] > yymax0) {
1521 break;
1522 }
1523 Double_t tanphi = (ycl[sLayer+3]-ycl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1524 if (TMath::Abs(tanphi) > kMaxPhi) {
1525 continue;
1526 }
1527 Double_t tantheta = (zcl[sLayer+3]-zcl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1528 if (TMath::Abs(tantheta) > kMaxTheta) {
1529 continue;
1530 }
1531 padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2() * 12.0);
1532
1533 // Expected position in 1 layer
1534 Double_t y1exp = ycl[sLayer+0] + (tanphi) * (xcl[sLayer+1]-xcl[sLayer+0]);
1535 Double_t z1exp = zcl[sLayer+0] + (tantheta) * (xcl[sLayer+1]-xcl[sLayer+0]);
1536 Float_t yymin1 = y1exp - kRoad0y - tanphi;
1537 Float_t yymax1 = y1exp + kRoad0y + tanphi;
1538 Int_t maxn1 = layer1;
1539
1540 for (Int_t icl1 = layer1.Find(yymin1); icl1 < maxn1; icl1++) {
1541
69b55c55 1542 AliTRDcluster *cl1 = layer1[icl1];
4551ea7c 1543 if (!cl1) {
1544 continue;
1545 }
69b55c55 1546 Int_t nusedCl = 0;
1547 if (cl3->IsUsed()) nusedCl++;
1548 if (cl0->IsUsed()) nusedCl++;
1549 if (cl1->IsUsed()) nusedCl++;
4551ea7c 1550 if (nusedCl > 1) {
1551 continue;
1552 }
69b55c55 1553 ycl[sLayer+1] = cl1->GetY();
1554 zcl[sLayer+1] = cl1->GetZ();
4551ea7c 1555 if (ycl[sLayer+1] > yymax1) {
1556 break;
1557 }
1558 if (TMath::Abs(ycl[sLayer+1]-y1exp) > kRoad0y+tanphi) {
1559 continue;
1560 }
1561 if (TMath::Abs(zcl[sLayer+1]-z1exp) > kRoad0z) {
1562 continue;
1563 }
1564 padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2() * 12.0);
1565
1566 Double_t y2exp = ycl[sLayer+0]+(tanphi) * (xcl[sLayer+2]-xcl[sLayer+0]) + (ycl[sLayer+1]-y1exp);
1567 Double_t z2exp = zcl[sLayer+0]+(tantheta) * (xcl[sLayer+2]-xcl[sLayer+0]);
1568 Int_t index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y,kRoad1z);
1569 if (index2 <= 0) {
1570 continue;
1571 }
1572 AliTRDcluster *cl2 = (AliTRDcluster *) GetCluster(index2);
1573 padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2() * 12.0);
1574 ycl[sLayer+2] = cl2->GetY();
1575 zcl[sLayer+2] = cl2->GetZ();
1576 if (TMath::Abs(cl2->GetZ()-z2exp) > kRoad0z) {
1577 continue;
1578 }
1579
69b55c55 1580 rieman.Reset();
1581 rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
1582 rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
1583 rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);
1584 rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
1585 rieman.Update();
4551ea7c 1586
1587 // Reset fitter
1588 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
69b55c55 1589 cseed[iLayer].Reset();
1590 }
4551ea7c 1591 chi2Z = 0.0;
1592 chi2R = 0.0;
1593
1594 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
ca1e1e5b 1595 cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1596 chi2Z += (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer])
1597 * (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer]);
1598 cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
1599 cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1600 chi2R += (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer])
1601 * (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer]);
1602 cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
69b55c55 1603 }
4551ea7c 1604 if (TMath::Sqrt(chi2R) > 1.0/iter) {
1605 continue;
1606 }
1607 if (TMath::Sqrt(chi2Z) > 7.0/iter) {
1608 continue;
69b55c55 1609 }
4551ea7c 1610
1611 Float_t minmax[2] = { -100.0, 100.0 };
1612 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1613 Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer] * 0.5
ca1e1e5b 1614 + 1.0 - cseed[sLayer+iLayer].GetZref(0);
4551ea7c 1615 if (max < minmax[1]) {
1616 minmax[1] = max;
1617 }
1618 Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer] * 0.5
ca1e1e5b 1619 - 1.0 - cseed[sLayer+iLayer].GetZref(0);
4551ea7c 1620 if (min > minmax[0]) {
1621 minmax[0] = min;
1622 }
1623 }
1624
69b55c55 1625 Bool_t isFake = kFALSE;
4551ea7c 1626 if (cl0->GetLabel(0) != cl3->GetLabel(0)) {
1627 isFake = kTRUE;
1628 }
1629 if (cl1->GetLabel(0) != cl3->GetLabel(0)) {
1630 isFake = kTRUE;
1631 }
1632 if (cl2->GetLabel(0) != cl3->GetLabel(0)) {
1633 isFake = kTRUE;
1634 }
1635
1636 if (AliTRDReconstructor::StreamLevel() > 0) {
1637 if ((!isFake) || ((icl3%10) == 0)) { // Debugging print
1638 TTreeSRedirector &cstream = *fDebugStreamer;
1639 cstream << "Seeds0"
1640 << "isFake=" << isFake
1641 << "Cl0.=" << cl0
1642 << "Cl1.=" << cl1
1643 << "Cl2.=" << cl2
1644 << "Cl3.=" << cl3
1645 << "Xref=" << xref
1646 << "X0=" << xcl[sLayer+0]
1647 << "X1=" << xcl[sLayer+1]
1648 << "X2=" << xcl[sLayer+2]
1649 << "X3=" << xcl[sLayer+3]
1650 << "Y2exp=" << y2exp
1651 << "Z2exp=" << z2exp
1652 << "Chi2R=" << chi2R
1653 << "Chi2Z=" << chi2Z
1654 << "Seed0.=" << &cseed[sLayer+0]
1655 << "Seed1.=" << &cseed[sLayer+1]
1656 << "Seed2.=" << &cseed[sLayer+2]
1657 << "Seed3.=" << &cseed[sLayer+3]
1658 << "Zmin=" << minmax[0]
1659 << "Zmax=" << minmax[1]
1660 << "\n";
d337ef8d 1661 }
69b55c55 1662 }
33744e5d 1663
4551ea7c 1664 ////////////////////////////////////////////////////////////////////////////////////
33744e5d 1665 //
4551ea7c 1666 // Fit seeding part
33744e5d 1667 //
4551ea7c 1668 ////////////////////////////////////////////////////////////////////////////////////
33744e5d 1669
69b55c55 1670 cl[sLayer+0] = cl0;
1671 cl[sLayer+1] = cl1;
1672 cl[sLayer+2] = cl2;
1673 cl[sLayer+3] = cl3;
4551ea7c 1674 Bool_t isOK = kTRUE;
1675
1676 for (Int_t jLayer = 0; jLayer < 4; jLayer++) {
1677
ca1e1e5b 1678 cseed[sLayer+jLayer].SetTilt(hL[sLayer+jLayer]);
1679 cseed[sLayer+jLayer].SetPadLength(padlength[sLayer+jLayer]);
1680 cseed[sLayer+jLayer].SetX0(xcl[sLayer+jLayer]);
4551ea7c 1681
1682 for (Int_t iter = 0; iter < 2; iter++) {
1683
69b55c55 1684 //
4551ea7c 1685 // In iteration 0 we try only one pad-row
1686 // If quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
69b55c55 1687 //
1688 AliTRDseed tseed = cseed[sLayer+jLayer];
4551ea7c 1689 Float_t roadz = padlength[sLayer+jLayer] * 0.5;
1690 if (iter > 0) {
1691 roadz = padlength[sLayer+jLayer];
1692 }
1693
1694 Float_t quality = 10000.0;
1695
1696 for (Int_t iTime = 2; iTime < 20; iTime++) {
1697
1698 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
1699 Double_t dxlayer = layer.GetX() - xcl[sLayer+jLayer];
1700 Double_t zexp = cl[sLayer+jLayer]->GetZ();
1701
1702 if (iter > 0) {
1703 // Try 2 pad-rows in second iteration
ca1e1e5b 1704 zexp = tseed.GetZref(0) + tseed.GetZref(1) * dxlayer;
4551ea7c 1705 if (zexp > cl[sLayer+jLayer]->GetZ()) {
1706 zexp = cl[sLayer+jLayer]->GetZ() + padlength[sLayer+jLayer]*0.5;
1707 }
1708 if (zexp < cl[sLayer+jLayer]->GetZ()) {
1709 zexp = cl[sLayer+jLayer]->GetZ() - padlength[sLayer+jLayer]*0.5;
1710 }
69b55c55 1711 }
4551ea7c 1712
ca1e1e5b 1713 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
4551ea7c 1714 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
1715 if (index <= 0) {
1716 continue;
1717 }
1718 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
1719
ca1e1e5b 1720 tseed.SetIndexes(iTime,index);
1721 tseed.SetClusters(iTime,cl); // Register cluster
1722 tseed.SetX(iTime,dxlayer); // Register cluster
1723 tseed.SetY(iTime,cl->GetY()); // Register cluster
1724 tseed.SetZ(iTime,cl->GetZ()); // Register cluster
4551ea7c 1725
1726 }
1727
69b55c55 1728 tseed.Update();
4551ea7c 1729
1730 // Count the number of clusters and distortions into quality
ca1e1e5b 1731 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1732 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
1733 + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
1734 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 1735 if ((iter == 0) && tseed.IsOK()) {
1736 cseed[sLayer+jLayer] = tseed;
1737 quality = tquality;
1738 if (tquality < 5) {
1739 break;
1740 }
1741 }
1742 if (tseed.IsOK() && (tquality < quality)) {
69b55c55 1743 cseed[sLayer+jLayer] = tseed;
69b55c55 1744 }
4551ea7c 1745
1746 } // Loop: iter
1747
1748 if (!cseed[sLayer+jLayer].IsOK()) {
69b55c55 1749 isOK = kFALSE;
1750 break;
4551ea7c 1751 }
1752
69b55c55 1753 cseed[sLayer+jLayer].CookLabels();
1754 cseed[sLayer+jLayer].UpdateUsed();
ca1e1e5b 1755 nusedCl += cseed[sLayer+jLayer].GetNUsed();
4551ea7c 1756 if (nusedCl > 25) {
69b55c55 1757 isOK = kFALSE;
1758 break;
4551ea7c 1759 }
1760
1761 } // Loop: jLayer
1762
1763 if (!isOK) {
1764 continue;
69b55c55 1765 }
4551ea7c 1766 nclusters = 0;
1767 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1768 if (cseed[sLayer+iLayer].IsOK()) {
ca1e1e5b 1769 nclusters += cseed[sLayer+iLayer].GetN2();
69b55c55 1770 }
1771 }
4551ea7c 1772
1773 // Iteration 0
69b55c55 1774 rieman.Reset();
4551ea7c 1775 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1776 rieman.AddPoint(xcl[sLayer+iLayer]
ca1e1e5b 1777 ,cseed[sLayer+iLayer].GetYfitR(0)
1778 ,cseed[sLayer+iLayer].GetZProb()
4551ea7c 1779 ,1
1780 ,10);
69b55c55 1781 }
1782 rieman.Update();
4551ea7c 1783
1784 chi2R = 0.0;
1785 chi2Z = 0.0;
1786
1787 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
ca1e1e5b 1788 cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1789 chi2R += (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0))
1790 * (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0));
1791 cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
1792 cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1793 chi2Z += (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz())
1794 * (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz());
1795 cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
69b55c55 1796 }
1797 Double_t curv = rieman.GetC();
4551ea7c 1798
69b55c55 1799 //
4551ea7c 1800 // Likelihoods
6c94f330 1801 //
ca1e1e5b 1802 Double_t sumda = TMath::Abs(cseed[sLayer+0].GetYfitR(1) - cseed[sLayer+0].GetYref(1))
1803 + TMath::Abs(cseed[sLayer+1].GetYfitR(1) - cseed[sLayer+1].GetYref(1))
1804 + TMath::Abs(cseed[sLayer+2].GetYfitR(1) - cseed[sLayer+2].GetYref(1))
1805 + TMath::Abs(cseed[sLayer+3].GetYfitR(1) - cseed[sLayer+3].GetYref(1));
4551ea7c 1806 Double_t likea = TMath::Exp(-sumda*10.6);
1807 Double_t likechi2 = 0.0000000001;
1808 if (chi2R < 0.5) {
1809 likechi2 += TMath::Exp(-TMath::Sqrt(chi2R) * 7.73);
1810 }
1811 Double_t likechi2z = TMath::Exp(-chi2Z * 0.088) / TMath::Exp(-chi2Z * 0.019);
1812 Double_t likeN = TMath::Exp(-(72 - nclusters) * 0.19);
1813 Double_t like = likea * likechi2 * likechi2z * likeN;
ca1e1e5b 1814 Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetYref(1) - 130.0*curv) * 1.9);
1815 Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetZref(1)
1816 - cseed[sLayer+0].GetZref(0) / xcl[sLayer+0]) * 5.9);
6c94f330 1817 Double_t likePrim = TMath::Max(likePrimY*likePrimZ,0.0005);
69b55c55 1818
4551ea7c 1819 seedquality[registered] = like;
1820 seedlayer[registered] = sLayer;
1821 if (TMath::Log(0.000000000000001 + like) < -15) {
1822 continue;
1823 }
69b55c55 1824 AliTRDseed seedb[6];
4551ea7c 1825 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
69b55c55 1826 seedb[iLayer] = cseed[iLayer];
1827 }
4551ea7c 1828
1829 ////////////////////////////////////////////////////////////////////////////////////
ad670fba 1830 //
4551ea7c 1831 // Full track fit part
33744e5d 1832 //
4551ea7c 1833 ////////////////////////////////////////////////////////////////////////////////////
1834
1835 Int_t nlayers = 0;
1836 Int_t nusedf = 0;
1837 Int_t findable = 0;
1838
69b55c55 1839 //
4551ea7c 1840 // Add new layers - avoid long extrapolation
69b55c55 1841 //
4551ea7c 1842 Int_t tLayer[2] = { 0, 0 };
1843 if (sLayer == 2) {
1844 tLayer[0] = 1;
1845 tLayer[1] = 0;
1846 }
1847 if (sLayer == 1) {
1848 tLayer[0] = 5;
1849 tLayer[1] = 0;
1850 }
1851 if (sLayer == 0) {
1852 tLayer[0] = 4;
1853 tLayer[1] = 5;
1854 }
1855
1856 for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
1857 Int_t jLayer = tLayer[iLayer]; // Set tracking layer
69b55c55 1858 cseed[jLayer].Reset();
ca1e1e5b 1859 cseed[jLayer].SetTilt(hL[jLayer]);
1860 cseed[jLayer].SetPadLength(padlength[jLayer]);
1861 cseed[jLayer].SetX0(xcl[jLayer]);
4551ea7c 1862 // Get pad length and rough cluster
ca1e1e5b 1863 Int_t indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].GetYref(0)
1864 ,cseed[jLayer].GetZref(0)
4551ea7c 1865 ,kRoad2y
1866 ,kRoad2z);
1867 if (indexdummy <= 0) {
1868 continue;
1869 }
1870 AliTRDcluster *cldummy = (AliTRDcluster *) GetCluster(indexdummy);
1871 padlength[jLayer] = TMath::Sqrt(cldummy->GetSigmaZ2() * 12.0);
69b55c55 1872 }
4551ea7c 1873 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
1874
1875 for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
1876
1877 Int_t jLayer = tLayer[iLayer]; // set tracking layer
1878 if ((jLayer == 0) && !(cseed[1].IsOK())) {
1879 continue; // break not allowed
1880 }
1881 if ((jLayer == 5) && !(cseed[4].IsOK())) {
1882 continue; // break not allowed
1883 }
ca1e1e5b 1884 Float_t zexp = cseed[jLayer].GetZref(0);
4551ea7c 1885 Double_t zroad = padlength[jLayer] * 0.5 + 1.0;
1886
1887 for (Int_t iter = 0; iter < 2; iter++) {
1888
69b55c55 1889 AliTRDseed tseed = cseed[jLayer];
4551ea7c 1890 Float_t quality = 10000.0;
1891
1892 for (Int_t iTime = 2; iTime < 20; iTime++) {
1893 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
1894 Double_t dxlayer = layer.GetX()-xcl[jLayer];
ca1e1e5b 1895 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
4551ea7c 1896 Float_t yroad = kRoad1y;
1897 Int_t index = layer.FindNearestCluster(yexp,zexp,yroad,zroad);
1898 if (index <= 0) {
1899 continue;
1900 }
1901 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
ca1e1e5b 1902 tseed.SetIndexes(iTime,index);
1903 tseed.SetClusters(iTime,cl); // Register cluster
1904 tseed.SetX(iTime,dxlayer); // Register cluster
1905 tseed.SetY(iTime,cl->GetY()); // Register cluster
1906 tseed.SetZ(iTime,cl->GetZ()); // Register cluster
4551ea7c 1907 }
1908
69b55c55 1909 tseed.Update();
4551ea7c 1910 if (tseed.IsOK()) {
ca1e1e5b 1911 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1912 Float_t tquality = (18.0 - tseed.GetN2())/2.0 + TMath::Abs(dangle) / 0.1
1913 + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
1914 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 1915 if (tquality < quality) {
1916 cseed[jLayer] = tseed;
1917 quality = tquality;
69b55c55 1918 }
1919 }
4551ea7c 1920
1921 zroad *= 2.0;
1922
1923 } // Loop: iter
1924
1925 if ( cseed[jLayer].IsOK()) {
69b55c55 1926 cseed[jLayer].CookLabels();
1927 cseed[jLayer].UpdateUsed();
ca1e1e5b 1928 nusedf += cseed[jLayer].GetNUsed();
4551ea7c 1929 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
69b55c55 1930 }
4551ea7c 1931
1932 } // Loop: iLayer
1933
1934 // Make copy
69b55c55 1935 AliTRDseed bseed[6];
4551ea7c 1936 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
69b55c55 1937 bseed[jLayer] = cseed[jLayer];
4551ea7c 1938 }
1939 Float_t lastquality = 10000.0;
1940 Float_t lastchi2 = 10000.0;
1941 Float_t chi2 = 1000.0;
1942
1943 for (Int_t iter = 0; iter < 4; iter++) {
1944
1945 // Sort tracklets according "quality", try to "improve" 4 worst
1946 Float_t sumquality = 0.0;
69b55c55 1947 Float_t squality[6];
1948 Int_t sortindexes[6];
4551ea7c 1949
1950 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1951 if (bseed[jLayer].IsOK()) {
69b55c55 1952 AliTRDseed &tseed = bseed[jLayer];
ca1e1e5b 1953 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
1954 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1955 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
1956 + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
1957 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 1958 squality[jLayer] = tquality;
1959 }
1960 else {
1961 squality[jLayer] = -1.0;
ad670fba 1962 }
6c94f330 1963 sumquality +=squality[jLayer];
69b55c55 1964 }
1965
4551ea7c 1966 if ((sumquality >= lastquality) ||
1967 (chi2 > lastchi2)) {
1968 break;
1969 }
69b55c55 1970 lastquality = sumquality;
1971 lastchi2 = chi2;
4551ea7c 1972 if (iter > 0) {
1973 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
69b55c55 1974 cseed[jLayer] = bseed[jLayer];
1975 }
1976 }
1977 TMath::Sort(6,squality,sortindexes,kFALSE);
4551ea7c 1978
1979 for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
1980
6c94f330 1981 Int_t bLayer = sortindexes[jLayer];
1982 AliTRDseed tseed = bseed[bLayer];
4551ea7c 1983
1984 for (Int_t iTime = 2; iTime < 20; iTime++) {
1985
1986 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
1987 Double_t dxlayer = layer.GetX() - xcl[bLayer];
ca1e1e5b 1988 Double_t zexp = tseed.GetZref(0);
1989 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
4551ea7c 1990 Float_t roadz = padlength[bLayer] + 1;
ca1e1e5b 1991 if (TMath::Abs(tseed.GetZProb() - zexp) > 0.5*padlength[bLayer]) {
4551ea7c 1992 roadz = padlength[bLayer] * 0.5;
1993 }
ca1e1e5b 1994 if (tseed.GetZfit(1)*tseed.GetZref(1) < 0.0) {
4551ea7c 1995 roadz = padlength[bLayer] * 0.5;
1996 }
ca1e1e5b 1997 if (TMath::Abs(tseed.GetZProb() - zexp) < 0.1*padlength[bLayer]) {
1998 zexp = tseed.GetZProb();
4551ea7c 1999 roadz = padlength[bLayer] * 0.5;
2000 }
2001
ca1e1e5b 2002 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer - zcor;
4551ea7c 2003 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
2004 if (index <= 0) {
2005 continue;
69b55c55 2006 }
4551ea7c 2007 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
2008
ca1e1e5b 2009 tseed.SetIndexes(iTime,index);
2010 tseed.SetClusters(iTime,cl); // Register cluster
2011 tseed.SetX(iTime,dxlayer); // Register cluster
2012 tseed.SetY(iTime,cl->GetY()); // Register cluster
2013 tseed.SetZ(iTime,cl->GetZ()); // Register cluster
4551ea7c 2014
2015 }
2016
69b55c55 2017 tseed.Update();
c6f438c0 2018 if (tseed.IsOK()) {
ca1e1e5b 2019 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
2020 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
2021 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0
2022 + TMath::Abs(dangle) / 0.1
2023 + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
2024 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 2025 if (tquality<squality[bLayer]) {
69b55c55 2026 bseed[bLayer] = tseed;
4551ea7c 2027 }
69b55c55 2028 }
4551ea7c 2029
2030 } // Loop: jLayer
2031
2032 chi2 = AliTRDseed::FitRiemanTilt(bseed,kTRUE);
2033
2034 } // Loop: iter
2035
2036 nclusters = 0;
2037 nlayers = 0;
2038 findable = 0;
2039 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ca1e1e5b 2040 if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) {
69b55c55 2041 findable++;
4551ea7c 2042 }
2043 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2044 nclusters += cseed[iLayer].GetN2();
69b55c55 2045 nlayers++;
2046 }
2047 }
4551ea7c 2048 if (nlayers < 3) {
2049 continue;
2050 }
69b55c55 2051 rieman.Reset();
4551ea7c 2052 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2053 if (cseed[iLayer].IsOK()) {
2054 rieman.AddPoint(xcl[iLayer]
ca1e1e5b 2055 ,cseed[iLayer].GetYfitR(0)
2056 ,cseed[iLayer].GetZProb()
4551ea7c 2057 ,1
2058 ,10);
2059 }
69b55c55 2060 }
2061 rieman.Update();
4551ea7c 2062
2063 chi2RF = 0.0;
2064 chi2ZF = 0.0;
2065 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2066 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2067 cseed[iLayer].SetYref(0,rieman.GetYat(xcl[iLayer]));
2068 chi2RF += (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0))
2069 * (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0));
2070 cseed[iLayer].SetYref(1,rieman.GetDYat(xcl[iLayer]));
2071 cseed[iLayer].SetZref(0,rieman.GetZat(xcl[iLayer]));
2072 chi2ZF += (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz())
2073 * (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz());
2074 cseed[iLayer].SetZref(1,rieman.GetDZat(xcl[iLayer]));
69b55c55 2075 }
2076 }
4551ea7c 2077 chi2RF /= TMath::Max((nlayers - 3.0),1.0);
2078 chi2ZF /= TMath::Max((nlayers - 3.0),1.0);
2079 curv = rieman.GetC();
2080
2081 Double_t xref2 = (xcl[2] + xcl[3]) * 0.5; // Middle of the chamber
2082 Double_t dzmf = rieman.GetDZat(xref2);
2083 Double_t zmf = rieman.GetZat(xref2);
69b55c55 2084
69b55c55 2085 //
4551ea7c 2086 // Fit hyperplane
69b55c55 2087 //
4551ea7c 2088 Int_t npointsT = 0;
69b55c55 2089 fitterTC.ClearPoints();
2090 fitterT2.ClearPoints();
2091 rieman2.Reset();
4551ea7c 2092
2093 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2094
2095 if (!cseed[iLayer].IsOK()) {
2096 continue;
2097 }
2098
2099 for (Int_t itime = 0; itime < 25; itime++) {
2100
ca1e1e5b 2101 if (!cseed[iLayer].IsUsable(itime)) {
4551ea7c 2102 continue;
2103 }
2104 // X relative to the middle chamber
ca1e1e5b 2105 Double_t x = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;
2106 Double_t y = cseed[iLayer].GetY(itime);
2107 Double_t z = cseed[iLayer].GetZ(itime);
69b55c55 2108 // ExB correction to the correction
4551ea7c 2109 // Tilted rieman
69b55c55 2110 Double_t uvt[6];
4551ea7c 2111 // Global x
ca1e1e5b 2112 Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();
4551ea7c 2113 Double_t t = 1.0 / (x2*x2 + y*y);
2114 uvt[1] = t; // t
2115 uvt[0] = 2.0 * x2 * uvt[1]; // u
2116 uvt[2] = 2.0 * hL[iLayer] * uvt[1];
2117 uvt[3] = 2.0 * hL[iLayer] * x * uvt[1];
2118 uvt[4] = 2.0 * (y + hL[iLayer]*z) * uvt[1];
2119 Double_t error = 2.0 * 0.2 * uvt[1];
69b55c55 2120 fitterT2.AddPoint(uvt,uvt[4],error);
4551ea7c 2121
69b55c55 2122 //
4551ea7c 2123 // Constrained rieman
69b55c55 2124 //
ca1e1e5b 2125 z = cseed[iLayer].GetZ(itime);
4551ea7c 2126 uvt[0] = 2.0 * x2 * t; // u
2127 uvt[1] = 2.0 * hL[iLayer] * x2 * uvt[1];
2128 uvt[2] = 2.0 * (y + hL[iLayer] * (z - GetZ())) * t;
69b55c55 2129 fitterTC.AddPoint(uvt,uvt[2],error);
69b55c55 2130 rieman2.AddPoint(x2,y,z,1,10);
2131 npointsT++;
4551ea7c 2132
69b55c55 2133 }
4551ea7c 2134
2135 } // Loop: iLayer
2136
69b55c55 2137 rieman2.Update();
2138 fitterTC.Eval();
2139 fitterT2.Eval();
2140 Double_t rpolz0 = fitterT2.GetParameter(3);
2141 Double_t rpolz1 = fitterT2.GetParameter(4);
4551ea7c 2142
69b55c55 2143 //
4551ea7c 2144 // Linear fitter - not possible to make boundaries
2145 // Do not accept non possible z and dzdx combinations
69b55c55 2146 //
4551ea7c 2147 Bool_t acceptablez = kTRUE;
2148 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2149 if (cseed[iLayer].IsOK()) {
2150 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
ca1e1e5b 2151 if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > padlength[iLayer] * 0.5 + 1.0) {
69b55c55 2152 acceptablez = kFALSE;
4551ea7c 2153 }
69b55c55 2154 }
2155 }
4551ea7c 2156 if (!acceptablez) {
69b55c55 2157 fitterT2.FixParameter(3,zmf);
2158 fitterT2.FixParameter(4,dzmf);
2159 fitterT2.Eval();
2160 fitterT2.ReleaseParameter(3);
2161 fitterT2.ReleaseParameter(4);
2162 rpolz0 = fitterT2.GetParameter(3);
2163 rpolz1 = fitterT2.GetParameter(4);
2164 }
4551ea7c 2165
2166 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
2167 Double_t chi2TC = fitterTC.GetChisquare() / Float_t(npointsT);
69b55c55 2168 Double_t polz1c = fitterTC.GetParameter(2);
4551ea7c 2169 Double_t polz0c = polz1c * xref2;
6c94f330 2170 Double_t aC = fitterTC.GetParameter(0);
2171 Double_t bC = fitterTC.GetParameter(1);
4551ea7c 2172 Double_t cC = aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
6c94f330 2173 Double_t aR = fitterT2.GetParameter(0);
2174 Double_t bR = fitterT2.GetParameter(1);
2175 Double_t dR = fitterT2.GetParameter(2);
4551ea7c 2176 Double_t cR = 1.0 + bR*bR - dR*aR;
2177 Double_t dca = 0.0;
2178 if (cR > 0.0) {
2179 dca = -dR / (TMath::Sqrt(1.0 + bR*bR - dR*aR) + TMath::Sqrt(1.0 + bR*bR));
2180 cR = aR / TMath::Sqrt(cR);
69b55c55 2181 }
4551ea7c 2182
2183 Double_t chi2ZT2 = 0.0;
2184 Double_t chi2ZTC = 0.0;
2185 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2186 if (cseed[iLayer].IsOK()) {
2187 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
2188 Double_t zTC = polz0c + polz1c * (xcl[iLayer] - xref2);
ca1e1e5b 2189 chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
2190 chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
69b55c55 2191 }
2192 }
4551ea7c 2193 chi2ZT2 /= TMath::Max((nlayers - 3.0),1.0);
2194 chi2ZTC /= TMath::Max((nlayers - 3.0),1.0);
2195
2196 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
2197 Float_t sumdaf = 0.0;
2198 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2199 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2200 sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))
2201 / cseed[iLayer].GetSigmaY2());
4551ea7c 2202 }
2203 }
2204 sumdaf /= Float_t (nlayers - 2.0);
2205
69b55c55 2206 //
4551ea7c 2207 // Likelihoods for full track
69b55c55 2208 //
4551ea7c 2209 Double_t likezf = TMath::Exp(-chi2ZF * 0.14);
2210 Double_t likechi2C = TMath::Exp(-chi2TC * 0.677);
2211 Double_t likechi2TR = TMath::Exp(-chi2TR * 0.78);
2212 Double_t likeaf = TMath::Exp(-sumdaf * 3.23);
2213 seedquality2[registered] = likezf * likechi2TR * likeaf;
33744e5d 2214
4551ea7c 2215 // Still needed ????
69b55c55 2216// Bool_t isGold = kFALSE;
2217//
6c94f330 2218// if (nlayers == 6 && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE; // gold
2219// if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE; // gold
69b55c55 2220// if (isGold &&nusedf<10){
2221// for (Int_t jLayer=0;jLayer<6;jLayer++){
6c94f330 2222// if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
69b55c55 2223// seed[index][jLayer].UseClusters(); //sign gold
2224// }
2225// }
4551ea7c 2226
2227 Int_t index0 = 0;
2228 if (!cseed[0].IsOK()) {
69b55c55 2229 index0 = 1;
4551ea7c 2230 if (!cseed[1].IsOK()) {
2231 index0 = 2;
2232 }
69b55c55 2233 }
ca1e1e5b 2234 seedparams[registered][0] = cseed[index0].GetX0();
2235 seedparams[registered][1] = cseed[index0].GetYref(0);
2236 seedparams[registered][2] = cseed[index0].GetZref(0);
c6f438c0 2237 seedparams[registered][5] = cR;
ca1e1e5b 2238 seedparams[registered][3] = cseed[index0].GetX0() * cR - TMath::Sin(TMath::ATan(cseed[0].GetYref(1)));
2239 seedparams[registered][4] = cseed[index0].GetZref(1)
2240 / TMath::Sqrt(1.0 + cseed[index0].GetYref(1) * cseed[index0].GetYref(1));
69b55c55 2241 seedparams[registered][6] = ns;
4551ea7c 2242
2243 Int_t labels[12];
2244 Int_t outlab[24];
2245 Int_t nlab = 0;
2246 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2247 if (!cseed[iLayer].IsOK()) {
2248 continue;
2249 }
ca1e1e5b 2250 if (cseed[iLayer].GetLabels(0) >= 0) {
2251 labels[nlab] = cseed[iLayer].GetLabels(0);
69b55c55 2252 nlab++;
2253 }
ca1e1e5b 2254 if (cseed[iLayer].GetLabels(1) >= 0) {
2255 labels[nlab] = cseed[iLayer].GetLabels(1);
69b55c55 2256 nlab++;
4551ea7c 2257 }
69b55c55 2258 }
2259 Freq(nlab,labels,outlab,kFALSE);
4551ea7c 2260 Int_t label = outlab[0];
2261 Int_t frequency = outlab[1];
2262 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ca1e1e5b 2263 cseed[iLayer].SetFreq(frequency);
2264 cseed[iLayer].SetC(cR);
2265 cseed[iLayer].SetCC(cC);
2266 cseed[iLayer].SetChi2(chi2TR);
2267 cseed[iLayer].SetChi2Z(chi2ZF);
69b55c55 2268 }
33744e5d 2269
2270 // Debugging print
4551ea7c 2271 if (1 || (!isFake)) {
69b55c55 2272 Float_t zvertex = GetZ();
4551ea7c 2273 TTreeSRedirector &cstream = *fDebugStreamer;
2274 if (AliTRDReconstructor::StreamLevel() > 0) {
2275 cstream << "Seeds1"
2276 << "isFake=" << isFake
2277 << "Vertex=" << zvertex
2278 << "Rieman2.=" << &rieman2
2279 << "Rieman.=" << &rieman
2280 << "Xref=" << xref
2281 << "X0=" << xcl[0]
2282 << "X1=" << xcl[1]
2283 << "X2=" << xcl[2]
2284 << "X3=" << xcl[3]
2285 << "X4=" << xcl[4]
2286 << "X5=" << xcl[5]
2287 << "Chi2R=" << chi2R
2288 << "Chi2Z=" << chi2Z
2289 << "Chi2RF=" << chi2RF // Chi2 of trackletes on full track
2290 << "Chi2ZF=" << chi2ZF // Chi2 z on tracklets on full track
2291 << "Chi2ZT2=" << chi2ZT2 // Chi2 z on tracklets on full track - rieman tilt
2292 << "Chi2ZTC=" << chi2ZTC // Chi2 z on tracklets on full track - rieman tilt const
2293 << "Chi2TR=" << chi2TR // Chi2 without vertex constrain
2294 << "Chi2TC=" << chi2TC // Chi2 with vertex constrain
2295 << "C=" << curv // Non constrained - no tilt correction
2296 << "DR=" << dR // DR parameter - tilt correction
2297 << "DCA=" << dca // DCA - tilt correction
2298 << "CR=" << cR // Non constrained curvature - tilt correction
2299 << "CC=" << cC // Constrained curvature
2300 << "Polz0=" << polz0c
2301 << "Polz1=" << polz1c
2302 << "RPolz0=" << rpolz0
2303 << "RPolz1=" << rpolz1
2304 << "Ncl=" << nclusters
2305 << "Nlayers=" << nlayers
2306 << "NUsedS=" << nusedCl
2307 << "NUsed=" << nusedf
2308 << "Findable=" << findable
2309 << "Like=" << like
2310 << "LikePrim=" << likePrim
2311 << "Likechi2C=" << likechi2C
2312 << "Likechi2TR=" << likechi2TR
2313 << "Likezf=" << likezf
2314 << "LikeF=" << seedquality2[registered]
2315 << "S0.=" << &cseed[0]
2316 << "S1.=" << &cseed[1]
2317 << "S2.=" << &cseed[2]
2318 << "S3.=" << &cseed[3]
2319 << "S4.=" << &cseed[4]
2320 << "S5.=" << &cseed[5]
2321 << "SB0.=" << &seedb[0]
2322 << "SB1.=" << &seedb[1]
2323 << "SB2.=" << &seedb[2]
2324 << "SB3.=" << &seedb[3]
2325 << "SB4.=" << &seedb[4]
2326 << "SB5.=" << &seedb[5]
2327 << "Label=" << label
2328 << "Freq=" << frequency
2329 << "sLayer=" << sLayer
2330 << "\n";
2331 }
69b55c55 2332 }
4551ea7c 2333
2334 if (registered<kMaxSeed - 1) {
69b55c55 2335 registered++;
2336 cseed = seed[registered];
2337 }
4551ea7c 2338
2339 } // End of loop over layer 1
2340
2341 } // End of loop over layer 0
2342
2343 } // End of loop over layer 3
2344
2345 } // End of loop over seeding time bins
2346
69b55c55 2347 //
4551ea7c 2348 // Choose best
69b55c55 2349 //
4551ea7c 2350
69b55c55 2351 TMath::Sort(registered,seedquality2,sort,kTRUE);
c6f438c0 2352 Bool_t signedseed[kMaxSeed];
4551ea7c 2353 for (Int_t i = 0; i < registered; i++) {
2354 signedseed[i] = kFALSE;
69b55c55 2355 }
4551ea7c 2356
2357 for (Int_t iter = 0; iter < 5; iter++) {
2358
2359 for (Int_t iseed = 0; iseed < registered; iseed++) {
2360
69b55c55 2361 Int_t index = sort[iseed];
4551ea7c 2362 if (signedseed[index]) {
2363 continue;
2364 }
69b55c55 2365 Int_t labelsall[1000];
4551ea7c 2366 Int_t nlabelsall = 0;
2367 Int_t naccepted = 0;;
2368 Int_t sLayer = seedlayer[index];
2369 Int_t ncl = 0;
2370 Int_t nused = 0;
2371 Int_t nlayers = 0;
69b55c55 2372 Int_t findable = 0;
4551ea7c 2373
2374 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2375
ca1e1e5b 2376 if (TMath::Abs(seed[index][jLayer].GetYref(0) / xcl[jLayer]) < 0.15) {
69b55c55 2377 findable++;
4551ea7c 2378 }
2379 if (seed[index][jLayer].IsOK()) {
69b55c55 2380 seed[index][jLayer].UpdateUsed();
ca1e1e5b 2381 ncl +=seed[index][jLayer].GetN2();
2382 nused +=seed[index][jLayer].GetNUsed();
69b55c55 2383 nlayers++;
4551ea7c 2384 // Cooking label
2385 for (Int_t itime = 0; itime < 25; itime++) {
ca1e1e5b 2386 if (seed[index][jLayer].IsUsable(itime)) {
69b55c55 2387 naccepted++;
4551ea7c 2388 for (Int_t ilab = 0; ilab < 3; ilab++) {
ca1e1e5b 2389 Int_t tindex = seed[index][jLayer].GetClusters(itime)->GetLabel(ilab);
4551ea7c 2390 if (tindex >= 0) {
69b55c55 2391 labelsall[nlabelsall] = tindex;
2392 nlabelsall++;
2393 }
2394 }
2395 }
2396 }
2397 }
4551ea7c 2398
69b55c55 2399 }
4551ea7c 2400
2401 if (nused > 30) {
2402 continue;
69b55c55 2403 }
4551ea7c 2404
2405 if (iter == 0) {
2406 if (nlayers < 6) {
2407 continue;
2408 }
2409 if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2410 continue; // Gold
2411 }
7ad19338 2412 }
4551ea7c 2413
2414 if (iter == 1) {
2415 if (nlayers < findable) {
2416 continue;
2417 }
2418 if (TMath::Log(0.000000001+seedquality2[index]) < -4.0) {
2419 continue;
2420 }
2421 }
2422
2423 if (iter == 2) {
2424 if ((nlayers == findable) ||
2425 (nlayers == 6)) {
2426 continue;
2427 }
2428 if (TMath::Log(0.000000001+seedquality2[index]) < -6.0) {
2429 continue;
2430 }
69b55c55 2431 }
4551ea7c 2432
2433 if (iter == 3) {
2434 if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2435 continue;
2436 }
69b55c55 2437 }
4551ea7c 2438
2439 if (iter == 4) {
2440 if (TMath::Log(0.000000001+seedquality2[index]) - nused/(nlayers-3.0) < -15.0) {
2441 continue;
2442 }
69b55c55 2443 }
4551ea7c 2444
69b55c55 2445 signedseed[index] = kTRUE;
4551ea7c 2446
2447 Int_t labels[1000];
2448 Int_t outlab[1000];
2449 Int_t nlab = 0;
2450 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2451 if (seed[index][iLayer].IsOK()) {
ca1e1e5b 2452 if (seed[index][iLayer].GetLabels(0) >= 0) {
2453 labels[nlab] = seed[index][iLayer].GetLabels(0);
69b55c55 2454 nlab++;
2455 }
ca1e1e5b 2456 if (seed[index][iLayer].GetLabels(1) >= 0) {
2457 labels[nlab] = seed[index][iLayer].GetLabels(1);
69b55c55 2458 nlab++;
4551ea7c 2459 }
2460 }
7ad19338 2461 }
69b55c55 2462 Freq(nlab,labels,outlab,kFALSE);
4551ea7c 2463 Int_t label = outlab[0];
2464 Int_t frequency = outlab[1];
69b55c55 2465 Freq(nlabelsall,labelsall,outlab,kFALSE);
4551ea7c 2466 Int_t label1 = outlab[0];
2467 Int_t label2 = outlab[2];
2468 Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
2469 Float_t ratio = Float_t(nused) / Float_t(ncl);
2470 if (ratio < 0.25) {
2471 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2472 if ((seed[index][jLayer].IsOK()) &&
ca1e1e5b 2473 (TMath::Abs(seed[index][jLayer].GetYfit(1) - seed[index][jLayer].GetYfit(1)) < 0.2)) {
4551ea7c 2474 seed[index][jLayer].UseClusters(); // Sign gold
2475 }
69b55c55 2476 }
7ad19338 2477 }
4551ea7c 2478
31fd97b2 2479 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 2480 TTreeSRedirector &cstream = *fDebugStreamer;
2481
69b55c55 2482 //
4551ea7c 2483 // Register seed
69b55c55 2484 //
4551ea7c 2485 AliTRDtrack *track = RegisterSeed(seed[index],seedparams[index]);
2486 AliTRDtrack dummy;
2487 if (!track) {
2488 track = &dummy;
2489 }
2490 else {
69b55c55 2491 AliESDtrack esdtrack;
4551ea7c 2492 esdtrack.UpdateTrackParams(track,AliESDtrack::kTRDout);
69b55c55 2493 esdtrack.SetLabel(label);
2494 esd->AddTrack(&esdtrack);
4551ea7c 2495 TTreeSRedirector &cstream = *fDebugStreamer;
2496 if (AliTRDReconstructor::StreamLevel() > 0) {
2497 cstream << "Tracks"
31fd97b2 2498 << "EventNrInFile=" << eventNrInFile
4551ea7c 2499 << "ESD.=" << &esdtrack
2500 << "trd.=" << track
2501 << "trdback.=" << track
2502 << "\n";
2503 }
7ad19338 2504 }
4551ea7c 2505
2506 if (AliTRDReconstructor::StreamLevel() > 0) {
2507 cstream << "Seeds2"
2508 << "Iter=" << iter
2509 << "Track.=" << track
2510 << "Like=" << seedquality[index]
2511 << "LikeF=" << seedquality2[index]
2512 << "S0.=" << &seed[index][0]
2513 << "S1.=" << &seed[index][1]
2514 << "S2.=" << &seed[index][2]
2515 << "S3.=" << &seed[index][3]
2516 << "S4.=" << &seed[index][4]
2517 << "S5.=" << &seed[index][5]
2518 << "Label=" << label
2519 << "Label1=" << label1
2520 << "Label2=" << label2
2521 << "FakeRatio=" << fakeratio
2522 << "Freq=" << frequency
2523 << "Ncl=" << ncl
2524 << "Nlayers=" << nlayers
2525 << "Findable=" << findable
2526 << "NUsed=" << nused
2527 << "sLayer=" << sLayer
31fd97b2 2528 << "EventNrInFile=" << eventNrInFile
4551ea7c 2529 << "\n";
2530 }
2531
2532 } // Loop: iseed
2533
2534 } // Loop: iter
2535
2536 } // End of loop over sectors
75bd7f81 2537
69b55c55 2538 delete [] pseed;
75bd7f81 2539
69b55c55 2540}
2541
5443e65e 2542//_____________________________________________________________________________
4551ea7c 2543Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *clusterTree) const
5443e65e 2544{
2545 //
a819a5f7 2546 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
2547 // from the file. The names of the cluster tree and branches
2548 // should match the ones used in AliTRDclusterizer::WriteClusters()
2549 //
75bd7f81 2550
4551ea7c 2551 Int_t nsize = Int_t(clusterTree->GetTotBytes() / (sizeof(AliTRDcluster)));
6c94f330 2552 TObjArray *clusterArray = new TObjArray(nsize+1000);
5443e65e 2553
4551ea7c 2554 TBranch *branch = clusterTree->GetBranch("TRDcluster");
c630aafd 2555 if (!branch) {
4551ea7c 2556 AliError("Can't get the branch !");
c630aafd 2557 return 1;
2558 }
029cd327 2559 branch->SetAddress(&clusterArray);
5443e65e 2560
a819a5f7 2561 // Loop through all entries in the tree
4551ea7c 2562 Int_t nEntries = (Int_t) clusterTree->GetEntries();
2563 Int_t nbytes = 0;
a819a5f7 2564 AliTRDcluster *c = 0;
a819a5f7 2565 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2566
2567 // Import the tree
4551ea7c 2568 nbytes += clusterTree->GetEvent(iEntry);
5443e65e 2569
a819a5f7 2570 // Get the number of points in the detector
029cd327 2571 Int_t nCluster = clusterArray->GetEntriesFast();
5443e65e 2572
a819a5f7 2573 // Loop through all TRD digits
2574 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
4551ea7c 2575 c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
4f1c04d3 2576 AliTRDcluster *co = c;
a819a5f7 2577 array->AddLast(co);
4f1c04d3 2578 clusterArray->RemoveAt(iCluster);
a819a5f7 2579 }
4551ea7c 2580
a819a5f7 2581 }
2582
029cd327 2583 delete clusterArray;
5443e65e 2584
c630aafd 2585 return 0;
75bd7f81 2586
a819a5f7 2587}
2588
75bd7f81 2589//_____________________________________________________________________________
4551ea7c 2590Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint &p) const
3551db50 2591{
2592 //
2593 // Get track space point with index i
2594 // Origin: C.Cheshkov
2595 //
2596
4551ea7c 2597 AliTRDcluster *cl = (AliTRDcluster *) fClusters->UncheckedAt(index);
2598 Int_t idet = cl->GetDetector();
2599 Int_t isector = fGeom->GetSector(idet);
2600 Int_t ichamber = fGeom->GetChamber(idet);
2601 Int_t iplan = fGeom->GetPlane(idet);
3551db50 2602 Double_t local[3];
4551ea7c 2603 local[0] = GetX(isector,iplan,cl->GetLocalTimeBin());
2604 local[1] = cl->GetY();
2605 local[2] = cl->GetZ();
3551db50 2606 Double_t global[3];
2607 fGeom->RotateBack(idet,local,global);
2608 p.SetXYZ(global[0],global[1],global[2]);
ae079791 2609 AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
3551db50 2610 switch (iplan) {
2611 case 0:
ae079791 2612 iLayer = AliGeomManager::kTRD1;
3551db50 2613 break;
2614 case 1:
ae079791 2615 iLayer = AliGeomManager::kTRD2;
3551db50 2616 break;
2617 case 2:
ae079791 2618 iLayer = AliGeomManager::kTRD3;
3551db50 2619 break;
2620 case 3:
ae079791 2621 iLayer = AliGeomManager::kTRD4;
3551db50 2622 break;
2623 case 4:
ae079791 2624 iLayer = AliGeomManager::kTRD5;
3551db50 2625 break;
2626 case 5:
ae079791 2627 iLayer = AliGeomManager::kTRD6;
3551db50 2628 break;
2629 };
4551ea7c 2630 Int_t modId = isector * fGeom->Ncham() + ichamber;
ae079791 2631 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,modId);
3551db50 2632 p.SetVolumeID(volid);
2633
2634 return kTRUE;
2635
2636}
2637
75bd7f81 2638//_____________________________________________________________________________
4551ea7c 2639void AliTRDtracker::CookLabel(AliKalmanTrack *pt, Float_t wrong) const
029cd327 2640{
2641 //
2642 // This cooks a label. Mmmmh, smells good...
2643 //
46d29e70 2644
4551ea7c 2645 Int_t label = 123456789;
2646 Int_t index;
2647 Int_t i;
2648 Int_t j;
2649 Int_t ncl = pt->GetNumberOfClusters();
2650
2651 const Int_t kRange = fTrSec[0]->GetOuterTimeBin() + 1;
5443e65e 2652
029cd327 2653 Bool_t labelAdded;
46d29e70 2654
6c94f330 2655 Int_t **s = new Int_t* [kRange];
4551ea7c 2656 for (i = 0; i < kRange; i++) {
d1b06c24 2657 s[i] = new Int_t[2];
2658 }
4551ea7c 2659 for (i = 0; i < kRange; i++) {
2660 s[i][0] = -1;
2661 s[i][1] = 0;
46d29e70 2662 }
2663
4551ea7c 2664 Int_t t0;
2665 Int_t t1;
2666 Int_t t2;
2667
2668 for (i = 0; i < ncl; i++) {
2669 index = pt->GetClusterIndex(i);
2670 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
6c94f330 2671 t0=c->GetLabel(0);
2672 t1=c->GetLabel(1);
2673 t2=c->GetLabel(2);
46d29e70 2674 }
2675
4551ea7c 2676 for (i = 0; i < ncl; i++) {
2677 index = pt->GetClusterIndex(i);
2678 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2679 for (Int_t k = 0; k < 3; k++) {
2680 label = c->GetLabel(k);
2681 labelAdded = kFALSE;
2682 j = 0;
46d29e70 2683 if (label >= 0) {
4551ea7c 2684 while ((!labelAdded) && (j < kRange)) {
2685 if ((s[j][0] == label) ||
2686 (s[j][1] == 0)) {
2687 s[j][0] = label;
2688 s[j][1] = s[j][1] + 1;
2689 labelAdded = kTRUE;
a9814c09 2690 }
2691 j++;
2692 }
46d29e70 2693 }
2694 }
2695 }
2696
4551ea7c 2697 Int_t max = 0;
6c94f330 2698 label = -123456789;
46d29e70 2699
4551ea7c 2700 for (i = 0; i < kRange; i++) {
2701 if (s[i][1] > max) {
2702 max = s[i][1];
2703 label = s[i][0];
46d29e70 2704 }
2705 }
5443e65e 2706
4551ea7c 2707 for (i = 0; i < kRange; i++) {
6c94f330 2708 delete []s[i];
5443e65e 2709 }
2710
6c94f330 2711 delete []s;
5443e65e 2712
4551ea7c 2713 if ((1.0 - Float_t(max)/ncl) > wrong) {
2714 label = -label;
2715 }
5443e65e 2716
2717 pt->SetLabel(label);
2718
46d29e70 2719}
2720
75bd7f81 2721//_____________________________________________________________________________
4551ea7c 2722void AliTRDtracker::UseClusters(const AliKalmanTrack *t, Int_t from) const
029cd327 2723{
2724 //
2725 // Use clusters, but don't abuse them!
2726 //
75bd7f81 2727
4551ea7c 2728 const Float_t kmaxchi2 = 18;
2729 const Float_t kmincl = 10;
2730
2731 AliTRDtrack *track = (AliTRDtrack *) t;
2732
2733 Int_t ncl = t->GetNumberOfClusters();
2734 for (Int_t i = from; i < ncl; i++) {
2735 Int_t index = t->GetClusterIndex(i);
2736 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
69b55c55 2737 Int_t iplane = fGeom->GetPlane(c->GetDetector());
27eaf44b 2738 if (track->GetTracklets(iplane).GetChi2() > kmaxchi2) {
4551ea7c 2739 continue;
2740 }
27eaf44b 2741 if (track->GetTracklets(iplane).GetN() < kmincl) {
4551ea7c 2742 continue;
2743 }
2744 if (!(c->IsUsed())) {
2745 c->Use();
2746 }
5443e65e 2747 }
5443e65e 2748
75bd7f81 2749}
5443e65e 2750
75bd7f81 2751//_____________________________________________________________________________
029cd327 2752Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
5443e65e 2753{
75bd7f81 2754 //
5443e65e 2755 // Parametrised "expected" error of the cluster reconstruction in Y
75bd7f81 2756 //
5443e65e 2757
2758 Double_t s = 0.08 * 0.08;
2759 return s;
75bd7f81 2760
5443e65e 2761}
2762
75bd7f81 2763//_____________________________________________________________________________
029cd327 2764Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
0a29d0f1 2765{
75bd7f81 2766 //
5443e65e 2767 // Parametrised "expected" error of the cluster reconstruction in Z
75bd7f81 2768 //
5443e65e 2769
4551ea7c 2770 Double_t s = 9.0 * 9.0 / 12.0;
5443e65e 2771 return s;
75bd7f81 2772
5443e65e 2773}
2774
75bd7f81 2775//_____________________________________________________________________________
029cd327 2776Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
5443e65e 2777{
2778 //
029cd327 2779 // Returns radial position which corresponds to time bin <localTB>
5443e65e 2780 // in tracking sector <sector> and plane <plane>
2781 //
2782
6c94f330 2783 Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
4551ea7c 2784 Int_t pl = fTrSec[sector]->GetLayerNumber(index);
2785
5443e65e 2786 return fTrSec[sector]->GetLayer(pl)->GetX();
2787
2788}
2789
75bd7f81 2790//_____________________________________________________________________________
2791AliTRDtracker::AliTRDpropagationLayer
2792 ::AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
2793 , Double_t radLength, Int_t tbIndex, Int_t plane)
ad670fba 2794 :fN(0)
2795 ,fSec(0)
2796 ,fClusters(NULL)
2797 ,fIndex(NULL)
2798 ,fX(x)
2799 ,fdX(dx)
2800 ,fRho(rho)
2801 ,fX0(radLength)
2802 ,fTimeBinIndex(tbIndex)
2803 ,fPlane(plane)
2804 ,fYmax(0)
2805 ,fYmaxSensitive(0)
2806 ,fHole(kFALSE)
2807 ,fHoleZc(0)
2808 ,fHoleZmax(0)
2809 ,fHoleYc(0)
2810 ,fHoleYmax(0)
2811 ,fHoleRho(0)
2812 ,fHoleX0(0)
5443e65e 2813{
0a29d0f1 2814 //
5443e65e 2815 // AliTRDpropagationLayer constructor
0a29d0f1 2816 //
46d29e70 2817
ad670fba 2818 for (Int_t i = 0; i < (Int_t) kZones; i++) {
2819 fZc[i] = 0;
2820 fZmax[i] = 0;
a819a5f7 2821 }
5443e65e 2822
ad670fba 2823 if (fTimeBinIndex >= 0) {
029cd327 2824 fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
ad670fba 2825 fIndex = new UInt_t[kMaxClusterPerTimeBin];
a819a5f7 2826 }
46d29e70 2827
ad670fba 2828 for (Int_t i = 0; i < 5; i++) {
2829 fIsHole[i] = kFALSE;
2830 }
5443e65e 2831
2832}
2833
75bd7f81 2834//_____________________________________________________________________________
2835void AliTRDtracker::AliTRDpropagationLayer
2836 ::SetHole(Double_t Zmax, Double_t Ymax, Double_t rho
2837 , Double_t radLength, Double_t Yc, Double_t Zc)
5443e65e 2838{
2839 //
2840 // Sets hole in the layer
2841 //
75bd7f81 2842
4551ea7c 2843 fHole = kTRUE;
2844 fHoleZc = Zc;
5443e65e 2845 fHoleZmax = Zmax;
4551ea7c 2846 fHoleYc = Yc;
5443e65e 2847 fHoleYmax = Ymax;
4551ea7c 2848 fHoleRho = rho;
2849 fHoleX0 = radLength;
75bd7f81 2850
5443e65e 2851}
46d29e70 2852
75bd7f81 2853//_____________________________________________________________________________
2854AliTRDtracker::AliTRDtrackingSector
ad670fba 2855 ::AliTRDtrackingSector(AliTRDgeometry *geo, Int_t gs)
2856 :fN(0)
2857 ,fGeom(geo)
2858 ,fGeomSector(gs)
5443e65e 2859{
2860 //
2861 // AliTRDtrackingSector Constructor
2862 //
75bd7f81 2863
7e448bcc 2864 AliTRDpadPlane *padPlane = 0;
2865 AliTRDpropagationLayer *ppl = 0;
a5cadd36 2866
ad670fba 2867 // Get holes description from geometry
3c625a9b 2868 Bool_t holes[AliTRDgeometry::kNcham];
ad670fba 2869 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3c625a9b 2870 holes[icham] = fGeom->IsHole(0,icham,gs);
3c625a9b 2871 }
3c625a9b 2872
ad670fba 2873 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
2874 fTimeBinIndex[i] = -1;
2875 }
5443e65e 2876
ad670fba 2877 Double_t x;
2878 Double_t dx;
2879 Double_t rho;
2880 Double_t radLength;
5443e65e 2881
ad670fba 2882 // Add layers for each of the planes
2883 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
a305677e 2884 //Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
5443e65e 2885
ad670fba 2886 const Int_t kNchambers = AliTRDgeometry::Ncham();
a305677e 2887 Int_t tbIndex;
ad670fba 2888 Double_t ymax = 0;
2889 Double_t ymaxsensitive = 0;
2890 Double_t *zc = new Double_t[kNchambers];
2891 Double_t *zmax = new Double_t[kNchambers];
3c625a9b 2892 Double_t *zmaxsensitive = new Double_t[kNchambers];
3551db50 2893
ad670fba 2894 for (Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
2895
2896 ymax = fGeom->GetChamberWidth(plane) / 2.0;
f162af62 2897 padPlane = fGeom->GetPadPlane(plane,0);
ad670fba 2898 ymaxsensitive = (padPlane->GetColSize(1) * padPlane->GetNcols() - 4.0) / 2.0;
2899
2900 for (Int_t ch = 0; ch < kNchambers; ch++) {
2901 zmax[ch] = fGeom->GetChamberLength(plane,ch) / 2.0;
2902 Float_t pad = padPlane->GetRowSize(1);
f162af62 2903 Float_t row0 = fGeom->GetRow0(plane,ch,0);
2904 Int_t nPads = fGeom->GetRowMax(plane,ch,0);
ad670fba 2905 zmaxsensitive[ch] = Float_t(nPads) * pad / 2.0;
2906 zc[ch] = -(pad * nPads) / 2.0 + row0;
5443e65e 2907 }
2908
ad670fba 2909 dx = AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
f162af62 2910 / AliTRDCommonParam::Instance()->GetSamplingFrequency();
ad670fba 2911 rho = 0.00295 * 0.85; //????
2912 radLength = 11.0;
5443e65e 2913
3551db50 2914 Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(plane);
a305677e 2915 //Double_t xbottom = x0 - dxDrift;
ad670fba 2916 //Double_t xtop = x0 + dxAmp;
2917
59393e34 2918 Int_t nTimeBins = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
ad670fba 2919 for (Int_t iTime = 0; iTime < nTimeBins; iTime++) {
2920
2921 Double_t xlayer = iTime * dx - dxAmp;
2922 //if (xlayer<0) xlayer = dxAmp / 2.0;
59393e34 2923 x = x0 - xlayer;
ad670fba 2924
2925 tbIndex = CookTimeBinIndex(plane,iTime);
2926 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex,plane);
3c625a9b 2927 ppl->SetYmax(ymax,ymaxsensitive);
ad670fba 2928 ppl->SetZ(zc,zmax,zmaxsensitive);
3c625a9b 2929 ppl->SetHoles(holes);
59393e34 2930 InsertLayer(ppl);
ad670fba 2931
5443e65e 2932 }
ad670fba 2933
5443e65e 2934 }
2935
5443e65e 2936 MapTimeBinLayers();
ad670fba 2937
029cd327 2938 delete [] zc;
2939 delete [] zmax;
4f1c04d3 2940 delete [] zmaxsensitive;
5443e65e 2941
2942}
2943
ad670fba 2944//_____________________________________________________________________________
2945AliTRDtracker::AliTRDtrackingSector
2946 ::AliTRDtrackingSector(const AliTRDtrackingSector &/*t*/)
2947 :fN(0)
2948 ,fGeom(0)
2949 ,fGeomSector(0)
2950{
2951 //
2952 // Copy constructor
2953 //
2954
2955}
2956
75bd7f81 2957//_____________________________________________________________________________
2958Int_t AliTRDtracker::AliTRDtrackingSector
2959 ::CookTimeBinIndex(Int_t plane, Int_t localTB) const
5443e65e 2960{
2961 //
6c94f330 2962 // depending on the digitization parameters calculates "global"
029cd327 2963 // time bin index for timebin <localTB> in plane <plane>
5443e65e 2964 //
59393e34 2965 //
75bd7f81 2966
59393e34 2967 Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
4551ea7c 2968 Int_t gtb = (plane+1) * tbPerPlane - localTB - 1;
2969 if (localTB < 0) {
2970 return -1;
2971 }
2972 if (gtb < 0) {
2973 return -1;
2974 }
75bd7f81 2975
5443e65e 2976 return gtb;
5443e65e 2977
75bd7f81 2978}
5443e65e 2979
75bd7f81 2980//_____________________________________________________________________________
2981void AliTRDtracker::AliTRDtrackingSector
2982 ::MapTimeBinLayers()
5443e65e 2983{
2984 //
2985 // For all sensitive time bins sets corresponding layer index
2986 // in the array fTimeBins
2987 //
2988
2989 Int_t index;
2990
4551ea7c 2991 for (Int_t i = 0; i < fN; i++) {
2992
5443e65e 2993 index = fLayers[i]->GetTimeBinIndex();
6c94f330 2994
4551ea7c 2995 if (index < 0) {
2996 continue;
2997 }
2998 if (index >= (Int_t) kMaxTimeBinIndex) {
2999 //AliWarning(Form("Index %d exceeds allowed maximum of %d!\n"
3000 // ,index,kMaxTimeBinIndex-1));
5443e65e 3001 continue;
3002 }
4551ea7c 3003
5443e65e 3004 fTimeBinIndex[index] = i;
4551ea7c 3005
5443e65e 3006 }
5443e65e 3007
75bd7f81 3008}
5443e65e 3009
75bd7f81 3010//_____________________________________________________________________________
3011Int_t AliTRDtracker::AliTRDtrackingSector
3012 ::GetLayerNumber(Double_t x) const
5443e65e 3013{
3014 //
3015 // Returns the number of time bin which in radial position is closest to <x>
3016 //
3017
4551ea7c 3018 if (x >= fLayers[fN-1]->GetX()) {
3019 return fN - 1;
3020 }
3021 if (x <= fLayers[ 0]->GetX()) {
3022 return 0;
3023 }
3024
3025 Int_t b = 0;
3026 Int_t e = fN - 1;
3027 Int_t m = (b + e) / 2;
ad670fba 3028
4551ea7c 3029 for ( ; b < e; m = (b + e) / 2) {
3030 if (x > fLayers[m]->GetX()) {
3031 b = m + 1;
3032 }
3033 else {
3034 e = m;
3035 }
5443e65e 3036 }
75bd7f81 3037
4551ea7c 3038 if (TMath::Abs(x - fLayers[m]->GetX()) > TMath::Abs(x - fLayers[m+1]->GetX())) {
3039 return m + 1;
3040 }
3041 else {
3042 return m;
3043 }
5443e65e 3044
3045}
3046
75bd7f81 3047//_____________________________________________________________________________
3048Int_t AliTRDtracker::AliTRDtrackingSector
3049 ::GetInnerTimeBin() const
5443e65e 3050{
3051 //
3052 // Returns number of the innermost SENSITIVE propagation layer
3053 //
3054
3055 return GetLayerNumber(0);
5443e65e 3056
75bd7f81 3057}
5443e65e 3058
75bd7f81 3059//_____________________________________________________________________________
3060Int_t AliTRDtracker::AliTRDtrackingSector
3061 ::GetOuterTimeBin() const
5443e65e 3062{
3063 //
3064 // Returns number of the outermost SENSITIVE time bin
3065 //
3066
3067 return GetLayerNumber(GetNumberOfTimeBins() - 1);
46d29e70 3068
75bd7f81 3069}
5443e65e 3070
75bd7f81 3071//_____________________________________________________________________________
3072Int_t AliTRDtracker::AliTRDtrackingSector
3073 ::GetNumberOfTimeBins() const
5443e65e 3074{
3075 //
3076 // Returns number of SENSITIVE time bins
3077 //
3078
4551ea7c 3079 Int_t tb;
3080 Int_t layer;
3081
3082 for (tb = kMaxTimeBinIndex - 1; tb >= 0; tb--) {
5443e65e 3083 layer = GetLayerNumber(tb);
4551ea7c 3084 if (layer >= 0) {
3085 break;
3086 }
5443e65e 3087 }
75bd7f81 3088
4551ea7c 3089 return tb + 1;
5443e65e 3090
75bd7f81 3091}
5443e65e 3092
75bd7f81 3093//_____________________________________________________________________________
3094void AliTRDtracker::AliTRDtrackingSector
4551ea7c 3095 ::InsertLayer(AliTRDpropagationLayer *pl)
5443e65e 3096{
3097 //
3098 // Insert layer <pl> in fLayers array.
3099 // Layers are sorted according to X coordinate.
75bd7f81 3100 //
5443e65e 3101
4551ea7c 3102 if (fN == ((Int_t) kMaxLayersPerSector)) {
3103 //AliWarning("Too many layers !\n");
3104 return;
3105 }
3106
3107 if (fN == 0) {
3108 fLayers[fN++] = pl;
ad670fba 3109 return;
3110 }
3111
4551ea7c 3112 Int_t i = Find(pl->GetX());
3113
3114 memmove(fLayers+i+1,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
3115
3116 fLayers[i] = pl;
3117 fN++;
5443e65e 3118
3119}
3120
75bd7f81 3121//_____________________________________________________________________________
3122Int_t AliTRDtracker::AliTRDtrackingSector
3123 ::Find(Double_t x) const
5443e65e 3124{
3125 //
3126 // Returns index of the propagation layer nearest to X
3127 //
3128
4551ea7c 3129 if (x <= fLayers[0]->GetX()) {
3130 return 0;
3131 }
3132
3133 if (x > fLayers[fN-1]->GetX()) {
3134 return fN;
3135 }
3136
3137 Int_t b = 0;
3138 Int_t e = fN-1;
3139 Int_t m = (b + e) / 2;
3140
3141 for (; b < e; m = (b + e) / 2) {
3142 if (x > fLayers[m]->GetX()) {
3143 b = m + 1;
3144 }
3145 else {
3146 e = m;
3147 }
5443e65e 3148 }
7ad19338 3149
75bd7f81 3150 return m;
7ad19338 3151
75bd7f81 3152}
7ad19338 3153
75bd7f81 3154//_____________________________________________________________________________
3155void AliTRDtracker::AliTRDpropagationLayer
4551ea7c 3156 ::SetZ(Double_t *center, Double_t *w, Double_t *wsensitive )
3c625a9b 3157{
3158 //
6c94f330 3159 // set centers and the width of sectors
75bd7f81 3160 //
3161
4551ea7c 3162 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3163 fZc[icham] = center[icham];
3164 fZmax[icham] = w[icham];
3c625a9b 3165 fZmaxSensitive[icham] = wsensitive[icham];
3c625a9b 3166 }
75bd7f81 3167
3c625a9b 3168}
5443e65e 3169
75bd7f81 3170//_____________________________________________________________________________
3c625a9b 3171void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
3172{
3173 //
6c94f330 3174 // set centers and the width of sectors
75bd7f81 3175 //
3176
3c625a9b 3177 fHole = kFALSE;
4551ea7c 3178
3179 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3c625a9b 3180 fIsHole[icham] = holes[icham];
4551ea7c 3181 if (holes[icham]) {
3182 fHole = kTRUE;
3183 }
3c625a9b 3184 }
5443e65e 3185
75bd7f81 3186}
5443e65e 3187
75bd7f81 3188//_____________________________________________________________________________
3189void AliTRDtracker::AliTRDpropagationLayer
4551ea7c 3190 ::InsertCluster(AliTRDcluster *c, UInt_t index)
75bd7f81 3191{
3192 //
3193 // Insert cluster in cluster array.
3194 // Clusters are sorted according to Y coordinate.
3195 //
5443e65e 3196
4551ea7c 3197 if (fTimeBinIndex < 0) {
3198 //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");
3199 return;
3200 }
3201
3202 if (fN == (Int_t) kMaxClusterPerTimeBin) {
3203 //AliWarning("Too many clusters !\n");
5443e65e 3204 return;
3205 }
ad670fba 3206
4551ea7c 3207 if (fN == 0) {
3208 fIndex[0] = index;
3209 fClusters[fN++] = c;
ad670fba 3210 return;
3211 }
4551ea7c 3212
3213 Int_t i = Find(c->GetY());
3214 memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
3215 memmove(fIndex +i+1,fIndex +i,(fN-i)*sizeof(UInt_t));
3216 fIndex[i] = index;
3217 fClusters[i] = c;
3218 fN++;
5443e65e 3219
75bd7f81 3220}
5443e65e 3221
75bd7f81 3222//_____________________________________________________________________________
3223Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const
3224{
3225 //
3226 // Returns index of the cluster nearest in Y
3227 //
5443e65e 3228
4551ea7c 3229 if (fN <= 0) {
3230 return 0;
3231 }
3232 if (y <= fClusters[0]->GetY()) {
3233 return 0;
3234 }
3235 if (y > fClusters[fN-1]->GetY()) {
3236 return fN;
3237 }
3238
3239 Int_t b = 0;
3240 Int_t e = fN - 1;
3241 Int_t m = (b + e) / 2;
3242
3243 for ( ; b < e; m = (b + e) / 2) {
3244 if (y > fClusters[m]->GetY()) {
3245 b = m + 1;
3246 }
3247 else {
3248 e = m;
3249 }
5443e65e 3250 }
75bd7f81 3251
5443e65e 3252 return m;
75bd7f81 3253
5443e65e 3254}
3255
75bd7f81 3256//_____________________________________________________________________________
3257Int_t AliTRDtracker::AliTRDpropagationLayer
3258 ::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad
3259 , Float_t maxroadz) const
7ad19338 3260{
3261 //
3262 // Returns index of the cluster nearest to the given y,z
3263 //
75bd7f81 3264
4551ea7c 3265 Int_t index = -1;
3266 Int_t maxn = fN;
69b55c55 3267 Float_t mindist = maxroad;
4551ea7c 3268
3269 for (Int_t i = Find(y-maxroad); i < maxn; i++) {
3270 AliTRDcluster *c = (AliTRDcluster *) (fClusters[i]);
69b55c55 3271 Float_t ycl = c->GetY();
4551ea7c 3272 if (ycl > (y + maxroad)) {
3273 break;
3274 }
3275 if (TMath::Abs(c->GetZ() - z) > maxroadz) {
3276 continue;
3277 }
3278 if (TMath::Abs(ycl - y) < mindist) {
3279 mindist = TMath::Abs(ycl - y);
3280 index = fIndex[i];
3281 }
6c94f330 3282 }
75bd7f81 3283
7ad19338 3284 return index;
7ad19338 3285
75bd7f81 3286}
7ad19338 3287
75bd7f81 3288//_____________________________________________________________________________
4551ea7c 3289Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c)
75bd7f81 3290{
3291 //
3292 // Returns correction factor for tilted pads geometry
3293 //
5443e65e 3294
4551ea7c 3295 Int_t det = c->GetDetector();
3296 Int_t plane = fGeom->GetPlane(det);
f162af62 3297 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(plane,0);
4551ea7c 3298 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
b8dc2353 3299
4551ea7c 3300 if (fNoTilt) {
3301 h01 = 0;
3302 }
75bd7f81 3303
fd621f36 3304 return h01;
5443e65e 3305
75bd7f81 3306}
c630aafd 3307
3308
75bd7f81 3309//_____________________________________________________________________________
3310Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
4551ea7c 3311 , AliTRDtrack *track
3312 , Int_t *clusters, AliTRDtracklet &tracklet)
4f1c04d3 3313{
6c94f330 3314 //
4f1c04d3 3315 //
75bd7f81 3316 // Try to find nearest clusters to the track in timebins from t0 to t1
4f1c04d3 3317 //
4551ea7c 3318 // Correction coeficients - depend on TRD parameters - to be changed accordingly
3319 //
3320
3321 Double_t x[100];
3322 Double_t yt[100];
3323 Double_t zt[100];
3324 Double_t xmean = 0.0; // Reference x
3325 Double_t dz[10][100];
3326 Double_t dy[10][100];
3327 Float_t zmean[100];
3328 Float_t nmean[100];
3329 Int_t clfound = 0;
3330 Int_t indexes[10][100]; // Indexes of the clusters in the road
3331 Int_t best[10][100]; // Index of best matching cluster
3332 AliTRDcluster *cl[10][100]; // Pointers to the clusters in the road
3333
3334 for (Int_t it = 0; it < 100; it++) {
3335 x[it] = 0.0;
3336 yt[it] = 0.0;
3337 zt[it] = 0.0;
3338 clusters[it] = -2;
3339 zmean[it] = 0.0;
3340 nmean[it] = 0.0;
3341 for (Int_t ih = 0; ih < 10;ih++) {
3342 indexes[ih][it] = -2; // Reset indexes1
3343 cl[ih][it] = 0;
3344 dz[ih][it] = -100.0;
3345 dy[ih][it] = -100.0;
3346 best[ih][it] = 0;
3347 }
3348 }
ad670fba 3349
4551ea7c 3350 Double_t x0 = track->GetX();
3351 Double_t sigmaz = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
3352 Int_t nall = 0;
3353 Int_t nfound = 0;
3354 Double_t h01 = 0.0;
3355 Int_t plane = -1;
3356 Int_t detector = -1;
3357 Float_t padlength = 0.0;
3358 AliTRDtrack track2(* track);
3359 Float_t snpy = track->GetSnp();
3360 Float_t tany = TMath::Sqrt(snpy*snpy / (1.0 - snpy*snpy));
3361 if (snpy < 0.0) {
3362 tany *= -1.0;
3363 }
ad670fba 3364
6c23ffed 3365 Double_t sy2 = ExpectedSigmaY2(x0,track->GetTgl(),track->GetSignedPt());
4551ea7c 3366 Double_t sz2 = ExpectedSigmaZ2(x0,track->GetTgl());
3367 Double_t road = 15.0 * TMath::Sqrt(track->GetSigmaY2() + sy2);
3368 if (road > 6.0) {
3369 road = 6.0;
3370 }
7e448bcc 3371 //road = 20.0;
4551ea7c 3372
3373 for (Int_t it = 0; it < t1-t0; it++) {
3374
3375 Double_t maxChi2[2] = { fgkMaxChi2, fgkMaxChi2 };
3376 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(it+t0));
3377 if (timeBin == 0) {
3378 continue; // No indexes1
6c94f330 3379 }
4551ea7c 3380
3381 Int_t maxn = timeBin;
3382 x[it] = timeBin.GetX();
7ad19338 3383 track2.PropagateTo(x[it]);
6c94f330 3384 yt[it] = track2.GetY();
3385 zt[it] = track2.GetZ();
7ad19338 3386
4551ea7c 3387 Double_t y = yt[it];
3388 Double_t z = zt[it];
3389 Double_t chi2 = 1000000.0;
4f1c04d3 3390 nall++;
4551ea7c 3391
4f1c04d3 3392 //
4551ea7c 3393 // Find 2 nearest cluster at given time bin
6c94f330 3394 //
7e448bcc 3395 int checkPoint[4] = {0,0,0,0};
3396 double minY = 123456789;
3397 double minD[2] = {1,1};
3398
4551ea7c 3399 for (Int_t i = timeBin.Find(y - road); i < maxn; i++) {
7e448bcc 3400 //for (Int_t i = 0; i < maxn; i++) {
4551ea7c 3401
3402 AliTRDcluster *c = (AliTRDcluster *) (timeBin[i]);
7ad19338 3403 h01 = GetTiltFactor(c);
4551ea7c 3404 if (plane < 0) {
c6f438c0 3405 Int_t det = c->GetDetector();
4551ea7c 3406 plane = fGeom->GetPlane(det);
3407 padlength = TMath::Sqrt(c->GetSigmaZ2() * 12.0);
3408 }
7e448bcc 3409
4551ea7c 3410 //if (c->GetLocalTimeBin()==0) continue;
3411 if (c->GetY() > (y + road)) {
3412 break;
3413 }
7e448bcc 3414
3415 fHDeltaX->Fill(c->GetX() - x[it]);
3416 //printf("%f\t%f\t%f \n", c->GetX(), x[it], c->GetX()-x[it]);
3417
3418 if (TMath::Abs(c->GetY()-y) < TMath::Abs(minY)) {
3419 minY = c->GetY()-y;
3420 minD[0] = c->GetY()-y;
3421 minD[1] = c->GetZ()-z;
3422 }
3423
3424 checkPoint[0]++;
3425
3426 fHMinZ->Fill(c->GetZ() - z);
3427 if ((c->GetZ() - z) * (c->GetZ() - z) > 2 * (12.0 * sz2)) {
4551ea7c 3428 continue;
7ad19338 3429 }
7e448bcc 3430 checkPoint[1]++;
7ad19338 3431
4551ea7c 3432 Double_t dist = TMath::Abs(c->GetZ() - z);
7e448bcc 3433 if (dist > (0.5 * padlength + 6.0 * sigmaz)) { // 0.5
4551ea7c 3434 continue; // 6 sigma boundary cut
3435 }
7e448bcc 3436 checkPoint[2]++;
3437
4551ea7c 3438 Double_t cost = 0.0;
3439 // Sigma boundary cost function
3440 if (dist> (0.5 * padlength - sigmaz)){
3441 cost = (dist - 0.5*padlength) / (2.0 * sigmaz);
3442 if (cost > -1) {
3443 cost = (cost + 1.0) * (cost + 1.0);
3444 }
3445 else {
3446 cost = 0.0;
3447 }
3448 }
3449 //Int_t label = TMath::Abs(track->GetLabel());
3450 //if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue;
3451 chi2 = track2.GetPredictedChi2(c,h01) + cost;
7ad19338 3452 clfound++;
6c94f330 3453
4551ea7c 3454 if (chi2 > maxChi2[1]) {
3455 continue;
3456 }
7e448bcc 3457 checkPoint[3]++;
3458
4551ea7c 3459 detector = c->GetDetector();
3460 // Store the clusters in the road
3461 for (Int_t ih = 2; ih < 9; ih++) {
3462 if (cl[ih][it] == 0) {
3463 cl[ih][it] = c;
3464 indexes[ih][it] = timeBin.GetIndex(i); // Index - 9 - reserved for outliers
7ad19338 3465 break;
3466 }
4f1c04d3 3467 }
4551ea7c 3468
3469 if (chi2 < maxChi2[0]) {
7ad19338 3470 maxChi2[1] = maxChi2[0];
3471 maxChi2[0] = chi2;
3472 indexes[1][it] = indexes[0][it];
3473 cl[1][it] = cl[0][it];
3474 indexes[0][it] = timeBin.GetIndex(i);
3475 cl[0][it] = c;
3476 continue;
3477 }
4551ea7c 3478 maxChi2[1] = chi2;
3479 cl[1][it] = c;
3480 indexes[1][it] = timeBin.GetIndex(i);
3481
3482 }
7e448bcc 3483
3484 for(int iCheckPoint = 0; iCheckPoint<4; iCheckPoint++)
3485 fHFindCl[iCheckPoint]->Fill(checkPoint[iCheckPoint]);
3486
3487 if (checkPoint[3]) {
6c23ffed 3488 if (track->GetSignedPt() > 0) fHMinYPos->Fill(minY);
7e448bcc 3489 else fHMinYNeg->Fill(minY);
3490
3491 fHMinD->Fill(minD[0], minD[1]);
3492 }
4551ea7c 3493
3494 if (cl[0][it]) {
7ad19338 3495 nfound++;
3496 xmean += x[it];
3497 }
4551ea7c 3498
4f1c04d3 3499 }
4551ea7c 3500
3501 if (nfound < 4) {
3502 return 0;
3503 }
3504 xmean /= Float_t(nfound); // Middle x
3505 track2.PropagateTo(xmean); // Propagate track to the center
3506
4f1c04d3 3507 //
4551ea7c 3508 // Choose one of the variants
7ad19338 3509 //
4551ea7c 3510 Int_t changes[10];
3511 Float_t sumz = 0.0;
3512 Float_t sum = 0.0;
3513 Double_t sumdy = 0.0;
3514 Double_t sumdy2 = 0.0;
3515 Double_t sumx = 0.0;
3516 Double_t sumxy = 0.0;
3517 Double_t sumx2 = 0.0;
3518 Double_t mpads = 0.0;
3519
3520 Int_t ngood[10];
3521 Int_t nbad[10];
3522
7ad19338 3523 Double_t meanz[10];
4551ea7c 3524 Double_t moffset[10]; // Mean offset
3525 Double_t mean[10]; // Mean value
3526 Double_t angle[10]; // Angle
3527
3528 Double_t smoffset[10]; // Sigma of mean offset
3529 Double_t smean[10]; // Sigma of mean value
3530 Double_t sangle[10]; // Sigma of angle
3531 Double_t smeanangle[10]; // Correlation
3532
7ad19338 3533 Double_t sigmas[10];
4551ea7c 3534 Double_t tchi2s[10]; // Chi2s for tracklet
ad670fba 3535
4551ea7c 3536 for (Int_t it = 0; it < 10; it++) {
ad670fba 3537
4551ea7c 3538 ngood[it] = 0;
3539 nbad[it] = 0;
3540
3541 meanz[it] = 0.0;
3542 moffset[it] = 0.0; // Mean offset
3543 mean[it] = 0.0; // Mean value
3544 angle[it] = 0.0; // Angle
3545
7e448bcc 3546 smoffset[it] = 1.0e5; // Sigma of mean offset
3547 smean[it] = 1.0e5; // Sigma of mean value
3548 sangle[it] = 1.0e5; // Sigma of angle
4551ea7c 3549 smeanangle[it] = 0.0; // Correlation
3550
7e448bcc 3551 sigmas[it] = 1.0e5;
3552 tchi2s[it] = 1.0e5; // Chi2s for tracklet
75bd7f81 3553
3554 }
3555
7ad19338 3556 //
4551ea7c 3557 // Calculate zmean
7ad19338 3558 //
4551ea7c 3559 for (Int_t it = 0; it < t1 - t0; it++) {
3560 if (!cl[0][it]) {
3561 continue;
3562 }
3563 for (Int_t dt = -3; dt <= 3; dt++) {
3564 if (it+dt < 0) {
3565 continue;
3566 }
3567 if (it+dt > t1-t0) {
3568 continue;
3569 }
3570 if (!cl[0][it+dt]) {
3571 continue;
3572 }
3573 zmean[it] += cl[0][it+dt]->GetZ();
3574 nmean[it] += 1.0;
7ad19338 3575 }
4551ea7c 3576 zmean[it] /= nmean[it];
7ad19338 3577 }
4551ea7c 3578
3579 for (Int_t it = 0; it < t1 - t0; it++) {
3580
3581 best[0][it] = 0;
3582
3583 for (Int_t ih = 0; ih < 10; ih++) {
3584 dz[ih][it] = -100.0;
3585 dy[ih][it] = -100.0;
3586 if (!cl[ih][it]) {
3587 continue;
3588 }
3589 Double_t xcluster = cl[ih][it]->GetX();
3590 Double_t ytrack;
3591 Double_t ztrack;
3592 track2.GetProlongation(xcluster,ytrack,ztrack );
3593 dz[ih][it] = cl[ih][it]->GetZ()- ztrack; // Calculate distance from track in z
3594 dy[ih][it] = cl[ih][it]->GetY() + dz[ih][it]*h01 - ytrack; // and in y
7ad19338 3595 }
4551ea7c 3596
3597 // Minimize changes
3598 if (!cl[0][it]) {
3599 continue;
3600 }
3601 if ((TMath::Abs(cl[0][it]->GetZ()-zmean[it]) > padlength * 0.8) &&
3602 (cl[1][it])) {
3603 if (TMath::Abs(cl[1][it]->GetZ()-zmean[it]) < padlength * 0.5) {
3604 best[0][it] = 1;
4f1c04d3 3605 }
4551ea7c 3606 }
3607
7ad19338 3608 }
4551ea7c 3609
7ad19338 3610 //
4551ea7c 3611 // Iterative choice of "best path"
6c94f330 3612 //
4551ea7c 3613 Int_t label = TMath::Abs(track->GetLabel());
3614 Int_t bestiter = 0;
3615
3616 for (Int_t iter = 0; iter < 9; iter++) {
3617
3618 changes[iter] = 0;
3619 sumz = 0;
3620 sum = 0;
3621 sumdy = 0;
3622 sumdy2 = 0;
3623 sumx = 0;
3624 sumx2 = 0;
3625 sumxy = 0;
3626 mpads = 0;
3627 ngood[iter] = 0;
3628 nbad[iter] = 0;
3629
3630 // Linear fit
3631 for (Int_t it = 0; it < t1 - t0; it++) {
3632
3633 if (!cl[best[iter][it]][it]) {
3634 continue;
3635 }
3636
3637 // Calculates pad-row changes
3638 Double_t zbefore = cl[best[iter][it]][it]->GetZ();
3639 Double_t zafter = cl[best[iter][it]][it]->GetZ();
3640 for (Int_t itd = it - 1; itd >= 0; itd--) {
7ad19338 3641 if (cl[best[iter][itd]][itd]) {
4551ea7c 3642 zbefore = cl[best[iter][itd]][itd]->GetZ();
7ad19338 3643 break;
3644 }
3645 }
4551ea7c 3646 for (Int_t itd = it + 1; itd < t1 - t0; itd++) {
7ad19338 3647 if (cl[best[iter][itd]][itd]) {
4551ea7c 3648 zafter = cl[best[iter][itd]][itd]->GetZ();
7ad19338 3649 break;
3650 }
3651 }
4551ea7c 3652 if ((TMath::Abs(cl[best[iter][it]][it]->GetZ()-zbefore) > 0.1) &&
3653 (TMath::Abs(cl[best[iter][it]][it]->GetZ()- zafter) > 0.1)) {
3654 changes[iter]++;
3655 }
3656
3657 Double_t dx = x[it]-xmean; // Distance to reference x
3658 sumz += cl[best[iter][it]][it]->GetZ();
4f1c04d3 3659 sum++;
4551ea7c 3660 sumdy += dy[best[iter][it]][it];
3661 sumdy2 += dy[best[iter][it]][it]*dy[best[iter][it]][it];
3662 sumx += dx;
3663 sumx2 += dx*dx;
7ad19338 3664 sumxy += dx*dy[best[iter][it]][it];
4551ea7c 3665 mpads += cl[best[iter][it]][it]->GetNPads();
3666 if ((cl[best[iter][it]][it]->GetLabel(0) == label) ||
3667 (cl[best[iter][it]][it]->GetLabel(1) == label) ||
3668 (cl[best[iter][it]][it]->GetLabel(2) == label)) {
7ad19338 3669 ngood[iter]++;
4f1c04d3 3670 }
4551ea7c 3671 else {
7ad19338 3672 nbad[iter]++;
4f1c04d3 3673 }
4551ea7c 3674
4f1c04d3 3675 }
4551ea7c 3676
7ad19338 3677 //
6c94f330 3678 // calculates line parameters
7ad19338 3679 //
4551ea7c 3680 Double_t det = sum*sumx2 - sumx*sumx;
3681 angle[iter] = (sum*sumxy - sumx*sumdy) / det;
3682 mean[iter] = (sumx2*sumdy - sumx*sumxy) / det;
3683 meanz[iter] = sumz / sum;
3684 moffset[iter] = sumdy / sum;
3685 mpads /= sum; // Mean number of pads
3686
3687 Double_t sigma2 = 0.0; // Normalized residuals - for line fit
3688 Double_t sigma1 = 0.0; // Normalized residuals - constant fit
3689
3690 for (Int_t it = 0; it < t1 - t0; it++) {
3691 if (!cl[best[iter][it]][it]) {
3692 continue;
3693 }
3694 Double_t dx = x[it] - xmean;
3695 Double_t ytr = mean[iter] + angle[iter] * dx;
3696 sigma2 += (dy[best[iter][it]][it] - ytr)
3697 * (dy[best[iter][it]][it] - ytr);
3698 sigma1 += (dy[best[iter][it]][it] - moffset[iter])
3699 * (dy[best[iter][it]][it] - moffset[iter]);
7ad19338 3700 sum++;
4f1c04d3 3701 }
4551ea7c 3702 sigma2 /= (sum - 2); // Normalized residuals
3703 sigma1 /= (sum - 1); // Normalized residuals
3704 smean[iter] = sigma2 * (sumx2 / det); // Estimated error2 of mean
3705 sangle[iter] = sigma2 * ( sum / det); // Estimated error2 of angle
3706 smeanangle[iter] = sigma2 * (-sumx / det); // Correlation
3707 sigmas[iter] = TMath::Sqrt(sigma1);
3708 smoffset[iter] = (sigma1 / sum) + 0.01*0.01; // Sigma of mean offset + unisochronity sigma
3709
6c94f330 3710 //
4551ea7c 3711 // Iterative choice of "better path"
6c94f330 3712 //
4551ea7c 3713 for (Int_t it = 0; it < t1 - t0; it++) {
3714
3715 if (!cl[best[iter][it]][it]) {
3716 continue;
3717 }
3718
3719 // Add unisochronity + angular effect contribution
3720 Double_t sigmatr2 = smoffset[iter] + 0.5*tany*tany;
3721 Double_t sweight = 1.0/sigmatr2 + 1.0/track->GetSigmaY2();
3722 Double_t weighty = (moffset[iter] / sigmatr2) / sweight; // Weighted mean
3723 Double_t sigmacl = TMath::Sqrt(sigma1*sigma1 + track->GetSigmaY2());
3724 Double_t mindist = 100000.0;
3725 Int_t ihbest = 0;
3726
3727 for (Int_t ih = 0; ih < 10; ih++) {
3728 if (!cl[ih][it]) {
3729 break;
3730 }
3731 Double_t dist2 = (dy[ih][it] - weighty) / sigmacl;
3732 dist2 *= dist2; // Chi2 distance
3733 if (dist2 < mindist) {
7ad19338 3734 mindist = dist2;
4551ea7c 3735 ihbest = ih;
7ad19338 3736 }
3737 }
4551ea7c 3738 best[iter+1][it] = ihbest;
4f1c04d3 3739 }
4551ea7c 3740
4f1c04d3 3741 //
4551ea7c 3742 // Update best hypothesy if better chi2 according tracklet position and angle
7ad19338 3743 //
69b55c55 3744 Double_t sy2 = smean[iter] + track->GetSigmaY2();
4551ea7c 3745 Double_t sa2 = sangle[iter] + track->GetSigmaSnp2(); // track->fCee;
3746 Double_t say = track->GetSigmaSnpY(); // track->fCey;
3747 //Double_t chi20 = mean[bestiter]*mean[bestiter ] / sy2+angle[bestiter]*angle[bestiter]/sa2;
3748 //Double_t chi21 = mean[iter]*mean[iter] / sy2+angle[iter]*angle[iter]/sa2;
6c94f330 3749
4551ea7c 3750 Double_t detchi = sy2*sa2 - say*say;
3751 Double_t invers[3] = {sa2/detchi,sy2/detchi,-say/detchi}; // Inverse value of covariance matrix
4f1c04d3 3752
4551ea7c 3753 Double_t chi20 = mean[bestiter] * mean[bestiter] * invers[0]
3754 + angle[bestiter] * angle[bestiter] * invers[1]
3755 + 2.0 * mean[bestiter] * angle[bestiter] * invers[2];
3756 Double_t chi21 = mean[iter] * mean[iter] * invers[0]
3757 + angle[iter] * angle[iter] * invers[1]
3758 + 2.0 * mean[iter] * angle[iter] * invers[2];
3759 tchi2s[iter] = chi21;
3760
3761 if ((changes[iter] <= changes[bestiter]) &&
3762 (chi21 < chi20)) {
3763 bestiter = iter;
7ad19338 3764 }
4551ea7c 3765
7ad19338 3766 }
4551ea7c 3767
7ad19338 3768 //
4551ea7c 3769 // Set clusters
7ad19338 3770 //
4551ea7c 3771 Double_t sigma2 = sigmas[0]; // Choose as sigma from 0 iteration
3772 Short_t maxpos = -1;
3773 Float_t maxcharge = 0.0;
3774 Short_t maxpos4 = -1;
3775 Float_t maxcharge4 = 0.0;
3776 Short_t maxpos5 = -1;
3777 Float_t maxcharge5 = 0.0;
8979685e 3778
7ad19338 3779 //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
3780 //if (tchi2s[bestiter]>25.) sigma2=1000.; // dont'accept
3781
4e009ce4 3782 Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
3783 ,-AliTracker::GetBz()*0.1);
4551ea7c 3784 Double_t expectederr = sigma2*sigma2 + 0.01*0.01;
3785 if (mpads > 3.5) {
3786 expectederr += (mpads - 3.5) * 0.04;
3787 }
3788 if (changes[bestiter] > 1) {
3789 expectederr += changes[bestiter] * 0.01;
3790 }
3791 expectederr += (0.03 * (tany-exB)*(tany-exB)) * 15.0;
3792 //if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.;
7ad19338 3793 //expectederr+=10000;
4551ea7c 3794
3795 for (Int_t it = 0; it < t1 - t0; it++) {
3796
3797 if (!cl[best[bestiter][it]][it]) {
3798 continue;
69b55c55 3799 }
4551ea7c 3800
3801 cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); // Set cluster error
3802 if (!cl[best[bestiter][it]][it]->IsUsed()) {
3803 cl[best[bestiter][it]][it]->SetY(cl[best[bestiter][it]][it]->GetY());
3804 //cl[best[bestiter][it]][it]->Use();
3805 }
3806
3807 // Time bins with maximal charge
3808 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge) {
69b55c55 3809 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3810 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 3811 }
6c94f330 3812
4551ea7c 3813 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
3814 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
69b55c55 3815 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3816 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 3817 }
3818 }
4551ea7c 3819
3820 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
3821 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
69b55c55 3822 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3823 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 3824 }
7ad19338 3825 }
4551ea7c 3826
3827 // Time bins with maximal charge
3828 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge) {
8979685e 3829 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3830 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 3831 }
6c94f330 3832
4551ea7c 3833 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
3834 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
8979685e 3835 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3836 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 3837 }
3838 }
4551ea7c 3839
3840 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
3841 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
8979685e 3842 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3843 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 3844 }
3845 }
4551ea7c 3846
7ad19338 3847 clusters[it+t0] = indexes[best[bestiter][it]][it];
4551ea7c 3848
3849 // Still needed ????
3850 //if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>4 &&
3851 //cl[best[bestiter][it]][it]->GetLocalTimeBin()<18) clusters[it+t0]
3852 // = indexes[best[bestiter][it]][it]; //Test
3853
7ad19338 3854 }
4551ea7c 3855
7ad19338 3856 //
4551ea7c 3857 // Set tracklet parameters
6c94f330 3858 //
4551ea7c 3859 Double_t trackleterr2 = smoffset[bestiter] + 0.01*0.01;
3860 if (mpads > 3.5) {
3861 trackleterr2 += (mpads - 3.5) * 0.04;
3862 }
3863 trackleterr2 += changes[bestiter] * 0.01;
3864 trackleterr2 *= TMath::Max(14.0 - nfound,1.0);
3865 trackleterr2 += 0.2 * (tany-exB)*(tany-exB);
3866
3867 // Set tracklet parameters
3868 tracklet.Set(xmean
3869 ,track2.GetY() + moffset[bestiter]
3870 ,meanz[bestiter]
3871 ,track2.GetAlpha()
3872 ,trackleterr2);
7ad19338 3873 tracklet.SetTilt(h01);
3874 tracklet.SetP0(mean[bestiter]);
3875 tracklet.SetP1(angle[bestiter]);
3876 tracklet.SetN(nfound);
3877 tracklet.SetNCross(changes[bestiter]);
3878 tracklet.SetPlane(plane);
3879 tracklet.SetSigma2(expectederr);
3880 tracklet.SetChi2(tchi2s[bestiter]);
8979685e 3881 tracklet.SetMaxPos(maxpos,maxpos4,maxpos5);
27eaf44b 3882 track->SetTracklets(plane,tracklet);
3883 track->SetNWrong(track->GetNWrong() + nbad[0]);
4551ea7c 3884
7ad19338 3885 //
3886 // Debuging part
3887 //
69b55c55 3888 TClonesArray array0("AliTRDcluster");
3889 TClonesArray array1("AliTRDcluster");
4551ea7c 3890 array0.ExpandCreateFast(t1 - t0 + 1);
3891 array1.ExpandCreateFast(t1 - t0 + 1);
3892 TTreeSRedirector &cstream = *fDebugStreamer;
7ad19338 3893 AliTRDcluster dummy;
3894 Double_t dy0[100];
8979685e 3895 Double_t dyb[100];
3896
4551ea7c 3897 for (Int_t it = 0; it < t1 - t0; it++) {
7ad19338 3898 dy0[it] = dy[0][it];
3899 dyb[it] = dy[best[bestiter][it]][it];
4551ea7c 3900 if (cl[0][it]) {
7ad19338 3901 new(array0[it]) AliTRDcluster(*cl[0][it]);
3902 }
4551ea7c 3903 else {
7ad19338 3904 new(array0[it]) AliTRDcluster(dummy);
3905 }
3906 if(cl[best[bestiter][it]][it]) {
3907 new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
3908 }
3909 else{
3910 new(array1[it]) AliTRDcluster(dummy);
3911 }
4f1c04d3 3912 }
7ad19338 3913 TGraph graph0(t1-t0,x,dy0);
3914 TGraph graph1(t1-t0,x,dyb);
3915 TGraph graphy(t1-t0,x,yt);
3916 TGraph graphz(t1-t0,x,zt);
4551ea7c 3917
3918 if (AliTRDReconstructor::StreamLevel() > 0) {
3919 cstream << "tracklet"
3920 << "track.=" << track // Track parameters
3921 << "tany=" << tany // Tangent of the local track angle
3922 << "xmean=" << xmean // Xmean - reference x of tracklet
3923 << "tilt=" << h01 // Tilt angle
3924 << "nall=" << nall // Number of foundable clusters
3925 << "nfound=" << nfound // Number of found clusters
3926 << "clfound=" << clfound // Total number of found clusters in road
3927 << "mpads=" << mpads // Mean number of pads per cluster
3928 << "plane=" << plane // Plane number
3929 << "detector=" << detector // Detector number
3930 << "road=" << road // The width of the used road
3931 << "graph0.=" << &graph0 // x - y = dy for closest cluster
3932 << "graph1.=" << &graph1 // x - y = dy for second closest cluster
3933 << "graphy.=" << &graphy // y position of the track
3934 << "graphz.=" << &graphz // z position of the track
3935 //<< "fCl.=" << &array0 // closest cluster
3936 //<< "fCl2.=" << &array1 // second closest cluster
3937 << "maxpos=" << maxpos // Maximal charge postion
3938 << "maxcharge=" << maxcharge // Maximal charge
3939 << "maxpos4=" << maxpos4 // Maximal charge postion - after bin 4
3940 << "maxcharge4=" << maxcharge4 // Maximal charge - after bin 4
3941 << "maxpos5=" << maxpos5 // Maximal charge postion - after bin 5
3942 << "maxcharge5=" << maxcharge5 // Maximal charge - after bin 5
3943 << "bestiter=" << bestiter // Best iteration number
3944 << "tracklet.=" << &tracklet // Corrspond to the best iteration
3945 << "tchi20=" << tchi2s[0] // Chi2 of cluster in the 0 iteration
3946 << "tchi2b=" << tchi2s[bestiter] // Chi2 of cluster in the best iteration
3947 << "sigmas0=" << sigmas[0] // Residuals sigma
3948 << "sigmasb=" << sigmas[bestiter] // Residulas sigma
3949 << "ngood0=" << ngood[0] // Number of good clusters in 0 iteration
3950 << "nbad0=" << nbad[0] // Number of bad clusters in 0 iteration
3951 << "ngoodb=" << ngood[bestiter] // in best iteration
3952 << "nbadb=" << nbad[bestiter] // in best iteration
3953 << "changes0=" << changes[0] // Changes of pardrows in iteration number 0
3954 << "changesb=" << changes[bestiter] // Changes of pardrows in best iteration
3955 << "moffset0=" << moffset[0] // Offset fixing angle in iter=0
3956 << "smoffset0=" << smoffset[0] // Sigma of offset fixing angle in iter=0
3957 << "moffsetb=" << moffset[bestiter] // Offset fixing angle in iter=best
3958 << "smoffsetb=" << smoffset[bestiter] // Sigma of offset fixing angle in iter=best
3959 << "mean0=" << mean[0] // Mean dy in iter=0;
3960 << "smean0=" << smean[0] // Sigma of mean dy in iter=0
3961 << "meanb=" << mean[bestiter] // Mean dy in iter=best
3962 << "smeanb=" << smean[bestiter] // Sigma of mean dy in iter=best
3963 << "angle0=" << angle[0] // Angle deviation in the iteration number 0
3964 << "sangle0=" << sangle[0] // Sigma of angular deviation in iteration number 0
3965 << "angleb=" << angle[bestiter] // Angle deviation in the best iteration
3966 << "sangleb=" << sangle[bestiter] // Sigma of angle deviation in the best iteration
3967 << "expectederr=" << expectederr // Expected error of cluster position
3968 << "\n";
3969 }
3970
4f1c04d3 3971 return nfound;
4f1c04d3 3972
75bd7f81 3973}
4f1c04d3 3974
75bd7f81 3975//_____________________________________________________________________________
3976Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist
3977 , Int_t *outlist, Bool_t down)
69b55c55 3978{
3979 //
4551ea7c 3980 // Sort eleements according occurancy
3981 // The size of output array has is 2*n
69b55c55 3982 //
75bd7f81 3983
a4d2a801 3984 if (n <= 0)
3985 return 0;
3986
4551ea7c 3987 Int_t *sindexS = new Int_t[n]; // Temporary array for sorting
3988 Int_t *sindexF = new Int_t[2*n];
3989 for (Int_t i = 0; i < n; i++) {
3990 sindexF[i] = 0;
3991 }
3992
3993 TMath::Sort(n,inlist,sindexS,down);
3994
3995 Int_t last = inlist[sindexS[0]];
3996 Int_t val = last;
3997 sindexF[0] = 1;
3998 sindexF[0+n] = last;
3999 Int_t countPos = 0;
4000
4001 // Find frequency
4002 for (Int_t i = 1; i < n; i++) {
69b55c55 4003 val = inlist[sindexS[i]];
4551ea7c 4004 if (last == val) {
4005 sindexF[countPos]++;
4006 }
4007 else {
69b55c55 4008 countPos++;
4009 sindexF[countPos+n] = val;
4010 sindexF[countPos]++;
4551ea7c 4011 last = val;
69b55c55 4012 }
4013 }
4551ea7c 4014 if (last == val) {
4015 countPos++;
4016 }
4017
4018 // Sort according frequency
4019 TMath::Sort(countPos,sindexF,sindexS,kTRUE);
4020
4021 for (Int_t i = 0; i < countPos; i++) {
69b55c55 4022 outlist[2*i ] = sindexF[sindexS[i]+n];
4023 outlist[2*i+1] = sindexF[sindexS[i]];
4024 }
4551ea7c 4025
69b55c55 4026 delete [] sindexS;
4027 delete [] sindexF;
4028
4029 return countPos;
75bd7f81 4030
69b55c55 4031}
4032
75bd7f81 4033//_____________________________________________________________________________
4551ea7c 4034AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params)
69b55c55 4035{
4036 //
75bd7f81 4037 // Register a seed
69b55c55 4038 //
4551ea7c 4039
4040 Double_t alpha = AliTRDgeometry::GetAlpha();
4041 Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
69b55c55 4042 Double_t c[15];
4551ea7c 4043
4044 c[ 0] = 0.2;
4045 c[ 1] = 0.0; c[ 2] = 2.0;
4046 c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
4047 c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0; c[ 9] = 0.1;
4048 c[10] = 0.0; c[11] = 0.0; c[12] = 0.0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
4049
4050 Int_t index = 0;
4051 AliTRDcluster *cl = 0;
4052
4053 for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
4054 if (seeds[ilayer].IsOK()) {
4055 for (Int_t itime = 22; itime > 0; itime--) {
ca1e1e5b 4056 if (seeds[ilayer].GetIndexes(itime) > 0) {
4057 index = seeds[ilayer].GetIndexes(itime);
4058 cl = seeds[ilayer].GetClusters(itime);
69b55c55 4059 break;
4060 }
4061 }
4062 }
4551ea7c 4063 if (index > 0) {
4064 break;
4065 }
69b55c55 4066 }
4551ea7c 4067 if (cl == 0) {
4068 return 0;
4069 }
4070
4071 AliTRDtrack *track = new AliTRDtrack(cl
4072 ,index
4073 ,&params[1]
4074 ,c
4075 ,params[0]
4076 ,params[6]*alpha+shift);
43bfc8af 4077 // SetCluster(cl, 0); // A. Bercuci
4078 track->PropagateTo(params[0]-5.0);
69b55c55 4079 track->ResetCovariance(1);
4551ea7c 4080
4081 Int_t rc = FollowBackProlongation(*track);
4082 if (rc < 30) {
69b55c55 4083 delete track;
4551ea7c 4084 track = 0;
4085 }
4086 else {
69b55c55 4087 track->CookdEdx();
43bfc8af 4088 track->CookdEdxTimBin();
4551ea7c 4089 CookLabel(track,0.9);
69b55c55 4090 }
69b55c55 4091
75bd7f81 4092 return track;
7e448bcc 4093}
4094
4095//////////////////////////////////////////////////////////////////////////////////////////
69b55c55 4096
7e448bcc 4097void AliTRDtracker::InitLogHists() {
4098
4099 fHBackfit = new TH1D("logTRD_backfit", "", 40, -0.5, 39.5);
4100 fHRefit = new TH1D("logTRD_refit", "", 40, -0.5, 39.5);
4101 fHClSearch = new TH1D("logTRD_clSearch", "", 60, -0.5, 59.5);
4102
4103 fHX = new TH1D("logTRD_X", ";x (cm)", 200, 50, 400);
4104 fHNCl = new TH1D("logTRD_ncl", "", 40, -0.5, 39.5);
4105 fHNClTrack = new TH1D("logTRD_nclTrack", "", 180, -0.5, 179.5);
4106
4107 fHMinYPos = new TH1D("logTRD_minYPos", ";#delta Y (cm)", 400, -6, 6);
4108 fHMinYNeg = new TH1D("logTRD_minYNeg", ";#delta Y (cm)", 400, -6, 6);
4109 fHMinZ = new TH1D("logTRD_minZ", ";#delta Z (cm)", 400, -20, 20);
4110 fHMinD = new TH2D("logTRD_minD", ";#delta Y (cm);#delta Z (cm)", 100, -6, 6, 100, -50, 50);
4111
4112 fHDeltaX = new TH1D("logTRD_deltaX", ";#delta X (cm)", 100, -5, 5);
4113 fHXCl = new TH1D("logTRD_xCl", ";cluster x position (cm)", 1000, 280, 380);
4114
4115 const char *nameFindCl[4] = {"logTRD_clY", "logTRD_clZ", "logTRD_clB", "logTRD_clG"};
4116
4117 for(int i=0; i<4; i++) {
4118 fHFindCl[i] = new TH1D(nameFindCl[i], "", 30, -0.5, 29.5);
4119 }
69b55c55 4120}
7e448bcc 4121
4122//////////////////////////////////////////////////////////////////////////////////////////
4123
4124void AliTRDtracker::SaveLogHists() {
4125
4126 TDirectory *sav = gDirectory;
4127 TFile *logFile = 0;
4128
4129 TSeqCollection *col = gROOT->GetListOfFiles();
4130 int N = col->GetEntries();
4131 for(int i=0; i<N; i++) {
4132 logFile = (TFile*)col->At(i);
4133 if (strstr(logFile->GetName(), "AliESDs.root")) break;
4134 }
4135
4136 logFile->cd();
4137 fHBackfit->Write(fHBackfit->GetName(), TObject::kOverwrite);
4138 fHRefit->Write(fHRefit->GetName(), TObject::kOverwrite);
4139 fHClSearch->Write(fHClSearch->GetName(), TObject::kOverwrite);
4140 fHX->Write(fHX->GetName(), TObject::kOverwrite);
4141 fHNCl->Write(fHNCl->GetName(), TObject::kOverwrite);
4142 fHNClTrack->Write(fHNClTrack->GetName(), TObject::kOverwrite);
4143
4144 fHMinYPos->Write(fHMinYPos->GetName(), TObject::kOverwrite);
4145 fHMinYNeg->Write(fHMinYNeg->GetName(), TObject::kOverwrite);
4146 fHMinD->Write(fHMinD->GetName(), TObject::kOverwrite);
4147 fHMinZ->Write(fHMinZ->GetName(), TObject::kOverwrite);
4148
4149 fHDeltaX->Write(fHDeltaX->GetName(), TObject::kOverwrite);
4150 fHXCl->Write(fHXCl->GetName(), TObject::kOverwrite);
4151
4152
4153 for(int i=0; i<4; i++)
4154 fHFindCl[i]->Write(fHFindCl[i]->GetName(), TObject::kOverwrite);
4155
4156 logFile->Flush();
4157
4158 sav->cd();
4159}
4160
4161//////////////////////////////////////////////////////////////////////////////////////////