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