Mods due to recently changed CDB interface
[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();
a819a5f7 167
4551ea7c 168 fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
46d29e70 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
5443e65e 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++) {
c630aafd 631
4551ea7c 632 //AliESDtrack *seed = event->GetTrack(i);
633 AliESDtrack *seed = event->GetTrack(index[i]);
634
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
705 if ((track->StatusForTOF() > 0) &&
706 (track->fNCross == 0) &&
707 (Float_t(track->fN)/Float_t(track->fNExpected) > 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());
8979685e 806 seed->SetTRDBudget(track->fBudget[0]);
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++;
1008 t.fNExpected++;
4551ea7c 1009 if (t.GetX() > 345.0) {
1010 t.fNExpectedLast++;
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) {
8979685e 1050 t.fNLast++;
4551ea7c 1051 t.fChi2Last += 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 }
1102 if (calibra->GetMItracking()) {
1103 calibra->Resettrack();
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;
1129 }
1130 if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) {
1131 break;
8979685e 1132 }
4551ea7c 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++;
1173 t.fNExpected++;
4551ea7c 1174 if (t.GetX() > 345.0) {
1175 t.fNExpectedLast++;
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) {
8979685e 1201 t.fNLast++;
4551ea7c 1202 t.fChi2Last += 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
77566f2a 1214 if (calibra->GetMItracking()) {
1215 calibra->UpdateHistograms(cl,&t);
1216 }
1217
4551ea7c 1218 // Reset material budget if 2 consecutive gold
1219 if (plane > 0) {
1220 if (t.fTracklets[plane].GetN() + t.fTracklets[plane-1].GetN() > 20) {
8979685e 1221 t.fBudget[2] = 0;
4551ea7c 1222 }
1223 }
1224
1225 }
1226
8979685e 1227 }
4551ea7c 1228
59393e34 1229 }
4551ea7c 1230
1231 ratio0 = ncl / Float_t(fTimeBinsPerPlane);
1232 Float_t ratio1 = Float_t(t.fN+1) / Float_t(t.fNExpected+1.0);
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.fN > 20)){
1240 t.MakeBackupTrack(); // Make backup of the track until is gold
59393e34 1241 }
7b580082 1242
8979685e 1243 }
46d29e70 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
46d29e70 1292 }
75bd7f81 1293
5443e65e 1294 return 1;
5443e65e 1295
59393e34 1296}
5443e65e 1297
46d29e70 1298//_____________________________________________________________________________
c630aafd 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
1354//_____________________________________________________________________________
b7a0917f 1355void AliTRDtracker::UnloadClusters()
5443e65e 1356{
0a29d0f1 1357 //
5443e65e 1358 // Clears the arrays of clusters and tracks. Resets sectors and timebins
0a29d0f1 1359 //
46d29e70 1360
4551ea7c 1361 Int_t i;
1362 Int_t nentr;
46d29e70 1363
5443e65e 1364 nentr = fClusters->GetEntriesFast();
4551ea7c 1365 for (i = 0; i < nentr; i++) {
1366 delete fClusters->RemoveAt(i);
1367 }
b7a0917f 1368 fNclusters = 0;
46d29e70 1369
5443e65e 1370 nentr = fSeeds->GetEntriesFast();
4551ea7c 1371 for (i = 0; i < nentr; i++) {
1372 delete fSeeds->RemoveAt(i);
1373 }
46d29e70 1374
5443e65e 1375 nentr = fTracks->GetEntriesFast();
4551ea7c 1376 for (i = 0; i < nentr; i++) {
1377 delete fTracks->RemoveAt(i);
1378 }
46d29e70 1379
5443e65e 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 }
46d29e70 1386
5443e65e 1387}
a819a5f7 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++) {
69b55c55 1609 cseed[sLayer+iLayer].fZref[0] = rieman.GetZat(xcl[sLayer+iLayer]);
4551ea7c 1610 chi2Z += (cseed[sLayer+iLayer].fZref[0]- zcl[sLayer+iLayer])
1611 * (cseed[sLayer+iLayer].fZref[0]- zcl[sLayer+iLayer]);
69b55c55 1612 cseed[sLayer+iLayer].fZref[1] = rieman.GetDZat(xcl[sLayer+iLayer]);
1613 cseed[sLayer+iLayer].fYref[0] = rieman.GetYat(xcl[sLayer+iLayer]);
4551ea7c 1614 chi2R += (cseed[sLayer+iLayer].fYref[0]- ycl[sLayer+iLayer])
1615 * (cseed[sLayer+iLayer].fYref[0]- ycl[sLayer+iLayer]);
69b55c55 1616 cseed[sLayer+iLayer].fYref[1] = rieman.GetDYat(xcl[sLayer+iLayer]);
1617 }
4551ea7c 1618 if (TMath::Sqrt(chi2R) > 1.0/iter) {
1619 continue;
69b55c55 1620 }
4551ea7c 1621 if (TMath::Sqrt(chi2Z) > 7.0/iter) {
1622 continue;
1623 }
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
1628 + 1.0 - cseed[sLayer+iLayer].fZref[0];
1629 if (max < minmax[1]) {
1630 minmax[1] = max;
1631 }
1632 Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer] * 0.5
1633 - 1.0 - cseed[sLayer+iLayer].fZref[0];
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
1692 cseed[sLayer+jLayer].fTilt = hL[sLayer+jLayer];
69b55c55 1693 cseed[sLayer+jLayer].fPadLength = padlength[sLayer+jLayer];
4551ea7c 1694 cseed[sLayer+jLayer].fX0 = xcl[sLayer+jLayer];
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
1718 zexp = tseed.fZref[0] + tseed.fZref[1]*dxlayer;
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 }
1725 }
1726
1727 Double_t yexp = tseed.fYref[0] + tseed.fYref[1]*dxlayer;
1728 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
1729 if (index <= 0) {
1730 continue;
69b55c55 1731 }
4551ea7c 1732 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
1733
69b55c55 1734 tseed.fIndexes[iTime] = index;
4551ea7c 1735 tseed.fClusters[iTime] = cl; // Register cluster
1736 tseed.fX[iTime] = dxlayer; // Register cluster
1737 tseed.fY[iTime] = cl->GetY(); // Register cluster
1738 tseed.fZ[iTime] = cl->GetZ(); // Register cluster
1739
1740 }
1741
69b55c55 1742 tseed.Update();
4551ea7c 1743
1744 // Count the number of clusters and distortions into quality
1745 Float_t dangle = tseed.fYfit[1] - tseed.fYref[1];
1746 Float_t tquality = (18.0 - tseed.fN2) / 2.0 + TMath::Abs(dangle) / 0.1
1747 + TMath::Abs(tseed.fYfit[0] - tseed.fYref[0]) / 0.2
1748 + 2.0 * TMath::Abs(tseed.fMeanz-tseed.fZref[0]) / padlength[jLayer];
1749 if ((iter == 0) && tseed.IsOK()) {
69b55c55 1750 cseed[sLayer+jLayer] = tseed;
4551ea7c 1751 quality = tquality;
1752 if (tquality < 5) {
1753 break;
1754 }
69b55c55 1755 }
4551ea7c 1756 if (tseed.IsOK() && (tquality < quality)) {
1757 cseed[sLayer+jLayer] = tseed;
1758 }
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();
4551ea7c 1769 nusedCl += cseed[sLayer+jLayer].fNUsed;
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()) {
1783 nclusters += cseed[sLayer+iLayer].fN2;
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]
1791 ,cseed[sLayer+iLayer].fYfitR[0]
1792 ,cseed[sLayer+iLayer].fZProb
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++) {
69b55c55 1802 cseed[sLayer+iLayer].fYref[0] = rieman.GetYat(xcl[sLayer+iLayer]);
4551ea7c 1803 chi2R += (cseed[sLayer+iLayer].fYref[0] - cseed[sLayer+iLayer].fYfitR[0])
1804 * (cseed[sLayer+iLayer].fYref[0] - cseed[sLayer+iLayer].fYfitR[0]);
69b55c55 1805 cseed[sLayer+iLayer].fYref[1] = rieman.GetDYat(xcl[sLayer+iLayer]);
1806 cseed[sLayer+iLayer].fZref[0] = rieman.GetZat(xcl[sLayer+iLayer]);
4551ea7c 1807 chi2Z += (cseed[sLayer+iLayer].fZref[0] - cseed[sLayer+iLayer].fMeanz)
1808 * (cseed[sLayer+iLayer].fZref[0]- cseed[sLayer+iLayer].fMeanz);
69b55c55 1809 cseed[sLayer+iLayer].fZref[1] = rieman.GetDZat(xcl[sLayer+iLayer]);
1810 }
1811 Double_t curv = rieman.GetC();
4551ea7c 1812
69b55c55 1813 //
4551ea7c 1814 // Likelihoods
6c94f330 1815 //
4551ea7c 1816 Double_t sumda = TMath::Abs(cseed[sLayer+0].fYfitR[1] - cseed[sLayer+0].fYref[1])
1817 + TMath::Abs(cseed[sLayer+1].fYfitR[1] - cseed[sLayer+1].fYref[1])
1818 + TMath::Abs(cseed[sLayer+2].fYfitR[1] - cseed[sLayer+2].fYref[1])
1819 + TMath::Abs(cseed[sLayer+3].fYfitR[1]- cseed[sLayer+3].fYref[1]);
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;
1828 Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].fYref[1] - 130.0*curv) * 1.9);
1829 Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].fZref[1]
1830 - cseed[sLayer+0].fZref[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 ////////////////////////////////////////////////////////////////////////////////////
69b55c55 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();
4551ea7c 1873 cseed[jLayer].fTilt = hL[jLayer];
69b55c55 1874 cseed[jLayer].fPadLength = padlength[jLayer];
4551ea7c 1875 cseed[jLayer].fX0 = xcl[jLayer];
1876 // Get pad length and rough cluster
1877 Int_t indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].fYref[0]
1878 ,cseed[jLayer].fZref[0]
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 }
1898 Float_t zexp = cseed[jLayer].fZref[0];
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];
1909 Double_t yexp = tseed.fYref[0] + tseed.fYref[1]*dxlayer;
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);
69b55c55 1916 tseed.fIndexes[iTime] = index;
4551ea7c 1917 tseed.fClusters[iTime] = cl; // Register cluster
1918 tseed.fX[iTime] = dxlayer; // Register cluster
1919 tseed.fY[iTime] = cl->GetY(); // Register cluster
1920 tseed.fZ[iTime] = cl->GetZ(); // Register cluster
1921 }
1922
69b55c55 1923 tseed.Update();
4551ea7c 1924 if (tseed.IsOK()) {
1925 Float_t dangle = tseed.fYfit[1] - tseed.fYref[1];
1926 Float_t tquality = (18.0 - tseed.fN2)/2.0 + TMath::Abs(dangle) / 0.1
1927 + TMath::Abs(tseed.fYfit[0] - tseed.fYref[0]) / 0.2
1928 + 2.0 * TMath::Abs(tseed.fMeanz-tseed.fZref[0]) / padlength[jLayer];
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();
4551ea7c 1942 nusedf += cseed[jLayer].fNUsed;
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];
4551ea7c 1967 Double_t zcor = tseed.fTilt * (tseed.fZProb - tseed.fZref[0]);
1968 Float_t dangle = tseed.fYfit[1] - tseed.fYref[1];
1969 Float_t tquality = (18.0 - tseed.fN2) / 2.0 + TMath::Abs(dangle) / 0.1
1970 + TMath::Abs(tseed.fYfit[0] - (tseed.fYref[0] - zcor)) / 0.2
1971 + 2.0 * TMath::Abs(tseed.fMeanz - tseed.fZref[0]) / padlength[jLayer];
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];
2002 Double_t zexp = tseed.fZref[0];
2003 Double_t zcor = tseed.fTilt * (tseed.fZProb - tseed.fZref[0]);
2004 Float_t roadz = padlength[bLayer] + 1;
2005 if (TMath::Abs(tseed.fZProb-zexp) > padlength[bLayer]*0.5) {
2006 roadz = padlength[bLayer] * 0.5;
2007 }
2008 if (tseed.fZfit[1]*tseed.fZref[1] < 0) {
2009 roadz = padlength[bLayer] * 0.5;
2010 }
2011 if (TMath::Abs(tseed.fZProb-zexp) < 0.1*padlength[bLayer]) {
2012 zexp = tseed.fZProb;
2013 roadz = padlength[bLayer] * 0.5;
2014 }
2015
2016 Double_t yexp = tseed.fYref[0] + tseed.fYref[1]*dxlayer-zcor;
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
69b55c55 2023 tseed.fIndexes[iTime] = index;
4551ea7c 2024 tseed.fClusters[iTime] = cl; // Register cluster
2025 tseed.fX[iTime] = dxlayer; // Register cluster
2026 tseed.fY[iTime] = cl->GetY(); // Register cluster
2027 tseed.fZ[iTime] = cl->GetZ(); // Register cluster
2028
2029 }
2030
69b55c55 2031 tseed.Update();
c6f438c0 2032 if (tseed.IsOK()) {
4551ea7c 2033 Float_t dangle = tseed.fYfit[1] - tseed.fYref[1];
2034 Double_t zcor = tseed.fTilt * (tseed.fZProb - tseed.fZref[0]);
2035 Float_t tquality = (18.0 - tseed.fN2) / 2.0
2036 + TMath::Abs(dangle) / 0.1
2037 + TMath::Abs(tseed.fYfit[0] - (tseed.fYref[0] - zcor)) / 0.2
2038 + 2.0 * TMath::Abs(tseed.fMeanz - tseed.fZref[0]) / padlength[jLayer];
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++) {
2054 if (TMath::Abs(cseed[iLayer].fYref[0] / cseed[iLayer].fX0) < 0.15) {
69b55c55 2055 findable++;
4551ea7c 2056 }
2057 if (cseed[iLayer].IsOK()) {
2058 nclusters += cseed[iLayer].fN2;
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]
2069 ,cseed[iLayer].fYfitR[0]
2070 ,cseed[iLayer].fZProb
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()) {
69b55c55 2081 cseed[iLayer].fYref[0] = rieman.GetYat(xcl[iLayer]);
4551ea7c 2082 chi2RF += (cseed[iLayer].fYref[0] - cseed[iLayer].fYfitR[0])
2083 * (cseed[iLayer].fYref[0] - cseed[iLayer].fYfitR[0]);
69b55c55 2084 cseed[iLayer].fYref[1] = rieman.GetDYat(xcl[iLayer]);
2085 cseed[iLayer].fZref[0] = rieman.GetZat(xcl[iLayer]);
4551ea7c 2086 chi2ZF += (cseed[iLayer].fZref[0] - cseed[iLayer].fMeanz)
2087 * (cseed[iLayer].fZref[0] - cseed[iLayer].fMeanz);
69b55c55 2088 cseed[iLayer].fZref[1] = rieman.GetDZat(xcl[iLayer]);
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
2115 if (!cseed[iLayer].fUsable[itime]) {
2116 continue;
2117 }
2118 // X relative to the middle chamber
2119 Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2;
2120 Double_t y = cseed[iLayer].fY[itime];
2121 Double_t z = cseed[iLayer].fZ[itime];
69b55c55 2122 // ExB correction to the correction
4551ea7c 2123 // Tilted rieman
69b55c55 2124 Double_t uvt[6];
4551ea7c 2125 // Global x
2126 Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0;
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 //
4551ea7c 2139 z = cseed[iLayer].fZ[itime];
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);
2165 if (TMath::Abs(cseed[iLayer].fZProb - 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);
2203 chi2ZT2 += TMath::Abs(cseed[iLayer].fMeanz - zT2);
2204 chi2ZTC += TMath::Abs(cseed[iLayer].fMeanz - 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()) {
2214 sumdaf += TMath::Abs((cseed[iLayer].fYfit[1] - cseed[iLayer].fYref[1])
2215 / cseed[iLayer].fSigmaY2);
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 }
2248 seedparams[registered][0] = cseed[index0].fX0;
2249 seedparams[registered][1] = cseed[index0].fYref[0];
2250 seedparams[registered][2] = cseed[index0].fZref[0];
c6f438c0 2251 seedparams[registered][5] = cR;
4551ea7c 2252 seedparams[registered][3] = cseed[index0].fX0 * cR - TMath::Sin(TMath::ATan(cseed[0].fYref[1]));
2253 seedparams[registered][4] = cseed[index0].fZref[1]
2254 / TMath::Sqrt(1.0 + cseed[index0].fYref[1] * cseed[index0].fYref[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 }
2264 if (cseed[iLayer].fLabels[0] >= 0) {
69b55c55 2265 labels[nlab] = cseed[iLayer].fLabels[0];
2266 nlab++;
2267 }
4551ea7c 2268 if (cseed[iLayer].fLabels[1] >= 0) {
69b55c55 2269 labels[nlab] = cseed[iLayer].fLabels[1];
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++) {
69b55c55 2277 cseed[iLayer].fFreq = frequency;
c6f438c0 2278 cseed[iLayer].fC = cR;
4551ea7c 2279 cseed[iLayer].fCC = cC;
69b55c55 2280 cseed[iLayer].fChi2 = chi2TR;
2281 cseed[iLayer].fChi2Z = chi2ZF;
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
2390 if (TMath::Abs(seed[index][jLayer].fYref[0] / xcl[jLayer]) < 0.15) {
69b55c55 2391 findable++;
4551ea7c 2392 }
2393 if (seed[index][jLayer].IsOK()) {
69b55c55 2394 seed[index][jLayer].UpdateUsed();
6c94f330 2395 ncl +=seed[index][jLayer].fN2;
2396 nused +=seed[index][jLayer].fNUsed;
69b55c55 2397 nlayers++;
4551ea7c 2398 // Cooking label
2399 for (Int_t itime = 0; itime < 25; itime++) {
2400 if (seed[index][jLayer].fUsable[itime]) {
69b55c55 2401 naccepted++;
4551ea7c 2402 for (Int_t ilab = 0; ilab < 3; ilab++) {
69b55c55 2403 Int_t tindex = seed[index][jLayer].fClusters[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 }
69b55c55 2435 }
4551ea7c 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 }
2457 }
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()) {
2466 if (seed[index][iLayer].fLabels[0] >= 0) {
69b55c55 2467 labels[nlab] = seed[index][iLayer].fLabels[0];
2468 nlab++;
2469 }
4551ea7c 2470 if (seed[index][iLayer].fLabels[1] >= 0) {
69b55c55 2471 labels[nlab] = seed[index][iLayer].fLabels[1];
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()) &&
2487 (TMath::Abs(seed[index][jLayer].fYfit[1] - seed[index][jLayer].fYfit[1]) < 0.2)) {
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 }
2518 }
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";
7ad19338 2544 }
4551ea7c 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
a819a5f7 2556//_____________________________________________________________________________
4551ea7c 2557Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *clusterTree) const
46d29e70 2558{
a819a5f7 2559 //
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());
4551ea7c 2752 if (track->fTracklets[iplane].GetChi2() > kmaxchi2) {
2753 continue;
2754 }
2755 if (track->fTracklets[iplane].GetN() < kmincl) {
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
5443e65e 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{
2828 //
2829 // AliTRDpropagationLayer constructor
2830 //
2831
ad670fba 2832 for (Int_t i = 0; i < (Int_t) kZones; i++) {
2833 fZc[i] = 0;
2834 fZmax[i] = 0;
5443e65e 2835 }
2836
ad670fba 2837 if (fTimeBinIndex >= 0) {
029cd327 2838 fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
ad670fba 2839 fIndex = new UInt_t[kMaxClusterPerTimeBin];
5443e65e 2840 }
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}
5443e65e 2866
75bd7f81 2867//_____________________________________________________________________________
2868AliTRDtracker::AliTRDtrackingSector
ad670fba 2869 ::AliTRDtrackingSector(AliTRDgeometry *geo, Int_t gs)
2870 :fN(0)
2871 ,fGeom(geo)
2872 ,fGeomSector(gs)
0a29d0f1 2873{
2874 //
5443e65e 2875 // AliTRDtrackingSector Constructor
0a29d0f1 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;
0a29d0f1 2962
46d29e70 2963}
2964
75bd7f81 2965//_____________________________________________________________________________
ad670fba 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
2978//_____________________________________________________________________________
75bd7f81 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) {
5443e65e 3017 continue;
3018 }
4551ea7c 3019 if (index >= (Int_t) kMaxTimeBinIndex) {
3020 //AliWarning(Form("Index %d exceeds allowed maximum of %d!\n"
3021 // ,index,kMaxTimeBinIndex-1));
3022 continue;
3023 }
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);
5443e65e 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