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