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