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