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