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