Replace AliTRDCalibra
[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
7e448bcc 29
a2cb5b3d 30#include <Riostream.h>
46d29e70 31#include <TBranch.h>
a1e17193 32#include <TFile.h>
4551ea7c 33#include <TGraph.h>
7e448bcc 34#include <TH1D.h>
35#include <TH2D.h>
a1e17193 36#include <TLinearFitter.h>
37#include <TObjArray.h>
38#include <TROOT.h>
39#include <TTree.h>
40#include <TTreeStream.h>
4551ea7c 41
42#include "AliESD.h"
43#include "AliAlignObj.h"
44#include "AliRieman.h"
45#include "AliTrackPointArray.h"
46d29e70 46
46d29e70 47#include "AliTRDgeometry.h"
a5cadd36 48#include "AliTRDpadPlane.h"
0b2318ec 49#include "AliTRDgeometry.h"
46d29e70 50#include "AliTRDcluster.h"
51#include "AliTRDtrack.h"
75bd7f81 52#include "AliTRDseed.h"
3551db50 53#include "AliTRDcalibDB.h"
54#include "AliTRDCommonParam.h"
46d29e70 55#include "AliTRDtracker.h"
7b580082 56#include "AliTRDReconstructor.h"
77566f2a 57#include "AliTRDCalibra.h"
7e448bcc 58ClassImp(AliTRDtracker)
46d29e70 59
7e448bcc 60const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5;
61const Float_t AliTRDtracker::fgkLabelFraction = 0.8; // ??
62const Double_t AliTRDtracker::fgkMaxChi2 = 12.0;
63const Double_t AliTRDtracker::fgkMaxSnp = 0.95; // Corresponds to tan = 3
64const Double_t AliTRDtracker::fgkMaxStep = 2.0; // Maximal step size in propagation
7ad19338 65
75bd7f81 66//_____________________________________________________________________________
4551ea7c 67AliTRDtracker::AliTRDtracker()
4e009ce4 68 :AliTracker()
69 ,fHBackfit(0x0)
70 ,fHClSearch(0x0)
71 ,fHRefit(0x0)
72 ,fHX(0x0)
73 ,fHNCl(0x0)
74 ,fHNClTrack(0x0)
75 ,fHMinYPos(0x0)
76 ,fHMinYNeg(0x0)
77 ,fHMinZ(0x0)
78 ,fHMinD(0x0)
79 ,fHDeltaX(0x0)
80 ,fHXCl(0x0)
81 ,fGeom(0)
82 ,fNclusters(0)
83 ,fClusters(0)
84 ,fNseeds(0)
85 ,fSeeds(0)
86 ,fNtracks(0)
87 ,fTracks(0)
88 ,fTimeBinsPerPlane(0)
89 ,fAddTRDseeds(kFALSE)
90 ,fNoTilt(kFALSE)
91 ,fDebugStreamer(0)
89f05372 92{
75bd7f81 93 //
b7a0917f 94 // Default constructor
75bd7f81 95 //
b7a0917f 96
4551ea7c 97 for (Int_t i = 0; i < kTrackingSectors; i++) {
98 fTrSec[i] = 0;
99 }
100 for (Int_t j = 0; j < 5; j++) {
101 for (Int_t k = 0; k < 18; k++) {
102 fHoles[j][k] = kFALSE;
103 }
104 }
75bd7f81 105
7e448bcc 106 InitLogHists();
4e009ce4 107
89f05372 108}
75bd7f81 109
110//_____________________________________________________________________________
6c94f330 111AliTRDtracker::AliTRDtracker(const AliTRDtracker &t)
112 :AliTracker(t)
4e009ce4 113 ,fHBackfit(0x0)
114 ,fHClSearch(0x0)
115 ,fHRefit(0x0)
116 ,fHX(0x0)
117 ,fHNCl(0x0)
118 ,fHNClTrack(0x0)
119 ,fHMinYPos(0x0)
120 ,fHMinYNeg(0x0)
121 ,fHMinZ(0x0)
122 ,fHMinD(0x0)
123 ,fHDeltaX(0x0)
124 ,fHXCl(0x0)
ad670fba 125 ,fGeom(0)
126 ,fNclusters(0)
127 ,fClusters(0)
128 ,fNseeds(0)
129 ,fSeeds(0)
130 ,fNtracks(0)
131 ,fTracks(0)
132 ,fTimeBinsPerPlane(0)
133 ,fAddTRDseeds(kFALSE)
134 ,fNoTilt(kFALSE)
135 ,fDebugStreamer(0)
6c94f330 136{
ad670fba 137 //
138 // Copy constructor
139 //
4e009ce4 140
ad670fba 141}
142
143//_____________________________________________________________________________
144AliTRDtracker::AliTRDtracker(const TFile *geomfile)
145 :AliTracker()
4e009ce4 146 ,fHBackfit(0x0)
147 ,fHClSearch(0x0)
148 ,fHRefit(0x0)
149 ,fHX(0x0)
150 ,fHNCl(0x0)
151 ,fHNClTrack(0x0)
152 ,fHMinYPos(0x0)
153 ,fHMinYNeg(0x0)
154 ,fHMinZ(0x0)
155 ,fHMinD(0x0)
156 ,fHDeltaX(0x0)
157 ,fHXCl(0x0)
ad670fba 158 ,fGeom(0)
159 ,fNclusters(0)
160 ,fClusters(new TObjArray(2000))
161 ,fNseeds(0)
162 ,fSeeds(new TObjArray(2000))
163 ,fNtracks(0)
164 ,fTracks(new TObjArray(1000))
165 ,fTimeBinsPerPlane(0)
166 ,fAddTRDseeds(kFALSE)
167 ,fNoTilt(kFALSE)
168 ,fDebugStreamer(0)
46d29e70 169{
5443e65e 170 //
171 // Main constructor
172 //
b8dc2353 173
4551ea7c 174 TDirectory *savedir = gDirectory;
175 TFile *in = (TFile *) geomfile;
176
5443e65e 177 if (!in->IsOpen()) {
4551ea7c 178 AliWarning("geometry file is not open!\n");
179 AliWarning("FULL TRD geometry and DEFAULT TRD parameter will be used\n");
5443e65e 180 }
181 else {
182 in->cd();
4551ea7c 183 fGeom = (AliTRDgeometry *) in->Get("TRDgeometry");
5443e65e 184 }
46d29e70 185
4551ea7c 186 if (!fGeom) {
187 AliWarning("Cannot find TRD geometry!\n");
0b2318ec 188 fGeom = new AliTRDgeometry();
c630aafd 189 }
c6f438c0 190 fGeom->ReadGeoMatrices();
c630aafd 191
5443e65e 192 savedir->cd();
46d29e70 193
4551ea7c 194 for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
439c63c8 195 Int_t trS = geomS;
4551ea7c 196 fTrSec[trS] = new AliTRDtrackingSector(fGeom,geomS);
197 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
198 fHoles[icham][trS] = fGeom->IsHole(0,icham,geomS);
3c625a9b 199 }
5443e65e 200 }
4551ea7c 201
3551db50 202 AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
7ad19338 203 Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
4551ea7c 204 if (tiltAngle < 0.1) {
b8dc2353 205 fNoTilt = kTRUE;
206 }
207
fc88348b 208 if (!AliTRDcalibDB::Instance()) {
209 AliFatal("Could not get calibration object");
210 }
211 fTimeBinsPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
a819a5f7 212
4551ea7c 213 fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
46d29e70 214
9c9d2487 215 savedir->cd();
7e448bcc 216
217 InitLogHists();
75bd7f81 218
5443e65e 219}
46d29e70 220
75bd7f81 221//_____________________________________________________________________________
5443e65e 222AliTRDtracker::~AliTRDtracker()
46d29e70 223{
029cd327 224 //
225 // Destructor of AliTRDtracker
226 //
227
89f05372 228 if (fClusters) {
229 fClusters->Delete();
230 delete fClusters;
231 }
4551ea7c 232
89f05372 233 if (fTracks) {
234 fTracks->Delete();
235 delete fTracks;
236 }
4551ea7c 237
89f05372 238 if (fSeeds) {
239 fSeeds->Delete();
240 delete fSeeds;
241 }
4551ea7c 242
5443e65e 243 delete fGeom;
0a29d0f1 244
4551ea7c 245 for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
029cd327 246 delete fTrSec[geomS];
5443e65e 247 }
4551ea7c 248
7ad19338 249 if (fDebugStreamer) {
7ad19338 250 delete fDebugStreamer;
251 }
9c9d2487 252
75bd7f81 253}
59393e34 254
75bd7f81 255//_____________________________________________________________________________
256Int_t AliTRDtracker::LocalToGlobalID(Int_t lid)
257{
59393e34 258 //
75bd7f81 259 // Transform internal TRD ID to global detector ID
59393e34 260 //
75bd7f81 261
4551ea7c 262 Int_t isector = fGeom->GetSector(lid);
263 Int_t ichamber = fGeom->GetChamber(lid);
264 Int_t iplan = fGeom->GetPlane(lid);
265
59393e34 266 AliAlignObj::ELayerID iLayer = AliAlignObj::kTRD1;
267 switch (iplan) {
268 case 0:
269 iLayer = AliAlignObj::kTRD1;
270 break;
271 case 1:
272 iLayer = AliAlignObj::kTRD2;
273 break;
274 case 2:
275 iLayer = AliAlignObj::kTRD3;
276 break;
277 case 3:
278 iLayer = AliAlignObj::kTRD4;
279 break;
280 case 4:
281 iLayer = AliAlignObj::kTRD5;
282 break;
283 case 5:
284 iLayer = AliAlignObj::kTRD6;
285 break;
286 };
4551ea7c 287
288 Int_t modId = isector * fGeom->Ncham() + ichamber;
59393e34 289 UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
75bd7f81 290
59393e34 291 return volid;
75bd7f81 292
59393e34 293}
294
75bd7f81 295//_____________________________________________________________________________
296Int_t AliTRDtracker::GlobalToLocalID(Int_t gid)
297{
59393e34 298 //
75bd7f81 299 // Transform global detector ID to local detector ID
59393e34 300 //
75bd7f81 301
4551ea7c 302 Int_t modId = 0;
303 AliAlignObj::ELayerID layerId = AliAlignObj::VolUIDToLayer(gid,modId);
304
305 Int_t isector = modId / fGeom->Ncham();
306 Int_t ichamber = modId % fGeom->Ncham();
307 Int_t iLayer = -1;
308
59393e34 309 switch (layerId) {
310 case AliAlignObj::kTRD1:
6c94f330 311 iLayer = 0;
59393e34 312 break;
313 case AliAlignObj::kTRD2:
6c94f330 314 iLayer = 1;
59393e34 315 break;
316 case AliAlignObj::kTRD3:
6c94f330 317 iLayer = 2;
59393e34 318 break;
319 case AliAlignObj::kTRD4:
6c94f330 320 iLayer = 3;
59393e34 321 break;
322 case AliAlignObj::kTRD5:
6c94f330 323 iLayer = 4;
59393e34 324 break;
325 case AliAlignObj::kTRD6:
6c94f330 326 iLayer = 5;
59393e34 327 break;
328 default:
6c94f330 329 iLayer =-1;
59393e34 330 }
4551ea7c 331
332 if (iLayer < 0) {
333 return -1;
334 }
335
59393e34 336 Int_t lid = fGeom->GetDetector(iLayer,ichamber,isector);
75bd7f81 337
59393e34 338 return lid;
59393e34 339
75bd7f81 340}
59393e34 341
75bd7f81 342//_____________________________________________________________________________
4551ea7c 343Bool_t AliTRDtracker::Transform(AliTRDcluster *cluster)
75bd7f81 344{
59393e34 345 //
75bd7f81 346 // Transform something ... whatever ...
c6f438c0 347 //
75bd7f81 348
33744e5d 349 // Magic constants for geo manager transformation
4551ea7c 350 const Double_t kX0shift = 2.52;
351 const Double_t kX0shift5 = 3.05;
33744e5d 352
6c94f330 353 //
33744e5d 354 // Apply alignment and calibration to transform cluster
6c94f330 355 //
356 Int_t detector = cluster->GetDetector();
4551ea7c 357 Int_t plane = fGeom->GetPlane(cluster->GetDetector());
358 Int_t chamber = fGeom->GetChamber(cluster->GetDetector());
359 Int_t sector = fGeom->GetSector(cluster->GetDetector());
360
361 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
362 Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.0); // Drift distance
c6f438c0 363
59393e34 364 //
59393e34 365 // ExB correction
366 //
367 Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
4e009ce4 368 Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift,-AliTracker::GetBz()*0.1);
4551ea7c 369
370 AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
371 AliTRDpadPlane *padPlane = commonParam->GetPadPlane(plane,chamber);
372 Double_t zshiftIdeal = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
373 Double_t localPos[3];
374 Double_t localPosTracker[3];
c6f438c0 375 localPos[0] = -cluster->GetX();
4551ea7c 376 localPos[1] = cluster->GetY() - driftX * exB;
377 localPos[2] = cluster->GetZ() - zshiftIdeal;
378
c6f438c0 379 cluster->SetY(cluster->GetY() - driftX*exB);
380 Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
381 cluster->SetX(xplane- cluster->GetX());
4551ea7c 382
383 TGeoHMatrix *matrix = fGeom->GetCorrectionMatrix(cluster->GetDetector());
384 if (!matrix) {
385 // No matrix found - if somebody used geometry with holes
c6f438c0 386 AliError("Invalid Geometry - Default Geometry used\n");
387 return kTRUE;
388 }
4551ea7c 389 matrix->LocalToMaster(localPos,localPosTracker);
390
391 if (AliTRDReconstructor::StreamLevel() > 1) {
392 (* fDebugStreamer) << "Transform"
393 << "Cl.=" << cluster
394 << "matrix.=" << matrix
395 << "Detector=" << detector
396 << "Sector=" << sector
397 << "Plane=" << plane
398 << "Chamber=" << chamber
399 << "lx0=" << localPosTracker[0]
400 << "ly0=" << localPosTracker[1]
401 << "lz0=" << localPosTracker[2]
402 << "\n";
c6f438c0 403 }
4551ea7c 404
405 if (plane == 5) {
c6f438c0 406 cluster->SetX(localPosTracker[0]+kX0shift5);
4551ea7c 407 }
408 else {
c6f438c0 409 cluster->SetX(localPosTracker[0]+kX0shift);
4551ea7c 410 }
411
c6f438c0 412 cluster->SetY(localPosTracker[1]);
413 cluster->SetZ(localPosTracker[2]);
75bd7f81 414
59393e34 415 return kTRUE;
75bd7f81 416
59393e34 417}
418
75bd7f81 419//_____________________________________________________________________________
4551ea7c 420// Bool_t AliTRDtracker::Transform(AliTRDcluster *cluster)
75bd7f81 421//{
c6f438c0 422// //
4551ea7c 423// // Is this still needed ????
c6f438c0 424// //
425// const Double_t kDriftCorrection = 1.01; // drift coeficient correction
426// const Double_t kTime0Cor = 0.32; // time0 correction
427// //
428// const Double_t kX0shift = 2.52;
429// const Double_t kX0shift5 = 3.05;
430
431// //
432// // apply alignment and calibration to transform cluster
433// //
434// //
435// Int_t detector = cluster->GetDetector();
436// Int_t plane = fGeom->GetPlane(cluster->GetDetector());
437// Int_t chamber = fGeom->GetChamber(cluster->GetDetector());
438// Int_t sector = fGeom->GetSector(cluster->GetDetector());
439
440// Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
441// Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.); // drift distance
442// //
443// // ExB correction
444// //
445// Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
4e009ce4 446// Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift,-AliTracker::GetBz()*0.1);
c6f438c0 447// //
448
449// AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
450// AliTRDpadPlane * padPlane = commonParam->GetPadPlane(plane,chamber);
451// Double_t zshiftIdeal = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
452// Double_t localPos[3], globalPos[3], localPosTracker[3], localPosTracker2[3];
453// localPos[2] = -cluster->GetX();
454// localPos[0] = cluster->GetY() - driftX*exB;
455// localPos[1] = cluster->GetZ() -zshiftIdeal;
456// TGeoHMatrix * matrix = fGeom->GetGeoMatrix(cluster->GetDetector());
457// matrix->LocalToMaster(localPos, globalPos);
458
459// Double_t sectorAngle = 20.*(sector%18)+10;
460// TGeoHMatrix rotSector;
461// rotSector.RotateZ(sectorAngle);
462// rotSector.LocalToMaster(globalPos, localPosTracker);
463// //
464// //
465// TGeoHMatrix matrix2(*matrix);
466// matrix2.MultiplyLeft(&rotSector);
467// matrix2.LocalToMaster(localPos,localPosTracker2);
468// //
469// //
470// //
471// cluster->SetY(cluster->GetY() - driftX*exB);
472// Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
473// cluster->SetX(xplane- kDriftCorrection*(cluster->GetX()-kTime0Cor));
474// (*fDebugStreamer)<<"Transform"<<
475// "Cl.="<<cluster<<
476// "matrix.="<<matrix<<
477// "matrix2.="<<&matrix2<<
478// "Detector="<<detector<<
479// "Sector="<<sector<<
480// "Plane="<<plane<<
481// "Chamber="<<chamber<<
482// "lx0="<<localPosTracker[0]<<
483// "ly0="<<localPosTracker[1]<<
484// "lz0="<<localPosTracker[2]<<
485// "lx2="<<localPosTracker2[0]<<
486// "ly2="<<localPosTracker2[1]<<
487// "lz2="<<localPosTracker2[2]<<
488// "\n";
489// //
490// if (plane==5)
491// cluster->SetX(localPosTracker[0]+kX0shift5);
492// else
493// cluster->SetX(localPosTracker[0]+kX0shift);
494
495// cluster->SetY(localPosTracker[1]);
496// cluster->SetZ(localPosTracker[2]);
497// return kTRUE;
498// }
499
75bd7f81 500//_____________________________________________________________________________
501Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track)
502{
9c9d2487 503 //
504 // Rotates the track when necessary
505 //
506
507 Double_t alpha = AliTRDgeometry::GetAlpha();
4551ea7c 508 Double_t y = track->GetY();
509 Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
9c9d2487 510
4551ea7c 511 // Is this still needed ????
c630aafd 512 //Int_t ns = AliTRDgeometry::kNsect;
9c9d2487 513 //Int_t s=Int_t(track->GetAlpha()/alpha)%ns;
514
4551ea7c 515 if (y > ymax) {
9c9d2487 516 //s = (s+1) % ns;
4551ea7c 517 if (!track->Rotate( alpha)) {
518 return kFALSE;
519 }
520 }
521 else if (y < -ymax) {
9c9d2487 522 //s = (s-1+ns) % ns;
4551ea7c 523 if (!track->Rotate(-alpha)) {
524 return kFALSE;
525 }
9c9d2487 526 }
527
528 return kTRUE;
9c9d2487 529
75bd7f81 530}
46e2d86c 531
75bd7f81 532//_____________________________________________________________________________
533AliTRDcluster *AliTRDtracker::GetCluster(AliTRDtrack *track, Int_t plane
534 , Int_t timebin, UInt_t &index)
535{
46e2d86c 536 //
75bd7f81 537 // Try to find cluster in the backup list
46e2d86c 538 //
75bd7f81 539
4551ea7c 540 AliTRDcluster *cl =0;
6c94f330 541 Int_t *indexes = track->GetBackupIndexes();
4551ea7c 542
543 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
544 if (indexes[i] == 0) {
545 break;
546 }
547 AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
548 if (!cli) {
549 break;
550 }
551 if (cli->GetLocalTimeBin() != timebin) {
552 continue;
553 }
46e2d86c 554 Int_t iplane = fGeom->GetPlane(cli->GetDetector());
4551ea7c 555 if (iplane == plane) {
556 cl = cli;
7ad19338 557 index = indexes[i];
46e2d86c 558 break;
559 }
560 }
75bd7f81 561
46e2d86c 562 return cl;
46e2d86c 563
75bd7f81 564}
3c625a9b 565
75bd7f81 566//_____________________________________________________________________________
4551ea7c 567Int_t AliTRDtracker::GetLastPlane(AliTRDtrack *track)
75bd7f81 568{
3c625a9b 569 //
75bd7f81 570 // Return last updated plane
571 //
572
4551ea7c 573 Int_t lastplane = 0;
574 Int_t *indexes = track->GetBackupIndexes();
575
576 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
577 AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
578 if (!cli) {
579 break;
580 }
3c625a9b 581 Int_t iplane = fGeom->GetPlane(cli->GetDetector());
4551ea7c 582 if (iplane > lastplane) {
3c625a9b 583 lastplane = iplane;
584 }
585 }
75bd7f81 586
3c625a9b 587 return lastplane;
75bd7f81 588
3c625a9b 589}
75bd7f81 590
591//_____________________________________________________________________________
4551ea7c 592Int_t AliTRDtracker::Clusters2Tracks(AliESD *event)
c630aafd 593{
594 //
595 // Finds tracks within the TRD. The ESD event is expected to contain seeds
596 // at the outer part of the TRD. The seeds
597 // are found within the TRD if fAddTRDseeds is TRUE.
598 // The tracks are propagated to the innermost time bin
599 // of the TRD and the ESD event is updated
600 //
601
4551ea7c 602 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
029cd327 603 Float_t foundMin = fgkMinClustersInTrack * timeBins;
4551ea7c 604 Int_t nseed = 0;
605 Int_t found = 0;
606 //Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
c630aafd 607
608 Int_t n = event->GetNumberOfTracks();
4551ea7c 609 for (Int_t i = 0; i < n; i++) {
610
611 AliESDtrack *seed = event->GetTrack(i);
612 ULong_t status = seed->GetStatus();
613 if ((status & AliESDtrack::kTRDout) == 0) {
614 continue;
615 }
616 if ((status & AliESDtrack::kTRDin) != 0) {
617 continue;
618 }
c630aafd 619 nseed++;
7ad19338 620
4551ea7c 621 AliTRDtrack *seed2 = new AliTRDtrack(*seed);
46e2d86c 622 //seed2->ResetCovariance();
4551ea7c 623 AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
624 AliTRDtrack &t = *pt;
7b580082 625 FollowProlongation(t);
c630aafd 626 if (t.GetNumberOfClusters() >= foundMin) {
627 UseClusters(&t);
4551ea7c 628 CookLabel(pt,1 - fgkLabelFraction);
629 //t.CookdEdx();
c630aafd 630 }
631 found++;
c630aafd 632
4551ea7c 633 Double_t xTPC = 250.0;
59393e34 634 if (PropagateToX(t,xTPC,fgkMaxStep)) {
c630aafd 635 seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
636 }
637 delete seed2;
638 delete pt;
4551ea7c 639
640 }
c630aafd 641
33744e5d 642 AliInfo(Form("Number of loaded seeds: %d",nseed));
643 AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
644 AliInfo(Form("Total number of found tracks: %d",found));
c630aafd 645
646 return 0;
75bd7f81 647
c630aafd 648}
5443e65e 649
5443e65e 650//_____________________________________________________________________________
4551ea7c 651Int_t AliTRDtracker::PropagateBack(AliESD *event)
75bd7f81 652{
c630aafd 653 //
654 // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
655 // backpropagated by the TPC tracker. Each seed is first propagated
656 // to the TRD, and then its prolongation is searched in the TRD.
657 // If sufficiently long continuation of the track is found in the TRD
658 // the track is updated, otherwise it's stored as originaly defined
659 // by the TPC tracker.
660 //
661
7e448bcc 662 Int_t found = 0; // number of tracks found
4551ea7c 663 Float_t foundMin = 20.0;
664 Int_t n = event->GetNumberOfTracks();
665
666 // Sort tracks
667 Float_t *quality = new Float_t[n];
668 Int_t *index = new Int_t[n];
669 for (Int_t i = 0; i < n; i++) {
670 AliESDtrack *seed = event->GetTrack(i);
4f1c04d3 671 Double_t covariance[15];
672 seed->GetExternalCovariance(covariance);
673 quality[i] = covariance[0]+covariance[2];
7e448bcc 674 //quality[i] = covariance[0];
4f1c04d3 675 }
676 TMath::Sort(n,quality,index,kFALSE);
4f1c04d3 677
4551ea7c 678 for (Int_t i = 0; i < n; i++) {
c630aafd 679
4551ea7c 680 //AliESDtrack *seed = event->GetTrack(i);
681 AliESDtrack *seed = event->GetTrack(index[i]);
7e448bcc 682 fHBackfit->Fill(0);
4551ea7c 683
684 ULong_t status = seed->GetStatus();
685 if ((status & AliESDtrack::kTPCout) == 0) {
7e448bcc 686 fHBackfit->Fill(1);
4551ea7c 687 continue;
688 }
7e448bcc 689
4551ea7c 690 if ((status & AliESDtrack::kTRDout) != 0) {
7e448bcc 691 fHBackfit->Fill(2);
4551ea7c 692 continue;
693 }
694
695 Int_t lbl = seed->GetLabel();
c630aafd 696 AliTRDtrack *track = new AliTRDtrack(*seed);
697 track->SetSeedLabel(lbl);
4551ea7c 698 seed->UpdateTrackParams(track,AliESDtrack::kTRDbackup); // Make backup
c630aafd 699 fNseeds++;
4551ea7c 700 Float_t p4 = track->GetC();
f6625211 701 Int_t expectedClr = FollowBackProlongation(*track);
7e448bcc 702
703 fHBackfit->Fill(3);
704 fHX->Fill(track->GetX());
705
706
707 // store the last measurement
708 /*
709 fHNClTrack->Fill(track->GetNumberOfClusters());
710 if (track->GetNumberOfClusters() >= foundMin) {
711
712 fHBackfit->Fill(4);
713 track->CookdEdx();
714 CookdEdxTimBin(*track);
715 CookLabel(track,1 - fgkLabelFraction);
716 if (track->GetBackupTrack()) {
717 //fHBackfit->Fill(5);
718 UseClusters(track->GetBackupTrack());
719 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
720 }
721 }
722 */
723
724 /**/
725 // inter-tracks competition ???
4551ea7c 726 if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
727 (TMath::Abs(track->GetPt()) > 0.8)) {
7e448bcc 728
729 fHBackfit->Fill(4);
730
4f1c04d3 731 //
4551ea7c 732 // Make backup for back propagation
4f1c04d3 733 //
7e448bcc 734
4f1c04d3 735 Int_t foundClr = track->GetNumberOfClusters();
736 if (foundClr >= foundMin) {
737 track->CookdEdx();
8979685e 738 CookdEdxTimBin(*track);
4551ea7c 739 CookLabel(track,1 - fgkLabelFraction);
740 if (track->GetBackupTrack()) {
741 UseClusters(track->GetBackupTrack());
742 }
743
744 // Sign only gold tracks
745 if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
7e448bcc 746 if ((seed->GetKinkIndex(0) == 0) &&
4551ea7c 747 (TMath::Abs(track->GetPt()) < 1.5)) {
748 UseClusters(track);
749 }
4f1c04d3 750 }
751 Bool_t isGold = kFALSE;
752
4551ea7c 753 // Full gold track
754 if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
755 //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
756 if (track->GetBackupTrack()) {
757 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
7e448bcc 758
4551ea7c 759 }
4f1c04d3 760 isGold = kTRUE;
7e448bcc 761 //fHBackfit->Fill()
4f1c04d3 762 }
4551ea7c 763
764 // Almost gold track
765 if ((!isGold) &&
766 (track->GetNCross() == 0) &&
767 (track->GetChi2() / track->GetNumberOfClusters() < 7)) {
768 //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
769 if (track->GetBackupTrack()) {
770 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
771 }
f4e9508c 772 isGold = kTRUE;
773 }
4551ea7c 774
775 if ((!isGold) &&
776 (track->GetBackupTrack())) {
7e448bcc 777 if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) &&
778 ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
4551ea7c 779 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
4f1c04d3 780 isGold = kTRUE;
781 }
782 }
4551ea7c 783
7e448bcc 784 if ((track->StatusForTOF() > 0) &&
785 (track->GetNCross() == 0) &&
27eaf44b 786 (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected()) > 0.4)) {
6c94f330 787 //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
7ad19338 788 }
4551ea7c 789
16d9fbba 790 }
c630aafd 791 }
7e448bcc 792 /**/
4551ea7c 793
7e448bcc 794 /**/
8979685e 795 // Debug part of tracking
4551ea7c 796 TTreeSRedirector &cstream = *fDebugStreamer;
31fd97b2 797 Int_t eventNrInFile = event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
4551ea7c 798 if (AliTRDReconstructor::StreamLevel() > 0) {
799 if (track->GetBackupTrack()) {
800 cstream << "Tracks"
31fd97b2 801 << "EventNrInFile=" << eventNrInFile
4551ea7c 802 << "ESD.=" << seed
803 << "trd.=" << track
804 << "trdback.=" << track->GetBackupTrack()
805 << "\n";
806 }
807 else {
808 cstream << "Tracks"
31fd97b2 809 << "EventNrInFile=" << eventNrInFile
4551ea7c 810 << "ESD.=" << seed
811 << "trd.=" << track
812 << "trdback.=" << track
813 << "\n";
d337ef8d 814 }
8979685e 815 }
7e448bcc 816 /**/
4551ea7c 817
818 // Propagation to the TOF (I.Belikov)
819 if (track->GetStop() == kFALSE) {
7e448bcc 820 fHBackfit->Fill(5);
821
4551ea7c 822 Double_t xtof = 371.0;
7e448bcc 823 Double_t xTOF0 = 370.0;
824
4551ea7c 825 Double_t c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
826 if (TMath::Abs(c2) >= 0.99) {
7e448bcc 827
828 fHBackfit->Fill(6);
c5a8e3df 829 delete track;
830 continue;
831 }
7e448bcc 832
59393e34 833 PropagateToX(*track,xTOF0,fgkMaxStep);
4551ea7c 834
835 // Energy losses taken to the account - check one more time
836 c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
837 if (TMath::Abs(c2) >= 0.99) {
7e448bcc 838
839 fHBackfit->Fill(7);
4f1c04d3 840 delete track;
841 continue;
842 }
7e448bcc 843
844 //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
845 // fHBackfit->Fill(7);
846 //delete track;
847 // continue;
848 //}
4f1c04d3 849
4551ea7c 850 Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
851 Double_t y;
852 track->GetYAt(xtof,GetBz(),y);
7e448bcc 853 if (y > ymax) {
4551ea7c 854 if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
7e448bcc 855 fHBackfit->Fill(8);
7ac6fa52 856 delete track;
7bed16a7 857 continue;
7ac6fa52 858 }
4551ea7c 859 }
860 else if (y < -ymax) {
7ac6fa52 861 if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
7e448bcc 862 fHBackfit->Fill(9);
7ac6fa52 863 delete track;
7bed16a7 864 continue;
7ac6fa52 865 }
3c625a9b 866 }
7e448bcc 867
3c625a9b 868 if (track->PropagateTo(xtof)) {
4551ea7c 869 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
7e448bcc 870 fHBackfit->Fill(10);
871
4551ea7c 872 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
873 for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
6d45eaef 874 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
875 }
876 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
eab5961e 877 }
4551ea7c 878 //seed->SetTRDtrack(new AliTRDtrack(*track));
879 if (track->GetNumberOfClusters() > foundMin) {
7e448bcc 880 fHBackfit->Fill(11);
4551ea7c 881 found++;
882 }
3c625a9b 883 }
4551ea7c 884
885 }
886 else {
7e448bcc 887
888 fHBackfit->Fill(12);
889
4551ea7c 890 if ((track->GetNumberOfClusters() > 15) &&
891 (track->GetNumberOfClusters() > 0.5*expectedClr)) {
7e448bcc 892
4551ea7c 893 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
7e448bcc 894 fHBackfit->Fill(13);
895
16d9fbba 896 //seed->SetStatus(AliESDtrack::kTRDStop);
4551ea7c 897 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
898 for (Int_t j = 0; j <AliESDtrack::kNSlice; j++) {
6d45eaef 899 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
900 }
901 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
eab5961e 902 }
7ad19338 903 //seed->SetTRDtrack(new AliTRDtrack(*track));
3c625a9b 904 found++;
905 }
4551ea7c 906
1e9bb598 907 }
4551ea7c 908
7ad19338 909 seed->SetTRDQuality(track->StatusForTOF());
27eaf44b 910 seed->SetTRDBudget(track->GetBudget(0));
8979685e 911
7e448bcc 912 fHBackfit->Fill(14);
d9b8978b 913 delete track;
c630aafd 914 }
ad670fba 915
33744e5d 916 AliInfo(Form("Number of seeds: %d",fNseeds));
917 AliInfo(Form("Number of back propagated TRD tracks: %d",found));
6c94f330 918
33744e5d 919 // New seeding
920 if (AliTRDReconstructor::SeedingOn()) {
4551ea7c 921 MakeSeedsMI(3,5,event);
33744e5d 922 }
923
924 fSeeds->Clear();
4551ea7c 925 fNseeds = 0;
7ad19338 926
4f1c04d3 927 delete [] index;
928 delete [] quality;
929
7e448bcc 930 SaveLogHists();
931
1e9bb598 932 return 0;
1e9bb598 933}
934
935//_____________________________________________________________________________
4551ea7c 936Int_t AliTRDtracker::RefitInward(AliESD *event)
1e9bb598 937{
938 //
939 // Refits tracks within the TRD. The ESD event is expected to contain seeds
940 // at the outer part of the TRD.
941 // The tracks are propagated to the innermost time bin
942 // of the TRD and the ESD event is updated
943 // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
944 //
945
4551ea7c 946 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
1e9bb598 947 Float_t foundMin = fgkMinClustersInTrack * timeBins;
4551ea7c 948 Int_t nseed = 0;
949 Int_t found = 0;
950 //Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
4f1c04d3 951 AliTRDtrack seed2;
6c94f330 952
1e9bb598 953 Int_t n = event->GetNumberOfTracks();
4551ea7c 954 for (Int_t i = 0; i < n; i++) {
955
956 AliESDtrack *seed = event->GetTrack(i);
957 new (&seed2) AliTRDtrack(*seed);
7e448bcc 958 fHRefit->Fill(0);
959
4551ea7c 960 if (seed2.GetX() < 270.0) {
961 seed->UpdateTrackParams(&seed2,AliESDtrack::kTRDbackup); // Backup TPC track - only update
7e448bcc 962 fHRefit->Fill(1);
f4e9508c 963 continue;
964 }
965
4551ea7c 966 ULong_t status = seed->GetStatus();
967 if ((status & AliESDtrack::kTRDout) == 0) {
7e448bcc 968 fHRefit->Fill(2);
0dd7d129 969 continue;
970 }
4551ea7c 971 if ((status & AliESDtrack::kTRDin) != 0) {
7e448bcc 972 fHRefit->Fill(3);
0dd7d129 973 continue;
974 }
7e448bcc 975
6c94f330 976 nseed++;
7e448bcc 977 fHRefit->Fill(4);
6c94f330 978
4551ea7c 979 seed2.ResetCovariance(50.0);
6c94f330 980
4f1c04d3 981 AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
4551ea7c 982 Int_t *indexes2 = seed2.GetIndexes();
983 for (Int_t i = 0; i < AliESDtrack::kNPlane;i++) {
984 for (Int_t j = 0; j < AliESDtrack::kNSlice;j++) {
6d45eaef 985 pt->SetPIDsignals(seed2.GetPIDsignals(i,j),i,j);
986 }
7ad19338 987 pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i);
988 }
eab5961e 989
4551ea7c 990 Int_t *indexes3 = pt->GetBackupIndexes();
991 for (Int_t i = 0; i < 200;i++) {
992 if (indexes2[i] == 0) {
993 break;
994 }
46e2d86c 995 indexes3[i] = indexes2[i];
4551ea7c 996 }
997
46e2d86c 998 //AliTRDtrack *pt = seed2;
4551ea7c 999 AliTRDtrack &t = *pt;
7b580082 1000 FollowProlongation(t);
1e9bb598 1001 if (t.GetNumberOfClusters() >= foundMin) {
4551ea7c 1002 //UseClusters(&t);
6c94f330 1003 //CookLabel(pt, 1-fgkLabelFraction);
7ad19338 1004 t.CookdEdx();
1005 CookdEdxTimBin(t);
1e9bb598 1006 }
1007 found++;
33744e5d 1008
4551ea7c 1009 Double_t xTPC = 250.0;
1010 if (PropagateToX(t,xTPC,fgkMaxStep)) {
7e448bcc 1011
4551ea7c 1012 seed->UpdateTrackParams(pt,AliESDtrack::kTRDrefit);
7e448bcc 1013 fHRefit->Fill(5);
1014
4551ea7c 1015 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
1016 for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
6d45eaef 1017 seed->SetTRDsignals(pt->GetPIDsignals(i,j),i,j);
1018 }
7ad19338 1019 seed->SetTRDTimBin(pt->GetPIDTimBin(i),i);
1020 }
4551ea7c 1021 }
1022 else {
1023 // If not prolongation to TPC - propagate without update
7e448bcc 1024 fHRefit->Fill(5);
4551ea7c 1025 AliTRDtrack *seed2 = new AliTRDtrack(*seed);
1026 seed2->ResetCovariance(5.0);
1027 AliTRDtrack *pt2 = new AliTRDtrack(*seed2,seed2->GetAlpha());
7bed16a7 1028 delete seed2;
59393e34 1029 if (PropagateToX(*pt2,xTPC,fgkMaxStep)) {
4551ea7c 1030 pt2->CookdEdx( );
eab5961e 1031 CookdEdxTimBin(*pt2);
4551ea7c 1032 seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit);
7e448bcc 1033 fHRefit->Fill(6);
1034
4551ea7c 1035 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
1036 for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
6d45eaef 1037 seed->SetTRDsignals(pt2->GetPIDsignals(i,j),i,j);
1038 }
7ad19338 1039 seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
1040 }
eab5961e 1041 }
7bed16a7 1042 delete pt2;
4551ea7c 1043 }
1044
1e9bb598 1045 delete pt;
4551ea7c 1046
eab5961e 1047 }
1e9bb598 1048
33744e5d 1049 AliInfo(Form("Number of loaded seeds: %d",nseed));
1050 AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
1e9bb598 1051
7e448bcc 1052 SaveLogHists();
c630aafd 1053 return 0;
c630aafd 1054}
1055
75bd7f81 1056//_____________________________________________________________________________
4551ea7c 1057Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t)
8979685e 1058{
75bd7f81 1059 //
8979685e 1060 // Starting from current position on track=t this function tries
1061 // to extrapolate the track up to timeBin=0 and to confirm prolongation
1062 // if a close cluster is found. Returns the number of clusters
1063 // expected to be found in sensitive layers
1064 // GeoManager used to estimate mean density
75bd7f81 1065 //
1066
4551ea7c 1067 Int_t sector;
1068 Int_t lastplane = GetLastPlane(&t);
3551db50 1069 Double_t radLength = 0.0;
4551ea7c 1070 Double_t rho = 0.0;
1071 Int_t expectedNumberOfClusters = 0;
1072
1073 for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
1074
1075 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
1076 Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
1077
8979685e 1078 //
4551ea7c 1079 // Propagate track close to the plane if neccessary
7b580082 1080 //
1081 Double_t currentx = fTrSec[0]->GetLayer(rowlast)->GetX();
4551ea7c 1082 if (currentx < (-fgkMaxStep + t.GetX())) {
1083 // Propagate closer to chamber - safety space fgkMaxStep
1084 if (!PropagateToX(t,currentx+fgkMaxStep,fgkMaxStep)) {
1085 break;
1086 }
1087 }
1088
1089 if (!AdjustSector(&t)) {
1090 break;
8979685e 1091 }
4551ea7c 1092
59393e34 1093 //
4551ea7c 1094 // Get material budget
8979685e 1095 //
4551ea7c 1096 Double_t xyz0[3];
1097 Double_t xyz1[3];
1098 Double_t param[7];
1099 Double_t x;
1100 Double_t y;
1101 Double_t z;
1102
1103 // Starting global position
1104 t.GetXYZ(xyz0);
1105 // End global position
7b580082 1106 x = fTrSec[0]->GetLayer(row0)->GetX();
4551ea7c 1107 if (!t.GetProlongation(x,y,z)) {
1108 break;
1109 }
1110 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1111 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1112 xyz1[2] = z;
7b580082 1113 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
4551ea7c 1114 rho = param[0];
1115 radLength = param[1]; // Get mean propagation parameters
1116
8979685e 1117 //
4551ea7c 1118 // Propagate and update
8979685e 1119 //
8979685e 1120 sector = t.GetSector();
4551ea7c 1121 //for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
1122 for (Int_t itime = 0 ; itime < GetTimeBinsPerPlane(); itime++) {
1123
1124 Int_t ilayer = GetGlobalTimeBin(0,iplane,itime);
8979685e 1125 expectedNumberOfClusters++;
27eaf44b 1126 t.SetNExpected(t.GetNExpected() + 1);
4551ea7c 1127 if (t.GetX() > 345.0) {
27eaf44b 1128 t.SetNExpectedLast(t.GetNExpectedLast() + 1);
4551ea7c 1129 }
1130 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
7e448bcc 1131 AliTRDcluster *cl = 0;
4551ea7c 1132 UInt_t index = 0;
1133 Double_t maxChi2 = fgkMaxChi2;
8979685e 1134 x = timeBin.GetX();
4551ea7c 1135
8979685e 1136 if (timeBin) {
4551ea7c 1137
1138 AliTRDcluster *cl0 = timeBin[0];
1139 if (!cl0) {
1140 // No clusters in given time bin
1141 continue;
1142 }
1143
1144 Int_t plane = fGeom->GetPlane(cl0->GetDetector());
1145 if (plane > lastplane) {
1146 continue;
1147 }
1148
8979685e 1149 Int_t timebin = cl0->GetLocalTimeBin();
4551ea7c 1150 AliTRDcluster *cl2 = GetCluster(&t,plane,timebin,index);
1151
8979685e 1152 if (cl2) {
4551ea7c 1153
1154 cl = cl2;
33744e5d 1155 //Double_t h01 = GetTiltFactor(cl); //I.B's fix
6c94f330 1156 //maxChi2=t.GetPredictedChi2(cl,h01);
4551ea7c 1157
7b580082 1158 }
8979685e 1159 if (cl) {
4551ea7c 1160
1161 //if (cl->GetNPads()<5)
59393e34 1162 Double_t dxsample = timeBin.GetdX();
1163 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
4551ea7c 1164 Double_t h01 = GetTiltFactor(cl);
1165 Int_t det = cl->GetDetector();
1166 Int_t plane = fGeom->GetPlane(det);
1167 if (t.GetX() > 345.0) {
27eaf44b 1168 t.SetNLast(t.GetNLast() + 1);
1169 t.SetChi2Last(t.GetChi2Last() + maxChi2);
8979685e 1170 }
4551ea7c 1171
59393e34 1172 Double_t xcluster = cl->GetX();
1173 t.PropagateTo(xcluster,radLength,rho);
6c94f330 1174
4551ea7c 1175 if (!AdjustSector(&t)) {
1176 break; //I.B's fix
1177 }
1178 maxChi2 = t.GetPredictedChi2(cl,h01);
6c94f330 1179
c2239ca8 1180 if (maxChi2<1e+10)
1181 if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1182 // ????
1183 }
4551ea7c 1184
6c94f330 1185 }
4551ea7c 1186
8979685e 1187 }
4551ea7c 1188
6c94f330 1189 }
4551ea7c 1190
8979685e 1191 }
8979685e 1192
75bd7f81 1193 return expectedNumberOfClusters;
69b55c55 1194
75bd7f81 1195}
7b580082 1196
75bd7f81 1197//_____________________________________________________________________________
4551ea7c 1198Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
8979685e 1199{
75bd7f81 1200 //
8979685e 1201 // Starting from current radial position of track <t> this function
1202 // extrapolates the track up to outer timebin and in the sensitive
1203 // layers confirms prolongation if a close cluster is found.
1204 // Returns the number of clusters expected to be found in sensitive layers
1205 // Use GEO manager for material Description
75bd7f81 1206 //
7e448bcc 1207 // return number of assigned clusters ?
1208 //
1209
59393e34 1210
4551ea7c 1211 Int_t sector;
1212 Int_t clusters[1000];
6c94f330 1213 Double_t radLength = 0.0;
4551ea7c 1214 Double_t rho = 0.0;
1215 Int_t expectedNumberOfClusters = 0;
1216 Float_t ratio0 = 0.0;
8979685e 1217 AliTRDtracklet tracklet;
33744e5d 1218
77566f2a 1219 // Calibration fill 2D
1220 AliTRDCalibra *calibra = AliTRDCalibra::Instance();
1221 if (!calibra) {
1222 AliInfo("Could not get Calibra instance\n");
1223 }
8ec526a4 1224 if (calibra->GetMITracking()) {
1225 calibra->ResetTrack();
77566f2a 1226 }
1227
33744e5d 1228 for (Int_t i = 0; i < 1000; i++) {
1229 clusters[i] = -1;
1230 }
1231
4551ea7c 1232 for (Int_t iplane = 0; iplane < AliESDtrack::kNPlane; iplane++) {
1233
7e448bcc 1234 int hb = iplane * 10;
1235 fHClSearch->Fill(hb);
1236
4551ea7c 1237 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
1238 Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
1239 Double_t currentx = fTrSec[0]->GetLayer(row0)->GetX();
1240 if (currentx < t.GetX()) {
7e448bcc 1241 fHClSearch->Fill(hb+1);
4551ea7c 1242 continue;
1243 }
1244
6c94f330 1245 //
4551ea7c 1246 // Propagate closer to chamber if neccessary
6c94f330 1247 //
4551ea7c 1248 if (currentx > (fgkMaxStep + t.GetX())) {
1249 if (!PropagateToX(t,currentx-fgkMaxStep,fgkMaxStep)) {
7e448bcc 1250 fHClSearch->Fill(hb+2);
4551ea7c 1251 break;
1252 }
1253 }
1254 if (!AdjustSector(&t)) {
7e448bcc 1255 fHClSearch->Fill(hb+3);
4551ea7c 1256 break;
1257 }
1258 if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) {
7e448bcc 1259 fHClSearch->Fill(hb+4);
4551ea7c 1260 break;
8979685e 1261 }
4551ea7c 1262
8979685e 1263 //
4551ea7c 1264 // Get material budget inside of chamber
59393e34 1265 //
4551ea7c 1266 Double_t xyz0[3];
1267 Double_t xyz1[3];
1268 Double_t param[7];
1269 Double_t x;
1270 Double_t y;
1271 Double_t z;
1272 // Starting global position
1273 t.GetXYZ(xyz0);
1274 // End global position
7b580082 1275 x = fTrSec[0]->GetLayer(rowlast)->GetX();
4551ea7c 1276 if (!t.GetProlongation(x,y,z)) {
7e448bcc 1277 fHClSearch->Fill(hb+5);
4551ea7c 1278 break;
1279 }
1280 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1281 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1282 xyz1[2] = z;
7b580082 1283 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
4551ea7c 1284 rho = param[0];
1285 radLength = param[1]; // Get mean propagation parameters
1286
8979685e 1287 //
7b580082 1288 // Find clusters
8979685e 1289 //
6c94f330 1290 sector = t.GetSector();
4551ea7c 1291 Float_t ncl = FindClusters(sector,row0,rowlast,&t,clusters,tracklet);
7e448bcc 1292 fHNCl->Fill(tracklet.GetN());
1293
4551ea7c 1294 if (tracklet.GetN() < GetTimeBinsPerPlane()/3) {
7e448bcc 1295 fHClSearch->Fill(hb+6);
4551ea7c 1296 continue;
1297 }
1298
8979685e 1299 //
7b580082 1300 // Propagate and update track
8979685e 1301 //
4551ea7c 1302 for (Int_t itime = GetTimeBinsPerPlane()-1; itime >= 0; itime--) {
1303
1304 Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
8979685e 1305 expectedNumberOfClusters++;
27eaf44b 1306 t.SetNExpected(t.GetNExpected() + 1);
4551ea7c 1307 if (t.GetX() > 345.0) {
27eaf44b 1308 t.SetNExpectedLast(t.GetNExpectedLast() + 1);
4551ea7c 1309 }
1310 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
7e448bcc 1311 AliTRDcluster *cl = 0;
4551ea7c 1312 UInt_t index = 0;
1313 Double_t maxChi2 = fgkMaxChi2;
8979685e 1314 x = timeBin.GetX();
4551ea7c 1315
8979685e 1316 if (timeBin) {
4551ea7c 1317
1318 if (clusters[ilayer] > 0) {
8979685e 1319 index = clusters[ilayer];
4551ea7c 1320 cl = (AliTRDcluster *)GetCluster(index);
33744e5d 1321 //Double_t h01 = GetTiltFactor(cl); // I.B's fix
1322 //maxChi2=t.GetPredictedChi2(cl,h01); //
8979685e 1323 }
1324
1325 if (cl) {
4551ea7c 1326
1327 //if (cl->GetNPads() < 5)
59393e34 1328 Double_t dxsample = timeBin.GetdX();
1329 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
4551ea7c 1330 Double_t h01 = GetTiltFactor(cl);
1331 Int_t det = cl->GetDetector();
1332 Int_t plane = fGeom->GetPlane(det);
1333 if (t.GetX() > 345.0) {
27eaf44b 1334 t.SetNLast(t.GetNLast() + 1);
1335 t.SetChi2Last(t.GetChi2Last() + maxChi2);
8979685e 1336 }
59393e34 1337 Double_t xcluster = cl->GetX();
1338 t.PropagateTo(xcluster,radLength,rho);
4551ea7c 1339 maxChi2 = t.GetPredictedChi2(cl,h01);
6c94f330 1340
c2239ca8 1341 if (maxChi2<1e+10)
1342 if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1343 if (!t.Update(cl,maxChi2,index,h01)) {
1344 // ????
1345 }
1346 }
4551ea7c 1347
8ec526a4 1348 if (calibra->GetMITracking()) {
77566f2a 1349 calibra->UpdateHistograms(cl,&t);
1350 }
1351
4551ea7c 1352 // Reset material budget if 2 consecutive gold
1353 if (plane > 0) {
27eaf44b 1354 if ((t.GetTracklets(plane).GetN() + t.GetTracklets(plane-1).GetN()) > 20) {
1355 t.SetBudget(2,0.0);
4551ea7c 1356 }
1357 }
1358
1359 }
1360
8979685e 1361 }
4551ea7c 1362
59393e34 1363 }
4551ea7c 1364
1365 ratio0 = ncl / Float_t(fTimeBinsPerPlane);
27eaf44b 1366 Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);
1367 if ((tracklet.GetChi2() < 18.0) &&
1368 (ratio0 > 0.8) &&
1369 (ratio1 > 0.6) &&
1370 (ratio0+ratio1 > 1.5) &&
1371 (t.GetNCross() == 0) &&
1372 (TMath::Abs(t.GetSnp()) < 0.85) &&
1373 (t.GetNumberOfClusters() > 20)){
7e448bcc 1374 //if (ratio0 > 0.8) {
4551ea7c 1375 t.MakeBackupTrack(); // Make backup of the track until is gold
59393e34 1376 }
7b580082 1377
8979685e 1378 }
46d29e70 1379
75bd7f81 1380 return expectedNumberOfClusters;
75bd7f81 1381}
1e9bb598 1382
75bd7f81 1383//_____________________________________________________________________________
4551ea7c 1384Int_t AliTRDtracker::PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep)
5443e65e 1385{
75bd7f81 1386 //
5443e65e 1387 // Starting from current radial position of track <t> this function
1388 // extrapolates the track up to radial position <xToGo>.
1389 // Returns 1 if track reaches the plane, and 0 otherwise
75bd7f81 1390 //
1391
59393e34 1392 const Double_t kEpsilon = 0.00001;
4551ea7c 1393 //Double_t tanmax = TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
1394 Double_t xpos = t.GetX();
1395 Double_t dir = (xpos<xToGo) ? 1.0 : -1.0;
1396
1397 while (((xToGo-xpos)*dir) > kEpsilon) {
1398
1399 Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos),maxStep);
1400
1401 Double_t xyz0[3];
1402 Double_t xyz1[3];
1403 Double_t param[7];
1404 Double_t x;
1405 Double_t y;
1406 Double_t z;
1407 // Starting global position
1408 t.GetXYZ(xyz0);
1409 x = xpos + step;
1410
1411 if (!t.GetProlongation(x,y,z)) {
1412 return 0; // No prolongation
1413 }
1414
1415 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1416 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1417 xyz1[2] = z;
1418
59393e34 1419 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
4551ea7c 1420 if (!t.PropagateTo(x,param[1],param[0])) {
1421 return 0;
1422 }
59393e34 1423 AdjustSector(&t);
1424 xpos = t.GetX();
4551ea7c 1425
46d29e70 1426 }
75bd7f81 1427
5443e65e 1428 return 1;
5443e65e 1429
59393e34 1430}
5443e65e 1431
46d29e70 1432//_____________________________________________________________________________
c630aafd 1433Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1434{
75bd7f81 1435 //
c630aafd 1436 // Fills clusters into TRD tracking_sectors
1437 // Note that the numbering scheme for the TRD tracking_sectors
1438 // differs from that of TRD sectors
75bd7f81 1439 //
1440
c630aafd 1441 if (ReadClusters(fClusters,cTree)) {
4551ea7c 1442 AliError("Problem with reading the clusters !");
c630aafd 1443 return 1;
1444 }
4551ea7c 1445 Int_t ncl = fClusters->GetEntriesFast();
1446 fNclusters = ncl;
1447 AliInfo(Form("Sorting %d clusters",ncl));
c630aafd 1448
1449 UInt_t index;
4551ea7c 1450 for (Int_t ichamber = 0; ichamber < 5; ichamber++) {
1451 for (Int_t isector = 0; isector < 18; isector++) {
1452 fHoles[ichamber][isector] = kTRUE;
3c625a9b 1453 }
4551ea7c 1454 }
ad670fba 1455
6c94f330 1456 while (ncl--) {
33744e5d 1457
4551ea7c 1458 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(ncl);
1459 Int_t detector = c->GetDetector();
1460 Int_t localTimeBin = c->GetLocalTimeBin();
1461 Int_t sector = fGeom->GetSector(detector);
1462 Int_t plane = fGeom->GetPlane(detector);
439c63c8 1463 Int_t trackingSector = sector;
4551ea7c 1464
7e448bcc 1465 //if (c->GetLabel(0) > 0) {
1466 if (c->GetQ() > 10) {
3c625a9b 1467 Int_t chamber = fGeom->GetChamber(detector);
4551ea7c 1468 fHoles[chamber][trackingSector] = kFALSE;
3c625a9b 1469 }
c630aafd 1470
6c94f330 1471 Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
4551ea7c 1472 if (gtb < 0) {
1473 continue;
1474 }
029cd327 1475 Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
c630aafd 1476
4551ea7c 1477 index = ncl;
33744e5d 1478
1479 // Apply pos correction
7e448bcc 1480 Transform(c);
1481 fHXCl->Fill(c->GetX());
33744e5d 1482
7e448bcc 1483 fTrSec[trackingSector]->GetLayer(layer)->SetX(c->GetX());
1484 fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
4551ea7c 1485 }
75bd7f81 1486
c630aafd 1487 return 0;
75bd7f81 1488
c630aafd 1489}
1490
1491//_____________________________________________________________________________
b7a0917f 1492void AliTRDtracker::UnloadClusters()
5443e65e 1493{
0a29d0f1 1494 //
5443e65e 1495 // Clears the arrays of clusters and tracks. Resets sectors and timebins
0a29d0f1 1496 //
46d29e70 1497
4551ea7c 1498 Int_t i;
1499 Int_t nentr;
46d29e70 1500
5443e65e 1501 nentr = fClusters->GetEntriesFast();
4551ea7c 1502 for (i = 0; i < nentr; i++) {
1503 delete fClusters->RemoveAt(i);
1504 }
b7a0917f 1505 fNclusters = 0;
46d29e70 1506
5443e65e 1507 nentr = fSeeds->GetEntriesFast();
4551ea7c 1508 for (i = 0; i < nentr; i++) {
1509 delete fSeeds->RemoveAt(i);
1510 }
46d29e70 1511
5443e65e 1512 nentr = fTracks->GetEntriesFast();
4551ea7c 1513 for (i = 0; i < nentr; i++) {
1514 delete fTracks->RemoveAt(i);
1515 }
46d29e70 1516
5443e65e 1517 Int_t nsec = AliTRDgeometry::kNsect;
5443e65e 1518 for (i = 0; i < nsec; i++) {
1519 for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1520 fTrSec[i]->GetLayer(pl)->Clear();
1521 }
1522 }
46d29e70 1523
5443e65e 1524}
a819a5f7 1525
75bd7f81 1526//_____________________________________________________________________________
4551ea7c 1527void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD *esd)
7ad19338 1528{
1529 //
1530 // Creates seeds using clusters between position inner plane and outer plane
1531 //
75bd7f81 1532
4551ea7c 1533 const Double_t kMaxTheta = 1.0;
1534 const Double_t kMaxPhi = 2.0;
1535
1536 const Double_t kRoad0y = 6.0; // Road for middle cluster
1537 const Double_t kRoad0z = 8.5; // Road for middle cluster
1538
1539 const Double_t kRoad1y = 2.0; // Road in y for seeded cluster
1540 const Double_t kRoad1z = 20.0; // Road in z for seeded cluster
1541
1542 const Double_t kRoad2y = 3.0; // Road in y for extrapolated cluster
1543 const Double_t kRoad2z = 20.0; // Road in z for extrapolated cluster
c6f438c0 1544 const Int_t kMaxSeed = 3000;
7ad19338 1545
4551ea7c 1546 Int_t maxSec = AliTRDgeometry::kNsect;
1547
1548 // Linear fitters in planes
1549 TLinearFitter fitterTC(2,"hyp2"); // Fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
1550 TLinearFitter fitterT2(4,"hyp4"); // Fitting with tilting pads - kz not fixed
69b55c55 1551 fitterTC.StoreData(kTRUE);
1552 fitterT2.StoreData(kTRUE);
4551ea7c 1553 AliRieman rieman(1000); // Rieman fitter
1554 AliRieman rieman2(1000); // Rieman fitter
1555
1556 // Find the maximal and minimal layer for the planes
7ad19338 1557 Int_t layers[6][2];
4551ea7c 1558 AliTRDpropagationLayer *reflayers[6];
1559 for (Int_t i = 0; i < 6; i++) {
1560 layers[i][0] = 10000;
1561 layers[i][1] = 0;
1562 }
1563 for (Int_t ns = 0; ns < maxSec; ns++) {
1564 for (Int_t ilayer = 0; ilayer < fTrSec[ns]->GetNumberOfLayers(); ilayer++) {
1565 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(ilayer));
1566 if (layer == 0) {
1567 continue;
1568 }
7ad19338 1569 Int_t det = layer[0]->GetDetector();
1570 Int_t plane = fGeom->GetPlane(det);
4551ea7c 1571 if (ilayer < layers[plane][0]) {
1572 layers[plane][0] = ilayer;
1573 }
1574 if (ilayer > layers[plane][1]) {
1575 layers[plane][1] = ilayer;
1576 }
7ad19338 1577 }
1578 }
4551ea7c 1579
3551db50 1580 AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
6c94f330 1581 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
4551ea7c 1582 Double_t hL[6]; // Tilting angle
1583 Double_t xcl[6]; // X - position of reference cluster
1584 Double_t ycl[6]; // Y - position of reference cluster
1585 Double_t zcl[6]; // Z - position of reference cluster
1586
1587 AliTRDcluster *cl[6] = { 0, 0, 0, 0, 0, 0 }; // Seeding clusters
1588 Float_t padlength[6] = { 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 }; // Current pad-length
1589
1590 Double_t chi2R = 0.0;
1591 Double_t chi2Z = 0.0;
1592 Double_t chi2RF = 0.0;
1593 Double_t chi2ZF = 0.0;
1594
1595 Int_t nclusters; // Total number of clusters
1596 for (Int_t i = 0; i < 6; i++) {
1597 hL[i] = h01;
1598 if (i%2==1) {
1599 hL[i]*=-1.0;
1600 }
1601 }
1602
1603 // Registered seed
c6f438c0 1604 AliTRDseed *pseed = new AliTRDseed[kMaxSeed*6];
1605 AliTRDseed *seed[kMaxSeed];
4551ea7c 1606 for (Int_t iseed = 0; iseed < kMaxSeed; iseed++) {
1607 seed[iseed]= &pseed[iseed*6];
1608 }
69b55c55 1609 AliTRDseed *cseed = seed[0];
4551ea7c 1610
1611 Double_t seedquality[kMaxSeed];
1612 Double_t seedquality2[kMaxSeed];
1613 Double_t seedparams[kMaxSeed][7];
1614 Int_t seedlayer[kMaxSeed];
1615 Int_t registered = 0;
1616 Int_t sort[kMaxSeed];
1617
1618 //
1619 // Seeding part
1620 //
1621 for (Int_t ns = 0; ns < maxSec; ns++) { // Loop over sectors
1622 //for (Int_t ns = 0; ns < 5; ns++) { // Loop over sectors
1623
1624 registered = 0; // Reset registerd seed counter
1625 cseed = seed[registered];
1626 Float_t iter = 0.0;
1627
1628 for (Int_t sLayer = 2; sLayer >= 0; sLayer--) {
1629 //for (Int_t dseed = 5; dseed < 15; dseed += 3) {
1630
1631 iter += 1.0;
1632 Int_t dseed = 5 + Int_t(iter) * 3;
1633
69b55c55 1634 // Initialize seeding layers
4551ea7c 1635 for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
69b55c55 1636 reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
1637 xcl[ilayer] = reflayers[ilayer]->GetX();
4551ea7c 1638 }
1639
1640 Double_t xref = (xcl[sLayer+1] + xcl[sLayer+2]) * 0.5;
1641 AliTRDpropagationLayer &layer0 = *reflayers[sLayer+0];
1642 AliTRDpropagationLayer &layer1 = *reflayers[sLayer+1];
1643 AliTRDpropagationLayer &layer2 = *reflayers[sLayer+2];
1644 AliTRDpropagationLayer &layer3 = *reflayers[sLayer+3];
1645
1646 Int_t maxn3 = layer3;
1647 for (Int_t icl3 = 0; icl3 < maxn3; icl3++) {
33744e5d 1648
69b55c55 1649 AliTRDcluster *cl3 = layer3[icl3];
4551ea7c 1650 if (!cl3) {
1651 continue;
1652 }
1653 padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2() * 12.0);
69b55c55 1654 ycl[sLayer+3] = cl3->GetY();
1655 zcl[sLayer+3] = cl3->GetZ();
4551ea7c 1656 Float_t yymin0 = ycl[sLayer+3] - 1.0 - kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1657 Float_t yymax0 = ycl[sLayer+3] + 1.0 + kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1658 Int_t maxn0 = layer0;
1659
1660 for (Int_t icl0 = layer0.Find(yymin0); icl0 < maxn0; icl0++) {
1661
69b55c55 1662 AliTRDcluster *cl0 = layer0[icl0];
4551ea7c 1663 if (!cl0) {
1664 continue;
1665 }
1666 if (cl3->IsUsed() && cl0->IsUsed()) {
1667 continue;
1668 }
69b55c55 1669 ycl[sLayer+0] = cl0->GetY();
1670 zcl[sLayer+0] = cl0->GetZ();
4551ea7c 1671 if (ycl[sLayer+0] > yymax0) {
1672 break;
1673 }
1674 Double_t tanphi = (ycl[sLayer+3]-ycl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1675 if (TMath::Abs(tanphi) > kMaxPhi) {
1676 continue;
1677 }
1678 Double_t tantheta = (zcl[sLayer+3]-zcl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1679 if (TMath::Abs(tantheta) > kMaxTheta) {
1680 continue;
1681 }
1682 padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2() * 12.0);
1683
1684 // Expected position in 1 layer
1685 Double_t y1exp = ycl[sLayer+0] + (tanphi) * (xcl[sLayer+1]-xcl[sLayer+0]);
1686 Double_t z1exp = zcl[sLayer+0] + (tantheta) * (xcl[sLayer+1]-xcl[sLayer+0]);
1687 Float_t yymin1 = y1exp - kRoad0y - tanphi;
1688 Float_t yymax1 = y1exp + kRoad0y + tanphi;
1689 Int_t maxn1 = layer1;
1690
1691 for (Int_t icl1 = layer1.Find(yymin1); icl1 < maxn1; icl1++) {
1692
69b55c55 1693 AliTRDcluster *cl1 = layer1[icl1];
4551ea7c 1694 if (!cl1) {
1695 continue;
1696 }
69b55c55 1697 Int_t nusedCl = 0;
1698 if (cl3->IsUsed()) nusedCl++;
1699 if (cl0->IsUsed()) nusedCl++;
1700 if (cl1->IsUsed()) nusedCl++;
4551ea7c 1701 if (nusedCl > 1) {
1702 continue;
1703 }
69b55c55 1704 ycl[sLayer+1] = cl1->GetY();
1705 zcl[sLayer+1] = cl1->GetZ();
4551ea7c 1706 if (ycl[sLayer+1] > yymax1) {
1707 break;
1708 }
1709 if (TMath::Abs(ycl[sLayer+1]-y1exp) > kRoad0y+tanphi) {
1710 continue;
1711 }
1712 if (TMath::Abs(zcl[sLayer+1]-z1exp) > kRoad0z) {
1713 continue;
1714 }
1715 padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2() * 12.0);
1716
1717 Double_t y2exp = ycl[sLayer+0]+(tanphi) * (xcl[sLayer+2]-xcl[sLayer+0]) + (ycl[sLayer+1]-y1exp);
1718 Double_t z2exp = zcl[sLayer+0]+(tantheta) * (xcl[sLayer+2]-xcl[sLayer+0]);
1719 Int_t index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y,kRoad1z);
1720 if (index2 <= 0) {
1721 continue;
1722 }
1723 AliTRDcluster *cl2 = (AliTRDcluster *) GetCluster(index2);
1724 padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2() * 12.0);
1725 ycl[sLayer+2] = cl2->GetY();
1726 zcl[sLayer+2] = cl2->GetZ();
1727 if (TMath::Abs(cl2->GetZ()-z2exp) > kRoad0z) {
1728 continue;
1729 }
1730
69b55c55 1731 rieman.Reset();
1732 rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
1733 rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
1734 rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);
1735 rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
1736 rieman.Update();
4551ea7c 1737
1738 // Reset fitter
1739 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
69b55c55 1740 cseed[iLayer].Reset();
1741 }
4551ea7c 1742 chi2Z = 0.0;
1743 chi2R = 0.0;
1744
1745 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
ca1e1e5b 1746 cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1747 chi2Z += (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer])
1748 * (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer]);
1749 cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
1750 cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1751 chi2R += (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer])
1752 * (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer]);
1753 cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
69b55c55 1754 }
4551ea7c 1755 if (TMath::Sqrt(chi2R) > 1.0/iter) {
1756 continue;
69b55c55 1757 }
4551ea7c 1758 if (TMath::Sqrt(chi2Z) > 7.0/iter) {
1759 continue;
1760 }
1761
1762 Float_t minmax[2] = { -100.0, 100.0 };
1763 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1764 Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer] * 0.5
ca1e1e5b 1765 + 1.0 - cseed[sLayer+iLayer].GetZref(0);
4551ea7c 1766 if (max < minmax[1]) {
1767 minmax[1] = max;
1768 }
1769 Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer] * 0.5
ca1e1e5b 1770 - 1.0 - cseed[sLayer+iLayer].GetZref(0);
4551ea7c 1771 if (min > minmax[0]) {
1772 minmax[0] = min;
1773 }
1774 }
1775
69b55c55 1776 Bool_t isFake = kFALSE;
4551ea7c 1777 if (cl0->GetLabel(0) != cl3->GetLabel(0)) {
1778 isFake = kTRUE;
1779 }
1780 if (cl1->GetLabel(0) != cl3->GetLabel(0)) {
1781 isFake = kTRUE;
1782 }
1783 if (cl2->GetLabel(0) != cl3->GetLabel(0)) {
1784 isFake = kTRUE;
1785 }
1786
1787 if (AliTRDReconstructor::StreamLevel() > 0) {
1788 if ((!isFake) || ((icl3%10) == 0)) { // Debugging print
1789 TTreeSRedirector &cstream = *fDebugStreamer;
1790 cstream << "Seeds0"
1791 << "isFake=" << isFake
1792 << "Cl0.=" << cl0
1793 << "Cl1.=" << cl1
1794 << "Cl2.=" << cl2
1795 << "Cl3.=" << cl3
1796 << "Xref=" << xref
1797 << "X0=" << xcl[sLayer+0]
1798 << "X1=" << xcl[sLayer+1]
1799 << "X2=" << xcl[sLayer+2]
1800 << "X3=" << xcl[sLayer+3]
1801 << "Y2exp=" << y2exp
1802 << "Z2exp=" << z2exp
1803 << "Chi2R=" << chi2R
1804 << "Chi2Z=" << chi2Z
1805 << "Seed0.=" << &cseed[sLayer+0]
1806 << "Seed1.=" << &cseed[sLayer+1]
1807 << "Seed2.=" << &cseed[sLayer+2]
1808 << "Seed3.=" << &cseed[sLayer+3]
1809 << "Zmin=" << minmax[0]
1810 << "Zmax=" << minmax[1]
1811 << "\n";
d337ef8d 1812 }
69b55c55 1813 }
33744e5d 1814
4551ea7c 1815 ////////////////////////////////////////////////////////////////////////////////////
33744e5d 1816 //
4551ea7c 1817 // Fit seeding part
33744e5d 1818 //
4551ea7c 1819 ////////////////////////////////////////////////////////////////////////////////////
33744e5d 1820
69b55c55 1821 cl[sLayer+0] = cl0;
1822 cl[sLayer+1] = cl1;
1823 cl[sLayer+2] = cl2;
1824 cl[sLayer+3] = cl3;
4551ea7c 1825 Bool_t isOK = kTRUE;
1826
1827 for (Int_t jLayer = 0; jLayer < 4; jLayer++) {
1828
ca1e1e5b 1829 cseed[sLayer+jLayer].SetTilt(hL[sLayer+jLayer]);
1830 cseed[sLayer+jLayer].SetPadLength(padlength[sLayer+jLayer]);
1831 cseed[sLayer+jLayer].SetX0(xcl[sLayer+jLayer]);
4551ea7c 1832
1833 for (Int_t iter = 0; iter < 2; iter++) {
1834
69b55c55 1835 //
4551ea7c 1836 // In iteration 0 we try only one pad-row
1837 // If quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
69b55c55 1838 //
1839 AliTRDseed tseed = cseed[sLayer+jLayer];
4551ea7c 1840 Float_t roadz = padlength[sLayer+jLayer] * 0.5;
1841 if (iter > 0) {
1842 roadz = padlength[sLayer+jLayer];
1843 }
1844
1845 Float_t quality = 10000.0;
1846
1847 for (Int_t iTime = 2; iTime < 20; iTime++) {
1848
1849 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
1850 Double_t dxlayer = layer.GetX() - xcl[sLayer+jLayer];
1851 Double_t zexp = cl[sLayer+jLayer]->GetZ();
1852
1853 if (iter > 0) {
1854 // Try 2 pad-rows in second iteration
ca1e1e5b 1855 zexp = tseed.GetZref(0) + tseed.GetZref(1) * dxlayer;
4551ea7c 1856 if (zexp > cl[sLayer+jLayer]->GetZ()) {
1857 zexp = cl[sLayer+jLayer]->GetZ() + padlength[sLayer+jLayer]*0.5;
1858 }
1859 if (zexp < cl[sLayer+jLayer]->GetZ()) {
1860 zexp = cl[sLayer+jLayer]->GetZ() - padlength[sLayer+jLayer]*0.5;
1861 }
1862 }
1863
ca1e1e5b 1864 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
4551ea7c 1865 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
1866 if (index <= 0) {
1867 continue;
69b55c55 1868 }
4551ea7c 1869 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
1870
ca1e1e5b 1871 tseed.SetIndexes(iTime,index);
1872 tseed.SetClusters(iTime,cl); // Register cluster
1873 tseed.SetX(iTime,dxlayer); // Register cluster
1874 tseed.SetY(iTime,cl->GetY()); // Register cluster
1875 tseed.SetZ(iTime,cl->GetZ()); // Register cluster
4551ea7c 1876
1877 }
1878
69b55c55 1879 tseed.Update();
4551ea7c 1880
1881 // Count the number of clusters and distortions into quality
ca1e1e5b 1882 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1883 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
1884 + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
1885 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 1886 if ((iter == 0) && tseed.IsOK()) {
69b55c55 1887 cseed[sLayer+jLayer] = tseed;
4551ea7c 1888 quality = tquality;
1889 if (tquality < 5) {
1890 break;
1891 }
69b55c55 1892 }
4551ea7c 1893 if (tseed.IsOK() && (tquality < quality)) {
1894 cseed[sLayer+jLayer] = tseed;
1895 }
1896
1897 } // Loop: iter
1898
1899 if (!cseed[sLayer+jLayer].IsOK()) {
69b55c55 1900 isOK = kFALSE;
1901 break;
4551ea7c 1902 }
1903
69b55c55 1904 cseed[sLayer+jLayer].CookLabels();
1905 cseed[sLayer+jLayer].UpdateUsed();
ca1e1e5b 1906 nusedCl += cseed[sLayer+jLayer].GetNUsed();
4551ea7c 1907 if (nusedCl > 25) {
69b55c55 1908 isOK = kFALSE;
1909 break;
4551ea7c 1910 }
1911
1912 } // Loop: jLayer
1913
1914 if (!isOK) {
1915 continue;
69b55c55 1916 }
4551ea7c 1917 nclusters = 0;
1918 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1919 if (cseed[sLayer+iLayer].IsOK()) {
ca1e1e5b 1920 nclusters += cseed[sLayer+iLayer].GetN2();
69b55c55 1921 }
1922 }
4551ea7c 1923
1924 // Iteration 0
69b55c55 1925 rieman.Reset();
4551ea7c 1926 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1927 rieman.AddPoint(xcl[sLayer+iLayer]
ca1e1e5b 1928 ,cseed[sLayer+iLayer].GetYfitR(0)
1929 ,cseed[sLayer+iLayer].GetZProb()
4551ea7c 1930 ,1
1931 ,10);
69b55c55 1932 }
1933 rieman.Update();
4551ea7c 1934
1935 chi2R = 0.0;
1936 chi2Z = 0.0;
1937
1938 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
ca1e1e5b 1939 cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1940 chi2R += (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0))
1941 * (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0));
1942 cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
1943 cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1944 chi2Z += (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz())
1945 * (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz());
1946 cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
69b55c55 1947 }
1948 Double_t curv = rieman.GetC();
4551ea7c 1949
69b55c55 1950 //
4551ea7c 1951 // Likelihoods
6c94f330 1952 //
ca1e1e5b 1953 Double_t sumda = TMath::Abs(cseed[sLayer+0].GetYfitR(1) - cseed[sLayer+0].GetYref(1))
1954 + TMath::Abs(cseed[sLayer+1].GetYfitR(1) - cseed[sLayer+1].GetYref(1))
1955 + TMath::Abs(cseed[sLayer+2].GetYfitR(1) - cseed[sLayer+2].GetYref(1))
1956 + TMath::Abs(cseed[sLayer+3].GetYfitR(1) - cseed[sLayer+3].GetYref(1));
4551ea7c 1957 Double_t likea = TMath::Exp(-sumda*10.6);
1958 Double_t likechi2 = 0.0000000001;
1959 if (chi2R < 0.5) {
1960 likechi2 += TMath::Exp(-TMath::Sqrt(chi2R) * 7.73);
1961 }
1962 Double_t likechi2z = TMath::Exp(-chi2Z * 0.088) / TMath::Exp(-chi2Z * 0.019);
1963 Double_t likeN = TMath::Exp(-(72 - nclusters) * 0.19);
1964 Double_t like = likea * likechi2 * likechi2z * likeN;
ca1e1e5b 1965 Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetYref(1) - 130.0*curv) * 1.9);
1966 Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetZref(1)
1967 - cseed[sLayer+0].GetZref(0) / xcl[sLayer+0]) * 5.9);
6c94f330 1968 Double_t likePrim = TMath::Max(likePrimY*likePrimZ,0.0005);
69b55c55 1969
4551ea7c 1970 seedquality[registered] = like;
1971 seedlayer[registered] = sLayer;
1972 if (TMath::Log(0.000000000000001 + like) < -15) {
1973 continue;
1974 }
69b55c55 1975 AliTRDseed seedb[6];
4551ea7c 1976 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
69b55c55 1977 seedb[iLayer] = cseed[iLayer];
1978 }
4551ea7c 1979
1980 ////////////////////////////////////////////////////////////////////////////////////
69b55c55 1981 //
4551ea7c 1982 // Full track fit part
33744e5d 1983 //
4551ea7c 1984 ////////////////////////////////////////////////////////////////////////////////////
1985
1986 Int_t nlayers = 0;
1987 Int_t nusedf = 0;
1988 Int_t findable = 0;
1989
69b55c55 1990 //
4551ea7c 1991 // Add new layers - avoid long extrapolation
69b55c55 1992 //
4551ea7c 1993 Int_t tLayer[2] = { 0, 0 };
1994 if (sLayer == 2) {
1995 tLayer[0] = 1;
1996 tLayer[1] = 0;
1997 }
1998 if (sLayer == 1) {
1999 tLayer[0] = 5;
2000 tLayer[1] = 0;
2001 }
2002 if (sLayer == 0) {
2003 tLayer[0] = 4;
2004 tLayer[1] = 5;
2005 }
2006
2007 for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
2008 Int_t jLayer = tLayer[iLayer]; // Set tracking layer
69b55c55 2009 cseed[jLayer].Reset();
ca1e1e5b 2010 cseed[jLayer].SetTilt(hL[jLayer]);
2011 cseed[jLayer].SetPadLength(padlength[jLayer]);
2012 cseed[jLayer].SetX0(xcl[jLayer]);
4551ea7c 2013 // Get pad length and rough cluster
ca1e1e5b 2014 Int_t indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].GetYref(0)
2015 ,cseed[jLayer].GetZref(0)
4551ea7c 2016 ,kRoad2y
2017 ,kRoad2z);
2018 if (indexdummy <= 0) {
2019 continue;
2020 }
2021 AliTRDcluster *cldummy = (AliTRDcluster *) GetCluster(indexdummy);
2022 padlength[jLayer] = TMath::Sqrt(cldummy->GetSigmaZ2() * 12.0);
69b55c55 2023 }
4551ea7c 2024 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
2025
2026 for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
2027
2028 Int_t jLayer = tLayer[iLayer]; // set tracking layer
2029 if ((jLayer == 0) && !(cseed[1].IsOK())) {
2030 continue; // break not allowed
2031 }
2032 if ((jLayer == 5) && !(cseed[4].IsOK())) {
2033 continue; // break not allowed
2034 }
ca1e1e5b 2035 Float_t zexp = cseed[jLayer].GetZref(0);
4551ea7c 2036 Double_t zroad = padlength[jLayer] * 0.5 + 1.0;
2037
2038 for (Int_t iter = 0; iter < 2; iter++) {
2039
69b55c55 2040 AliTRDseed tseed = cseed[jLayer];
4551ea7c 2041 Float_t quality = 10000.0;
2042
2043 for (Int_t iTime = 2; iTime < 20; iTime++) {
2044 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
2045 Double_t dxlayer = layer.GetX()-xcl[jLayer];
ca1e1e5b 2046 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
4551ea7c 2047 Float_t yroad = kRoad1y;
2048 Int_t index = layer.FindNearestCluster(yexp,zexp,yroad,zroad);
2049 if (index <= 0) {
2050 continue;
2051 }
2052 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
ca1e1e5b 2053 tseed.SetIndexes(iTime,index);
2054 tseed.SetClusters(iTime,cl); // Register cluster
2055 tseed.SetX(iTime,dxlayer); // Register cluster
2056 tseed.SetY(iTime,cl->GetY()); // Register cluster
2057 tseed.SetZ(iTime,cl->GetZ()); // Register cluster
4551ea7c 2058 }
2059
69b55c55 2060 tseed.Update();
4551ea7c 2061 if (tseed.IsOK()) {
ca1e1e5b 2062 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
2063 Float_t tquality = (18.0 - tseed.GetN2())/2.0 + TMath::Abs(dangle) / 0.1
2064 + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
2065 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 2066 if (tquality < quality) {
2067 cseed[jLayer] = tseed;
2068 quality = tquality;
69b55c55 2069 }
2070 }
4551ea7c 2071
2072 zroad *= 2.0;
2073
2074 } // Loop: iter
2075
2076 if ( cseed[jLayer].IsOK()) {
69b55c55 2077 cseed[jLayer].CookLabels();
2078 cseed[jLayer].UpdateUsed();
ca1e1e5b 2079 nusedf += cseed[jLayer].GetNUsed();
4551ea7c 2080 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
69b55c55 2081 }
4551ea7c 2082
2083 } // Loop: iLayer
2084
2085 // Make copy
69b55c55 2086 AliTRDseed bseed[6];
4551ea7c 2087 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
69b55c55 2088 bseed[jLayer] = cseed[jLayer];
4551ea7c 2089 }
2090 Float_t lastquality = 10000.0;
2091 Float_t lastchi2 = 10000.0;
2092 Float_t chi2 = 1000.0;
2093
2094 for (Int_t iter = 0; iter < 4; iter++) {
2095
2096 // Sort tracklets according "quality", try to "improve" 4 worst
2097 Float_t sumquality = 0.0;
69b55c55 2098 Float_t squality[6];
2099 Int_t sortindexes[6];
4551ea7c 2100
2101 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2102 if (bseed[jLayer].IsOK()) {
69b55c55 2103 AliTRDseed &tseed = bseed[jLayer];
ca1e1e5b 2104 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
2105 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
2106 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
2107 + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
2108 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 2109 squality[jLayer] = tquality;
2110 }
2111 else {
2112 squality[jLayer] = -1.0;
ad670fba 2113 }
6c94f330 2114 sumquality +=squality[jLayer];
69b55c55 2115 }
2116
4551ea7c 2117 if ((sumquality >= lastquality) ||
2118 (chi2 > lastchi2)) {
2119 break;
2120 }
69b55c55 2121 lastquality = sumquality;
2122 lastchi2 = chi2;
4551ea7c 2123 if (iter > 0) {
2124 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
69b55c55 2125 cseed[jLayer] = bseed[jLayer];
2126 }
2127 }
2128 TMath::Sort(6,squality,sortindexes,kFALSE);
4551ea7c 2129
2130 for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
2131
6c94f330 2132 Int_t bLayer = sortindexes[jLayer];
2133 AliTRDseed tseed = bseed[bLayer];
4551ea7c 2134
2135 for (Int_t iTime = 2; iTime < 20; iTime++) {
2136
2137 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
2138 Double_t dxlayer = layer.GetX() - xcl[bLayer];
ca1e1e5b 2139 Double_t zexp = tseed.GetZref(0);
2140 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
4551ea7c 2141 Float_t roadz = padlength[bLayer] + 1;
ca1e1e5b 2142 if (TMath::Abs(tseed.GetZProb() - zexp) > 0.5*padlength[bLayer]) {
4551ea7c 2143 roadz = padlength[bLayer] * 0.5;
2144 }
ca1e1e5b 2145 if (tseed.GetZfit(1)*tseed.GetZref(1) < 0.0) {
4551ea7c 2146 roadz = padlength[bLayer] * 0.5;
2147 }
ca1e1e5b 2148 if (TMath::Abs(tseed.GetZProb() - zexp) < 0.1*padlength[bLayer]) {
2149 zexp = tseed.GetZProb();
4551ea7c 2150 roadz = padlength[bLayer] * 0.5;
2151 }
2152
ca1e1e5b 2153 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer - zcor;
4551ea7c 2154 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
2155 if (index <= 0) {
2156 continue;
69b55c55 2157 }
4551ea7c 2158 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
2159
ca1e1e5b 2160 tseed.SetIndexes(iTime,index);
2161 tseed.SetClusters(iTime,cl); // Register cluster
2162 tseed.SetX(iTime,dxlayer); // Register cluster
2163 tseed.SetY(iTime,cl->GetY()); // Register cluster
2164 tseed.SetZ(iTime,cl->GetZ()); // Register cluster
4551ea7c 2165
2166 }
2167
69b55c55 2168 tseed.Update();
c6f438c0 2169 if (tseed.IsOK()) {
ca1e1e5b 2170 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
2171 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
2172 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0
2173 + TMath::Abs(dangle) / 0.1
2174 + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
2175 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 2176 if (tquality<squality[bLayer]) {
69b55c55 2177 bseed[bLayer] = tseed;
4551ea7c 2178 }
69b55c55 2179 }
4551ea7c 2180
2181 } // Loop: jLayer
2182
2183 chi2 = AliTRDseed::FitRiemanTilt(bseed,kTRUE);
2184
2185 } // Loop: iter
2186
2187 nclusters = 0;
2188 nlayers = 0;
2189 findable = 0;
2190 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ca1e1e5b 2191 if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) {
69b55c55 2192 findable++;
4551ea7c 2193 }
2194 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2195 nclusters += cseed[iLayer].GetN2();
69b55c55 2196 nlayers++;
2197 }
2198 }
4551ea7c 2199 if (nlayers < 3) {
2200 continue;
2201 }
69b55c55 2202 rieman.Reset();
4551ea7c 2203 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2204 if (cseed[iLayer].IsOK()) {
2205 rieman.AddPoint(xcl[iLayer]
ca1e1e5b 2206 ,cseed[iLayer].GetYfitR(0)
2207 ,cseed[iLayer].GetZProb()
4551ea7c 2208 ,1
2209 ,10);
2210 }
69b55c55 2211 }
2212 rieman.Update();
4551ea7c 2213
2214 chi2RF = 0.0;
2215 chi2ZF = 0.0;
2216 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2217 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2218 cseed[iLayer].SetYref(0,rieman.GetYat(xcl[iLayer]));
2219 chi2RF += (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0))
2220 * (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0));
2221 cseed[iLayer].SetYref(1,rieman.GetDYat(xcl[iLayer]));
2222 cseed[iLayer].SetZref(0,rieman.GetZat(xcl[iLayer]));
2223 chi2ZF += (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz())
2224 * (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz());
2225 cseed[iLayer].SetZref(1,rieman.GetDZat(xcl[iLayer]));
69b55c55 2226 }
2227 }
4551ea7c 2228 chi2RF /= TMath::Max((nlayers - 3.0),1.0);
2229 chi2ZF /= TMath::Max((nlayers - 3.0),1.0);
2230 curv = rieman.GetC();
2231
2232 Double_t xref2 = (xcl[2] + xcl[3]) * 0.5; // Middle of the chamber
2233 Double_t dzmf = rieman.GetDZat(xref2);
2234 Double_t zmf = rieman.GetZat(xref2);
69b55c55 2235
69b55c55 2236 //
4551ea7c 2237 // Fit hyperplane
69b55c55 2238 //
4551ea7c 2239 Int_t npointsT = 0;
69b55c55 2240 fitterTC.ClearPoints();
2241 fitterT2.ClearPoints();
2242 rieman2.Reset();
4551ea7c 2243
2244 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2245
2246 if (!cseed[iLayer].IsOK()) {
2247 continue;
2248 }
2249
2250 for (Int_t itime = 0; itime < 25; itime++) {
2251
ca1e1e5b 2252 if (!cseed[iLayer].IsUsable(itime)) {
4551ea7c 2253 continue;
2254 }
2255 // X relative to the middle chamber
ca1e1e5b 2256 Double_t x = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;
2257 Double_t y = cseed[iLayer].GetY(itime);
2258 Double_t z = cseed[iLayer].GetZ(itime);
69b55c55 2259 // ExB correction to the correction
4551ea7c 2260 // Tilted rieman
69b55c55 2261 Double_t uvt[6];
4551ea7c 2262 // Global x
ca1e1e5b 2263 Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();
4551ea7c 2264 Double_t t = 1.0 / (x2*x2 + y*y);
2265 uvt[1] = t; // t
2266 uvt[0] = 2.0 * x2 * uvt[1]; // u
2267 uvt[2] = 2.0 * hL[iLayer] * uvt[1];
2268 uvt[3] = 2.0 * hL[iLayer] * x * uvt[1];
2269 uvt[4] = 2.0 * (y + hL[iLayer]*z) * uvt[1];
2270 Double_t error = 2.0 * 0.2 * uvt[1];
69b55c55 2271 fitterT2.AddPoint(uvt,uvt[4],error);
4551ea7c 2272
69b55c55 2273 //
4551ea7c 2274 // Constrained rieman
69b55c55 2275 //
ca1e1e5b 2276 z = cseed[iLayer].GetZ(itime);
4551ea7c 2277 uvt[0] = 2.0 * x2 * t; // u
2278 uvt[1] = 2.0 * hL[iLayer] * x2 * uvt[1];
2279 uvt[2] = 2.0 * (y + hL[iLayer] * (z - GetZ())) * t;
69b55c55 2280 fitterTC.AddPoint(uvt,uvt[2],error);
69b55c55 2281 rieman2.AddPoint(x2,y,z,1,10);
2282 npointsT++;
4551ea7c 2283
69b55c55 2284 }
4551ea7c 2285
2286 } // Loop: iLayer
2287
69b55c55 2288 rieman2.Update();
2289 fitterTC.Eval();
2290 fitterT2.Eval();
2291 Double_t rpolz0 = fitterT2.GetParameter(3);
2292 Double_t rpolz1 = fitterT2.GetParameter(4);
4551ea7c 2293
69b55c55 2294 //
4551ea7c 2295 // Linear fitter - not possible to make boundaries
2296 // Do not accept non possible z and dzdx combinations
69b55c55 2297 //
4551ea7c 2298 Bool_t acceptablez = kTRUE;
2299 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2300 if (cseed[iLayer].IsOK()) {
2301 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
ca1e1e5b 2302 if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > padlength[iLayer] * 0.5 + 1.0) {
69b55c55 2303 acceptablez = kFALSE;
4551ea7c 2304 }
69b55c55 2305 }
2306 }
4551ea7c 2307 if (!acceptablez) {
69b55c55 2308 fitterT2.FixParameter(3,zmf);
2309 fitterT2.FixParameter(4,dzmf);
2310 fitterT2.Eval();
2311 fitterT2.ReleaseParameter(3);
2312 fitterT2.ReleaseParameter(4);
2313 rpolz0 = fitterT2.GetParameter(3);
2314 rpolz1 = fitterT2.GetParameter(4);
2315 }
4551ea7c 2316
2317 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
2318 Double_t chi2TC = fitterTC.GetChisquare() / Float_t(npointsT);
69b55c55 2319 Double_t polz1c = fitterTC.GetParameter(2);
4551ea7c 2320 Double_t polz0c = polz1c * xref2;
6c94f330 2321 Double_t aC = fitterTC.GetParameter(0);
2322 Double_t bC = fitterTC.GetParameter(1);
4551ea7c 2323 Double_t cC = aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
6c94f330 2324 Double_t aR = fitterT2.GetParameter(0);
2325 Double_t bR = fitterT2.GetParameter(1);
2326 Double_t dR = fitterT2.GetParameter(2);
4551ea7c 2327 Double_t cR = 1.0 + bR*bR - dR*aR;
2328 Double_t dca = 0.0;
2329 if (cR > 0.0) {
2330 dca = -dR / (TMath::Sqrt(1.0 + bR*bR - dR*aR) + TMath::Sqrt(1.0 + bR*bR));
2331 cR = aR / TMath::Sqrt(cR);
69b55c55 2332 }
4551ea7c 2333
2334 Double_t chi2ZT2 = 0.0;
2335 Double_t chi2ZTC = 0.0;
2336 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2337 if (cseed[iLayer].IsOK()) {
2338 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
2339 Double_t zTC = polz0c + polz1c * (xcl[iLayer] - xref2);
ca1e1e5b 2340 chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
2341 chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
69b55c55 2342 }
2343 }
4551ea7c 2344 chi2ZT2 /= TMath::Max((nlayers - 3.0),1.0);
2345 chi2ZTC /= TMath::Max((nlayers - 3.0),1.0);
2346
2347 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
2348 Float_t sumdaf = 0.0;
2349 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2350 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2351 sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))
2352 / cseed[iLayer].GetSigmaY2());
4551ea7c 2353 }
2354 }
2355 sumdaf /= Float_t (nlayers - 2.0);
2356
69b55c55 2357 //
4551ea7c 2358 // Likelihoods for full track
69b55c55 2359 //
4551ea7c 2360 Double_t likezf = TMath::Exp(-chi2ZF * 0.14);
2361 Double_t likechi2C = TMath::Exp(-chi2TC * 0.677);
2362 Double_t likechi2TR = TMath::Exp(-chi2TR * 0.78);
2363 Double_t likeaf = TMath::Exp(-sumdaf * 3.23);
2364 seedquality2[registered] = likezf * likechi2TR * likeaf;
33744e5d 2365
4551ea7c 2366 // Still needed ????
69b55c55 2367// Bool_t isGold = kFALSE;
2368//
6c94f330 2369// if (nlayers == 6 && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE; // gold
2370// if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE; // gold
69b55c55 2371// if (isGold &&nusedf<10){
2372// for (Int_t jLayer=0;jLayer<6;jLayer++){
6c94f330 2373// if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
69b55c55 2374// seed[index][jLayer].UseClusters(); //sign gold
2375// }
2376// }
4551ea7c 2377
2378 Int_t index0 = 0;
2379 if (!cseed[0].IsOK()) {
69b55c55 2380 index0 = 1;
4551ea7c 2381 if (!cseed[1].IsOK()) {
2382 index0 = 2;
2383 }
69b55c55 2384 }
ca1e1e5b 2385 seedparams[registered][0] = cseed[index0].GetX0();
2386 seedparams[registered][1] = cseed[index0].GetYref(0);
2387 seedparams[registered][2] = cseed[index0].GetZref(0);
c6f438c0 2388 seedparams[registered][5] = cR;
ca1e1e5b 2389 seedparams[registered][3] = cseed[index0].GetX0() * cR - TMath::Sin(TMath::ATan(cseed[0].GetYref(1)));
2390 seedparams[registered][4] = cseed[index0].GetZref(1)
2391 / TMath::Sqrt(1.0 + cseed[index0].GetYref(1) * cseed[index0].GetYref(1));
69b55c55 2392 seedparams[registered][6] = ns;
4551ea7c 2393
2394 Int_t labels[12];
2395 Int_t outlab[24];
2396 Int_t nlab = 0;
2397 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2398 if (!cseed[iLayer].IsOK()) {
2399 continue;
2400 }
ca1e1e5b 2401 if (cseed[iLayer].GetLabels(0) >= 0) {
2402 labels[nlab] = cseed[iLayer].GetLabels(0);
69b55c55 2403 nlab++;
2404 }
ca1e1e5b 2405 if (cseed[iLayer].GetLabels(1) >= 0) {
2406 labels[nlab] = cseed[iLayer].GetLabels(1);
69b55c55 2407 nlab++;
4551ea7c 2408 }
69b55c55 2409 }
2410 Freq(nlab,labels,outlab,kFALSE);
4551ea7c 2411 Int_t label = outlab[0];
2412 Int_t frequency = outlab[1];
2413 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ca1e1e5b 2414 cseed[iLayer].SetFreq(frequency);
2415 cseed[iLayer].SetC(cR);
2416 cseed[iLayer].SetCC(cC);
2417 cseed[iLayer].SetChi2(chi2TR);
2418 cseed[iLayer].SetChi2Z(chi2ZF);
69b55c55 2419 }
33744e5d 2420
2421 // Debugging print
4551ea7c 2422 if (1 || (!isFake)) {
69b55c55 2423 Float_t zvertex = GetZ();
4551ea7c 2424 TTreeSRedirector &cstream = *fDebugStreamer;
2425 if (AliTRDReconstructor::StreamLevel() > 0) {
2426 cstream << "Seeds1"
2427 << "isFake=" << isFake
2428 << "Vertex=" << zvertex
2429 << "Rieman2.=" << &rieman2
2430 << "Rieman.=" << &rieman
2431 << "Xref=" << xref
2432 << "X0=" << xcl[0]
2433 << "X1=" << xcl[1]
2434 << "X2=" << xcl[2]
2435 << "X3=" << xcl[3]
2436 << "X4=" << xcl[4]
2437 << "X5=" << xcl[5]
2438 << "Chi2R=" << chi2R
2439 << "Chi2Z=" << chi2Z
2440 << "Chi2RF=" << chi2RF // Chi2 of trackletes on full track
2441 << "Chi2ZF=" << chi2ZF // Chi2 z on tracklets on full track
2442 << "Chi2ZT2=" << chi2ZT2 // Chi2 z on tracklets on full track - rieman tilt
2443 << "Chi2ZTC=" << chi2ZTC // Chi2 z on tracklets on full track - rieman tilt const
2444 << "Chi2TR=" << chi2TR // Chi2 without vertex constrain
2445 << "Chi2TC=" << chi2TC // Chi2 with vertex constrain
2446 << "C=" << curv // Non constrained - no tilt correction
2447 << "DR=" << dR // DR parameter - tilt correction
2448 << "DCA=" << dca // DCA - tilt correction
2449 << "CR=" << cR // Non constrained curvature - tilt correction
2450 << "CC=" << cC // Constrained curvature
2451 << "Polz0=" << polz0c
2452 << "Polz1=" << polz1c
2453 << "RPolz0=" << rpolz0
2454 << "RPolz1=" << rpolz1
2455 << "Ncl=" << nclusters
2456 << "Nlayers=" << nlayers
2457 << "NUsedS=" << nusedCl
2458 << "NUsed=" << nusedf
2459 << "Findable=" << findable
2460 << "Like=" << like
2461 << "LikePrim=" << likePrim
2462 << "Likechi2C=" << likechi2C
2463 << "Likechi2TR=" << likechi2TR
2464 << "Likezf=" << likezf
2465 << "LikeF=" << seedquality2[registered]
2466 << "S0.=" << &cseed[0]
2467 << "S1.=" << &cseed[1]
2468 << "S2.=" << &cseed[2]
2469 << "S3.=" << &cseed[3]
2470 << "S4.=" << &cseed[4]
2471 << "S5.=" << &cseed[5]
2472 << "SB0.=" << &seedb[0]
2473 << "SB1.=" << &seedb[1]
2474 << "SB2.=" << &seedb[2]
2475 << "SB3.=" << &seedb[3]
2476 << "SB4.=" << &seedb[4]
2477 << "SB5.=" << &seedb[5]
2478 << "Label=" << label
2479 << "Freq=" << frequency
2480 << "sLayer=" << sLayer
2481 << "\n";
2482 }
69b55c55 2483 }
4551ea7c 2484
2485 if (registered<kMaxSeed - 1) {
69b55c55 2486 registered++;
2487 cseed = seed[registered];
2488 }
4551ea7c 2489
2490 } // End of loop over layer 1
2491
2492 } // End of loop over layer 0
2493
2494 } // End of loop over layer 3
2495
2496 } // End of loop over seeding time bins
2497
69b55c55 2498 //
4551ea7c 2499 // Choose best
69b55c55 2500 //
4551ea7c 2501
69b55c55 2502 TMath::Sort(registered,seedquality2,sort,kTRUE);
c6f438c0 2503 Bool_t signedseed[kMaxSeed];
4551ea7c 2504 for (Int_t i = 0; i < registered; i++) {
2505 signedseed[i] = kFALSE;
69b55c55 2506 }
4551ea7c 2507
2508 for (Int_t iter = 0; iter < 5; iter++) {
2509
2510 for (Int_t iseed = 0; iseed < registered; iseed++) {
2511
69b55c55 2512 Int_t index = sort[iseed];
4551ea7c 2513 if (signedseed[index]) {
2514 continue;
2515 }
69b55c55 2516 Int_t labelsall[1000];
4551ea7c 2517 Int_t nlabelsall = 0;
2518 Int_t naccepted = 0;;
2519 Int_t sLayer = seedlayer[index];
2520 Int_t ncl = 0;
2521 Int_t nused = 0;
2522 Int_t nlayers = 0;
69b55c55 2523 Int_t findable = 0;
4551ea7c 2524
2525 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2526
ca1e1e5b 2527 if (TMath::Abs(seed[index][jLayer].GetYref(0) / xcl[jLayer]) < 0.15) {
69b55c55 2528 findable++;
4551ea7c 2529 }
2530 if (seed[index][jLayer].IsOK()) {
69b55c55 2531 seed[index][jLayer].UpdateUsed();
ca1e1e5b 2532 ncl +=seed[index][jLayer].GetN2();
2533 nused +=seed[index][jLayer].GetNUsed();
69b55c55 2534 nlayers++;
4551ea7c 2535 // Cooking label
2536 for (Int_t itime = 0; itime < 25; itime++) {
ca1e1e5b 2537 if (seed[index][jLayer].IsUsable(itime)) {
69b55c55 2538 naccepted++;
4551ea7c 2539 for (Int_t ilab = 0; ilab < 3; ilab++) {
ca1e1e5b 2540 Int_t tindex = seed[index][jLayer].GetClusters(itime)->GetLabel(ilab);
4551ea7c 2541 if (tindex >= 0) {
69b55c55 2542 labelsall[nlabelsall] = tindex;
2543 nlabelsall++;
2544 }
2545 }
2546 }
2547 }
2548 }
4551ea7c 2549
69b55c55 2550 }
4551ea7c 2551
2552 if (nused > 30) {
2553 continue;
69b55c55 2554 }
4551ea7c 2555
2556 if (iter == 0) {
2557 if (nlayers < 6) {
2558 continue;
2559 }
2560 if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2561 continue; // Gold
2562 }
7ad19338 2563 }
4551ea7c 2564
2565 if (iter == 1) {
2566 if (nlayers < findable) {
2567 continue;
2568 }
2569 if (TMath::Log(0.000000001+seedquality2[index]) < -4.0) {
2570 continue;
2571 }
69b55c55 2572 }
4551ea7c 2573
2574 if (iter == 2) {
2575 if ((nlayers == findable) ||
2576 (nlayers == 6)) {
2577 continue;
2578 }
2579 if (TMath::Log(0.000000001+seedquality2[index]) < -6.0) {
2580 continue;
2581 }
69b55c55 2582 }
4551ea7c 2583
2584 if (iter == 3) {
2585 if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2586 continue;
2587 }
69b55c55 2588 }
4551ea7c 2589
2590 if (iter == 4) {
2591 if (TMath::Log(0.000000001+seedquality2[index]) - nused/(nlayers-3.0) < -15.0) {
2592 continue;
2593 }
2594 }
2595
69b55c55 2596 signedseed[index] = kTRUE;
4551ea7c 2597
2598 Int_t labels[1000];
2599 Int_t outlab[1000];
2600 Int_t nlab = 0;
2601 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2602 if (seed[index][iLayer].IsOK()) {
ca1e1e5b 2603 if (seed[index][iLayer].GetLabels(0) >= 0) {
2604 labels[nlab] = seed[index][iLayer].GetLabels(0);
69b55c55 2605 nlab++;
2606 }
ca1e1e5b 2607 if (seed[index][iLayer].GetLabels(1) >= 0) {
2608 labels[nlab] = seed[index][iLayer].GetLabels(1);
69b55c55 2609 nlab++;
4551ea7c 2610 }
2611 }
7ad19338 2612 }
69b55c55 2613 Freq(nlab,labels,outlab,kFALSE);
4551ea7c 2614 Int_t label = outlab[0];
2615 Int_t frequency = outlab[1];
69b55c55 2616 Freq(nlabelsall,labelsall,outlab,kFALSE);
4551ea7c 2617 Int_t label1 = outlab[0];
2618 Int_t label2 = outlab[2];
2619 Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
2620 Float_t ratio = Float_t(nused) / Float_t(ncl);
2621 if (ratio < 0.25) {
2622 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2623 if ((seed[index][jLayer].IsOK()) &&
ca1e1e5b 2624 (TMath::Abs(seed[index][jLayer].GetYfit(1) - seed[index][jLayer].GetYfit(1)) < 0.2)) {
4551ea7c 2625 seed[index][jLayer].UseClusters(); // Sign gold
2626 }
69b55c55 2627 }
7ad19338 2628 }
4551ea7c 2629
31fd97b2 2630 Int_t eventNrInFile = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
4551ea7c 2631 TTreeSRedirector &cstream = *fDebugStreamer;
2632
69b55c55 2633 //
4551ea7c 2634 // Register seed
69b55c55 2635 //
4551ea7c 2636 AliTRDtrack *track = RegisterSeed(seed[index],seedparams[index]);
2637 AliTRDtrack dummy;
2638 if (!track) {
2639 track = &dummy;
2640 }
2641 else {
69b55c55 2642 AliESDtrack esdtrack;
4551ea7c 2643 esdtrack.UpdateTrackParams(track,AliESDtrack::kTRDout);
69b55c55 2644 esdtrack.SetLabel(label);
2645 esd->AddTrack(&esdtrack);
4551ea7c 2646 TTreeSRedirector &cstream = *fDebugStreamer;
2647 if (AliTRDReconstructor::StreamLevel() > 0) {
2648 cstream << "Tracks"
31fd97b2 2649 << "EventNrInFile=" << eventNrInFile
4551ea7c 2650 << "ESD.=" << &esdtrack
2651 << "trd.=" << track
2652 << "trdback.=" << track
2653 << "\n";
2654 }
2655 }
2656
2657 if (AliTRDReconstructor::StreamLevel() > 0) {
2658 cstream << "Seeds2"
2659 << "Iter=" << iter
2660 << "Track.=" << track
2661 << "Like=" << seedquality[index]
2662 << "LikeF=" << seedquality2[index]
2663 << "S0.=" << &seed[index][0]
2664 << "S1.=" << &seed[index][1]
2665 << "S2.=" << &seed[index][2]
2666 << "S3.=" << &seed[index][3]
2667 << "S4.=" << &seed[index][4]
2668 << "S5.=" << &seed[index][5]
2669 << "Label=" << label
2670 << "Label1=" << label1
2671 << "Label2=" << label2
2672 << "FakeRatio=" << fakeratio
2673 << "Freq=" << frequency
2674 << "Ncl=" << ncl
2675 << "Nlayers=" << nlayers
2676 << "Findable=" << findable
2677 << "NUsed=" << nused
2678 << "sLayer=" << sLayer
31fd97b2 2679 << "EventNrInFile=" << eventNrInFile
4551ea7c 2680 << "\n";
7ad19338 2681 }
4551ea7c 2682
2683 } // Loop: iseed
2684
2685 } // Loop: iter
2686
2687 } // End of loop over sectors
75bd7f81 2688
69b55c55 2689 delete [] pseed;
75bd7f81 2690
69b55c55 2691}
2692
a819a5f7 2693//_____________________________________________________________________________
4551ea7c 2694Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *clusterTree) const
46d29e70 2695{
a819a5f7 2696 //
2697 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
2698 // from the file. The names of the cluster tree and branches
2699 // should match the ones used in AliTRDclusterizer::WriteClusters()
2700 //
75bd7f81 2701
4551ea7c 2702 Int_t nsize = Int_t(clusterTree->GetTotBytes() / (sizeof(AliTRDcluster)));
6c94f330 2703 TObjArray *clusterArray = new TObjArray(nsize+1000);
5443e65e 2704
4551ea7c 2705 TBranch *branch = clusterTree->GetBranch("TRDcluster");
c630aafd 2706 if (!branch) {
4551ea7c 2707 AliError("Can't get the branch !");
c630aafd 2708 return 1;
2709 }
029cd327 2710 branch->SetAddress(&clusterArray);
5443e65e 2711
a819a5f7 2712 // Loop through all entries in the tree
4551ea7c 2713 Int_t nEntries = (Int_t) clusterTree->GetEntries();
2714 Int_t nbytes = 0;
a819a5f7 2715 AliTRDcluster *c = 0;
a819a5f7 2716 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2717
2718 // Import the tree
4551ea7c 2719 nbytes += clusterTree->GetEvent(iEntry);
5443e65e 2720
a819a5f7 2721 // Get the number of points in the detector
029cd327 2722 Int_t nCluster = clusterArray->GetEntriesFast();
5443e65e 2723
a819a5f7 2724 // Loop through all TRD digits
2725 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
4551ea7c 2726 c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
4f1c04d3 2727 AliTRDcluster *co = c;
a819a5f7 2728 array->AddLast(co);
4f1c04d3 2729 clusterArray->RemoveAt(iCluster);
a819a5f7 2730 }
4551ea7c 2731
a819a5f7 2732 }
2733
029cd327 2734 delete clusterArray;
5443e65e 2735
c630aafd 2736 return 0;
75bd7f81 2737
a819a5f7 2738}
2739
75bd7f81 2740//_____________________________________________________________________________
4551ea7c 2741Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint &p) const
3551db50 2742{
2743 //
2744 // Get track space point with index i
2745 // Origin: C.Cheshkov
2746 //
2747
4551ea7c 2748 AliTRDcluster *cl = (AliTRDcluster *) fClusters->UncheckedAt(index);
2749 Int_t idet = cl->GetDetector();
2750 Int_t isector = fGeom->GetSector(idet);
2751 Int_t ichamber = fGeom->GetChamber(idet);
2752 Int_t iplan = fGeom->GetPlane(idet);
3551db50 2753 Double_t local[3];
4551ea7c 2754 local[0] = GetX(isector,iplan,cl->GetLocalTimeBin());
2755 local[1] = cl->GetY();
2756 local[2] = cl->GetZ();
3551db50 2757 Double_t global[3];
2758 fGeom->RotateBack(idet,local,global);
2759 p.SetXYZ(global[0],global[1],global[2]);
2760 AliAlignObj::ELayerID iLayer = AliAlignObj::kTRD1;
2761 switch (iplan) {
2762 case 0:
2763 iLayer = AliAlignObj::kTRD1;
2764 break;
2765 case 1:
2766 iLayer = AliAlignObj::kTRD2;
2767 break;
2768 case 2:
2769 iLayer = AliAlignObj::kTRD3;
2770 break;
2771 case 3:
2772 iLayer = AliAlignObj::kTRD4;
2773 break;
2774 case 4:
2775 iLayer = AliAlignObj::kTRD5;
2776 break;
2777 case 5:
2778 iLayer = AliAlignObj::kTRD6;
2779 break;
2780 };
4551ea7c 2781 Int_t modId = isector * fGeom->Ncham() + ichamber;
3551db50 2782 UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
2783 p.SetVolumeID(volid);
2784
2785 return kTRUE;
2786
2787}
2788
75bd7f81 2789//_____________________________________________________________________________
4551ea7c 2790void AliTRDtracker::CookLabel(AliKalmanTrack *pt, Float_t wrong) const
029cd327 2791{
2792 //
2793 // This cooks a label. Mmmmh, smells good...
2794 //
46d29e70 2795
4551ea7c 2796 Int_t label = 123456789;
2797 Int_t index;
2798 Int_t i;
2799 Int_t j;
2800 Int_t ncl = pt->GetNumberOfClusters();
2801
2802 const Int_t kRange = fTrSec[0]->GetOuterTimeBin() + 1;
5443e65e 2803
029cd327 2804 Bool_t labelAdded;
46d29e70 2805
6c94f330 2806 Int_t **s = new Int_t* [kRange];
4551ea7c 2807 for (i = 0; i < kRange; i++) {
d1b06c24 2808 s[i] = new Int_t[2];
2809 }
4551ea7c 2810 for (i = 0; i < kRange; i++) {
2811 s[i][0] = -1;
2812 s[i][1] = 0;
46d29e70 2813 }
2814
4551ea7c 2815 Int_t t0;
2816 Int_t t1;
2817 Int_t t2;
2818
2819 for (i = 0; i < ncl; i++) {
2820 index = pt->GetClusterIndex(i);
2821 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
6c94f330 2822 t0=c->GetLabel(0);
2823 t1=c->GetLabel(1);
2824 t2=c->GetLabel(2);
46d29e70 2825 }
2826
4551ea7c 2827 for (i = 0; i < ncl; i++) {
2828 index = pt->GetClusterIndex(i);
2829 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2830 for (Int_t k = 0; k < 3; k++) {
2831 label = c->GetLabel(k);
2832 labelAdded = kFALSE;
2833 j = 0;
46d29e70 2834 if (label >= 0) {
4551ea7c 2835 while ((!labelAdded) && (j < kRange)) {
2836 if ((s[j][0] == label) ||
2837 (s[j][1] == 0)) {
2838 s[j][0] = label;
2839 s[j][1] = s[j][1] + 1;
2840 labelAdded = kTRUE;
a9814c09 2841 }
2842 j++;
2843 }
46d29e70 2844 }
2845 }
2846 }
2847
4551ea7c 2848 Int_t max = 0;
6c94f330 2849 label = -123456789;
46d29e70 2850
4551ea7c 2851 for (i = 0; i < kRange; i++) {
2852 if (s[i][1] > max) {
2853 max = s[i][1];
2854 label = s[i][0];
46d29e70 2855 }
2856 }
5443e65e 2857
4551ea7c 2858 for (i = 0; i < kRange; i++) {
6c94f330 2859 delete []s[i];
5443e65e 2860 }
2861
6c94f330 2862 delete []s;
5443e65e 2863
4551ea7c 2864 if ((1.0 - Float_t(max)/ncl) > wrong) {
2865 label = -label;
2866 }
5443e65e 2867
2868 pt->SetLabel(label);
2869
46d29e70 2870}
2871
75bd7f81 2872//_____________________________________________________________________________
4551ea7c 2873void AliTRDtracker::UseClusters(const AliKalmanTrack *t, Int_t from) const
029cd327 2874{
2875 //
2876 // Use clusters, but don't abuse them!
2877 //
75bd7f81 2878
4551ea7c 2879 const Float_t kmaxchi2 = 18;
2880 const Float_t kmincl = 10;
2881
2882 AliTRDtrack *track = (AliTRDtrack *) t;
2883
2884 Int_t ncl = t->GetNumberOfClusters();
2885 for (Int_t i = from; i < ncl; i++) {
2886 Int_t index = t->GetClusterIndex(i);
2887 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
69b55c55 2888 Int_t iplane = fGeom->GetPlane(c->GetDetector());
27eaf44b 2889 if (track->GetTracklets(iplane).GetChi2() > kmaxchi2) {
4551ea7c 2890 continue;
2891 }
27eaf44b 2892 if (track->GetTracklets(iplane).GetN() < kmincl) {
4551ea7c 2893 continue;
2894 }
2895 if (!(c->IsUsed())) {
2896 c->Use();
2897 }
5443e65e 2898 }
5443e65e 2899
75bd7f81 2900}
5443e65e 2901
75bd7f81 2902//_____________________________________________________________________________
029cd327 2903Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
5443e65e 2904{
75bd7f81 2905 //
5443e65e 2906 // Parametrised "expected" error of the cluster reconstruction in Y
75bd7f81 2907 //
5443e65e 2908
2909 Double_t s = 0.08 * 0.08;
2910 return s;
75bd7f81 2911
5443e65e 2912}
2913
75bd7f81 2914//_____________________________________________________________________________
029cd327 2915Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
5443e65e 2916{
75bd7f81 2917 //
5443e65e 2918 // Parametrised "expected" error of the cluster reconstruction in Z
75bd7f81 2919 //
5443e65e 2920
4551ea7c 2921 Double_t s = 9.0 * 9.0 / 12.0;
5443e65e 2922 return s;
75bd7f81 2923
5443e65e 2924}
2925
75bd7f81 2926//_____________________________________________________________________________
029cd327 2927Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
5443e65e 2928{
2929 //
029cd327 2930 // Returns radial position which corresponds to time bin <localTB>
5443e65e 2931 // in tracking sector <sector> and plane <plane>
2932 //
2933
6c94f330 2934 Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
4551ea7c 2935 Int_t pl = fTrSec[sector]->GetLayerNumber(index);
2936
5443e65e 2937 return fTrSec[sector]->GetLayer(pl)->GetX();
2938
2939}
2940
75bd7f81 2941//_____________________________________________________________________________
2942AliTRDtracker::AliTRDpropagationLayer
2943 ::AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
2944 , Double_t radLength, Int_t tbIndex, Int_t plane)
ad670fba 2945 :fN(0)
2946 ,fSec(0)
2947 ,fClusters(NULL)
2948 ,fIndex(NULL)
2949 ,fX(x)
2950 ,fdX(dx)
2951 ,fRho(rho)
2952 ,fX0(radLength)
2953 ,fTimeBinIndex(tbIndex)
2954 ,fPlane(plane)
2955 ,fYmax(0)
2956 ,fYmaxSensitive(0)
2957 ,fHole(kFALSE)
2958 ,fHoleZc(0)
2959 ,fHoleZmax(0)
2960 ,fHoleYc(0)
2961 ,fHoleYmax(0)
2962 ,fHoleRho(0)
2963 ,fHoleX0(0)
5443e65e 2964{
2965 //
2966 // AliTRDpropagationLayer constructor
2967 //
2968
ad670fba 2969 for (Int_t i = 0; i < (Int_t) kZones; i++) {
2970 fZc[i] = 0;
2971 fZmax[i] = 0;
5443e65e 2972 }
2973
ad670fba 2974 if (fTimeBinIndex >= 0) {
029cd327 2975 fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
ad670fba 2976 fIndex = new UInt_t[kMaxClusterPerTimeBin];
5443e65e 2977 }
2978
ad670fba 2979 for (Int_t i = 0; i < 5; i++) {
2980 fIsHole[i] = kFALSE;
2981 }
5443e65e 2982
2983}
2984
75bd7f81 2985//_____________________________________________________________________________
2986void AliTRDtracker::AliTRDpropagationLayer
2987 ::SetHole(Double_t Zmax, Double_t Ymax, Double_t rho
2988 , Double_t radLength, Double_t Yc, Double_t Zc)
5443e65e 2989{
2990 //
2991 // Sets hole in the layer
2992 //
75bd7f81 2993
4551ea7c 2994 fHole = kTRUE;
2995 fHoleZc = Zc;
5443e65e 2996 fHoleZmax = Zmax;
4551ea7c 2997 fHoleYc = Yc;
5443e65e 2998 fHoleYmax = Ymax;
4551ea7c 2999 fHoleRho = rho;
3000 fHoleX0 = radLength;
75bd7f81 3001
5443e65e 3002}
5443e65e 3003
75bd7f81 3004//_____________________________________________________________________________
3005AliTRDtracker::AliTRDtrackingSector
ad670fba 3006 ::AliTRDtrackingSector(AliTRDgeometry *geo, Int_t gs)
3007 :fN(0)
3008 ,fGeom(geo)
3009 ,fGeomSector(gs)
0a29d0f1 3010{
3011 //
5443e65e 3012 // AliTRDtrackingSector Constructor
0a29d0f1 3013 //
75bd7f81 3014
7e448bcc 3015 AliTRDpadPlane *padPlane = 0;
3016 AliTRDpropagationLayer *ppl = 0;
a5cadd36 3017
ad670fba 3018 // Get holes description from geometry
3c625a9b 3019 Bool_t holes[AliTRDgeometry::kNcham];
ad670fba 3020 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3c625a9b 3021 holes[icham] = fGeom->IsHole(0,icham,gs);
3c625a9b 3022 }
3c625a9b 3023
ad670fba 3024 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
3025 fTimeBinIndex[i] = -1;
3026 }
5443e65e 3027
ad670fba 3028 Double_t x;
3029 Double_t dx;
3030 Double_t rho;
3031 Double_t radLength;
5443e65e 3032
ad670fba 3033 // Add layers for each of the planes
3034 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
a305677e 3035 //Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
5443e65e 3036
ad670fba 3037 const Int_t kNchambers = AliTRDgeometry::Ncham();
a305677e 3038 Int_t tbIndex;
ad670fba 3039 Double_t ymax = 0;
3040 Double_t ymaxsensitive = 0;
3041 Double_t *zc = new Double_t[kNchambers];
3042 Double_t *zmax = new Double_t[kNchambers];
3c625a9b 3043 Double_t *zmaxsensitive = new Double_t[kNchambers];
5443e65e 3044
ad670fba 3045 AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
3046 if (!commonParam) {
030b4415 3047 AliErrorGeneral("AliTRDtrackingSector::Ctor"
3048 ,"Could not get common parameters\n");
3551db50 3049 return;
3050 }
3051
ad670fba 3052 for (Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
3053
3054 ymax = fGeom->GetChamberWidth(plane) / 2.0;
3055 padPlane = commonParam->GetPadPlane(plane,0);
3056 ymaxsensitive = (padPlane->GetColSize(1) * padPlane->GetNcols() - 4.0) / 2.0;