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