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