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