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