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