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