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