]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDtracker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bug fixed in the $ field
[u/mrichter/AliRoot.git] / TRD / AliTRDtracker.cxx
CommitLineData
46d29e70 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
0fa7dfa7 16/* $Id$ */
bbf92647 17
029cd327 18///////////////////////////////////////////////////////////////////////////////
19// //
c6f438c0 20// The standard TRD tracker //
21// Based on Kalman filltering approach //
029cd327 22// //
4551ea7c 23// Authors: //
24// M. Ivanov (Marian.Ivanov@cern.ch) //
25// Y. Belikov (Jouri.Belikov@cern.ch) //
26// //
029cd327 27///////////////////////////////////////////////////////////////////////////////
28
7e448bcc 29
a2cb5b3d 30#include <Riostream.h>
46d29e70 31#include <TBranch.h>
a1e17193 32#include <TFile.h>
4551ea7c 33#include <TGraph.h>
7e448bcc 34#include <TH1D.h>
35#include <TH2D.h>
a1e17193 36#include <TLinearFitter.h>
37#include <TObjArray.h>
38#include <TROOT.h>
39#include <TTree.h>
40#include <TTreeStream.h>
4551ea7c 41
af885e0f 42#include "AliESDEvent.h"
f3f19a6e 43#include "AliESDtrack.h"
4551ea7c 44#include "AliAlignObj.h"
45#include "AliRieman.h"
46#include "AliTrackPointArray.h"
46d29e70 47
46d29e70 48#include "AliTRDgeometry.h"
a5cadd36 49#include "AliTRDpadPlane.h"
0b2318ec 50#include "AliTRDgeometry.h"
46d29e70 51#include "AliTRDcluster.h"
52#include "AliTRDtrack.h"
75bd7f81 53#include "AliTRDseed.h"
3551db50 54#include "AliTRDcalibDB.h"
55#include "AliTRDCommonParam.h"
46d29e70 56#include "AliTRDtracker.h"
7b580082 57#include "AliTRDReconstructor.h"
082a2283 58#include "AliTRDCalibraFillHisto.h"
59
7e448bcc 60ClassImp(AliTRDtracker)
46d29e70 61
b08a2e7e 62const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5; //
63const Float_t AliTRDtracker::fgkLabelFraction = 0.8; //
64const Double_t AliTRDtracker::fgkMaxChi2 = 12.0; //
7e448bcc 65const Double_t AliTRDtracker::fgkMaxSnp = 0.95; // Corresponds to tan = 3
66const Double_t AliTRDtracker::fgkMaxStep = 2.0; // Maximal step size in propagation
7ad19338 67
75bd7f81 68//_____________________________________________________________________________
4551ea7c 69AliTRDtracker::AliTRDtracker()
4e009ce4 70 :AliTracker()
71 ,fHBackfit(0x0)
72 ,fHClSearch(0x0)
73 ,fHRefit(0x0)
74 ,fHX(0x0)
75 ,fHNCl(0x0)
76 ,fHNClTrack(0x0)
77 ,fHMinYPos(0x0)
78 ,fHMinYNeg(0x0)
79 ,fHMinZ(0x0)
80 ,fHMinD(0x0)
81 ,fHDeltaX(0x0)
82 ,fHXCl(0x0)
83 ,fGeom(0)
84 ,fNclusters(0)
85 ,fClusters(0)
86 ,fNseeds(0)
87 ,fSeeds(0)
88 ,fNtracks(0)
89 ,fTracks(0)
90 ,fTimeBinsPerPlane(0)
91 ,fAddTRDseeds(kFALSE)
92 ,fNoTilt(kFALSE)
93 ,fDebugStreamer(0)
89f05372 94{
75bd7f81 95 //
b7a0917f 96 // Default constructor
75bd7f81 97 //
b7a0917f 98
4551ea7c 99 for (Int_t i = 0; i < kTrackingSectors; i++) {
100 fTrSec[i] = 0;
101 }
102 for (Int_t j = 0; j < 5; j++) {
103 for (Int_t k = 0; k < 18; k++) {
104 fHoles[j][k] = kFALSE;
105 }
106 }
75bd7f81 107
7e448bcc 108 InitLogHists();
4e009ce4 109
89f05372 110}
75bd7f81 111
112//_____________________________________________________________________________
6c94f330 113AliTRDtracker::AliTRDtracker(const AliTRDtracker &t)
114 :AliTracker(t)
4e009ce4 115 ,fHBackfit(0x0)
116 ,fHClSearch(0x0)
117 ,fHRefit(0x0)
118 ,fHX(0x0)
119 ,fHNCl(0x0)
120 ,fHNClTrack(0x0)
121 ,fHMinYPos(0x0)
122 ,fHMinYNeg(0x0)
123 ,fHMinZ(0x0)
124 ,fHMinD(0x0)
125 ,fHDeltaX(0x0)
126 ,fHXCl(0x0)
ad670fba 127 ,fGeom(0)
128 ,fNclusters(0)
129 ,fClusters(0)
130 ,fNseeds(0)
131 ,fSeeds(0)
132 ,fNtracks(0)
133 ,fTracks(0)
134 ,fTimeBinsPerPlane(0)
135 ,fAddTRDseeds(kFALSE)
136 ,fNoTilt(kFALSE)
137 ,fDebugStreamer(0)
6c94f330 138{
ad670fba 139 //
140 // Copy constructor
141 //
4e009ce4 142
ad670fba 143}
144
145//_____________________________________________________________________________
eccb46a3 146AliTRDtracker::AliTRDtracker(const TFile */*geomfile*/)
ad670fba 147 :AliTracker()
4e009ce4 148 ,fHBackfit(0x0)
149 ,fHClSearch(0x0)
150 ,fHRefit(0x0)
151 ,fHX(0x0)
152 ,fHNCl(0x0)
153 ,fHNClTrack(0x0)
154 ,fHMinYPos(0x0)
155 ,fHMinYNeg(0x0)
156 ,fHMinZ(0x0)
157 ,fHMinD(0x0)
158 ,fHDeltaX(0x0)
159 ,fHXCl(0x0)
ad670fba 160 ,fGeom(0)
161 ,fNclusters(0)
162 ,fClusters(new TObjArray(2000))
163 ,fNseeds(0)
164 ,fSeeds(new TObjArray(2000))
165 ,fNtracks(0)
166 ,fTracks(new TObjArray(1000))
167 ,fTimeBinsPerPlane(0)
168 ,fAddTRDseeds(kFALSE)
169 ,fNoTilt(kFALSE)
170 ,fDebugStreamer(0)
46d29e70 171{
5443e65e 172 //
173 // Main constructor
174 //
b8dc2353 175
4551ea7c 176 TDirectory *savedir = gDirectory;
46d29e70 177
f162af62 178 fGeom = new AliTRDgeometry();
c6f438c0 179 fGeom->ReadGeoMatrices();
c630aafd 180
4551ea7c 181 for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
439c63c8 182 Int_t trS = geomS;
4551ea7c 183 fTrSec[trS] = new AliTRDtrackingSector(fGeom,geomS);
184 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
eccb46a3 185 // Could also go ...
4551ea7c 186 fHoles[icham][trS] = fGeom->IsHole(0,icham,geomS);
3c625a9b 187 }
5443e65e 188 }
4551ea7c 189
f162af62 190 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0);
7ad19338 191 Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
4551ea7c 192 if (tiltAngle < 0.1) {
b8dc2353 193 fNoTilt = kTRUE;
194 }
195
fc88348b 196 if (!AliTRDcalibDB::Instance()) {
197 AliFatal("Could not get calibration object");
198 }
199 fTimeBinsPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
46d29e70 200
4551ea7c 201 fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
0a29d0f1 202
9c9d2487 203 savedir->cd();
7e448bcc 204
205 InitLogHists();
75bd7f81 206
5443e65e 207}
46d29e70 208
75bd7f81 209//_____________________________________________________________________________
5443e65e 210AliTRDtracker::~AliTRDtracker()
46d29e70 211{
029cd327 212 //
213 // Destructor of AliTRDtracker
214 //
215
89f05372 216 if (fClusters) {
217 fClusters->Delete();
218 delete fClusters;
219 }
4551ea7c 220
89f05372 221 if (fTracks) {
222 fTracks->Delete();
223 delete fTracks;
224 }
4551ea7c 225
89f05372 226 if (fSeeds) {
227 fSeeds->Delete();
228 delete fSeeds;
229 }
4551ea7c 230
f162af62 231 if (fGeom) {
232 delete fGeom;
233 }
0a29d0f1 234
4551ea7c 235 for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
029cd327 236 delete fTrSec[geomS];
5443e65e 237 }
4551ea7c 238
7ad19338 239 if (fDebugStreamer) {
7ad19338 240 delete fDebugStreamer;
241 }
9c9d2487 242
75bd7f81 243}
59393e34 244
75bd7f81 245//_____________________________________________________________________________
246Int_t AliTRDtracker::LocalToGlobalID(Int_t lid)
247{
59393e34 248 //
75bd7f81 249 // Transform internal TRD ID to global detector ID
59393e34 250 //
75bd7f81 251
4551ea7c 252 Int_t isector = fGeom->GetSector(lid);
253 Int_t ichamber = fGeom->GetChamber(lid);
254 Int_t iplan = fGeom->GetPlane(lid);
255
ae079791 256 AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
59393e34 257 switch (iplan) {
258 case 0:
ae079791 259 iLayer = AliGeomManager::kTRD1;
59393e34 260 break;
261 case 1:
ae079791 262 iLayer = AliGeomManager::kTRD2;
59393e34 263 break;
264 case 2:
ae079791 265 iLayer = AliGeomManager::kTRD3;
59393e34 266 break;
267 case 3:
ae079791 268 iLayer = AliGeomManager::kTRD4;
59393e34 269 break;
270 case 4:
ae079791 271 iLayer = AliGeomManager::kTRD5;
59393e34 272 break;
273 case 5:
ae079791 274 iLayer = AliGeomManager::kTRD6;
59393e34 275 break;
276 };
4551ea7c 277
278 Int_t modId = isector * fGeom->Ncham() + ichamber;
ae079791 279 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,modId);
75bd7f81 280
59393e34 281 return volid;
75bd7f81 282
59393e34 283}
284
75bd7f81 285//_____________________________________________________________________________
286Int_t AliTRDtracker::GlobalToLocalID(Int_t gid)
287{
59393e34 288 //
75bd7f81 289 // Transform global detector ID to local detector ID
59393e34 290 //
75bd7f81 291
4551ea7c 292 Int_t modId = 0;
ae079791 293 AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(gid,modId);
4551ea7c 294
295 Int_t isector = modId / fGeom->Ncham();
296 Int_t ichamber = modId % fGeom->Ncham();
297 Int_t iLayer = -1;
298
59393e34 299 switch (layerId) {
ae079791 300 case AliGeomManager::kTRD1:
6c94f330 301 iLayer = 0;
59393e34 302 break;
ae079791 303 case AliGeomManager::kTRD2:
6c94f330 304 iLayer = 1;
59393e34 305 break;
ae079791 306 case AliGeomManager::kTRD3:
6c94f330 307 iLayer = 2;
59393e34 308 break;
ae079791 309 case AliGeomManager::kTRD4:
6c94f330 310 iLayer = 3;
59393e34 311 break;
ae079791 312 case AliGeomManager::kTRD5:
6c94f330 313 iLayer = 4;
59393e34 314 break;
ae079791 315 case AliGeomManager::kTRD6:
6c94f330 316 iLayer = 5;
59393e34 317 break;
318 default:
6c94f330 319 iLayer =-1;
59393e34 320 }
4551ea7c 321
322 if (iLayer < 0) {
323 return -1;
324 }
325
59393e34 326 Int_t lid = fGeom->GetDetector(iLayer,ichamber,isector);
75bd7f81 327
59393e34 328 return lid;
59393e34 329
75bd7f81 330}
59393e34 331
75bd7f81 332//_____________________________________________________________________________
4551ea7c 333Bool_t AliTRDtracker::Transform(AliTRDcluster *cluster)
75bd7f81 334{
59393e34 335 //
eb37e0af 336 // Transform from cluster system to tracking system
c6f438c0 337 //
75bd7f81 338
33744e5d 339 // Magic constants for geo manager transformation
4551ea7c 340 const Double_t kX0shift = 2.52;
33744e5d 341
6c94f330 342 //
33744e5d 343 // Apply alignment and calibration to transform cluster
6c94f330 344 //
345 Int_t detector = cluster->GetDetector();
4551ea7c 346 Int_t plane = fGeom->GetPlane(cluster->GetDetector());
347 Int_t chamber = fGeom->GetChamber(cluster->GetDetector());
348 Int_t sector = fGeom->GetSector(cluster->GetDetector());
349
350 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
351 Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.0); // Drift distance
c6f438c0 352
59393e34 353 //
59393e34 354 // ExB correction
355 //
356 Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
4e009ce4 357 Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift,-AliTracker::GetBz()*0.1);
4551ea7c 358
f162af62 359 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(plane,chamber);
4551ea7c 360 Double_t zshiftIdeal = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
361 Double_t localPos[3];
362 Double_t localPosTracker[3];
c6f438c0 363 localPos[0] = -cluster->GetX();
4551ea7c 364 localPos[1] = cluster->GetY() - driftX * exB;
365 localPos[2] = cluster->GetZ() - zshiftIdeal;
366
c6f438c0 367 cluster->SetY(cluster->GetY() - driftX*exB);
368 Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
2fa01832 369 cluster->SetX(xplane - cluster->GetX());
4551ea7c 370
2fa01832 371 TGeoHMatrix *matrix = fGeom->GetCorrectionMatrix(cluster->GetDetector());
4551ea7c 372 if (!matrix) {
373 // No matrix found - if somebody used geometry with holes
c6f438c0 374 AliError("Invalid Geometry - Default Geometry used\n");
375 return kTRUE;
376 }
4551ea7c 377 matrix->LocalToMaster(localPos,localPosTracker);
378
379 if (AliTRDReconstructor::StreamLevel() > 1) {
380 (* fDebugStreamer) << "Transform"
381 << "Cl.=" << cluster
382 << "matrix.=" << matrix
383 << "Detector=" << detector
384 << "Sector=" << sector
385 << "Plane=" << plane
386 << "Chamber=" << chamber
387 << "lx0=" << localPosTracker[0]
388 << "ly0=" << localPosTracker[1]
389 << "lz0=" << localPosTracker[2]
390 << "\n";
c6f438c0 391 }
4551ea7c 392
eb37e0af 393 cluster->SetX(localPosTracker[0]+kX0shift);
c6f438c0 394 cluster->SetY(localPosTracker[1]);
395 cluster->SetZ(localPosTracker[2]);
75bd7f81 396
59393e34 397 return kTRUE;
75bd7f81 398
59393e34 399}
400
75bd7f81 401//_____________________________________________________________________________
4551ea7c 402// Bool_t AliTRDtracker::Transform(AliTRDcluster *cluster)
75bd7f81 403//{
c6f438c0 404// //
4551ea7c 405// // Is this still needed ????
c6f438c0 406// //
407// const Double_t kDriftCorrection = 1.01; // drift coeficient correction
408// const Double_t kTime0Cor = 0.32; // time0 correction
409// //
410// const Double_t kX0shift = 2.52;
411// const Double_t kX0shift5 = 3.05;
412
413// //
414// // apply alignment and calibration to transform cluster
415// //
416// //
417// Int_t detector = cluster->GetDetector();
418// Int_t plane = fGeom->GetPlane(cluster->GetDetector());
419// Int_t chamber = fGeom->GetChamber(cluster->GetDetector());
420// Int_t sector = fGeom->GetSector(cluster->GetDetector());
421
422// Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
423// Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.); // drift distance
424// //
425// // ExB correction
426// //
427// Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
4e009ce4 428// Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift,-AliTracker::GetBz()*0.1);
c6f438c0 429// //
430
431// AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
432// AliTRDpadPlane * padPlane = commonParam->GetPadPlane(plane,chamber);
433// Double_t zshiftIdeal = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
434// Double_t localPos[3], globalPos[3], localPosTracker[3], localPosTracker2[3];
435// localPos[2] = -cluster->GetX();
436// localPos[0] = cluster->GetY() - driftX*exB;
437// localPos[1] = cluster->GetZ() -zshiftIdeal;
438// TGeoHMatrix * matrix = fGeom->GetGeoMatrix(cluster->GetDetector());
439// matrix->LocalToMaster(localPos, globalPos);
440
441// Double_t sectorAngle = 20.*(sector%18)+10;
442// TGeoHMatrix rotSector;
443// rotSector.RotateZ(sectorAngle);
444// rotSector.LocalToMaster(globalPos, localPosTracker);
445// //
446// //
447// TGeoHMatrix matrix2(*matrix);
448// matrix2.MultiplyLeft(&rotSector);
449// matrix2.LocalToMaster(localPos,localPosTracker2);
450// //
451// //
452// //
453// cluster->SetY(cluster->GetY() - driftX*exB);
454// Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
455// cluster->SetX(xplane- kDriftCorrection*(cluster->GetX()-kTime0Cor));
456// (*fDebugStreamer)<<"Transform"<<
457// "Cl.="<<cluster<<
458// "matrix.="<<matrix<<
459// "matrix2.="<<&matrix2<<
460// "Detector="<<detector<<
461// "Sector="<<sector<<
462// "Plane="<<plane<<
463// "Chamber="<<chamber<<
464// "lx0="<<localPosTracker[0]<<
465// "ly0="<<localPosTracker[1]<<
466// "lz0="<<localPosTracker[2]<<
467// "lx2="<<localPosTracker2[0]<<
468// "ly2="<<localPosTracker2[1]<<
469// "lz2="<<localPosTracker2[2]<<
470// "\n";
471// //
472// if (plane==5)
473// cluster->SetX(localPosTracker[0]+kX0shift5);
474// else
475// cluster->SetX(localPosTracker[0]+kX0shift);
476
477// cluster->SetY(localPosTracker[1]);
478// cluster->SetZ(localPosTracker[2]);
479// return kTRUE;
480// }
481
75bd7f81 482//_____________________________________________________________________________
483Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track)
484{
9c9d2487 485 //
486 // Rotates the track when necessary
487 //
488
489 Double_t alpha = AliTRDgeometry::GetAlpha();
4551ea7c 490 Double_t y = track->GetY();
491 Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
9c9d2487 492
4551ea7c 493 // Is this still needed ????
c630aafd 494 //Int_t ns = AliTRDgeometry::kNsect;
9c9d2487 495 //Int_t s=Int_t(track->GetAlpha()/alpha)%ns;
496
4551ea7c 497 if (y > ymax) {
9c9d2487 498 //s = (s+1) % ns;
4551ea7c 499 if (!track->Rotate( alpha)) {
500 return kFALSE;
501 }
502 }
503 else if (y < -ymax) {
9c9d2487 504 //s = (s-1+ns) % ns;
4551ea7c 505 if (!track->Rotate(-alpha)) {
506 return kFALSE;
507 }
9c9d2487 508 }
509
510 return kTRUE;
9c9d2487 511
75bd7f81 512}
46e2d86c 513
75bd7f81 514//_____________________________________________________________________________
515AliTRDcluster *AliTRDtracker::GetCluster(AliTRDtrack *track, Int_t plane
516 , Int_t timebin, UInt_t &index)
517{
46e2d86c 518 //
75bd7f81 519 // Try to find cluster in the backup list
46e2d86c 520 //
75bd7f81 521
4551ea7c 522 AliTRDcluster *cl =0;
6c94f330 523 Int_t *indexes = track->GetBackupIndexes();
4551ea7c 524
525 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
526 if (indexes[i] == 0) {
527 break;
528 }
529 AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
530 if (!cli) {
531 break;
532 }
533 if (cli->GetLocalTimeBin() != timebin) {
534 continue;
535 }
46e2d86c 536 Int_t iplane = fGeom->GetPlane(cli->GetDetector());
4551ea7c 537 if (iplane == plane) {
538 cl = cli;
7ad19338 539 index = indexes[i];
46e2d86c 540 break;
541 }
542 }
75bd7f81 543
46e2d86c 544 return cl;
46e2d86c 545
75bd7f81 546}
3c625a9b 547
75bd7f81 548//_____________________________________________________________________________
4551ea7c 549Int_t AliTRDtracker::GetLastPlane(AliTRDtrack *track)
75bd7f81 550{
3c625a9b 551 //
75bd7f81 552 // Return last updated plane
553 //
554
4551ea7c 555 Int_t lastplane = 0;
556 Int_t *indexes = track->GetBackupIndexes();
557
558 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
559 AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
560 if (!cli) {
561 break;
562 }
3c625a9b 563 Int_t iplane = fGeom->GetPlane(cli->GetDetector());
4551ea7c 564 if (iplane > lastplane) {
3c625a9b 565 lastplane = iplane;
566 }
567 }
75bd7f81 568
3c625a9b 569 return lastplane;
75bd7f81 570
3c625a9b 571}
75bd7f81 572
573//_____________________________________________________________________________
af885e0f 574Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *event)
c630aafd 575{
576 //
577 // Finds tracks within the TRD. The ESD event is expected to contain seeds
578 // at the outer part of the TRD. The seeds
579 // are found within the TRD if fAddTRDseeds is TRUE.
580 // The tracks are propagated to the innermost time bin
581 // of the TRD and the ESD event is updated
582 //
583
4551ea7c 584 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
029cd327 585 Float_t foundMin = fgkMinClustersInTrack * timeBins;
4551ea7c 586 Int_t nseed = 0;
587 Int_t found = 0;
588 //Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
c630aafd 589
590 Int_t n = event->GetNumberOfTracks();
4551ea7c 591 for (Int_t i = 0; i < n; i++) {
592
593 AliESDtrack *seed = event->GetTrack(i);
594 ULong_t status = seed->GetStatus();
595 if ((status & AliESDtrack::kTRDout) == 0) {
596 continue;
597 }
598 if ((status & AliESDtrack::kTRDin) != 0) {
599 continue;
600 }
c630aafd 601 nseed++;
7ad19338 602
4551ea7c 603 AliTRDtrack *seed2 = new AliTRDtrack(*seed);
46e2d86c 604 //seed2->ResetCovariance();
4551ea7c 605 AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
606 AliTRDtrack &t = *pt;
7b580082 607 FollowProlongation(t);
c630aafd 608 if (t.GetNumberOfClusters() >= foundMin) {
609 UseClusters(&t);
4551ea7c 610 CookLabel(pt,1 - fgkLabelFraction);
611 //t.CookdEdx();
c630aafd 612 }
613 found++;
c630aafd 614
4551ea7c 615 Double_t xTPC = 250.0;
59393e34 616 if (PropagateToX(t,xTPC,fgkMaxStep)) {
c630aafd 617 seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
618 }
619 delete seed2;
620 delete pt;
4551ea7c 621
622 }
c630aafd 623
33744e5d 624 AliInfo(Form("Number of loaded seeds: %d",nseed));
625 AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
626 AliInfo(Form("Total number of found tracks: %d",found));
c630aafd 627
628 return 0;
75bd7f81 629
c630aafd 630}
5443e65e 631
c630aafd 632//_____________________________________________________________________________
af885e0f 633Int_t AliTRDtracker::PropagateBack(AliESDEvent *event)
75bd7f81 634{
c630aafd 635 //
636 // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
637 // backpropagated by the TPC tracker. Each seed is first propagated
638 // to the TRD, and then its prolongation is searched in the TRD.
639 // If sufficiently long continuation of the track is found in the TRD
640 // the track is updated, otherwise it's stored as originaly defined
641 // by the TPC tracker.
642 //
643
7e448bcc 644 Int_t found = 0; // number of tracks found
4551ea7c 645 Float_t foundMin = 20.0;
646 Int_t n = event->GetNumberOfTracks();
647
648 // Sort tracks
649 Float_t *quality = new Float_t[n];
650 Int_t *index = new Int_t[n];
651 for (Int_t i = 0; i < n; i++) {
652 AliESDtrack *seed = event->GetTrack(i);
4f1c04d3 653 Double_t covariance[15];
654 seed->GetExternalCovariance(covariance);
655 quality[i] = covariance[0]+covariance[2];
7e448bcc 656 //quality[i] = covariance[0];
4f1c04d3 657 }
658 TMath::Sort(n,quality,index,kFALSE);
4f1c04d3 659
4551ea7c 660 for (Int_t i = 0; i < n; i++) {
661
662 //AliESDtrack *seed = event->GetTrack(i);
663 AliESDtrack *seed = event->GetTrack(index[i]);
7e448bcc 664 fHBackfit->Fill(0);
c630aafd 665
4551ea7c 666 ULong_t status = seed->GetStatus();
667 if ((status & AliESDtrack::kTPCout) == 0) {
7e448bcc 668 fHBackfit->Fill(1);
4551ea7c 669 continue;
670 }
7e448bcc 671
4551ea7c 672 if ((status & AliESDtrack::kTRDout) != 0) {
7e448bcc 673 fHBackfit->Fill(2);
4551ea7c 674 continue;
675 }
676
677 Int_t lbl = seed->GetLabel();
c630aafd 678 AliTRDtrack *track = new AliTRDtrack(*seed);
679 track->SetSeedLabel(lbl);
4551ea7c 680 seed->UpdateTrackParams(track,AliESDtrack::kTRDbackup); // Make backup
c630aafd 681 fNseeds++;
4551ea7c 682 Float_t p4 = track->GetC();
f6625211 683 Int_t expectedClr = FollowBackProlongation(*track);
7e448bcc 684
685 fHBackfit->Fill(3);
686 fHX->Fill(track->GetX());
687
688
689 // store the last measurement
690 /*
691 fHNClTrack->Fill(track->GetNumberOfClusters());
692 if (track->GetNumberOfClusters() >= foundMin) {
693
694 fHBackfit->Fill(4);
695 track->CookdEdx();
696 CookdEdxTimBin(*track);
697 CookLabel(track,1 - fgkLabelFraction);
698 if (track->GetBackupTrack()) {
699 //fHBackfit->Fill(5);
700 UseClusters(track->GetBackupTrack());
701 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
702 }
703 }
704 */
705
706 /**/
707 // inter-tracks competition ???
4551ea7c 708 if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
709 (TMath::Abs(track->GetPt()) > 0.8)) {
7e448bcc 710
711 fHBackfit->Fill(4);
712
4f1c04d3 713 //
4551ea7c 714 // Make backup for back propagation
4f1c04d3 715 //
7e448bcc 716
4f1c04d3 717 Int_t foundClr = track->GetNumberOfClusters();
718 if (foundClr >= foundMin) {
719 track->CookdEdx();
43bfc8af 720 track->CookdEdxTimBin(); // A.Bercuci 25.07.07
4551ea7c 721 CookLabel(track,1 - fgkLabelFraction);
722 if (track->GetBackupTrack()) {
723 UseClusters(track->GetBackupTrack());
724 }
725
726 // Sign only gold tracks
727 if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
7e448bcc 728 if ((seed->GetKinkIndex(0) == 0) &&
4551ea7c 729 (TMath::Abs(track->GetPt()) < 1.5)) {
730 UseClusters(track);
731 }
4f1c04d3 732 }
733 Bool_t isGold = kFALSE;
734
4551ea7c 735 // Full gold track
736 if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
737 //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
738 if (track->GetBackupTrack()) {
739 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
7e448bcc 740
4551ea7c 741 }
4f1c04d3 742 isGold = kTRUE;
7e448bcc 743 //fHBackfit->Fill()
4f1c04d3 744 }
4551ea7c 745
746 // Almost gold track
747 if ((!isGold) &&
748 (track->GetNCross() == 0) &&
749 (track->GetChi2() / track->GetNumberOfClusters() < 7)) {
750 //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
751 if (track->GetBackupTrack()) {
752 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
753 }
f4e9508c 754 isGold = kTRUE;
755 }
4551ea7c 756
757 if ((!isGold) &&
758 (track->GetBackupTrack())) {
7e448bcc 759 if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) &&
760 ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
4551ea7c 761 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
4f1c04d3 762 isGold = kTRUE;
763 }
764 }
4551ea7c 765
7e448bcc 766 if ((track->StatusForTOF() > 0) &&
767 (track->GetNCross() == 0) &&
27eaf44b 768 (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected()) > 0.4)) {
6c94f330 769 //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
7ad19338 770 }
4551ea7c 771
16d9fbba 772 }
c630aafd 773 }
7e448bcc 774 /**/
4551ea7c 775
7e448bcc 776 /**/
8979685e 777 // Debug part of tracking
4551ea7c 778 TTreeSRedirector &cstream = *fDebugStreamer;
31fd97b2 779 Int_t eventNrInFile = event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
4551ea7c 780 if (AliTRDReconstructor::StreamLevel() > 0) {
781 if (track->GetBackupTrack()) {
782 cstream << "Tracks"
31fd97b2 783 << "EventNrInFile=" << eventNrInFile
4551ea7c 784 << "ESD.=" << seed
785 << "trd.=" << track
786 << "trdback.=" << track->GetBackupTrack()
787 << "\n";
788 }
789 else {
790 cstream << "Tracks"
31fd97b2 791 << "EventNrInFile=" << eventNrInFile
4551ea7c 792 << "ESD.=" << seed
793 << "trd.=" << track
794 << "trdback.=" << track
795 << "\n";
d337ef8d 796 }
8979685e 797 }
7e448bcc 798 /**/
4551ea7c 799
800 // Propagation to the TOF (I.Belikov)
801 if (track->GetStop() == kFALSE) {
7e448bcc 802 fHBackfit->Fill(5);
803
4551ea7c 804 Double_t xtof = 371.0;
7e448bcc 805 Double_t xTOF0 = 370.0;
806
4551ea7c 807 Double_t c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
808 if (TMath::Abs(c2) >= 0.99) {
7e448bcc 809
810 fHBackfit->Fill(6);
c5a8e3df 811 delete track;
812 continue;
813 }
7e448bcc 814
59393e34 815 PropagateToX(*track,xTOF0,fgkMaxStep);
4551ea7c 816
817 // Energy losses taken to the account - check one more time
818 c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
819 if (TMath::Abs(c2) >= 0.99) {
7e448bcc 820
821 fHBackfit->Fill(7);
4f1c04d3 822 delete track;
823 continue;
824 }
7e448bcc 825
826 //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
827 // fHBackfit->Fill(7);
828 //delete track;
829 // continue;
830 //}
4f1c04d3 831
4551ea7c 832 Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
833 Double_t y;
834 track->GetYAt(xtof,GetBz(),y);
7e448bcc 835 if (y > ymax) {
4551ea7c 836 if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
7e448bcc 837 fHBackfit->Fill(8);
7ac6fa52 838 delete track;
7bed16a7 839 continue;
7ac6fa52 840 }
4551ea7c 841 }
842 else if (y < -ymax) {
7ac6fa52 843 if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
7e448bcc 844 fHBackfit->Fill(9);
7ac6fa52 845 delete track;
7bed16a7 846 continue;
7ac6fa52 847 }
3c625a9b 848 }
7e448bcc 849
3c625a9b 850 if (track->PropagateTo(xtof)) {
4551ea7c 851 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
7e448bcc 852 fHBackfit->Fill(10);
853
4551ea7c 854 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
855 for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
6d45eaef 856 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
857 }
858 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
eab5961e 859 }
4551ea7c 860 //seed->SetTRDtrack(new AliTRDtrack(*track));
861 if (track->GetNumberOfClusters() > foundMin) {
7e448bcc 862 fHBackfit->Fill(11);
4551ea7c 863 found++;
864 }
3c625a9b 865 }
4551ea7c 866
867 }
868 else {
7e448bcc 869
870 fHBackfit->Fill(12);
871
4551ea7c 872 if ((track->GetNumberOfClusters() > 15) &&
873 (track->GetNumberOfClusters() > 0.5*expectedClr)) {
7e448bcc 874
4551ea7c 875 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
7e448bcc 876 fHBackfit->Fill(13);
877
16d9fbba 878 //seed->SetStatus(AliESDtrack::kTRDStop);
4551ea7c 879 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
880 for (Int_t j = 0; j <AliESDtrack::kNSlice; j++) {
6d45eaef 881 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
882 }
883 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
eab5961e 884 }
7ad19338 885 //seed->SetTRDtrack(new AliTRDtrack(*track));
3c625a9b 886 found++;
887 }
4551ea7c 888
1e9bb598 889 }
4551ea7c 890
7ad19338 891 seed->SetTRDQuality(track->StatusForTOF());
27eaf44b 892 seed->SetTRDBudget(track->GetBudget(0));
8979685e 893
7e448bcc 894 fHBackfit->Fill(14);
d9b8978b 895 delete track;
c630aafd 896 }
ad670fba 897
33744e5d 898 AliInfo(Form("Number of seeds: %d",fNseeds));
899 AliInfo(Form("Number of back propagated TRD tracks: %d",found));
6c94f330 900
33744e5d 901 // New seeding
902 if (AliTRDReconstructor::SeedingOn()) {
4551ea7c 903 MakeSeedsMI(3,5,event);
33744e5d 904 }
905
906 fSeeds->Clear();
4551ea7c 907 fNseeds = 0;
7ad19338 908
4f1c04d3 909 delete [] index;
910 delete [] quality;
911
7e448bcc 912 SaveLogHists();
913
1e9bb598 914 return 0;
1e9bb598 915}
916
917//_____________________________________________________________________________
af885e0f 918Int_t AliTRDtracker::RefitInward(AliESDEvent *event)
1e9bb598 919{
920 //
921 // Refits tracks within the TRD. The ESD event is expected to contain seeds
922 // at the outer part of the TRD.
923 // The tracks are propagated to the innermost time bin
924 // of the TRD and the ESD event is updated
925 // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
926 //
927
4551ea7c 928 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
1e9bb598 929 Float_t foundMin = fgkMinClustersInTrack * timeBins;
4551ea7c 930 Int_t nseed = 0;
931 Int_t found = 0;
932 //Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
4f1c04d3 933 AliTRDtrack seed2;
6c94f330 934
1e9bb598 935 Int_t n = event->GetNumberOfTracks();
4551ea7c 936 for (Int_t i = 0; i < n; i++) {
937
938 AliESDtrack *seed = event->GetTrack(i);
939 new (&seed2) AliTRDtrack(*seed);
7e448bcc 940 fHRefit->Fill(0);
941
4551ea7c 942 if (seed2.GetX() < 270.0) {
943 seed->UpdateTrackParams(&seed2,AliESDtrack::kTRDbackup); // Backup TPC track - only update
7e448bcc 944 fHRefit->Fill(1);
f4e9508c 945 continue;
946 }
947
4551ea7c 948 ULong_t status = seed->GetStatus();
949 if ((status & AliESDtrack::kTRDout) == 0) {
7e448bcc 950 fHRefit->Fill(2);
0dd7d129 951 continue;
952 }
4551ea7c 953 if ((status & AliESDtrack::kTRDin) != 0) {
7e448bcc 954 fHRefit->Fill(3);
0dd7d129 955 continue;
956 }
7e448bcc 957
6c94f330 958 nseed++;
7e448bcc 959 fHRefit->Fill(4);
6c94f330 960
4551ea7c 961 seed2.ResetCovariance(50.0);
6c94f330 962
4f1c04d3 963 AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
4551ea7c 964 Int_t *indexes2 = seed2.GetIndexes();
965 for (Int_t i = 0; i < AliESDtrack::kNPlane;i++) {
966 for (Int_t j = 0; j < AliESDtrack::kNSlice;j++) {
6d45eaef 967 pt->SetPIDsignals(seed2.GetPIDsignals(i,j),i,j);
968 }
7ad19338 969 pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i);
970 }
eab5961e 971
4551ea7c 972 Int_t *indexes3 = pt->GetBackupIndexes();
973 for (Int_t i = 0; i < 200;i++) {
974 if (indexes2[i] == 0) {
975 break;
976 }
46e2d86c 977 indexes3[i] = indexes2[i];
4551ea7c 978 }
979
46e2d86c 980 //AliTRDtrack *pt = seed2;
43bfc8af 981 // A.Bercuci 25.07.07
982 //AliTRDtrack &t = *pt;
983 FollowProlongation(*pt);
984 if (pt->GetNumberOfClusters() >= foundMin) {
4551ea7c 985 //UseClusters(&t);
6c94f330 986 //CookLabel(pt, 1-fgkLabelFraction);
43bfc8af 987 pt->CookdEdx();
988 pt->CookdEdxTimBin();
989 pt->CookPID(seed);
990 //pt->Calibrate(); // slot for calibration
1e9bb598 991 }
992 found++;
33744e5d 993
4551ea7c 994 Double_t xTPC = 250.0;
43bfc8af 995 if (PropagateToX(*pt,xTPC,fgkMaxStep)) {
7e448bcc 996
4551ea7c 997 seed->UpdateTrackParams(pt,AliESDtrack::kTRDrefit);
7e448bcc 998 fHRefit->Fill(5);
999
4551ea7c 1000 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
1001 for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
6d45eaef 1002 seed->SetTRDsignals(pt->GetPIDsignals(i,j),i,j);
1003 }
7ad19338 1004 seed->SetTRDTimBin(pt->GetPIDTimBin(i),i);
1005 }
4551ea7c 1006 }
1007 else {
1008 // If not prolongation to TPC - propagate without update
7e448bcc 1009 fHRefit->Fill(5);
4551ea7c 1010 AliTRDtrack *seed2 = new AliTRDtrack(*seed);
1011 seed2->ResetCovariance(5.0);
1012 AliTRDtrack *pt2 = new AliTRDtrack(*seed2,seed2->GetAlpha());
7bed16a7 1013 delete seed2;
59393e34 1014 if (PropagateToX(*pt2,xTPC,fgkMaxStep)) {
4551ea7c 1015 pt2->CookdEdx( );
43bfc8af 1016 pt2->CookdEdxTimBin(); // A.Bercuci 25.07.07
1017 seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit);
1018 fHRefit->Fill(6);
7e448bcc 1019
4551ea7c 1020 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
1021 for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
6d45eaef 1022 seed->SetTRDsignals(pt2->GetPIDsignals(i,j),i,j);
43bfc8af 1023 }
7ad19338 1024 seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
1025 }
eab5961e 1026 }
43bfc8af 1027 // A.Bercuci 25.07.07
1028 // Add TRD track to ESDfriendTrack - maybe this tracks are not useful for post-processing - TODO make decission
1029 if (AliTRDReconstructor::StreamLevel()>0) seed->AddCalibObject(new AliTRDtrack(*pt2/*, kTRUE*/));
7bed16a7 1030 delete pt2;
4551ea7c 1031 }
1032
1e9bb598 1033 delete pt;
43bfc8af 1034 // A.Bercuci 25.07.07
1035 // Add TRD track to ESDfriendTrack
1036 if (AliTRDReconstructor::StreamLevel()>0) seed->AddCalibObject(new AliTRDtrack(*pt/*, kTRUE*/));
4551ea7c 1037
eab5961e 1038 }
1e9bb598 1039
33744e5d 1040 AliInfo(Form("Number of loaded seeds: %d",nseed));
1041 AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
1e9bb598 1042
7e448bcc 1043 SaveLogHists();
c630aafd 1044 return 0;
c630aafd 1045}
1046
75bd7f81 1047//_____________________________________________________________________________
4551ea7c 1048Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t)
8979685e 1049{
43bfc8af 1050 //
1051 // Starting from current position on track=t this function tries
1052 // to extrapolate the track up to timeBin=0 and to confirm prolongation
1053 // if a close cluster is found. Returns the number of clusters
1054 // expected to be found in sensitive layers
1055 // GeoManager used to estimate mean density
1056 //
4551ea7c 1057
43bfc8af 1058 Int_t sector;
1059 Int_t lastplane = GetLastPlane(&t);
1060 Double_t radLength = 0.0;
1061 Double_t rho = 0.0;
1062 Int_t expectedNumberOfClusters = 0;
4551ea7c 1063
43bfc8af 1064 for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
4551ea7c 1065
43bfc8af 1066 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
1067 Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
4551ea7c 1068
43bfc8af 1069 //
1070 // Propagate track close to the plane if neccessary
1071 //
1072 Double_t currentx = fTrSec[0]->GetLayer(rowlast)->GetX();
1073 if (currentx < (-fgkMaxStep + t.GetX())) {
1074 // Propagate closer to chamber - safety space fgkMaxStep
1075 if (!PropagateToX(t,currentx+fgkMaxStep,fgkMaxStep)) {
1076 break;
1077 }
1078 }
4551ea7c 1079
43bfc8af 1080 if (!AdjustSector(&t)) {
1081 break;
1082 }
4551ea7c 1083
43bfc8af 1084 //
1085 // Get material budget
1086 //
1087 Double_t xyz0[3];
1088 Double_t xyz1[3];
1089 Double_t param[7];
1090 Double_t x;
1091 Double_t y;
1092 Double_t z;
1093
1094 // Starting global position
1095 t.GetXYZ(xyz0);
1096 // End global position
1097 x = fTrSec[0]->GetLayer(row0)->GetX();
1098 if (!t.GetProlongation(x,y,z)) {
1099 break;
1100 }
1101 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1102 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1103 xyz1[2] = z;
1104 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1105 rho = param[0];
1106 radLength = param[1]; // Get mean propagation parameters
1107 // A.Bercuci 25.07.07
1108 // Flags for marking the track momentum and direction. The track is
1109 // marked only if it has at least 1 cluster picked up in the current
1110 // chamber.
1111 Bool_t kUPDATED = kFALSE, kMARKED = kFALSE;
4551ea7c 1112
43bfc8af 1113 //
1114 // Propagate and update
1115 //
1116 sector = t.GetSector();
1117 //for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
1118 for (Int_t itime = 0 ; itime < GetTimeBinsPerPlane(); itime++) {
1119 // A.Bercuci 25.07.07
1120 // Mark track kinematics
1121 if(itime > 10 && kUPDATED && !kMARKED){
1122 t.SetTrackSegmentDirMom(iplane);
1123 kMARKED = kTRUE;
1124 }
1125
1126 Int_t ilayer = GetGlobalTimeBin(0,iplane,itime);
1127 expectedNumberOfClusters++;
1128 t.SetNExpected(t.GetNExpected() + 1);
1129 if (t.GetX() > 345.0) {
1130 t.SetNExpectedLast(t.GetNExpectedLast() + 1);
1131 }
1132 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
1133 AliTRDcluster *cl = 0;
1134 UInt_t index = 0;
1135 Double_t maxChi2 = fgkMaxChi2;
1136 x = timeBin.GetX();
1137
1138 if (timeBin) {
1139
1140 AliTRDcluster *cl0 = timeBin[0];
1141 if (!cl0) {
1142 // No clusters in given time bin
1143 continue;
1144 }
1145
1146 Int_t plane = fGeom->GetPlane(cl0->GetDetector());
1147 if (plane > lastplane) {
1148 continue;
1149 }
1150
1151 Int_t timebin = cl0->GetLocalTimeBin();
1152 AliTRDcluster *cl2 = GetCluster(&t,plane,timebin,index);
1153
1154 if (cl2) {
1155 cl = cl2;
1156 //Double_t h01 = GetTiltFactor(cl); //I.B's fix
1157 //maxChi2=t.GetPredictedChi2(cl,h01);
1158 }
1159 if (cl) {
1160 //if (cl->GetNPads()<5)
1161 Double_t dxsample = timeBin.GetdX();
1162 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
1163 Double_t h01 = GetTiltFactor(cl);
1164 Int_t det = cl->GetDetector();
1165 Int_t plane = fGeom->GetPlane(det);
1166 if (t.GetX() > 345.0) {
1167 t.SetNLast(t.GetNLast() + 1);
1168 t.SetChi2Last(t.GetChi2Last() + maxChi2);
1169 }
1170
1171 Double_t xcluster = cl->GetX();
1172 t.PropagateTo(xcluster,radLength,rho);
1173
1174 if (!AdjustSector(&t)) {
1175 break; //I.B's fix
1176 }
1177 maxChi2 = t.GetPredictedChi2(cl,h01);
1178
1179 if (maxChi2<1e+10)
1180 if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1181 // ????
1182 } else {
1183 // A.Bercuci 25.07.07
1184 //SetCluster(cl, GetNumberOfClusters()-1);
1185 kUPDATED = kTRUE;
1186 }
1187 }
1188 }
1189 }
4551ea7c 1190 }
1191
43bfc8af 1192 return expectedNumberOfClusters;
75bd7f81 1193}
7b580082 1194
75bd7f81 1195//_____________________________________________________________________________
4551ea7c 1196Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
8979685e 1197{
75bd7f81 1198 //
8979685e 1199 // Starting from current radial position of track <t> this function
1200 // extrapolates the track up to outer timebin and in the sensitive
1201 // layers confirms prolongation if a close cluster is found.
1202 // Returns the number of clusters expected to be found in sensitive layers
1203 // Use GEO manager for material Description
75bd7f81 1204 //
7e448bcc 1205 // return number of assigned clusters ?
1206 //
1207
59393e34 1208
4551ea7c 1209 Int_t sector;
1210 Int_t clusters[1000];
6c94f330 1211 Double_t radLength = 0.0;
4551ea7c 1212 Double_t rho = 0.0;
1213 Int_t expectedNumberOfClusters = 0;
1214 Float_t ratio0 = 0.0;
8979685e 1215 AliTRDtracklet tracklet;
33744e5d 1216
77566f2a 1217 // Calibration fill 2D
082a2283 1218 AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
77566f2a 1219 if (!calibra) {
1220 AliInfo("Could not get Calibra instance\n");
1221 }
8ec526a4 1222 if (calibra->GetMITracking()) {
1223 calibra->ResetTrack();
77566f2a 1224 }
1225
33744e5d 1226 for (Int_t i = 0; i < 1000; i++) {
1227 clusters[i] = -1;
1228 }
1229
4551ea7c 1230 for (Int_t iplane = 0; iplane < AliESDtrack::kNPlane; iplane++) {
1231
7e448bcc 1232 int hb = iplane * 10;
1233 fHClSearch->Fill(hb);
1234
4551ea7c 1235 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
1236 Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
1237 Double_t currentx = fTrSec[0]->GetLayer(row0)->GetX();
1238 if (currentx < t.GetX()) {
7e448bcc 1239 fHClSearch->Fill(hb+1);
4551ea7c 1240 continue;
1241 }
1242
6c94f330 1243 //
4551ea7c 1244 // Propagate closer to chamber if neccessary
6c94f330 1245 //
4551ea7c 1246 if (currentx > (fgkMaxStep + t.GetX())) {
1247 if (!PropagateToX(t,currentx-fgkMaxStep,fgkMaxStep)) {
7e448bcc 1248 fHClSearch->Fill(hb+2);
4551ea7c 1249 break;
1250 }
1251 }
1252 if (!AdjustSector(&t)) {
7e448bcc 1253 fHClSearch->Fill(hb+3);
4551ea7c 1254 break;
8979685e 1255 }
4551ea7c 1256 if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) {
7e448bcc 1257 fHClSearch->Fill(hb+4);
4551ea7c 1258 break;
1259 }
1260
8979685e 1261 //
4551ea7c 1262 // Get material budget inside of chamber
59393e34 1263 //
4551ea7c 1264 Double_t xyz0[3];
1265 Double_t xyz1[3];
1266 Double_t param[7];
1267 Double_t x;
1268 Double_t y;
1269 Double_t z;
1270 // Starting global position
1271 t.GetXYZ(xyz0);
1272 // End global position
7b580082 1273 x = fTrSec[0]->GetLayer(rowlast)->GetX();
4551ea7c 1274 if (!t.GetProlongation(x,y,z)) {
7e448bcc 1275 fHClSearch->Fill(hb+5);
4551ea7c 1276 break;
1277 }
1278 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1279 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1280 xyz1[2] = z;
7b580082 1281 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
4551ea7c 1282 rho = param[0];
1283 radLength = param[1]; // Get mean propagation parameters
1284
8979685e 1285 //
7b580082 1286 // Find clusters
8979685e 1287 //
6c94f330 1288 sector = t.GetSector();
4551ea7c 1289 Float_t ncl = FindClusters(sector,row0,rowlast,&t,clusters,tracklet);
7e448bcc 1290 fHNCl->Fill(tracklet.GetN());
1291
4551ea7c 1292 if (tracklet.GetN() < GetTimeBinsPerPlane()/3) {
7e448bcc 1293 fHClSearch->Fill(hb+6);
4551ea7c 1294 continue;
1295 }
1296
8979685e 1297 //
7b580082 1298 // Propagate and update track
8979685e 1299 //
4551ea7c 1300 for (Int_t itime = GetTimeBinsPerPlane()-1; itime >= 0; itime--) {
1301
1302 Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
8979685e 1303 expectedNumberOfClusters++;
27eaf44b 1304 t.SetNExpected(t.GetNExpected() + 1);
4551ea7c 1305 if (t.GetX() > 345.0) {
27eaf44b 1306 t.SetNExpectedLast(t.GetNExpectedLast() + 1);
4551ea7c 1307 }
1308 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
7e448bcc 1309 AliTRDcluster *cl = 0;
4551ea7c 1310 UInt_t index = 0;
1311 Double_t maxChi2 = fgkMaxChi2;
8979685e 1312 x = timeBin.GetX();
4551ea7c 1313
8979685e 1314 if (timeBin) {
4551ea7c 1315
1316 if (clusters[ilayer] > 0) {
8979685e 1317 index = clusters[ilayer];
4551ea7c 1318 cl = (AliTRDcluster *)GetCluster(index);
33744e5d 1319 //Double_t h01 = GetTiltFactor(cl); // I.B's fix
1320 //maxChi2=t.GetPredictedChi2(cl,h01); //
8979685e 1321 }
1322
1323 if (cl) {
4551ea7c 1324
1325 //if (cl->GetNPads() < 5)
59393e34 1326 Double_t dxsample = timeBin.GetdX();
1327 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
4551ea7c 1328 Double_t h01 = GetTiltFactor(cl);
1329 Int_t det = cl->GetDetector();
1330 Int_t plane = fGeom->GetPlane(det);
1331 if (t.GetX() > 345.0) {
27eaf44b 1332 t.SetNLast(t.GetNLast() + 1);
1333 t.SetChi2Last(t.GetChi2Last() + maxChi2);
8979685e 1334 }
59393e34 1335 Double_t xcluster = cl->GetX();
1336 t.PropagateTo(xcluster,radLength,rho);
4551ea7c 1337 maxChi2 = t.GetPredictedChi2(cl,h01);
6c94f330 1338
c2239ca8 1339 if (maxChi2<1e+10)
1340 if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1341 if (!t.Update(cl,maxChi2,index,h01)) {
1342 // ????
1343 }
43bfc8af 1344 } // else SetCluster(cl, GetNumberOfClusters()-1); // A.Bercuci 25.07.07
1345
4551ea7c 1346
8ec526a4 1347 if (calibra->GetMITracking()) {
77566f2a 1348 calibra->UpdateHistograms(cl,&t);
1349 }
1350
4551ea7c 1351 // Reset material budget if 2 consecutive gold
1352 if (plane > 0) {
27eaf44b 1353 if ((t.GetTracklets(plane).GetN() + t.GetTracklets(plane-1).GetN()) > 20) {
1354 t.SetBudget(2,0.0);
4551ea7c 1355 }
1356 }
1357
1358 }
1359
8979685e 1360 }
4551ea7c 1361
59393e34 1362 }
4551ea7c 1363
1364 ratio0 = ncl / Float_t(fTimeBinsPerPlane);
27eaf44b 1365 Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);
1366 if ((tracklet.GetChi2() < 18.0) &&
1367 (ratio0 > 0.8) &&
1368 (ratio1 > 0.6) &&
1369 (ratio0+ratio1 > 1.5) &&
1370 (t.GetNCross() == 0) &&
1371 (TMath::Abs(t.GetSnp()) < 0.85) &&
1372 (t.GetNumberOfClusters() > 20)){
7e448bcc 1373 //if (ratio0 > 0.8) {
4551ea7c 1374 t.MakeBackupTrack(); // Make backup of the track until is gold
59393e34 1375 }
7b580082 1376
8979685e 1377 }
5443e65e 1378
75bd7f81 1379 return expectedNumberOfClusters;
75bd7f81 1380}
1e9bb598 1381
75bd7f81 1382//_____________________________________________________________________________
4551ea7c 1383Int_t AliTRDtracker::PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep)
5443e65e 1384{
75bd7f81 1385 //
5443e65e 1386 // Starting from current radial position of track <t> this function
1387 // extrapolates the track up to radial position <xToGo>.
1388 // Returns 1 if track reaches the plane, and 0 otherwise
75bd7f81 1389 //
1390
59393e34 1391 const Double_t kEpsilon = 0.00001;
4551ea7c 1392 //Double_t tanmax = TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
1393 Double_t xpos = t.GetX();
1394 Double_t dir = (xpos<xToGo) ? 1.0 : -1.0;
1395
1396 while (((xToGo-xpos)*dir) > kEpsilon) {
1397
1398 Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos),maxStep);
1399
1400 Double_t xyz0[3];
1401 Double_t xyz1[3];
1402 Double_t param[7];
1403 Double_t x;
1404 Double_t y;
1405 Double_t z;
1406 // Starting global position
1407 t.GetXYZ(xyz0);
1408 x = xpos + step;
1409
1410 if (!t.GetProlongation(x,y,z)) {
1411 return 0; // No prolongation
1412 }
1413
1414 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1415 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1416 xyz1[2] = z;
1417
59393e34 1418 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
4551ea7c 1419 if (!t.PropagateTo(x,param[1],param[0])) {
1420 return 0;
1421 }
59393e34 1422 AdjustSector(&t);
1423 xpos = t.GetX();
4551ea7c 1424
5443e65e 1425 }
75bd7f81 1426
5443e65e 1427 return 1;
5443e65e 1428
59393e34 1429}
5443e65e 1430
c630aafd 1431//_____________________________________________________________________________
1432Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1433{
75bd7f81 1434 //
c630aafd 1435 // Fills clusters into TRD tracking_sectors
1436 // Note that the numbering scheme for the TRD tracking_sectors
1437 // differs from that of TRD sectors
75bd7f81 1438 //
1439
c630aafd 1440 if (ReadClusters(fClusters,cTree)) {
4551ea7c 1441 AliError("Problem with reading the clusters !");
c630aafd 1442 return 1;
1443 }
4551ea7c 1444 Int_t ncl = fClusters->GetEntriesFast();
1445 fNclusters = ncl;
1446 AliInfo(Form("Sorting %d clusters",ncl));
c630aafd 1447
1448 UInt_t index;
4551ea7c 1449 for (Int_t ichamber = 0; ichamber < 5; ichamber++) {
1450 for (Int_t isector = 0; isector < 18; isector++) {
1451 fHoles[ichamber][isector] = kTRUE;
3c625a9b 1452 }
4551ea7c 1453 }
ad670fba 1454
6c94f330 1455 while (ncl--) {
33744e5d 1456
4551ea7c 1457 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(ncl);
1458 Int_t detector = c->GetDetector();
1459 Int_t localTimeBin = c->GetLocalTimeBin();
1460 Int_t sector = fGeom->GetSector(detector);
1461 Int_t plane = fGeom->GetPlane(detector);
439c63c8 1462 Int_t trackingSector = sector;
4551ea7c 1463
7e448bcc 1464 //if (c->GetLabel(0) > 0) {
1465 if (c->GetQ() > 10) {
3c625a9b 1466 Int_t chamber = fGeom->GetChamber(detector);
4551ea7c 1467 fHoles[chamber][trackingSector] = kFALSE;
3c625a9b 1468 }
c630aafd 1469
6c94f330 1470 Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
4551ea7c 1471 if (gtb < 0) {
1472 continue;
1473 }
029cd327 1474 Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
c630aafd 1475
4551ea7c 1476 index = ncl;
33744e5d 1477
1478 // Apply pos correction
7e448bcc 1479 Transform(c);
1480 fHXCl->Fill(c->GetX());
33744e5d 1481
7e448bcc 1482 fTrSec[trackingSector]->GetLayer(layer)->SetX(c->GetX());
1483 fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
4551ea7c 1484 }
75bd7f81 1485
c630aafd 1486 return 0;
75bd7f81 1487
c630aafd 1488}
1489
5443e65e 1490//_____________________________________________________________________________
b7a0917f 1491void AliTRDtracker::UnloadClusters()
5443e65e 1492{
1493 //
1494 // Clears the arrays of clusters and tracks. Resets sectors and timebins
1495 //
1496
4551ea7c 1497 Int_t i;
1498 Int_t nentr;
5443e65e 1499
1500 nentr = fClusters->GetEntriesFast();
4551ea7c 1501 for (i = 0; i < nentr; i++) {
1502 delete fClusters->RemoveAt(i);
1503 }
b7a0917f 1504 fNclusters = 0;
5443e65e 1505
1506 nentr = fSeeds->GetEntriesFast();
4551ea7c 1507 for (i = 0; i < nentr; i++) {
1508 delete fSeeds->RemoveAt(i);
1509 }
5443e65e 1510
1511 nentr = fTracks->GetEntriesFast();
4551ea7c 1512 for (i = 0; i < nentr; i++) {
1513 delete fTracks->RemoveAt(i);
1514 }
5443e65e 1515
1516 Int_t nsec = AliTRDgeometry::kNsect;
5443e65e 1517 for (i = 0; i < nsec; i++) {
1518 for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1519 fTrSec[i]->GetLayer(pl)->Clear();
1520 }
1521 }
1522
1523}
1524
75bd7f81 1525//_____________________________________________________________________________
af885e0f 1526void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESDEvent *esd)
7ad19338 1527{
1528 //
1529 // Creates seeds using clusters between position inner plane and outer plane
1530 //
75bd7f81 1531
4551ea7c 1532 const Double_t kMaxTheta = 1.0;
1533 const Double_t kMaxPhi = 2.0;
1534
1535 const Double_t kRoad0y = 6.0; // Road for middle cluster
1536 const Double_t kRoad0z = 8.5; // Road for middle cluster
1537
1538 const Double_t kRoad1y = 2.0; // Road in y for seeded cluster
1539 const Double_t kRoad1z = 20.0; // Road in z for seeded cluster
1540
1541 const Double_t kRoad2y = 3.0; // Road in y for extrapolated cluster
1542 const Double_t kRoad2z = 20.0; // Road in z for extrapolated cluster
c6f438c0 1543 const Int_t kMaxSeed = 3000;
7ad19338 1544
4551ea7c 1545 Int_t maxSec = AliTRDgeometry::kNsect;
1546
1547 // Linear fitters in planes
1548 TLinearFitter fitterTC(2,"hyp2"); // Fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
1549 TLinearFitter fitterT2(4,"hyp4"); // Fitting with tilting pads - kz not fixed
69b55c55 1550 fitterTC.StoreData(kTRUE);
1551 fitterT2.StoreData(kTRUE);
4551ea7c 1552 AliRieman rieman(1000); // Rieman fitter
1553 AliRieman rieman2(1000); // Rieman fitter
1554
1555 // Find the maximal and minimal layer for the planes
7ad19338 1556 Int_t layers[6][2];
4551ea7c 1557 AliTRDpropagationLayer *reflayers[6];
1558 for (Int_t i = 0; i < 6; i++) {
1559 layers[i][0] = 10000;
1560 layers[i][1] = 0;
1561 }
1562 for (Int_t ns = 0; ns < maxSec; ns++) {
1563 for (Int_t ilayer = 0; ilayer < fTrSec[ns]->GetNumberOfLayers(); ilayer++) {
1564 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(ilayer));
1565 if (layer == 0) {
1566 continue;
1567 }
7ad19338 1568 Int_t det = layer[0]->GetDetector();
1569 Int_t plane = fGeom->GetPlane(det);
4551ea7c 1570 if (ilayer < layers[plane][0]) {
1571 layers[plane][0] = ilayer;
1572 }
1573 if (ilayer > layers[plane][1]) {
1574 layers[plane][1] = ilayer;
1575 }
7ad19338 1576 }
1577 }
4551ea7c 1578
f162af62 1579 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0);
6c94f330 1580 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
4551ea7c 1581 Double_t hL[6]; // Tilting angle
1582 Double_t xcl[6]; // X - position of reference cluster
1583 Double_t ycl[6]; // Y - position of reference cluster
1584 Double_t zcl[6]; // Z - position of reference cluster
1585
1586 AliTRDcluster *cl[6] = { 0, 0, 0, 0, 0, 0 }; // Seeding clusters
1587 Float_t padlength[6] = { 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 }; // Current pad-length
1588
1589 Double_t chi2R = 0.0;
1590 Double_t chi2Z = 0.0;
1591 Double_t chi2RF = 0.0;
1592 Double_t chi2ZF = 0.0;
1593
1594 Int_t nclusters; // Total number of clusters
1595 for (Int_t i = 0; i < 6; i++) {
1596 hL[i] = h01;
1597 if (i%2==1) {
1598 hL[i]*=-1.0;
1599 }
1600 }
1601
1602 // Registered seed
c6f438c0 1603 AliTRDseed *pseed = new AliTRDseed[kMaxSeed*6];
1604 AliTRDseed *seed[kMaxSeed];
4551ea7c 1605 for (Int_t iseed = 0; iseed < kMaxSeed; iseed++) {
1606 seed[iseed]= &pseed[iseed*6];
1607 }
69b55c55 1608 AliTRDseed *cseed = seed[0];
4551ea7c 1609
1610 Double_t seedquality[kMaxSeed];
1611 Double_t seedquality2[kMaxSeed];
1612 Double_t seedparams[kMaxSeed][7];
1613 Int_t seedlayer[kMaxSeed];
1614 Int_t registered = 0;
1615 Int_t sort[kMaxSeed];
1616
1617 //
1618 // Seeding part
1619 //
1620 for (Int_t ns = 0; ns < maxSec; ns++) { // Loop over sectors
1621 //for (Int_t ns = 0; ns < 5; ns++) { // Loop over sectors
1622
1623 registered = 0; // Reset registerd seed counter
1624 cseed = seed[registered];
1625 Float_t iter = 0.0;
1626
1627 for (Int_t sLayer = 2; sLayer >= 0; sLayer--) {
1628 //for (Int_t dseed = 5; dseed < 15; dseed += 3) {
1629
1630 iter += 1.0;
1631 Int_t dseed = 5 + Int_t(iter) * 3;
1632
69b55c55 1633 // Initialize seeding layers
4551ea7c 1634 for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
69b55c55 1635 reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
1636 xcl[ilayer] = reflayers[ilayer]->GetX();
4551ea7c 1637 }
1638
1639 Double_t xref = (xcl[sLayer+1] + xcl[sLayer+2]) * 0.5;
1640 AliTRDpropagationLayer &layer0 = *reflayers[sLayer+0];
1641 AliTRDpropagationLayer &layer1 = *reflayers[sLayer+1];
1642 AliTRDpropagationLayer &layer2 = *reflayers[sLayer+2];
1643 AliTRDpropagationLayer &layer3 = *reflayers[sLayer+3];
1644
1645 Int_t maxn3 = layer3;
1646 for (Int_t icl3 = 0; icl3 < maxn3; icl3++) {
33744e5d 1647
69b55c55 1648 AliTRDcluster *cl3 = layer3[icl3];
4551ea7c 1649 if (!cl3) {
1650 continue;
1651 }
1652 padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2() * 12.0);
69b55c55 1653 ycl[sLayer+3] = cl3->GetY();
1654 zcl[sLayer+3] = cl3->GetZ();
4551ea7c 1655 Float_t yymin0 = ycl[sLayer+3] - 1.0 - kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1656 Float_t yymax0 = ycl[sLayer+3] + 1.0 + kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1657 Int_t maxn0 = layer0;
1658
1659 for (Int_t icl0 = layer0.Find(yymin0); icl0 < maxn0; icl0++) {
1660
69b55c55 1661 AliTRDcluster *cl0 = layer0[icl0];
4551ea7c 1662 if (!cl0) {
1663 continue;
1664 }
1665 if (cl3->IsUsed() && cl0->IsUsed()) {
1666 continue;
1667 }
69b55c55 1668 ycl[sLayer+0] = cl0->GetY();
1669 zcl[sLayer+0] = cl0->GetZ();
4551ea7c 1670 if (ycl[sLayer+0] > yymax0) {
1671 break;
1672 }
1673 Double_t tanphi = (ycl[sLayer+3]-ycl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1674 if (TMath::Abs(tanphi) > kMaxPhi) {
1675 continue;
1676 }
1677 Double_t tantheta = (zcl[sLayer+3]-zcl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1678 if (TMath::Abs(tantheta) > kMaxTheta) {
1679 continue;
1680 }
1681 padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2() * 12.0);
1682
1683 // Expected position in 1 layer
1684 Double_t y1exp = ycl[sLayer+0] + (tanphi) * (xcl[sLayer+1]-xcl[sLayer+0]);
1685 Double_t z1exp = zcl[sLayer+0] + (tantheta) * (xcl[sLayer+1]-xcl[sLayer+0]);
1686 Float_t yymin1 = y1exp - kRoad0y - tanphi;
1687 Float_t yymax1 = y1exp + kRoad0y + tanphi;
1688 Int_t maxn1 = layer1;
1689
1690 for (Int_t icl1 = layer1.Find(yymin1); icl1 < maxn1; icl1++) {
1691
69b55c55 1692 AliTRDcluster *cl1 = layer1[icl1];
4551ea7c 1693 if (!cl1) {
1694 continue;
1695 }
69b55c55 1696 Int_t nusedCl = 0;
1697 if (cl3->IsUsed()) nusedCl++;
1698 if (cl0->IsUsed()) nusedCl++;
1699 if (cl1->IsUsed()) nusedCl++;
4551ea7c 1700 if (nusedCl > 1) {
1701 continue;
1702 }
69b55c55 1703 ycl[sLayer+1] = cl1->GetY();
1704 zcl[sLayer+1] = cl1->GetZ();
4551ea7c 1705 if (ycl[sLayer+1] > yymax1) {
1706 break;
1707 }
1708 if (TMath::Abs(ycl[sLayer+1]-y1exp) > kRoad0y+tanphi) {
1709 continue;
1710 }
1711 if (TMath::Abs(zcl[sLayer+1]-z1exp) > kRoad0z) {
1712 continue;
1713 }
1714 padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2() * 12.0);
1715
1716 Double_t y2exp = ycl[sLayer+0]+(tanphi) * (xcl[sLayer+2]-xcl[sLayer+0]) + (ycl[sLayer+1]-y1exp);
1717 Double_t z2exp = zcl[sLayer+0]+(tantheta) * (xcl[sLayer+2]-xcl[sLayer+0]);
1718 Int_t index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y,kRoad1z);
1719 if (index2 <= 0) {
1720 continue;
1721 }
1722 AliTRDcluster *cl2 = (AliTRDcluster *) GetCluster(index2);
1723 padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2() * 12.0);
1724 ycl[sLayer+2] = cl2->GetY();
1725 zcl[sLayer+2] = cl2->GetZ();
1726 if (TMath::Abs(cl2->GetZ()-z2exp) > kRoad0z) {
1727 continue;
1728 }
1729
69b55c55 1730 rieman.Reset();
1731 rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
1732 rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
1733 rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);
1734 rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
1735 rieman.Update();
4551ea7c 1736
1737 // Reset fitter
1738 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
69b55c55 1739 cseed[iLayer].Reset();
1740 }
4551ea7c 1741 chi2Z = 0.0;
1742 chi2R = 0.0;
1743
1744 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
ca1e1e5b 1745 cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1746 chi2Z += (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer])
1747 * (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer]);
1748 cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
1749 cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1750 chi2R += (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer])
1751 * (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer]);
1752 cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
69b55c55 1753 }
4551ea7c 1754 if (TMath::Sqrt(chi2R) > 1.0/iter) {
1755 continue;
1756 }
1757 if (TMath::Sqrt(chi2Z) > 7.0/iter) {
1758 continue;
69b55c55 1759 }
4551ea7c 1760
1761 Float_t minmax[2] = { -100.0, 100.0 };
1762 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1763 Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer] * 0.5
ca1e1e5b 1764 + 1.0 - cseed[sLayer+iLayer].GetZref(0);
4551ea7c 1765 if (max < minmax[1]) {
1766 minmax[1] = max;
1767 }
1768 Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer] * 0.5
ca1e1e5b 1769 - 1.0 - cseed[sLayer+iLayer].GetZref(0);
4551ea7c 1770 if (min > minmax[0]) {
1771 minmax[0] = min;
1772 }
1773 }
1774
69b55c55 1775 Bool_t isFake = kFALSE;
4551ea7c 1776 if (cl0->GetLabel(0) != cl3->GetLabel(0)) {
1777 isFake = kTRUE;
1778 }
1779 if (cl1->GetLabel(0) != cl3->GetLabel(0)) {
1780 isFake = kTRUE;
1781 }
1782 if (cl2->GetLabel(0) != cl3->GetLabel(0)) {
1783 isFake = kTRUE;
1784 }
1785
1786 if (AliTRDReconstructor::StreamLevel() > 0) {
1787 if ((!isFake) || ((icl3%10) == 0)) { // Debugging print
1788 TTreeSRedirector &cstream = *fDebugStreamer;
1789 cstream << "Seeds0"
1790 << "isFake=" << isFake
1791 << "Cl0.=" << cl0
1792 << "Cl1.=" << cl1
1793 << "Cl2.=" << cl2
1794 << "Cl3.=" << cl3
1795 << "Xref=" << xref
1796 << "X0=" << xcl[sLayer+0]
1797 << "X1=" << xcl[sLayer+1]
1798 << "X2=" << xcl[sLayer+2]
1799 << "X3=" << xcl[sLayer+3]
1800 << "Y2exp=" << y2exp
1801 << "Z2exp=" << z2exp
1802 << "Chi2R=" << chi2R
1803 << "Chi2Z=" << chi2Z
1804 << "Seed0.=" << &cseed[sLayer+0]
1805 << "Seed1.=" << &cseed[sLayer+1]
1806 << "Seed2.=" << &cseed[sLayer+2]
1807 << "Seed3.=" << &cseed[sLayer+3]
1808 << "Zmin=" << minmax[0]
1809 << "Zmax=" << minmax[1]
1810 << "\n";
d337ef8d 1811 }
69b55c55 1812 }
33744e5d 1813
4551ea7c 1814 ////////////////////////////////////////////////////////////////////////////////////
33744e5d 1815 //
4551ea7c 1816 // Fit seeding part
33744e5d 1817 //
4551ea7c 1818 ////////////////////////////////////////////////////////////////////////////////////
33744e5d 1819
69b55c55 1820 cl[sLayer+0] = cl0;
1821 cl[sLayer+1] = cl1;
1822 cl[sLayer+2] = cl2;
1823 cl[sLayer+3] = cl3;
4551ea7c 1824 Bool_t isOK = kTRUE;
1825
1826 for (Int_t jLayer = 0; jLayer < 4; jLayer++) {
1827
ca1e1e5b 1828 cseed[sLayer+jLayer].SetTilt(hL[sLayer+jLayer]);
1829 cseed[sLayer+jLayer].SetPadLength(padlength[sLayer+jLayer]);
1830 cseed[sLayer+jLayer].SetX0(xcl[sLayer+jLayer]);
4551ea7c 1831
1832 for (Int_t iter = 0; iter < 2; iter++) {
1833
69b55c55 1834 //
4551ea7c 1835 // In iteration 0 we try only one pad-row
1836 // If quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
69b55c55 1837 //
1838 AliTRDseed tseed = cseed[sLayer+jLayer];
4551ea7c 1839 Float_t roadz = padlength[sLayer+jLayer] * 0.5;
1840 if (iter > 0) {
1841 roadz = padlength[sLayer+jLayer];
1842 }
1843
1844 Float_t quality = 10000.0;
1845
1846 for (Int_t iTime = 2; iTime < 20; iTime++) {
1847
1848 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
1849 Double_t dxlayer = layer.GetX() - xcl[sLayer+jLayer];
1850 Double_t zexp = cl[sLayer+jLayer]->GetZ();
1851
1852 if (iter > 0) {
1853 // Try 2 pad-rows in second iteration
ca1e1e5b 1854 zexp = tseed.GetZref(0) + tseed.GetZref(1) * dxlayer;
4551ea7c 1855 if (zexp > cl[sLayer+jLayer]->GetZ()) {
1856 zexp = cl[sLayer+jLayer]->GetZ() + padlength[sLayer+jLayer]*0.5;
1857 }
1858 if (zexp < cl[sLayer+jLayer]->GetZ()) {
1859 zexp = cl[sLayer+jLayer]->GetZ() - padlength[sLayer+jLayer]*0.5;
1860 }
69b55c55 1861 }
4551ea7c 1862
ca1e1e5b 1863 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
4551ea7c 1864 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
1865 if (index <= 0) {
1866 continue;
1867 }
1868 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
1869
ca1e1e5b 1870 tseed.SetIndexes(iTime,index);
1871 tseed.SetClusters(iTime,cl); // Register cluster
1872 tseed.SetX(iTime,dxlayer); // Register cluster
1873 tseed.SetY(iTime,cl->GetY()); // Register cluster
1874 tseed.SetZ(iTime,cl->GetZ()); // Register cluster
4551ea7c 1875
1876 }
1877
69b55c55 1878 tseed.Update();
4551ea7c 1879
1880 // Count the number of clusters and distortions into quality
ca1e1e5b 1881 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1882 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
1883 + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
1884 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 1885 if ((iter == 0) && tseed.IsOK()) {
1886 cseed[sLayer+jLayer] = tseed;
1887 quality = tquality;
1888 if (tquality < 5) {
1889 break;
1890 }
1891 }
1892 if (tseed.IsOK() && (tquality < quality)) {
69b55c55 1893 cseed[sLayer+jLayer] = tseed;
69b55c55 1894 }
4551ea7c 1895
1896 } // Loop: iter
1897
1898 if (!cseed[sLayer+jLayer].IsOK()) {
69b55c55 1899 isOK = kFALSE;
1900 break;
4551ea7c 1901 }
1902
69b55c55 1903 cseed[sLayer+jLayer].CookLabels();
1904 cseed[sLayer+jLayer].UpdateUsed();
ca1e1e5b 1905 nusedCl += cseed[sLayer+jLayer].GetNUsed();
4551ea7c 1906 if (nusedCl > 25) {
69b55c55 1907 isOK = kFALSE;
1908 break;
4551ea7c 1909 }
1910
1911 } // Loop: jLayer
1912
1913 if (!isOK) {
1914 continue;
69b55c55 1915 }
4551ea7c 1916 nclusters = 0;
1917 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1918 if (cseed[sLayer+iLayer].IsOK()) {
ca1e1e5b 1919 nclusters += cseed[sLayer+iLayer].GetN2();
69b55c55 1920 }
1921 }
4551ea7c 1922
1923 // Iteration 0
69b55c55 1924 rieman.Reset();
4551ea7c 1925 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1926 rieman.AddPoint(xcl[sLayer+iLayer]
ca1e1e5b 1927 ,cseed[sLayer+iLayer].GetYfitR(0)
1928 ,cseed[sLayer+iLayer].GetZProb()
4551ea7c 1929 ,1
1930 ,10);
69b55c55 1931 }
1932 rieman.Update();
4551ea7c 1933
1934 chi2R = 0.0;
1935 chi2Z = 0.0;
1936
1937 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
ca1e1e5b 1938 cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1939 chi2R += (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0))
1940 * (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0));
1941 cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
1942 cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1943 chi2Z += (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz())
1944 * (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz());
1945 cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
69b55c55 1946 }
1947 Double_t curv = rieman.GetC();
4551ea7c 1948
69b55c55 1949 //
4551ea7c 1950 // Likelihoods
6c94f330 1951 //
ca1e1e5b 1952 Double_t sumda = TMath::Abs(cseed[sLayer+0].GetYfitR(1) - cseed[sLayer+0].GetYref(1))
1953 + TMath::Abs(cseed[sLayer+1].GetYfitR(1) - cseed[sLayer+1].GetYref(1))
1954 + TMath::Abs(cseed[sLayer+2].GetYfitR(1) - cseed[sLayer+2].GetYref(1))
1955 + TMath::Abs(cseed[sLayer+3].GetYfitR(1) - cseed[sLayer+3].GetYref(1));
4551ea7c 1956 Double_t likea = TMath::Exp(-sumda*10.6);
1957 Double_t likechi2 = 0.0000000001;
1958 if (chi2R < 0.5) {
1959 likechi2 += TMath::Exp(-TMath::Sqrt(chi2R) * 7.73);
1960 }
1961 Double_t likechi2z = TMath::Exp(-chi2Z * 0.088) / TMath::Exp(-chi2Z * 0.019);
1962 Double_t likeN = TMath::Exp(-(72 - nclusters) * 0.19);
1963 Double_t like = likea * likechi2 * likechi2z * likeN;
ca1e1e5b 1964 Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetYref(1) - 130.0*curv) * 1.9);
1965 Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetZref(1)
1966 - cseed[sLayer+0].GetZref(0) / xcl[sLayer+0]) * 5.9);
6c94f330 1967 Double_t likePrim = TMath::Max(likePrimY*likePrimZ,0.0005);
69b55c55 1968
4551ea7c 1969 seedquality[registered] = like;
1970 seedlayer[registered] = sLayer;
1971 if (TMath::Log(0.000000000000001 + like) < -15) {
1972 continue;
1973 }
69b55c55 1974 AliTRDseed seedb[6];
4551ea7c 1975 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
69b55c55 1976 seedb[iLayer] = cseed[iLayer];
1977 }
4551ea7c 1978
1979 ////////////////////////////////////////////////////////////////////////////////////
ad670fba 1980 //
4551ea7c 1981 // Full track fit part
33744e5d 1982 //
4551ea7c 1983 ////////////////////////////////////////////////////////////////////////////////////
1984
1985 Int_t nlayers = 0;
1986 Int_t nusedf = 0;
1987 Int_t findable = 0;
1988
69b55c55 1989 //
4551ea7c 1990 // Add new layers - avoid long extrapolation
69b55c55 1991 //
4551ea7c 1992 Int_t tLayer[2] = { 0, 0 };
1993 if (sLayer == 2) {
1994 tLayer[0] = 1;
1995 tLayer[1] = 0;
1996 }
1997 if (sLayer == 1) {
1998 tLayer[0] = 5;
1999 tLayer[1] = 0;
2000 }
2001 if (sLayer == 0) {
2002 tLayer[0] = 4;
2003 tLayer[1] = 5;
2004 }
2005
2006 for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
2007 Int_t jLayer = tLayer[iLayer]; // Set tracking layer
69b55c55 2008 cseed[jLayer].Reset();
ca1e1e5b 2009 cseed[jLayer].SetTilt(hL[jLayer]);
2010 cseed[jLayer].SetPadLength(padlength[jLayer]);
2011 cseed[jLayer].SetX0(xcl[jLayer]);
4551ea7c 2012 // Get pad length and rough cluster
ca1e1e5b 2013 Int_t indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].GetYref(0)
2014 ,cseed[jLayer].GetZref(0)
4551ea7c 2015 ,kRoad2y
2016 ,kRoad2z);
2017 if (indexdummy <= 0) {
2018 continue;
2019 }
2020 AliTRDcluster *cldummy = (AliTRDcluster *) GetCluster(indexdummy);
2021 padlength[jLayer] = TMath::Sqrt(cldummy->GetSigmaZ2() * 12.0);
69b55c55 2022 }
4551ea7c 2023 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
2024
2025 for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
2026
2027 Int_t jLayer = tLayer[iLayer]; // set tracking layer
2028 if ((jLayer == 0) && !(cseed[1].IsOK())) {
2029 continue; // break not allowed
2030 }
2031 if ((jLayer == 5) && !(cseed[4].IsOK())) {
2032 continue; // break not allowed
2033 }
ca1e1e5b 2034 Float_t zexp = cseed[jLayer].GetZref(0);
4551ea7c 2035 Double_t zroad = padlength[jLayer] * 0.5 + 1.0;
2036
2037 for (Int_t iter = 0; iter < 2; iter++) {
2038
69b55c55 2039 AliTRDseed tseed = cseed[jLayer];
4551ea7c 2040 Float_t quality = 10000.0;
2041
2042 for (Int_t iTime = 2; iTime < 20; iTime++) {
2043 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
2044 Double_t dxlayer = layer.GetX()-xcl[jLayer];
ca1e1e5b 2045 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
4551ea7c 2046 Float_t yroad = kRoad1y;
2047 Int_t index = layer.FindNearestCluster(yexp,zexp,yroad,zroad);
2048 if (index <= 0) {
2049 continue;
2050 }
2051 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
ca1e1e5b 2052 tseed.SetIndexes(iTime,index);
2053 tseed.SetClusters(iTime,cl); // Register cluster
2054 tseed.SetX(iTime,dxlayer); // Register cluster
2055 tseed.SetY(iTime,cl->GetY()); // Register cluster
2056 tseed.SetZ(iTime,cl->GetZ()); // Register cluster
4551ea7c 2057 }
2058
69b55c55 2059 tseed.Update();
4551ea7c 2060 if (tseed.IsOK()) {
ca1e1e5b 2061 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
2062 Float_t tquality = (18.0 - tseed.GetN2())/2.0 + TMath::Abs(dangle) / 0.1
2063 + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
2064 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 2065 if (tquality < quality) {
2066 cseed[jLayer] = tseed;
2067 quality = tquality;
69b55c55 2068 }
2069 }
4551ea7c 2070
2071 zroad *= 2.0;
2072
2073 } // Loop: iter
2074
2075 if ( cseed[jLayer].IsOK()) {
69b55c55 2076 cseed[jLayer].CookLabels();
2077 cseed[jLayer].UpdateUsed();
ca1e1e5b 2078 nusedf += cseed[jLayer].GetNUsed();
4551ea7c 2079 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
69b55c55 2080 }
4551ea7c 2081
2082 } // Loop: iLayer
2083
2084 // Make copy
69b55c55 2085 AliTRDseed bseed[6];
4551ea7c 2086 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
69b55c55 2087 bseed[jLayer] = cseed[jLayer];
4551ea7c 2088 }
2089 Float_t lastquality = 10000.0;
2090 Float_t lastchi2 = 10000.0;
2091 Float_t chi2 = 1000.0;
2092
2093 for (Int_t iter = 0; iter < 4; iter++) {
2094
2095 // Sort tracklets according "quality", try to "improve" 4 worst
2096 Float_t sumquality = 0.0;
69b55c55 2097 Float_t squality[6];
2098 Int_t sortindexes[6];
4551ea7c 2099
2100 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2101 if (bseed[jLayer].IsOK()) {
69b55c55 2102 AliTRDseed &tseed = bseed[jLayer];
ca1e1e5b 2103 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
2104 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
2105 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
2106 + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
2107 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 2108 squality[jLayer] = tquality;
2109 }
2110 else {
2111 squality[jLayer] = -1.0;
ad670fba 2112 }
6c94f330 2113 sumquality +=squality[jLayer];
69b55c55 2114 }
2115
4551ea7c 2116 if ((sumquality >= lastquality) ||
2117 (chi2 > lastchi2)) {
2118 break;
2119 }
69b55c55 2120 lastquality = sumquality;
2121 lastchi2 = chi2;
4551ea7c 2122 if (iter > 0) {
2123 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
69b55c55 2124 cseed[jLayer] = bseed[jLayer];
2125 }
2126 }
2127 TMath::Sort(6,squality,sortindexes,kFALSE);
4551ea7c 2128
2129 for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
2130
6c94f330 2131 Int_t bLayer = sortindexes[jLayer];
2132 AliTRDseed tseed = bseed[bLayer];
4551ea7c 2133
2134 for (Int_t iTime = 2; iTime < 20; iTime++) {
2135
2136 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
2137 Double_t dxlayer = layer.GetX() - xcl[bLayer];
ca1e1e5b 2138 Double_t zexp = tseed.GetZref(0);
2139 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
4551ea7c 2140 Float_t roadz = padlength[bLayer] + 1;
ca1e1e5b 2141 if (TMath::Abs(tseed.GetZProb() - zexp) > 0.5*padlength[bLayer]) {
4551ea7c 2142 roadz = padlength[bLayer] * 0.5;
2143 }
ca1e1e5b 2144 if (tseed.GetZfit(1)*tseed.GetZref(1) < 0.0) {
4551ea7c 2145 roadz = padlength[bLayer] * 0.5;
2146 }
ca1e1e5b 2147 if (TMath::Abs(tseed.GetZProb() - zexp) < 0.1*padlength[bLayer]) {
2148 zexp = tseed.GetZProb();
4551ea7c 2149 roadz = padlength[bLayer] * 0.5;
2150 }
2151
ca1e1e5b 2152 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer - zcor;
4551ea7c 2153 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
2154 if (index <= 0) {
2155 continue;
69b55c55 2156 }
4551ea7c 2157 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
2158
ca1e1e5b 2159 tseed.SetIndexes(iTime,index);
2160 tseed.SetClusters(iTime,cl); // Register cluster
2161 tseed.SetX(iTime,dxlayer); // Register cluster
2162 tseed.SetY(iTime,cl->GetY()); // Register cluster
2163 tseed.SetZ(iTime,cl->GetZ()); // Register cluster
4551ea7c 2164
2165 }
2166
69b55c55 2167 tseed.Update();
c6f438c0 2168 if (tseed.IsOK()) {
ca1e1e5b 2169 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
2170 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
2171 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0
2172 + TMath::Abs(dangle) / 0.1
2173 + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
2174 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 2175 if (tquality<squality[bLayer]) {
69b55c55 2176 bseed[bLayer] = tseed;
4551ea7c 2177 }
69b55c55 2178 }
4551ea7c 2179
2180 } // Loop: jLayer
2181
2182 chi2 = AliTRDseed::FitRiemanTilt(bseed,kTRUE);
2183
2184 } // Loop: iter
2185
2186 nclusters = 0;
2187 nlayers = 0;
2188 findable = 0;
2189 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ca1e1e5b 2190 if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) {
69b55c55 2191 findable++;
4551ea7c 2192 }
2193 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2194 nclusters += cseed[iLayer].GetN2();
69b55c55 2195 nlayers++;
2196 }
2197 }
4551ea7c 2198 if (nlayers < 3) {
2199 continue;
2200 }
69b55c55 2201 rieman.Reset();
4551ea7c 2202 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2203 if (cseed[iLayer].IsOK()) {
2204 rieman.AddPoint(xcl[iLayer]
ca1e1e5b 2205 ,cseed[iLayer].GetYfitR(0)
2206 ,cseed[iLayer].GetZProb()
4551ea7c 2207 ,1
2208 ,10);
2209 }
69b55c55 2210 }
2211 rieman.Update();
4551ea7c 2212
2213 chi2RF = 0.0;
2214 chi2ZF = 0.0;
2215 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2216 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2217 cseed[iLayer].SetYref(0,rieman.GetYat(xcl[iLayer]));
2218 chi2RF += (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0))
2219 * (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0));
2220 cseed[iLayer].SetYref(1,rieman.GetDYat(xcl[iLayer]));
2221 cseed[iLayer].SetZref(0,rieman.GetZat(xcl[iLayer]));
2222 chi2ZF += (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz())
2223 * (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz());
2224 cseed[iLayer].SetZref(1,rieman.GetDZat(xcl[iLayer]));
69b55c55 2225 }
2226 }
4551ea7c 2227 chi2RF /= TMath::Max((nlayers - 3.0),1.0);
2228 chi2ZF /= TMath::Max((nlayers - 3.0),1.0);
2229 curv = rieman.GetC();
2230
2231 Double_t xref2 = (xcl[2] + xcl[3]) * 0.5; // Middle of the chamber
2232 Double_t dzmf = rieman.GetDZat(xref2);
2233 Double_t zmf = rieman.GetZat(xref2);
69b55c55 2234
69b55c55 2235 //
4551ea7c 2236 // Fit hyperplane
69b55c55 2237 //
4551ea7c 2238 Int_t npointsT = 0;
69b55c55 2239 fitterTC.ClearPoints();
2240 fitterT2.ClearPoints();
2241 rieman2.Reset();
4551ea7c 2242
2243 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2244
2245 if (!cseed[iLayer].IsOK()) {
2246 continue;
2247 }
2248
2249 for (Int_t itime = 0; itime < 25; itime++) {
2250
ca1e1e5b 2251 if (!cseed[iLayer].IsUsable(itime)) {
4551ea7c 2252 continue;
2253 }
2254 // X relative to the middle chamber
ca1e1e5b 2255 Double_t x = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;
2256 Double_t y = cseed[iLayer].GetY(itime);
2257 Double_t z = cseed[iLayer].GetZ(itime);
69b55c55 2258 // ExB correction to the correction
4551ea7c 2259 // Tilted rieman
69b55c55 2260 Double_t uvt[6];
4551ea7c 2261 // Global x
ca1e1e5b 2262 Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();
4551ea7c 2263 Double_t t = 1.0 / (x2*x2 + y*y);
2264 uvt[1] = t; // t
2265 uvt[0] = 2.0 * x2 * uvt[1]; // u
2266 uvt[2] = 2.0 * hL[iLayer] * uvt[1];
2267 uvt[3] = 2.0 * hL[iLayer] * x * uvt[1];
2268 uvt[4] = 2.0 * (y + hL[iLayer]*z) * uvt[1];
2269 Double_t error = 2.0 * 0.2 * uvt[1];
69b55c55 2270 fitterT2.AddPoint(uvt,uvt[4],error);
4551ea7c 2271
69b55c55 2272 //
4551ea7c 2273 // Constrained rieman
69b55c55 2274 //
ca1e1e5b 2275 z = cseed[iLayer].GetZ(itime);
4551ea7c 2276 uvt[0] = 2.0 * x2 * t; // u
2277 uvt[1] = 2.0 * hL[iLayer] * x2 * uvt[1];
2278 uvt[2] = 2.0 * (y + hL[iLayer] * (z - GetZ())) * t;
69b55c55 2279 fitterTC.AddPoint(uvt,uvt[2],error);
69b55c55 2280 rieman2.AddPoint(x2,y,z,1,10);
2281 npointsT++;
4551ea7c 2282
69b55c55 2283 }
4551ea7c 2284
2285 } // Loop: iLayer
2286
69b55c55 2287 rieman2.Update();
2288 fitterTC.Eval();
2289 fitterT2.Eval();
2290 Double_t rpolz0 = fitterT2.GetParameter(3);
2291 Double_t rpolz1 = fitterT2.GetParameter(4);
4551ea7c 2292
69b55c55 2293 //
4551ea7c 2294 // Linear fitter - not possible to make boundaries
2295 // Do not accept non possible z and dzdx combinations
69b55c55 2296 //
4551ea7c 2297 Bool_t acceptablez = kTRUE;
2298 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2299 if (cseed[iLayer].IsOK()) {
2300 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
ca1e1e5b 2301 if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > padlength[iLayer] * 0.5 + 1.0) {
69b55c55 2302 acceptablez = kFALSE;
4551ea7c 2303 }
69b55c55 2304 }
2305 }
4551ea7c 2306 if (!acceptablez) {
69b55c55 2307 fitterT2.FixParameter(3,zmf);
2308 fitterT2.FixParameter(4,dzmf);
2309 fitterT2.Eval();
2310 fitterT2.ReleaseParameter(3);
2311 fitterT2.ReleaseParameter(4);
2312 rpolz0 = fitterT2.GetParameter(3);
2313 rpolz1 = fitterT2.GetParameter(4);
2314 }
4551ea7c 2315
2316 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
2317 Double_t chi2TC = fitterTC.GetChisquare() / Float_t(npointsT);
69b55c55 2318 Double_t polz1c = fitterTC.GetParameter(2);
4551ea7c 2319 Double_t polz0c = polz1c * xref2;
6c94f330 2320 Double_t aC = fitterTC.GetParameter(0);
2321 Double_t bC = fitterTC.GetParameter(1);
4551ea7c 2322 Double_t cC = aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
6c94f330 2323 Double_t aR = fitterT2.GetParameter(0);
2324 Double_t bR = fitterT2.GetParameter(1);
2325 Double_t dR = fitterT2.GetParameter(2);
4551ea7c 2326 Double_t cR = 1.0 + bR*bR - dR*aR;
2327 Double_t dca = 0.0;
2328 if (cR > 0.0) {
2329 dca = -dR / (TMath::Sqrt(1.0 + bR*bR - dR*aR) + TMath::Sqrt(1.0 + bR*bR));
2330 cR = aR / TMath::Sqrt(cR);
69b55c55 2331 }
4551ea7c 2332
2333 Double_t chi2ZT2 = 0.0;
2334 Double_t chi2ZTC = 0.0;
2335 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2336 if (cseed[iLayer].IsOK()) {
2337 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
2338 Double_t zTC = polz0c + polz1c * (xcl[iLayer] - xref2);
ca1e1e5b 2339 chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
2340 chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
69b55c55 2341 }
2342 }
4551ea7c 2343 chi2ZT2 /= TMath::Max((nlayers - 3.0),1.0);
2344 chi2ZTC /= TMath::Max((nlayers - 3.0),1.0);
2345
2346 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
2347 Float_t sumdaf = 0.0;
2348 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2349 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2350 sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))
2351 / cseed[iLayer].GetSigmaY2());
4551ea7c 2352 }
2353 }
2354 sumdaf /= Float_t (nlayers - 2.0);
2355
69b55c55 2356 //
4551ea7c 2357 // Likelihoods for full track
69b55c55 2358 //
4551ea7c 2359 Double_t likezf = TMath::Exp(-chi2ZF * 0.14);
2360 Double_t likechi2C = TMath::Exp(-chi2TC * 0.677);
2361 Double_t likechi2TR = TMath::Exp(-chi2TR * 0.78);
2362 Double_t likeaf = TMath::Exp(-sumdaf * 3.23);
2363 seedquality2[registered] = likezf * likechi2TR * likeaf;
33744e5d 2364
4551ea7c 2365 // Still needed ????
69b55c55 2366// Bool_t isGold = kFALSE;
2367//
6c94f330 2368// if (nlayers == 6 && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE; // gold
2369// if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE; // gold
69b55c55 2370// if (isGold &&nusedf<10){
2371// for (Int_t jLayer=0;jLayer<6;jLayer++){
6c94f330 2372// if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
69b55c55 2373// seed[index][jLayer].UseClusters(); //sign gold
2374// }
2375// }
4551ea7c 2376
2377 Int_t index0 = 0;
2378 if (!cseed[0].IsOK()) {
69b55c55 2379 index0 = 1;
4551ea7c 2380 if (!cseed[1].IsOK()) {
2381 index0 = 2;
2382 }
69b55c55 2383 }
ca1e1e5b 2384 seedparams[registered][0] = cseed[index0].GetX0();
2385 seedparams[registered][1] = cseed[index0].GetYref(0);
2386 seedparams[registered][2] = cseed[index0].GetZref(0);
c6f438c0 2387 seedparams[registered][5] = cR;
ca1e1e5b 2388 seedparams[registered][3] = cseed[index0].GetX0() * cR - TMath::Sin(TMath::ATan(cseed[0].GetYref(1)));
2389 seedparams[registered][4] = cseed[index0].GetZref(1)
2390 / TMath::Sqrt(1.0 + cseed[index0].GetYref(1) * cseed[index0].GetYref(1));
69b55c55 2391 seedparams[registered][6] = ns;
4551ea7c 2392
2393 Int_t labels[12];
2394 Int_t outlab[24];
2395 Int_t nlab = 0;
2396 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2397 if (!cseed[iLayer].IsOK()) {
2398 continue;
2399 }
ca1e1e5b 2400 if (cseed[iLayer].GetLabels(0) >= 0) {
2401 labels[nlab] = cseed[iLayer].GetLabels(0);
69b55c55 2402 nlab++;
2403 }
ca1e1e5b 2404 if (cseed[iLayer].GetLabels(1) >= 0) {
2405 labels[nlab] = cseed[iLayer].GetLabels(1);
69b55c55 2406 nlab++;
4551ea7c 2407 }
69b55c55 2408 }
2409 Freq(nlab,labels,outlab,kFALSE);
4551ea7c 2410 Int_t label = outlab[0];
2411 Int_t frequency = outlab[1];
2412 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ca1e1e5b 2413 cseed[iLayer].SetFreq(frequency);
2414 cseed[iLayer].SetC(cR);
2415 cseed[iLayer].SetCC(cC);
2416 cseed[iLayer].SetChi2(chi2TR);
2417 cseed[iLayer].SetChi2Z(chi2ZF);
69b55c55 2418 }
33744e5d 2419
2420 // Debugging print
4551ea7c 2421 if (1 || (!isFake)) {
69b55c55 2422 Float_t zvertex = GetZ();
4551ea7c 2423 TTreeSRedirector &cstream = *fDebugStreamer;
2424 if (AliTRDReconstructor::StreamLevel() > 0) {
2425 cstream << "Seeds1"
2426 << "isFake=" << isFake
2427 << "Vertex=" << zvertex
2428 << "Rieman2.=" << &rieman2
2429 << "Rieman.=" << &rieman
2430 << "Xref=" << xref
2431 << "X0=" << xcl[0]
2432 << "X1=" << xcl[1]
2433 << "X2=" << xcl[2]
2434 << "X3=" << xcl[3]
2435 << "X4=" << xcl[4]
2436 << "X5=" << xcl[5]
2437 << "Chi2R=" << chi2R
2438 << "Chi2Z=" << chi2Z
2439 << "Chi2RF=" << chi2RF // Chi2 of trackletes on full track
2440 << "Chi2ZF=" << chi2ZF // Chi2 z on tracklets on full track
2441 << "Chi2ZT2=" << chi2ZT2 // Chi2 z on tracklets on full track - rieman tilt
2442 << "Chi2ZTC=" << chi2ZTC // Chi2 z on tracklets on full track - rieman tilt const
2443 << "Chi2TR=" << chi2TR // Chi2 without vertex constrain
2444 << "Chi2TC=" << chi2TC // Chi2 with vertex constrain
2445 << "C=" << curv // Non constrained - no tilt correction
2446 << "DR=" << dR // DR parameter - tilt correction
2447 << "DCA=" << dca // DCA - tilt correction
2448 << "CR=" << cR // Non constrained curvature - tilt correction
2449 << "CC=" << cC // Constrained curvature
2450 << "Polz0=" << polz0c
2451 << "Polz1=" << polz1c
2452 << "RPolz0=" << rpolz0
2453 << "RPolz1=" << rpolz1
2454 << "Ncl=" << nclusters
2455 << "Nlayers=" << nlayers
2456 << "NUsedS=" << nusedCl
2457 << "NUsed=" << nusedf
2458 << "Findable=" << findable
2459 << "Like=" << like
2460 << "LikePrim=" << likePrim
2461 << "Likechi2C=" << likechi2C
2462 << "Likechi2TR=" << likechi2TR
2463 << "Likezf=" << likezf
2464 << "LikeF=" << seedquality2[registered]
2465 << "S0.=" << &cseed[0]
2466 << "S1.=" << &cseed[1]
2467 << "S2.=" << &cseed[2]
2468 << "S3.=" << &cseed[3]
2469 << "S4.=" << &cseed[4]
2470 << "S5.=" << &cseed[5]
2471 << "SB0.=" << &seedb[0]
2472 << "SB1.=" << &seedb[1]
2473 << "SB2.=" << &seedb[2]
2474 << "SB3.=" << &seedb[3]
2475 << "SB4.=" << &seedb[4]
2476 << "SB5.=" << &seedb[5]
2477 << "Label=" << label
2478 << "Freq=" << frequency
2479 << "sLayer=" << sLayer
2480 << "\n";
2481 }
69b55c55 2482 }
4551ea7c 2483
2484 if (registered<kMaxSeed - 1) {
69b55c55 2485 registered++;
2486 cseed = seed[registered];
2487 }
4551ea7c 2488
2489 } // End of loop over layer 1
2490
2491 } // End of loop over layer 0
2492
2493 } // End of loop over layer 3
2494
2495 } // End of loop over seeding time bins
2496
69b55c55 2497 //
4551ea7c 2498 // Choose best
69b55c55 2499 //
4551ea7c 2500
69b55c55 2501 TMath::Sort(registered,seedquality2,sort,kTRUE);
c6f438c0 2502 Bool_t signedseed[kMaxSeed];
4551ea7c 2503 for (Int_t i = 0; i < registered; i++) {
2504 signedseed[i] = kFALSE;
69b55c55 2505 }
4551ea7c 2506
2507 for (Int_t iter = 0; iter < 5; iter++) {
2508
2509 for (Int_t iseed = 0; iseed < registered; iseed++) {
2510
69b55c55 2511 Int_t index = sort[iseed];
4551ea7c 2512 if (signedseed[index]) {
2513 continue;
2514 }
69b55c55 2515 Int_t labelsall[1000];
4551ea7c 2516 Int_t nlabelsall = 0;
2517 Int_t naccepted = 0;;
2518 Int_t sLayer = seedlayer[index];
2519 Int_t ncl = 0;
2520 Int_t nused = 0;
2521 Int_t nlayers = 0;
69b55c55 2522 Int_t findable = 0;
4551ea7c 2523
2524 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2525
ca1e1e5b 2526 if (TMath::Abs(seed[index][jLayer].GetYref(0) / xcl[jLayer]) < 0.15) {
69b55c55 2527 findable++;
4551ea7c 2528 }
2529 if (seed[index][jLayer].IsOK()) {
69b55c55 2530 seed[index][jLayer].UpdateUsed();
ca1e1e5b 2531 ncl +=seed[index][jLayer].GetN2();
2532 nused +=seed[index][jLayer].GetNUsed();
69b55c55 2533 nlayers++;
4551ea7c 2534 // Cooking label
2535 for (Int_t itime = 0; itime < 25; itime++) {
ca1e1e5b 2536 if (seed[index][jLayer].IsUsable(itime)) {
69b55c55 2537 naccepted++;
4551ea7c 2538 for (Int_t ilab = 0; ilab < 3; ilab++) {
ca1e1e5b 2539 Int_t tindex = seed[index][jLayer].GetClusters(itime)->GetLabel(ilab);
4551ea7c 2540 if (tindex >= 0) {
69b55c55 2541 labelsall[nlabelsall] = tindex;
2542 nlabelsall++;
2543 }
2544 }
2545 }
2546 }
2547 }
4551ea7c 2548
69b55c55 2549 }
4551ea7c 2550
2551 if (nused > 30) {
2552 continue;
69b55c55 2553 }
4551ea7c 2554
2555 if (iter == 0) {
2556 if (nlayers < 6) {
2557 continue;
2558 }
2559 if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2560 continue; // Gold
2561 }
7ad19338 2562 }
4551ea7c 2563
2564 if (iter == 1) {
2565 if (nlayers < findable) {
2566 continue;
2567 }
2568 if (TMath::Log(0.000000001+seedquality2[index]) < -4.0) {
2569 continue;
2570 }
2571 }
2572
2573 if (iter == 2) {
2574 if ((nlayers == findable) ||
2575 (nlayers == 6)) {
2576 continue;
2577 }
2578 if (TMath::Log(0.000000001+seedquality2[index]) < -6.0) {
2579 continue;
2580 }
69b55c55 2581 }
4551ea7c 2582
2583 if (iter == 3) {
2584 if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2585 continue;
2586 }
69b55c55 2587 }
4551ea7c 2588
2589 if (iter == 4) {
2590 if (TMath::Log(0.000000001+seedquality2[index]) - nused/(nlayers-3.0) < -15.0) {
2591 continue;
2592 }
69b55c55 2593 }
4551ea7c 2594
69b55c55 2595 signedseed[index] = kTRUE;
4551ea7c 2596
2597 Int_t labels[1000];
2598 Int_t outlab[1000];
2599 Int_t nlab = 0;
2600 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2601 if (seed[index][iLayer].IsOK()) {
ca1e1e5b 2602 if (seed[index][iLayer].GetLabels(0) >= 0) {
2603 labels[nlab] = seed[index][iLayer].GetLabels(0);
69b55c55 2604 nlab++;
2605 }
ca1e1e5b 2606 if (seed[index][iLayer].GetLabels(1) >= 0) {
2607 labels[nlab] = seed[index][iLayer].GetLabels(1);
69b55c55 2608 nlab++;
4551ea7c 2609 }
2610 }
7ad19338 2611 }
69b55c55 2612 Freq(nlab,labels,outlab,kFALSE);
4551ea7c 2613 Int_t label = outlab[0];
2614 Int_t frequency = outlab[1];
69b55c55 2615 Freq(nlabelsall,labelsall,outlab,kFALSE);
4551ea7c 2616 Int_t label1 = outlab[0];
2617 Int_t label2 = outlab[2];
2618 Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
2619 Float_t ratio = Float_t(nused) / Float_t(ncl);
2620 if (ratio < 0.25) {
2621 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2622 if ((seed[index][jLayer].IsOK()) &&
ca1e1e5b 2623 (TMath::Abs(seed[index][jLayer].GetYfit(1) - seed[index][jLayer].GetYfit(1)) < 0.2)) {
4551ea7c 2624 seed[index][jLayer].UseClusters(); // Sign gold
2625 }
69b55c55 2626 }
7ad19338 2627 }
4551ea7c 2628
31fd97b2 2629 Int_t eventNrInFile = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
4551ea7c 2630 TTreeSRedirector &cstream = *fDebugStreamer;
2631
69b55c55 2632 //
4551ea7c 2633 // Register seed
69b55c55 2634 //
4551ea7c 2635 AliTRDtrack *track = RegisterSeed(seed[index],seedparams[index]);
2636 AliTRDtrack dummy;
2637 if (!track) {
2638 track = &dummy;
2639 }
2640 else {
69b55c55 2641 AliESDtrack esdtrack;
4551ea7c 2642 esdtrack.UpdateTrackParams(track,AliESDtrack::kTRDout);
69b55c55 2643 esdtrack.SetLabel(label);
2644 esd->AddTrack(&esdtrack);
4551ea7c 2645 TTreeSRedirector &cstream = *fDebugStreamer;
2646 if (AliTRDReconstructor::StreamLevel() > 0) {
2647 cstream << "Tracks"
31fd97b2 2648 << "EventNrInFile=" << eventNrInFile
4551ea7c 2649 << "ESD.=" << &esdtrack
2650 << "trd.=" << track
2651 << "trdback.=" << track
2652 << "\n";
2653 }
7ad19338 2654 }
4551ea7c 2655
2656 if (AliTRDReconstructor::StreamLevel() > 0) {
2657 cstream << "Seeds2"
2658 << "Iter=" << iter
2659 << "Track.=" << track
2660 << "Like=" << seedquality[index]
2661 << "LikeF=" << seedquality2[index]
2662 << "S0.=" << &seed[index][0]
2663 << "S1.=" << &seed[index][1]
2664 << "S2.=" << &seed[index][2]
2665 << "S3.=" << &seed[index][3]
2666 << "S4.=" << &seed[index][4]
2667 << "S5.=" << &seed[index][5]
2668 << "Label=" << label
2669 << "Label1=" << label1
2670 << "Label2=" << label2
2671 << "FakeRatio=" << fakeratio
2672 << "Freq=" << frequency
2673 << "Ncl=" << ncl
2674 << "Nlayers=" << nlayers
2675 << "Findable=" << findable
2676 << "NUsed=" << nused
2677 << "sLayer=" << sLayer
31fd97b2 2678 << "EventNrInFile=" << eventNrInFile
4551ea7c 2679 << "\n";
2680 }
2681
2682 } // Loop: iseed
2683
2684 } // Loop: iter
2685
2686 } // End of loop over sectors
75bd7f81 2687
69b55c55 2688 delete [] pseed;
75bd7f81 2689
69b55c55 2690}
2691
5443e65e 2692//_____________________________________________________________________________
4551ea7c 2693Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *clusterTree) const
5443e65e 2694{
2695 //
a819a5f7 2696 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
2697 // from the file. The names of the cluster tree and branches
2698 // should match the ones used in AliTRDclusterizer::WriteClusters()
2699 //
75bd7f81 2700
4551ea7c 2701 Int_t nsize = Int_t(clusterTree->GetTotBytes() / (sizeof(AliTRDcluster)));
6c94f330 2702 TObjArray *clusterArray = new TObjArray(nsize+1000);
5443e65e 2703
4551ea7c 2704 TBranch *branch = clusterTree->GetBranch("TRDcluster");
c630aafd 2705 if (!branch) {
4551ea7c 2706 AliError("Can't get the branch !");
c630aafd 2707 return 1;
2708 }
029cd327 2709 branch->SetAddress(&clusterArray);
5443e65e 2710
a819a5f7 2711 // Loop through all entries in the tree
4551ea7c 2712 Int_t nEntries = (Int_t) clusterTree->GetEntries();
2713 Int_t nbytes = 0;
a819a5f7 2714 AliTRDcluster *c = 0;
a819a5f7 2715 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2716
2717 // Import the tree
4551ea7c 2718 nbytes += clusterTree->GetEvent(iEntry);
5443e65e 2719
a819a5f7 2720 // Get the number of points in the detector
029cd327 2721 Int_t nCluster = clusterArray->GetEntriesFast();
5443e65e 2722
a819a5f7 2723 // Loop through all TRD digits
2724 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
4551ea7c 2725 c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
4f1c04d3 2726 AliTRDcluster *co = c;
a819a5f7 2727 array->AddLast(co);
4f1c04d3 2728 clusterArray->RemoveAt(iCluster);
a819a5f7 2729 }
4551ea7c 2730
a819a5f7 2731 }
2732
029cd327 2733 delete clusterArray;
5443e65e 2734
c630aafd 2735 return 0;
75bd7f81 2736
a819a5f7 2737}
2738
75bd7f81 2739//_____________________________________________________________________________
4551ea7c 2740Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint &p) const
3551db50 2741{
2742 //
2743 // Get track space point with index i
2744 // Origin: C.Cheshkov
2745 //
2746
4551ea7c 2747 AliTRDcluster *cl = (AliTRDcluster *) fClusters->UncheckedAt(index);
2748 Int_t idet = cl->GetDetector();
2749 Int_t isector = fGeom->GetSector(idet);
2750 Int_t ichamber = fGeom->GetChamber(idet);
2751 Int_t iplan = fGeom->GetPlane(idet);
3551db50 2752 Double_t local[3];
4551ea7c 2753 local[0] = GetX(isector,iplan,cl->GetLocalTimeBin());
2754 local[1] = cl->GetY();
2755 local[2] = cl->GetZ();
3551db50 2756 Double_t global[3];
2757 fGeom->RotateBack(idet,local,global);
2758 p.SetXYZ(global[0],global[1],global[2]);
ae079791 2759 AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
3551db50 2760 switch (iplan) {
2761 case 0:
ae079791 2762 iLayer = AliGeomManager::kTRD1;
3551db50 2763 break;
2764 case 1:
ae079791 2765 iLayer = AliGeomManager::kTRD2;
3551db50 2766 break;
2767 case 2:
ae079791 2768 iLayer = AliGeomManager::kTRD3;
3551db50 2769 break;
2770 case 3:
ae079791 2771 iLayer = AliGeomManager::kTRD4;
3551db50 2772 break;
2773 case 4:
ae079791 2774 iLayer = AliGeomManager::kTRD5;
3551db50 2775 break;
2776 case 5:
ae079791 2777 iLayer = AliGeomManager::kTRD6;
3551db50 2778 break;
2779 };
4551ea7c 2780 Int_t modId = isector * fGeom->Ncham() + ichamber;
ae079791 2781 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,modId);
3551db50 2782 p.SetVolumeID(volid);
2783
2784 return kTRUE;
2785
2786}
2787
75bd7f81 2788//_____________________________________________________________________________
4551ea7c 2789void AliTRDtracker::CookLabel(AliKalmanTrack *pt, Float_t wrong) const
029cd327 2790{
2791 //
2792 // This cooks a label. Mmmmh, smells good...
2793 //
46d29e70 2794
4551ea7c 2795 Int_t label = 123456789;
2796 Int_t index;
2797 Int_t i;
2798 Int_t j;
2799 Int_t ncl = pt->GetNumberOfClusters();
2800
2801 const Int_t kRange = fTrSec[0]->GetOuterTimeBin() + 1;
5443e65e 2802
029cd327 2803 Bool_t labelAdded;
46d29e70 2804
6c94f330 2805 Int_t **s = new Int_t* [kRange];
4551ea7c 2806 for (i = 0; i < kRange; i++) {
d1b06c24 2807 s[i] = new Int_t[2];
2808 }
4551ea7c 2809 for (i = 0; i < kRange; i++) {
2810 s[i][0] = -1;
2811 s[i][1] = 0;
46d29e70 2812 }
2813
4551ea7c 2814 Int_t t0;
2815 Int_t t1;
2816 Int_t t2;
2817
2818 for (i = 0; i < ncl; i++) {
2819 index = pt->GetClusterIndex(i);
2820 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
6c94f330 2821 t0=c->GetLabel(0);
2822 t1=c->GetLabel(1);
2823 t2=c->GetLabel(2);
46d29e70 2824 }
2825
4551ea7c 2826 for (i = 0; i < ncl; i++) {
2827 index = pt->GetClusterIndex(i);
2828 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2829 for (Int_t k = 0; k < 3; k++) {
2830 label = c->GetLabel(k);
2831 labelAdded = kFALSE;
2832 j = 0;
46d29e70 2833 if (label >= 0) {
4551ea7c 2834 while ((!labelAdded) && (j < kRange)) {
2835 if ((s[j][0] == label) ||
2836 (s[j][1] == 0)) {
2837 s[j][0] = label;
2838 s[j][1] = s[j][1] + 1;
2839 labelAdded = kTRUE;
a9814c09 2840 }
2841 j++;
2842 }
46d29e70 2843 }
2844 }
2845 }
2846
4551ea7c 2847 Int_t max = 0;
6c94f330 2848 label = -123456789;
46d29e70 2849
4551ea7c 2850 for (i = 0; i < kRange; i++) {
2851 if (s[i][1] > max) {
2852 max = s[i][1];
2853 label = s[i][0];
46d29e70 2854 }
2855 }
5443e65e 2856
4551ea7c 2857 for (i = 0; i < kRange; i++) {
6c94f330 2858 delete []s[i];
5443e65e 2859 }
2860
6c94f330 2861 delete []s;
5443e65e 2862
4551ea7c 2863 if ((1.0 - Float_t(max)/ncl) > wrong) {
2864 label = -label;
2865 }
5443e65e 2866
2867 pt->SetLabel(label);
2868
46d29e70 2869}
2870
75bd7f81 2871//_____________________________________________________________________________
4551ea7c 2872void AliTRDtracker::UseClusters(const AliKalmanTrack *t, Int_t from) const
029cd327 2873{
2874 //
2875 // Use clusters, but don't abuse them!
2876 //
75bd7f81 2877
4551ea7c 2878 const Float_t kmaxchi2 = 18;
2879 const Float_t kmincl = 10;
2880
2881 AliTRDtrack *track = (AliTRDtrack *) t;
2882
2883 Int_t ncl = t->GetNumberOfClusters();
2884 for (Int_t i = from; i < ncl; i++) {
2885 Int_t index = t->GetClusterIndex(i);
2886 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
69b55c55 2887 Int_t iplane = fGeom->GetPlane(c->GetDetector());
27eaf44b 2888 if (track->GetTracklets(iplane).GetChi2() > kmaxchi2) {
4551ea7c 2889 continue;
2890 }
27eaf44b 2891 if (track->GetTracklets(iplane).GetN() < kmincl) {
4551ea7c 2892 continue;
2893 }
2894 if (!(c->IsUsed())) {
2895 c->Use();
2896 }
5443e65e 2897 }
5443e65e 2898
75bd7f81 2899}
5443e65e 2900
75bd7f81 2901//_____________________________________________________________________________
029cd327 2902Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
5443e65e 2903{
75bd7f81 2904 //
5443e65e 2905 // Parametrised "expected" error of the cluster reconstruction in Y
75bd7f81 2906 //
5443e65e 2907
2908 Double_t s = 0.08 * 0.08;
2909 return s;
75bd7f81 2910
5443e65e 2911}
2912
75bd7f81 2913//_____________________________________________________________________________
029cd327 2914Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
0a29d0f1 2915{
75bd7f81 2916 //
5443e65e 2917 // Parametrised "expected" error of the cluster reconstruction in Z
75bd7f81 2918 //
5443e65e 2919
4551ea7c 2920 Double_t s = 9.0 * 9.0 / 12.0;
5443e65e 2921 return s;
75bd7f81 2922
5443e65e 2923}
2924
75bd7f81 2925//_____________________________________________________________________________
029cd327 2926Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
5443e65e 2927{
2928 //
029cd327 2929 // Returns radial position which corresponds to time bin <localTB>
5443e65e 2930 // in tracking sector <sector> and plane <plane>
2931 //
2932
6c94f330 2933 Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
4551ea7c 2934 Int_t pl = fTrSec[sector]->GetLayerNumber(index);
2935
5443e65e 2936 return fTrSec[sector]->GetLayer(pl)->GetX();
2937
2938}
2939
75bd7f81 2940//_____________________________________________________________________________
2941AliTRDtracker::AliTRDpropagationLayer
2942 ::AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
2943 , Double_t radLength, Int_t tbIndex, Int_t plane)
ad670fba 2944 :fN(0)
2945 ,fSec(0)
2946 ,fClusters(NULL)
2947 ,fIndex(NULL)
2948 ,fX(x)
2949 ,fdX(dx)
2950 ,fRho(rho)
2951 ,fX0(radLength)
2952 ,fTimeBinIndex(tbIndex)
2953 ,fPlane(plane)
2954 ,fYmax(0)
2955 ,fYmaxSensitive(0)
2956 ,fHole(kFALSE)
2957 ,fHoleZc(0)
2958 ,fHoleZmax(0)
2959 ,fHoleYc(0)
2960 ,fHoleYmax(0)
2961 ,fHoleRho(0)
2962 ,fHoleX0(0)
5443e65e 2963{
0a29d0f1 2964 //
5443e65e 2965 // AliTRDpropagationLayer constructor
0a29d0f1 2966 //
46d29e70 2967
ad670fba 2968 for (Int_t i = 0; i < (Int_t) kZones; i++) {
2969 fZc[i] = 0;
2970 fZmax[i] = 0;
a819a5f7 2971 }
5443e65e 2972
ad670fba 2973 if (fTimeBinIndex >= 0) {
029cd327 2974 fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
ad670fba 2975 fIndex = new UInt_t[kMaxClusterPerTimeBin];
a819a5f7 2976 }
46d29e70 2977
ad670fba 2978 for (Int_t i = 0; i < 5; i++) {
2979 fIsHole[i] = kFALSE;
2980 }
5443e65e 2981
2982}
2983
75bd7f81 2984//_____________________________________________________________________________
2985void AliTRDtracker::AliTRDpropagationLayer
2986 ::SetHole(Double_t Zmax, Double_t Ymax, Double_t rho
2987 , Double_t radLength, Double_t Yc, Double_t Zc)
5443e65e 2988{
2989 //
2990 // Sets hole in the layer
2991 //
75bd7f81 2992
4551ea7c 2993 fHole = kTRUE;
2994 fHoleZc = Zc;
5443e65e 2995 fHoleZmax = Zmax;
4551ea7c 2996 fHoleYc = Yc;
5443e65e 2997 fHoleYmax = Ymax;
4551ea7c 2998 fHoleRho = rho;
2999 fHoleX0 = radLength;
75bd7f81 3000
5443e65e 3001}
46d29e70 3002
75bd7f81 3003//_____________________________________________________________________________
3004AliTRDtracker::AliTRDtrackingSector
ad670fba 3005 ::AliTRDtrackingSector(AliTRDgeometry *geo, Int_t gs)
3006 :fN(0)
3007 ,fGeom(geo)
3008 ,fGeomSector(gs)
5443e65e 3009{
3010 //
3011 // AliTRDtrackingSector Constructor
3012 //
75bd7f81 3013
7e448bcc 3014 AliTRDpadPlane *padPlane = 0;
3015 AliTRDpropagationLayer *ppl = 0;
a5cadd36 3016
ad670fba 3017 // Get holes description from geometry
3c625a9b 3018 Bool_t holes[AliTRDgeometry::kNcham];
ad670fba 3019 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3c625a9b 3020 holes[icham] = fGeom->IsHole(0,icham,gs);
3c625a9b 3021 }
3c625a9b 3022
ad670fba 3023 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
3024 fTimeBinIndex[i] = -1;
3025 }
5443e65e 3026
ad670fba 3027 Double_t x;
3028 Double_t dx;
3029 Double_t rho;
3030 Double_t radLength;
5443e65e 3031
ad670fba 3032 // Add layers for each of the planes
3033 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
a305677e 3034 //Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
5443e65e 3035
ad670fba 3036 const Int_t kNchambers = AliTRDgeometry::Ncham();
a305677e 3037 Int_t tbIndex;
ad670fba 3038 Double_t ymax = 0;
3039 Double_t ymaxsensitive = 0;
3040 Double_t *zc = new Double_t[kNchambers];
3041 Double_t *zmax = new Double_t[kNchambers];
3c625a9b 3042 Double_t *zmaxsensitive = new Double_t[kNchambers];
3551db50 3043
ad670fba 3044 for (Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
3045
3046 ymax = fGeom->GetChamberWidth(plane) / 2.0;
f162af62 3047 padPlane = fGeom->GetPadPlane(plane,0);
ad670fba 3048 ymaxsensitive = (padPlane->GetColSize(1) * padPlane->GetNcols() - 4.0) / 2.0;
3049
3050 for (Int_t ch = 0; ch < kNchambers; ch++) {
3051 zmax[ch] = fGeom->GetChamberLength(plane,ch) / 2.0;
3052 Float_t pad = padPlane->GetRowSize(1);
f162af62 3053 Float_t row0 = fGeom->GetRow0(plane,ch,0);
3054 Int_t nPads = fGeom->GetRowMax(plane,ch,0);
ad670fba 3055 zmaxsensitive[ch] = Float_t(nPads) * pad / 2.0;
3056 zc[ch] = -(pad * nPads) / 2.0 + row0;
5443e65e 3057 }
3058
ad670fba 3059 dx = AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
f162af62 3060 / AliTRDCommonParam::Instance()->GetSamplingFrequency();
ad670fba 3061 rho = 0.00295 * 0.85; //????
3062 radLength = 11.0;
5443e65e 3063
3551db50 3064 Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(plane);
a305677e 3065 //Double_t xbottom = x0 - dxDrift;
ad670fba 3066 //Double_t xtop = x0 + dxAmp;
3067
59393e34 3068 Int_t nTimeBins = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
ad670fba 3069 for (Int_t iTime = 0; iTime < nTimeBins; iTime++) {
3070
3071 Double_t xlayer = iTime * dx - dxAmp;
3072 //if (xlayer<0) xlayer = dxAmp / 2.0;
59393e34 3073 x = x0 - xlayer;
ad670fba 3074
3075 tbIndex = CookTimeBinIndex(plane,iTime);
3076 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex,plane);
3c625a9b 3077 ppl->SetYmax(ymax,ymaxsensitive);
ad670fba 3078 ppl->SetZ(zc,zmax,zmaxsensitive);
3c625a9b 3079 ppl->SetHoles(holes);
59393e34 3080 InsertLayer(ppl);
ad670fba 3081
5443e65e 3082 }
ad670fba 3083
5443e65e 3084 }
3085
5443e65e 3086 MapTimeBinLayers();
ad670fba 3087
029cd327 3088 delete [] zc;
3089 delete [] zmax;
4f1c04d3 3090 delete [] zmaxsensitive;
5443e65e 3091
3092}
3093
ad670fba 3094//_____________________________________________________________________________
3095AliTRDtracker::AliTRDtrackingSector
3096 ::AliTRDtrackingSector(const AliTRDtrackingSector &/*t*/)
3097 :fN(0)
3098 ,fGeom(0)
3099 ,fGeomSector(0)
3100{
3101 //
3102 // Copy constructor
3103 //
3104
3105}
3106
75bd7f81 3107//_____________________________________________________________________________
3108Int_t AliTRDtracker::AliTRDtrackingSector
3109 ::CookTimeBinIndex(Int_t plane, Int_t localTB) const
5443e65e 3110{
3111 //
6c94f330 3112 // depending on the digitization parameters calculates "global"
029cd327 3113 // time bin index for timebin <localTB> in plane <plane>
5443e65e 3114 //
59393e34 3115 //
75bd7f81 3116
59393e34 3117 Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
4551ea7c 3118 Int_t gtb = (plane+1) * tbPerPlane - localTB - 1;
3119 if (localTB < 0) {
3120 return -1;
3121 }
3122 if (gtb < 0) {
3123 return -1;
3124 }
75bd7f81 3125
5443e65e 3126 return gtb;
5443e65e 3127
75bd7f81 3128}
5443e65e 3129
75bd7f81 3130//_____________________________________________________________________________
3131void AliTRDtracker::AliTRDtrackingSector
3132 ::MapTimeBinLayers()
5443e65e 3133{
3134 //
3135 // For all sensitive time bins sets corresponding layer index
3136 // in the array fTimeBins
3137 //
3138
3139 Int_t index;
3140
4551ea7c 3141 for (Int_t i = 0; i < fN; i++) {
3142
5443e65e 3143 index = fLayers[i]->GetTimeBinIndex();
6c94f330 3144
4551ea7c 3145 if (index < 0) {
3146 continue;
3147 }
3148 if (index >= (Int_t) kMaxTimeBinIndex) {
3149 //AliWarning(Form("Index %d exceeds allowed maximum of %d!\n"
3150 // ,index,kMaxTimeBinIndex-1));
5443e65e 3151 continue;
3152 }
4551ea7c 3153
5443e65e 3154 fTimeBinIndex[index] = i;
4551ea7c 3155
5443e65e 3156 }
5443e65e 3157
75bd7f81 3158}
5443e65e 3159
75bd7f81 3160//_____________________________________________________________________________
3161Int_t AliTRDtracker::AliTRDtrackingSector
3162 ::GetLayerNumber(Double_t x) const
5443e65e 3163{
3164 //
3165 // Returns the number of time bin which in radial position is closest to <x>
3166 //
3167
4551ea7c 3168 if (x >= fLayers[fN-1]->GetX()) {
3169 return fN - 1;
3170 }
3171 if (x <= fLayers[ 0]->GetX()) {
3172 return 0;
3173 }
3174
3175 Int_t b = 0;
3176 Int_t e = fN - 1;
3177 Int_t m = (b + e) / 2;
ad670fba 3178
4551ea7c 3179 for ( ; b < e; m = (b + e) / 2) {
3180 if (x > fLayers[m]->GetX()) {
3181 b = m + 1;
3182 }
3183 else {
3184 e = m;
3185 }
5443e65e 3186 }
75bd7f81 3187
4551ea7c 3188 if (TMath::Abs(x - fLayers[m]->GetX()) > TMath::Abs(x - fLayers[m+1]->GetX())) {
3189 return m + 1;
3190 }
3191 else {
3192 return m;
3193 }
5443e65e 3194
3195}
3196
75bd7f81 3197//_____________________________________________________________________________
3198Int_t AliTRDtracker::AliTRDtrackingSector
3199 ::GetInnerTimeBin() const
5443e65e 3200{
3201 //
3202 // Returns number of the innermost SENSITIVE propagation layer
3203 //
3204
3205 return GetLayerNumber(0);
5443e65e 3206
75bd7f81 3207}
5443e65e 3208
75bd7f81 3209//_____________________________________________________________________________
3210Int_t AliTRDtracker::AliTRDtrackingSector
3211 ::GetOuterTimeBin() const
5443e65e 3212{
3213 //
3214 // Returns number of the outermost SENSITIVE time bin
3215 //
3216
3217 return GetLayerNumber(GetNumberOfTimeBins() - 1);
46d29e70 3218
75bd7f81 3219}
5443e65e 3220
75bd7f81 3221//_____________________________________________________________________________
3222Int_t AliTRDtracker::AliTRDtrackingSector
3223 ::GetNumberOfTimeBins() const
5443e65e 3224{
3225 //
3226 // Returns number of SENSITIVE time bins
3227 //
3228
4551ea7c 3229 Int_t tb;
3230 Int_t layer;
3231
3232 for (tb = kMaxTimeBinIndex - 1; tb >= 0; tb--) {
5443e65e 3233 layer = GetLayerNumber(tb);
4551ea7c 3234 if (layer >= 0) {
3235 break;
3236 }
5443e65e 3237 }
75bd7f81 3238
4551ea7c 3239 return tb + 1;
5443e65e 3240
75bd7f81 3241}
5443e65e 3242
75bd7f81 3243//_____________________________________________________________________________
3244void AliTRDtracker::AliTRDtrackingSector
4551ea7c 3245 ::InsertLayer(AliTRDpropagationLayer *pl)
5443e65e 3246{
3247 //
3248 // Insert layer <pl> in fLayers array.
3249 // Layers are sorted according to X coordinate.
75bd7f81 3250 //
5443e65e 3251
4551ea7c 3252 if (fN == ((Int_t) kMaxLayersPerSector)) {
3253 //AliWarning("Too many layers !\n");
3254 return;
3255 }
3256
3257 if (fN == 0) {
3258 fLayers[fN++] = pl;
ad670fba 3259 return;
3260 }
3261
4551ea7c 3262 Int_t i = Find(pl->GetX());
3263
3264 memmove(fLayers+i+1,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
3265
3266 fLayers[i] = pl;
3267 fN++;
5443e65e 3268
3269}
3270
75bd7f81 3271//_____________________________________________________________________________
3272Int_t AliTRDtracker::AliTRDtrackingSector
3273 ::Find(Double_t x) const
5443e65e 3274{
3275 //
3276 // Returns index of the propagation layer nearest to X
3277 //
3278
4551ea7c 3279 if (x <= fLayers[0]->GetX()) {
3280 return 0;
3281 }
3282
3283 if (x > fLayers[fN-1]->GetX()) {
3284 return fN;
3285 }
3286
3287 Int_t b = 0;
3288 Int_t e = fN-1;
3289 Int_t m = (b + e) / 2;
3290
3291 for (; b < e; m = (b + e) / 2) {
3292 if (x > fLayers[m]->GetX()) {
3293 b = m + 1;
3294 }
3295 else {
3296 e = m;
3297 }
5443e65e 3298 }
7ad19338 3299
75bd7f81 3300 return m;
7ad19338 3301
75bd7f81 3302}
7ad19338 3303
75bd7f81 3304//_____________________________________________________________________________
3305void AliTRDtracker::AliTRDpropagationLayer
4551ea7c 3306 ::SetZ(Double_t *center, Double_t *w, Double_t *wsensitive )
3c625a9b 3307{
3308 //
6c94f330 3309 // set centers and the width of sectors
75bd7f81 3310 //
3311
4551ea7c 3312 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3313 fZc[icham] = center[icham];
3314 fZmax[icham] = w[icham];
3c625a9b 3315 fZmaxSensitive[icham] = wsensitive[icham];
3c625a9b 3316 }
75bd7f81 3317
3c625a9b 3318}
5443e65e 3319
75bd7f81 3320//_____________________________________________________________________________
3c625a9b 3321void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
3322{
3323 //
6c94f330 3324 // set centers and the width of sectors
75bd7f81 3325 //
3326
3c625a9b 3327 fHole = kFALSE;
4551ea7c 3328
3329 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3c625a9b 3330 fIsHole[icham] = holes[icham];
4551ea7c 3331 if (holes[icham]) {
3332 fHole = kTRUE;
3333 }
3c625a9b 3334 }
5443e65e 3335
75bd7f81 3336}
5443e65e 3337
75bd7f81 3338//_____________________________________________________________________________
3339void AliTRDtracker::AliTRDpropagationLayer
4551ea7c 3340 ::InsertCluster(AliTRDcluster *c, UInt_t index)
75bd7f81 3341{
3342 //
3343 // Insert cluster in cluster array.
3344 // Clusters are sorted according to Y coordinate.
3345 //
5443e65e 3346
4551ea7c 3347 if (fTimeBinIndex < 0) {
3348 //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");
3349 return;
3350 }
3351
3352 if (fN == (Int_t) kMaxClusterPerTimeBin) {
3353 //AliWarning("Too many clusters !\n");
5443e65e 3354 return;
3355 }
ad670fba 3356
4551ea7c 3357 if (fN == 0) {
3358 fIndex[0] = index;
3359 fClusters[fN++] = c;
ad670fba 3360 return;
3361 }
4551ea7c 3362
3363 Int_t i = Find(c->GetY());
3364 memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
3365 memmove(fIndex +i+1,fIndex +i,(fN-i)*sizeof(UInt_t));
3366 fIndex[i] = index;
3367 fClusters[i] = c;
3368 fN++;
5443e65e 3369
75bd7f81 3370}
5443e65e 3371
75bd7f81 3372//_____________________________________________________________________________
3373Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const
3374{
3375 //
3376 // Returns index of the cluster nearest in Y
3377 //
5443e65e 3378
4551ea7c 3379 if (fN <= 0) {
3380 return 0;
3381 }
3382 if (y <= fClusters[0]->GetY()) {
3383 return 0;
3384 }
3385 if (y > fClusters[fN-1]->GetY()) {
3386 return fN;
3387 }
3388
3389 Int_t b = 0;
3390 Int_t e = fN - 1;
3391 Int_t m = (b + e) / 2;
3392
3393 for ( ; b < e; m = (b + e) / 2) {
3394 if (y > fClusters[m]->GetY()) {
3395 b = m + 1;
3396 }
3397 else {
3398 e = m;
3399 }
5443e65e 3400 }
75bd7f81 3401
5443e65e 3402 return m;
75bd7f81 3403
5443e65e 3404}
3405
75bd7f81 3406//_____________________________________________________________________________
3407Int_t AliTRDtracker::AliTRDpropagationLayer
3408 ::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad
3409 , Float_t maxroadz) const
7ad19338 3410{
3411 //
3412 // Returns index of the cluster nearest to the given y,z
3413 //
75bd7f81 3414
4551ea7c 3415 Int_t index = -1;
3416 Int_t maxn = fN;
69b55c55 3417 Float_t mindist = maxroad;
4551ea7c 3418
3419 for (Int_t i = Find(y-maxroad); i < maxn; i++) {
3420 AliTRDcluster *c = (AliTRDcluster *) (fClusters[i]);
69b55c55 3421 Float_t ycl = c->GetY();
4551ea7c 3422 if (ycl > (y + maxroad)) {
3423 break;
3424 }
3425 if (TMath::Abs(c->GetZ() - z) > maxroadz) {
3426 continue;
3427 }
3428 if (TMath::Abs(ycl - y) < mindist) {
3429 mindist = TMath::Abs(ycl - y);
3430 index = fIndex[i];
3431 }
6c94f330 3432 }
75bd7f81 3433
7ad19338 3434 return index;
7ad19338 3435
75bd7f81 3436}
7ad19338 3437
75bd7f81 3438//_____________________________________________________________________________
4551ea7c 3439Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c)
75bd7f81 3440{
3441 //
3442 // Returns correction factor for tilted pads geometry
3443 //
5443e65e 3444
4551ea7c 3445 Int_t det = c->GetDetector();
3446 Int_t plane = fGeom->GetPlane(det);
f162af62 3447 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(plane,0);
4551ea7c 3448 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
b8dc2353 3449
4551ea7c 3450 if (fNoTilt) {
3451 h01 = 0;
3452 }
75bd7f81 3453
fd621f36 3454 return h01;
5443e65e 3455
75bd7f81 3456}
c630aafd 3457
3458
75bd7f81 3459//_____________________________________________________________________________
3460Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
4551ea7c 3461 , AliTRDtrack *track
3462 , Int_t *clusters, AliTRDtracklet &tracklet)
4f1c04d3 3463{
6c94f330 3464 //
4f1c04d3 3465 //
75bd7f81 3466 // Try to find nearest clusters to the track in timebins from t0 to t1
4f1c04d3 3467 //
4551ea7c 3468 // Correction coeficients - depend on TRD parameters - to be changed accordingly
3469 //
3470
3471 Double_t x[100];
3472 Double_t yt[100];
3473 Double_t zt[100];
3474 Double_t xmean = 0.0; // Reference x
3475 Double_t dz[10][100];
3476 Double_t dy[10][100];
3477 Float_t zmean[100];
3478 Float_t nmean[100];
3479 Int_t clfound = 0;
3480 Int_t indexes[10][100]; // Indexes of the clusters in the road
3481 Int_t best[10][100]; // Index of best matching cluster
3482 AliTRDcluster *cl[10][100]; // Pointers to the clusters in the road
3483
3484 for (Int_t it = 0; it < 100; it++) {
3485 x[it] = 0.0;
3486 yt[it] = 0.0;
3487 zt[it] = 0.0;
3488 clusters[it] = -2;
3489 zmean[it] = 0.0;
3490 nmean[it] = 0.0;
3491 for (Int_t ih = 0; ih < 10;ih++) {
3492 indexes[ih][it] = -2; // Reset indexes1
3493 cl[ih][it] = 0;
3494 dz[ih][it] = -100.0;
3495 dy[ih][it] = -100.0;
3496 best[ih][it] = 0;
3497 }
3498 }
ad670fba 3499
4551ea7c 3500 Double_t x0 = track->GetX();
3501 Double_t sigmaz = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
3502 Int_t nall = 0;
3503 Int_t nfound = 0;
3504 Double_t h01 = 0.0;
3505 Int_t plane = -1;
3506 Int_t detector = -1;
3507 Float_t padlength = 0.0;
3508 AliTRDtrack track2(* track);
3509 Float_t snpy = track->GetSnp();
3510 Float_t tany = TMath::Sqrt(snpy*snpy / (1.0 - snpy*snpy));
3511 if (snpy < 0.0) {
3512 tany *= -1.0;
3513 }
ad670fba 3514
4551ea7c 3515 Double_t sy2 = ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt());
3516 Double_t sz2 = ExpectedSigmaZ2(x0,track->GetTgl());
3517 Double_t road = 15.0 * TMath::Sqrt(track->GetSigmaY2() + sy2);
3518 if (road > 6.0) {
3519 road = 6.0;
3520 }
7e448bcc 3521 //road = 20.0;
4551ea7c 3522
3523 for (Int_t it = 0; it < t1-t0; it++) {
3524
3525 Double_t maxChi2[2] = { fgkMaxChi2, fgkMaxChi2 };
3526 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(it+t0));
3527 if (timeBin == 0) {
3528 continue; // No indexes1
6c94f330 3529 }
4551ea7c 3530
3531 Int_t maxn = timeBin;
3532 x[it] = timeBin.GetX();
7ad19338 3533 track2.PropagateTo(x[it]);
6c94f330 3534 yt[it] = track2.GetY();
3535 zt[it] = track2.GetZ();
7ad19338 3536
4551ea7c 3537 Double_t y = yt[it];
3538 Double_t z = zt[it];
3539 Double_t chi2 = 1000000.0;
4f1c04d3 3540 nall++;
4551ea7c 3541
4f1c04d3 3542 //
4551ea7c 3543 // Find 2 nearest cluster at given time bin
6c94f330 3544 //
7e448bcc 3545 int checkPoint[4] = {0,0,0,0};
3546 double minY = 123456789;
3547 double minD[2] = {1,1};
3548
4551ea7c 3549 for (Int_t i = timeBin.Find(y - road); i < maxn; i++) {
7e448bcc 3550 //for (Int_t i = 0; i < maxn; i++) {
4551ea7c 3551
3552 AliTRDcluster *c = (AliTRDcluster *) (timeBin[i]);
7ad19338 3553 h01 = GetTiltFactor(c);
4551ea7c 3554 if (plane < 0) {
c6f438c0 3555 Int_t det = c->GetDetector();
4551ea7c 3556 plane = fGeom->GetPlane(det);
3557 padlength = TMath::Sqrt(c->GetSigmaZ2() * 12.0);
3558 }
7e448bcc 3559
4551ea7c 3560 //if (c->GetLocalTimeBin()==0) continue;
3561 if (c->GetY() > (y + road)) {
3562 break;
3563 }
7e448bcc 3564
3565 fHDeltaX->Fill(c->GetX() - x[it]);
3566 //printf("%f\t%f\t%f \n", c->GetX(), x[it], c->GetX()-x[it]);
3567
3568 if (TMath::Abs(c->GetY()-y) < TMath::Abs(minY)) {
3569 minY = c->GetY()-y;
3570 minD[0] = c->GetY()-y;
3571 minD[1] = c->GetZ()-z;
3572 }
3573
3574 checkPoint[0]++;
3575
3576 fHMinZ->Fill(c->GetZ() - z);
3577 if ((c->GetZ() - z) * (c->GetZ() - z) > 2 * (12.0 * sz2)) {
4551ea7c 3578 continue;
7ad19338 3579 }
7e448bcc 3580 checkPoint[1]++;
7ad19338 3581
4551ea7c 3582 Double_t dist = TMath::Abs(c->GetZ() - z);
7e448bcc 3583 if (dist > (0.5 * padlength + 6.0 * sigmaz)) { // 0.5
4551ea7c 3584 continue; // 6 sigma boundary cut
3585 }
7e448bcc 3586 checkPoint[2]++;
3587
4551ea7c 3588 Double_t cost = 0.0;
3589 // Sigma boundary cost function
3590 if (dist> (0.5 * padlength - sigmaz)){
3591 cost = (dist - 0.5*padlength) / (2.0 * sigmaz);
3592 if (cost > -1) {
3593 cost = (cost + 1.0) * (cost + 1.0);
3594 }
3595 else {
3596 cost = 0.0;
3597 }
3598 }
3599 //Int_t label = TMath::Abs(track->GetLabel());
3600 //if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue;
3601 chi2 = track2.GetPredictedChi2(c,h01) + cost;
7ad19338 3602 clfound++;
6c94f330 3603
4551ea7c 3604 if (chi2 > maxChi2[1]) {
3605 continue;
3606 }
7e448bcc 3607 checkPoint[3]++;
3608
4551ea7c 3609 detector = c->GetDetector();
3610 // Store the clusters in the road
3611 for (Int_t ih = 2; ih < 9; ih++) {
3612 if (cl[ih][it] == 0) {
3613 cl[ih][it] = c;
3614 indexes[ih][it] = timeBin.GetIndex(i); // Index - 9 - reserved for outliers
7ad19338 3615 break;
3616 }
4f1c04d3 3617 }
4551ea7c 3618
3619 if (chi2 < maxChi2[0]) {
7ad19338 3620 maxChi2[1] = maxChi2[0];
3621 maxChi2[0] = chi2;
3622 indexes[1][it] = indexes[0][it];
3623 cl[1][it] = cl[0][it];
3624 indexes[0][it] = timeBin.GetIndex(i);
3625 cl[0][it] = c;
3626 continue;
3627 }
4551ea7c 3628 maxChi2[1] = chi2;
3629 cl[1][it] = c;
3630 indexes[1][it] = timeBin.GetIndex(i);
3631
3632 }
7e448bcc 3633
3634 for(int iCheckPoint = 0; iCheckPoint<4; iCheckPoint++)
3635 fHFindCl[iCheckPoint]->Fill(checkPoint[iCheckPoint]);
3636
3637 if (checkPoint[3]) {
3638 if (track->GetPt() > 0) fHMinYPos->Fill(minY);
3639 else fHMinYNeg->Fill(minY);
3640
3641 fHMinD->Fill(minD[0], minD[1]);
3642 }
4551ea7c 3643
3644 if (cl[0][it]) {
7ad19338 3645 nfound++;
3646 xmean += x[it];
3647 }
4551ea7c 3648
4f1c04d3 3649 }
4551ea7c 3650
3651 if (nfound < 4) {
3652 return 0;
3653 }
3654 xmean /= Float_t(nfound); // Middle x
3655 track2.PropagateTo(xmean); // Propagate track to the center
3656
4f1c04d3 3657 //
4551ea7c 3658 // Choose one of the variants
7ad19338 3659 //
4551ea7c 3660 Int_t changes[10];
3661 Float_t sumz = 0.0;
3662 Float_t sum = 0.0;
3663 Double_t sumdy = 0.0;
3664 Double_t sumdy2 = 0.0;
3665 Double_t sumx = 0.0;
3666 Double_t sumxy = 0.0;
3667 Double_t sumx2 = 0.0;
3668 Double_t mpads = 0.0;
3669
3670 Int_t ngood[10];
3671 Int_t nbad[10];
3672
7ad19338 3673 Double_t meanz[10];
4551ea7c 3674 Double_t moffset[10]; // Mean offset
3675 Double_t mean[10]; // Mean value
3676 Double_t angle[10]; // Angle
3677
3678 Double_t smoffset[10]; // Sigma of mean offset
3679 Double_t smean[10]; // Sigma of mean value
3680 Double_t sangle[10]; // Sigma of angle
3681 Double_t smeanangle[10]; // Correlation
3682
7ad19338 3683 Double_t sigmas[10];
4551ea7c 3684 Double_t tchi2s[10]; // Chi2s for tracklet
ad670fba 3685
4551ea7c 3686 for (Int_t it = 0; it < 10; it++) {
ad670fba 3687
4551ea7c 3688 ngood[it] = 0;
3689 nbad[it] = 0;
3690
3691 meanz[it] = 0.0;
3692 moffset[it] = 0.0; // Mean offset
3693 mean[it] = 0.0; // Mean value
3694 angle[it] = 0.0; // Angle
3695
7e448bcc 3696 smoffset[it] = 1.0e5; // Sigma of mean offset
3697 smean[it] = 1.0e5; // Sigma of mean value
3698 sangle[it] = 1.0e5; // Sigma of angle
4551ea7c 3699 smeanangle[it] = 0.0; // Correlation
3700
7e448bcc 3701 sigmas[it] = 1.0e5;
3702 tchi2s[it] = 1.0e5; // Chi2s for tracklet
75bd7f81 3703
3704 }
3705
7ad19338 3706 //
4551ea7c 3707 // Calculate zmean
7ad19338 3708 //
4551ea7c 3709 for (Int_t it = 0; it < t1 - t0; it++) {
3710 if (!cl[0][it]) {
3711 continue;
3712 }
3713 for (Int_t dt = -3; dt <= 3; dt++) {
3714 if (it+dt < 0) {
3715 continue;
3716 }
3717 if (it+dt > t1-t0) {
3718 continue;
3719 }
3720 if (!cl[0][it+dt]) {
3721 continue;
3722 }
3723 zmean[it] += cl[0][it+dt]->GetZ();
3724 nmean[it] += 1.0;
7ad19338 3725 }
4551ea7c 3726 zmean[it] /= nmean[it];
7ad19338 3727 }
4551ea7c 3728
3729 for (Int_t it = 0; it < t1 - t0; it++) {
3730
3731 best[0][it] = 0;
3732
3733 for (Int_t ih = 0; ih < 10; ih++) {
3734 dz[ih][it] = -100.0;
3735 dy[ih][it] = -100.0;
3736 if (!cl[ih][it]) {
3737 continue;
3738 }
3739 Double_t xcluster = cl[ih][it]->GetX();
3740 Double_t ytrack;
3741 Double_t ztrack;
3742 track2.GetProlongation(xcluster,ytrack,ztrack );
3743 dz[ih][it] = cl[ih][it]->GetZ()- ztrack; // Calculate distance from track in z
3744 dy[ih][it] = cl[ih][it]->GetY() + dz[ih][it]*h01 - ytrack; // and in y
7ad19338 3745 }
4551ea7c 3746
3747 // Minimize changes
3748 if (!cl[0][it]) {
3749 continue;
3750 }
3751 if ((TMath::Abs(cl[0][it]->GetZ()-zmean[it]) > padlength * 0.8) &&
3752 (cl[1][it])) {
3753 if (TMath::Abs(cl[1][it]->GetZ()-zmean[it]) < padlength * 0.5) {
3754 best[0][it] = 1;
4f1c04d3 3755 }
4551ea7c 3756 }
3757
7ad19338 3758 }
4551ea7c 3759
7ad19338 3760 //
4551ea7c 3761 // Iterative choice of "best path"
6c94f330 3762 //
4551ea7c 3763 Int_t label = TMath::Abs(track->GetLabel());
3764 Int_t bestiter = 0;
3765
3766 for (Int_t iter = 0; iter < 9; iter++) {
3767
3768 changes[iter] = 0;
3769 sumz = 0;
3770 sum = 0;
3771 sumdy = 0;
3772 sumdy2 = 0;
3773 sumx = 0;
3774 sumx2 = 0;
3775 sumxy = 0;
3776 mpads = 0;
3777 ngood[iter] = 0;
3778 nbad[iter] = 0;
3779
3780 // Linear fit
3781 for (Int_t it = 0; it < t1 - t0; it++) {
3782
3783 if (!cl[best[iter][it]][it]) {
3784 continue;
3785 }
3786
3787 // Calculates pad-row changes
3788 Double_t zbefore = cl[best[iter][it]][it]->GetZ();
3789 Double_t zafter = cl[best[iter][it]][it]->GetZ();
3790 for (Int_t itd = it - 1; itd >= 0; itd--) {
7ad19338 3791 if (cl[best[iter][itd]][itd]) {
4551ea7c 3792 zbefore = cl[best[iter][itd]][itd]->GetZ();
7ad19338 3793 break;
3794 }
3795 }
4551ea7c 3796 for (Int_t itd = it + 1; itd < t1 - t0; itd++) {
7ad19338 3797 if (cl[best[iter][itd]][itd]) {
4551ea7c 3798 zafter = cl[best[iter][itd]][itd]->GetZ();
7ad19338 3799 break;
3800 }
3801 }
4551ea7c 3802 if ((TMath::Abs(cl[best[iter][it]][it]->GetZ()-zbefore) > 0.1) &&
3803 (TMath::Abs(cl[best[iter][it]][it]->GetZ()- zafter) > 0.1)) {
3804 changes[iter]++;
3805 }
3806
3807 Double_t dx = x[it]-xmean; // Distance to reference x
3808 sumz += cl[best[iter][it]][it]->GetZ();
4f1c04d3 3809 sum++;
4551ea7c 3810 sumdy += dy[best[iter][it]][it];
3811 sumdy2 += dy[best[iter][it]][it]*dy[best[iter][it]][it];
3812 sumx += dx;
3813 sumx2 += dx*dx;
7ad19338 3814 sumxy += dx*dy[best[iter][it]][it];
4551ea7c 3815 mpads += cl[best[iter][it]][it]->GetNPads();
3816 if ((cl[best[iter][it]][it]->GetLabel(0) == label) ||
3817 (cl[best[iter][it]][it]->GetLabel(1) == label) ||
3818 (cl[best[iter][it]][it]->GetLabel(2) == label)) {
7ad19338 3819 ngood[iter]++;
4f1c04d3 3820 }
4551ea7c 3821 else {
7ad19338 3822 nbad[iter]++;
4f1c04d3 3823 }
4551ea7c 3824
4f1c04d3 3825 }
4551ea7c 3826
7ad19338 3827 //
6c94f330 3828 // calculates line parameters
7ad19338 3829 //
4551ea7c 3830 Double_t det = sum*sumx2 - sumx*sumx;
3831 angle[iter] = (sum*sumxy - sumx*sumdy) / det;
3832 mean[iter] = (sumx2*sumdy - sumx*sumxy) / det;
3833 meanz[iter] = sumz / sum;
3834 moffset[iter] = sumdy / sum;
3835 mpads /= sum; // Mean number of pads
3836
3837 Double_t sigma2 = 0.0; // Normalized residuals - for line fit
3838 Double_t sigma1 = 0.0; // Normalized residuals - constant fit
3839
3840 for (Int_t it = 0; it < t1 - t0; it++) {
3841 if (!cl[best[iter][it]][it]) {
3842 continue;
3843 }
3844 Double_t dx = x[it] - xmean;
3845 Double_t ytr = mean[iter] + angle[iter] * dx;
3846 sigma2 += (dy[best[iter][it]][it] - ytr)
3847 * (dy[best[iter][it]][it] - ytr);
3848 sigma1 += (dy[best[iter][it]][it] - moffset[iter])
3849 * (dy[best[iter][it]][it] - moffset[iter]);
7ad19338 3850 sum++;
4f1c04d3 3851 }
4551ea7c 3852 sigma2 /= (sum - 2); // Normalized residuals
3853 sigma1 /= (sum - 1); // Normalized residuals
3854 smean[iter] = sigma2 * (sumx2 / det); // Estimated error2 of mean
3855 sangle[iter] = sigma2 * ( sum / det); // Estimated error2 of angle
3856 smeanangle[iter] = sigma2 * (-sumx / det); // Correlation
3857 sigmas[iter] = TMath::Sqrt(sigma1);
3858 smoffset[iter] = (sigma1 / sum) + 0.01*0.01; // Sigma of mean offset + unisochronity sigma
3859
6c94f330 3860 //
4551ea7c 3861 // Iterative choice of "better path"
6c94f330 3862 //
4551ea7c 3863 for (Int_t it = 0; it < t1 - t0; it++) {
3864
3865 if (!cl[best[iter][it]][it]) {
3866 continue;
3867 }
3868
3869 // Add unisochronity + angular effect contribution
3870 Double_t sigmatr2 = smoffset[iter] + 0.5*tany*tany;
3871 Double_t sweight = 1.0/sigmatr2 + 1.0/track->GetSigmaY2();
3872 Double_t weighty = (moffset[iter] / sigmatr2) / sweight; // Weighted mean
3873 Double_t sigmacl = TMath::Sqrt(sigma1*sigma1 + track->GetSigmaY2());
3874 Double_t mindist = 100000.0;
3875 Int_t ihbest = 0;
3876
3877 for (Int_t ih = 0; ih < 10; ih++) {
3878 if (!cl[ih][it]) {
3879 break;
3880 }
3881 Double_t dist2 = (dy[ih][it] - weighty) / sigmacl;
3882 dist2 *= dist2; // Chi2 distance
3883 if (dist2 < mindist) {
7ad19338 3884 mindist = dist2;
4551ea7c 3885 ihbest = ih;
7ad19338 3886 }
3887 }
4551ea7c 3888 best[iter+1][it] = ihbest;
4f1c04d3 3889 }
4551ea7c 3890
4f1c04d3 3891 //
4551ea7c 3892 // Update best hypothesy if better chi2 according tracklet position and angle
7ad19338 3893 //
69b55c55 3894 Double_t sy2 = smean[iter] + track->GetSigmaY2();
4551ea7c 3895 Double_t sa2 = sangle[iter] + track->GetSigmaSnp2(); // track->fCee;
3896 Double_t say = track->GetSigmaSnpY(); // track->fCey;
3897 //Double_t chi20 = mean[bestiter]*mean[bestiter ] / sy2+angle[bestiter]*angle[bestiter]/sa2;
3898 //Double_t chi21 = mean[iter]*mean[iter] / sy2+angle[iter]*angle[iter]/sa2;
6c94f330 3899
4551ea7c 3900 Double_t detchi = sy2*sa2 - say*say;
3901 Double_t invers[3] = {sa2/detchi,sy2/detchi,-say/detchi}; // Inverse value of covariance matrix
4f1c04d3 3902
4551ea7c 3903 Double_t chi20 = mean[bestiter] * mean[bestiter] * invers[0]
3904 + angle[bestiter] * angle[bestiter] * invers[1]
3905 + 2.0 * mean[bestiter] * angle[bestiter] * invers[2];
3906 Double_t chi21 = mean[iter] * mean[iter] * invers[0]
3907 + angle[iter] * angle[iter] * invers[1]
3908 + 2.0 * mean[iter] * angle[iter] * invers[2];
3909 tchi2s[iter] = chi21;
3910
3911 if ((changes[iter] <= changes[bestiter]) &&
3912 (chi21 < chi20)) {
3913 bestiter = iter;
7ad19338 3914 }
4551ea7c 3915
7ad19338 3916 }
4551ea7c 3917
7ad19338 3918 //
4551ea7c 3919 // Set clusters
7ad19338 3920 //
4551ea7c 3921 Double_t sigma2 = sigmas[0]; // Choose as sigma from 0 iteration
3922 Short_t maxpos = -1;
3923 Float_t maxcharge = 0.0;
3924 Short_t maxpos4 = -1;
3925 Float_t maxcharge4 = 0.0;
3926 Short_t maxpos5 = -1;
3927 Float_t maxcharge5 = 0.0;
8979685e 3928
7ad19338 3929 //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
3930 //if (tchi2s[bestiter]>25.) sigma2=1000.; // dont'accept
3931
4e009ce4 3932 Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
3933 ,-AliTracker::GetBz()*0.1);
4551ea7c 3934 Double_t expectederr = sigma2*sigma2 + 0.01*0.01;
3935 if (mpads > 3.5) {
3936 expectederr += (mpads - 3.5) * 0.04;
3937 }
3938 if (changes[bestiter] > 1) {
3939 expectederr += changes[bestiter] * 0.01;
3940 }
3941 expectederr += (0.03 * (tany-exB)*(tany-exB)) * 15.0;
3942 //if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.;
7ad19338 3943 //expectederr+=10000;
4551ea7c 3944
3945 for (Int_t it = 0; it < t1 - t0; it++) {
3946
3947 if (!cl[best[bestiter][it]][it]) {
3948 continue;
69b55c55 3949 }
4551ea7c 3950
3951 cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); // Set cluster error
3952 if (!cl[best[bestiter][it]][it]->IsUsed()) {
3953 cl[best[bestiter][it]][it]->SetY(cl[best[bestiter][it]][it]->GetY());
3954 //cl[best[bestiter][it]][it]->Use();
3955 }
3956
3957 // Time bins with maximal charge
3958 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge) {
69b55c55 3959 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3960 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 3961 }
6c94f330 3962
4551ea7c 3963 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
3964 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
69b55c55 3965 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3966 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 3967 }
3968 }
4551ea7c 3969
3970 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
3971 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
69b55c55 3972 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3973 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 3974 }
7ad19338 3975 }
4551ea7c 3976
3977 // Time bins with maximal charge
3978 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge) {
8979685e 3979 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3980 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 3981 }
6c94f330 3982
4551ea7c 3983 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
3984 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
8979685e 3985 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3986 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 3987 }
3988 }
4551ea7c 3989
3990 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
3991 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
8979685e 3992 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3993 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 3994 }
3995 }
4551ea7c 3996
7ad19338 3997 clusters[it+t0] = indexes[best[bestiter][it]][it];
4551ea7c 3998
3999 // Still needed ????
4000 //if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>4 &&
4001 //cl[best[bestiter][it]][it]->GetLocalTimeBin()<18) clusters[it+t0]
4002 // = indexes[best[bestiter][it]][it]; //Test
4003
7ad19338 4004 }
4551ea7c 4005
7ad19338 4006 //
4551ea7c 4007 // Set tracklet parameters
6c94f330 4008 //
4551ea7c 4009 Double_t trackleterr2 = smoffset[bestiter] + 0.01*0.01;
4010 if (mpads > 3.5) {
4011 trackleterr2 += (mpads - 3.5) * 0.04;
4012 }
4013 trackleterr2 += changes[bestiter] * 0.01;
4014 trackleterr2 *= TMath::Max(14.0 - nfound,1.0);
4015 trackleterr2 += 0.2 * (tany-exB)*(tany-exB);
4016
4017 // Set tracklet parameters
4018 tracklet.Set(xmean
4019 ,track2.GetY() + moffset[bestiter]
4020 ,meanz[bestiter]
4021 ,track2.GetAlpha()
4022 ,trackleterr2);
7ad19338 4023 tracklet.SetTilt(h01);
4024 tracklet.SetP0(mean[bestiter]);
4025 tracklet.SetP1(angle[bestiter]);
4026 tracklet.SetN(nfound);
4027 tracklet.SetNCross(changes[bestiter]);
4028 tracklet.SetPlane(plane);
4029 tracklet.SetSigma2(expectederr);
4030 tracklet.SetChi2(tchi2s[bestiter]);
8979685e 4031 tracklet.SetMaxPos(maxpos,maxpos4,maxpos5);
27eaf44b 4032 track->SetTracklets(plane,tracklet);
4033 track->SetNWrong(track->GetNWrong() + nbad[0]);
4551ea7c 4034
7ad19338 4035 //
4036 // Debuging part
4037 //
69b55c55 4038 TClonesArray array0("AliTRDcluster");
4039 TClonesArray array1("AliTRDcluster");
4551ea7c 4040 array0.ExpandCreateFast(t1 - t0 + 1);
4041 array1.ExpandCreateFast(t1 - t0 + 1);
4042 TTreeSRedirector &cstream = *fDebugStreamer;
7ad19338 4043 AliTRDcluster dummy;
4044 Double_t dy0[100];
8979685e 4045 Double_t dyb[100];
4046
4551ea7c 4047 for (Int_t it = 0; it < t1 - t0; it++) {
7ad19338 4048 dy0[it] = dy[0][it];
4049 dyb[it] = dy[best[bestiter][it]][it];
4551ea7c 4050 if (cl[0][it]) {
7ad19338 4051 new(array0[it]) AliTRDcluster(*cl[0][it]);
4052 }
4551ea7c 4053 else {
7ad19338 4054 new(array0[it]) AliTRDcluster(dummy);
4055 }
4056 if(cl[best[bestiter][it]][it]) {
4057 new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
4058 }
4059 else{
4060 new(array1[it]) AliTRDcluster(dummy);
4061 }
4f1c04d3 4062 }
7ad19338 4063 TGraph graph0(t1-t0,x,dy0);
4064 TGraph graph1(t1-t0,x,dyb);
4065 TGraph graphy(t1-t0,x,yt);
4066 TGraph graphz(t1-t0,x,zt);
4551ea7c 4067
4068 if (AliTRDReconstructor::StreamLevel() > 0) {
4069 cstream << "tracklet"
4070 << "track.=" << track // Track parameters
4071 << "tany=" << tany // Tangent of the local track angle
4072 << "xmean=" << xmean // Xmean - reference x of tracklet
4073 << "tilt=" << h01 // Tilt angle
4074 << "nall=" << nall // Number of foundable clusters
4075 << "nfound=" << nfound // Number of found clusters
4076 << "clfound=" << clfound // Total number of found clusters in road
4077 << "mpads=" << mpads // Mean number of pads per cluster
4078 << "plane=" << plane // Plane number
4079 << "detector=" << detector // Detector number
4080 << "road=" << road // The width of the used road
4081 << "graph0.=" << &graph0 // x - y = dy for closest cluster
4082 << "graph1.=" << &graph1 // x - y = dy for second closest cluster
4083 << "graphy.=" << &graphy // y position of the track
4084 << "graphz.=" << &graphz // z position of the track
4085 //<< "fCl.=" << &array0 // closest cluster
4086 //<< "fCl2.=" << &array1 // second closest cluster
4087 << "maxpos=" << maxpos // Maximal charge postion
4088 << "maxcharge=" << maxcharge // Maximal charge
4089 << "maxpos4=" << maxpos4 // Maximal charge postion - after bin 4
4090 << "maxcharge4=" << maxcharge4 // Maximal charge - after bin 4
4091 << "maxpos5=" << maxpos5 // Maximal charge postion - after bin 5
4092 << "maxcharge5=" << maxcharge5 // Maximal charge - after bin 5
4093 << "bestiter=" << bestiter // Best iteration number
4094 << "tracklet.=" << &tracklet // Corrspond to the best iteration
4095 << "tchi20=" << tchi2s[0] // Chi2 of cluster in the 0 iteration
4096 << "tchi2b=" << tchi2s[bestiter] // Chi2 of cluster in the best iteration
4097 << "sigmas0=" << sigmas[0] // Residuals sigma
4098 << "sigmasb=" << sigmas[bestiter] // Residulas sigma
4099 << "ngood0=" << ngood[0] // Number of good clusters in 0 iteration
4100 << "nbad0=" << nbad[0] // Number of bad clusters in 0 iteration
4101 << "ngoodb=" << ngood[bestiter] // in best iteration
4102 << "nbadb=" << nbad[bestiter] // in best iteration
4103 << "changes0=" << changes[0] // Changes of pardrows in iteration number 0
4104 << "changesb=" << changes[bestiter] // Changes of pardrows in best iteration
4105 << "moffset0=" << moffset[0] // Offset fixing angle in iter=0
4106 << "smoffset0=" << smoffset[0] // Sigma of offset fixing angle in iter=0
4107 << "moffsetb=" << moffset[bestiter] // Offset fixing angle in iter=best
4108 << "smoffsetb=" << smoffset[bestiter] // Sigma of offset fixing angle in iter=best
4109 << "mean0=" << mean[0] // Mean dy in iter=0;
4110 << "smean0=" << smean[0] // Sigma of mean dy in iter=0
4111 << "meanb=" << mean[bestiter] // Mean dy in iter=best
4112 << "smeanb=" << smean[bestiter] // Sigma of mean dy in iter=best
4113 << "angle0=" << angle[0] // Angle deviation in the iteration number 0
4114 << "sangle0=" << sangle[0] // Sigma of angular deviation in iteration number 0
4115 << "angleb=" << angle[bestiter] // Angle deviation in the best iteration
4116 << "sangleb=" << sangle[bestiter] // Sigma of angle deviation in the best iteration
4117 << "expectederr=" << expectederr // Expected error of cluster position
4118 << "\n";
4119 }
4120
4f1c04d3 4121 return nfound;
4f1c04d3 4122
75bd7f81 4123}
4f1c04d3 4124
75bd7f81 4125//_____________________________________________________________________________
4126Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist
4127 , Int_t *outlist, Bool_t down)
69b55c55 4128{
4129 //
4551ea7c 4130 // Sort eleements according occurancy
4131 // The size of output array has is 2*n
69b55c55 4132 //
75bd7f81 4133
4551ea7c 4134 Int_t *sindexS = new Int_t[n]; // Temporary array for sorting
4135 Int_t *sindexF = new Int_t[2*n];
4136 for (Int_t i = 0; i < n; i++) {
4137 sindexF[i] = 0;
4138 }
4139
4140 TMath::Sort(n,inlist,sindexS,down);
4141
4142 Int_t last = inlist[sindexS[0]];
4143 Int_t val = last;
4144 sindexF[0] = 1;
4145 sindexF[0+n] = last;
4146 Int_t countPos = 0;
4147
4148 // Find frequency
4149 for (Int_t i = 1; i < n; i++) {
69b55c55 4150 val = inlist[sindexS[i]];
4551ea7c 4151 if (last == val) {
4152 sindexF[countPos]++;
4153 }
4154 else {
69b55c55 4155 countPos++;
4156 sindexF[countPos+n] = val;
4157 sindexF[countPos]++;
4551ea7c 4158 last = val;
69b55c55 4159 }
4160 }
4551ea7c 4161 if (last == val) {
4162 countPos++;
4163 }
4164
4165 // Sort according frequency
4166 TMath::Sort(countPos,sindexF,sindexS,kTRUE);
4167
4168 for (Int_t i = 0; i < countPos; i++) {
69b55c55 4169 outlist[2*i ] = sindexF[sindexS[i]+n];
4170 outlist[2*i+1] = sindexF[sindexS[i]];
4171 }
4551ea7c 4172
69b55c55 4173 delete [] sindexS;
4174 delete [] sindexF;
4175
4176 return countPos;
75bd7f81 4177
69b55c55 4178}
4179
75bd7f81 4180//_____________________________________________________________________________
4551ea7c 4181AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params)
69b55c55 4182{
4183 //
75bd7f81 4184 // Register a seed
69b55c55 4185 //
4551ea7c 4186
4187 Double_t alpha = AliTRDgeometry::GetAlpha();
4188 Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
69b55c55 4189 Double_t c[15];
4551ea7c 4190
4191 c[ 0] = 0.2;
4192 c[ 1] = 0.0; c[ 2] = 2.0;
4193 c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
4194 c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0; c[ 9] = 0.1;
4195 c[10] = 0.0; c[11] = 0.0; c[12] = 0.0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
4196
4197 Int_t index = 0;
4198 AliTRDcluster *cl = 0;
4199
4200 for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
4201 if (seeds[ilayer].IsOK()) {
4202 for (Int_t itime = 22; itime > 0; itime--) {
ca1e1e5b 4203 if (seeds[ilayer].GetIndexes(itime) > 0) {
4204 index = seeds[ilayer].GetIndexes(itime);
4205 cl = seeds[ilayer].GetClusters(itime);
69b55c55 4206 break;
4207 }
4208 }
4209 }
4551ea7c 4210 if (index > 0) {
4211 break;
4212 }
69b55c55 4213 }
4551ea7c 4214 if (cl == 0) {
4215 return 0;
4216 }
4217
4218 AliTRDtrack *track = new AliTRDtrack(cl
4219 ,index
4220 ,&params[1]
4221 ,c
4222 ,params[0]
4223 ,params[6]*alpha+shift);
43bfc8af 4224 // SetCluster(cl, 0); // A. Bercuci
4225 track->PropagateTo(params[0]-5.0);
69b55c55 4226 track->ResetCovariance(1);
4551ea7c 4227
4228 Int_t rc = FollowBackProlongation(*track);
4229 if (rc < 30) {
69b55c55 4230 delete track;
4551ea7c 4231 track = 0;
4232 }
4233 else {
69b55c55 4234 track->CookdEdx();
43bfc8af 4235 track->CookdEdxTimBin();
4551ea7c 4236 CookLabel(track,0.9);
69b55c55 4237 }
69b55c55 4238
75bd7f81 4239 return track;
7e448bcc 4240}
4241
4242//////////////////////////////////////////////////////////////////////////////////////////
69b55c55 4243
7e448bcc 4244void AliTRDtracker::InitLogHists() {
4245
4246 fHBackfit = new TH1D("logTRD_backfit", "", 40, -0.5, 39.5);
4247 fHRefit = new TH1D("logTRD_refit", "", 40, -0.5, 39.5);
4248 fHClSearch = new TH1D("logTRD_clSearch", "", 60, -0.5, 59.5);
4249
4250 fHX = new TH1D("logTRD_X", ";x (cm)", 200, 50, 400);
4251 fHNCl = new TH1D("logTRD_ncl", "", 40, -0.5, 39.5);
4252 fHNClTrack = new TH1D("logTRD_nclTrack", "", 180, -0.5, 179.5);
4253
4254 fHMinYPos = new TH1D("logTRD_minYPos", ";#delta Y (cm)", 400, -6, 6);
4255 fHMinYNeg = new TH1D("logTRD_minYNeg", ";#delta Y (cm)", 400, -6, 6);
4256 fHMinZ = new TH1D("logTRD_minZ", ";#delta Z (cm)", 400, -20, 20);
4257 fHMinD = new TH2D("logTRD_minD", ";#delta Y (cm);#delta Z (cm)", 100, -6, 6, 100, -50, 50);
4258
4259 fHDeltaX = new TH1D("logTRD_deltaX", ";#delta X (cm)", 100, -5, 5);
4260 fHXCl = new TH1D("logTRD_xCl", ";cluster x position (cm)", 1000, 280, 380);
4261
4262 const char *nameFindCl[4] = {"logTRD_clY", "logTRD_clZ", "logTRD_clB", "logTRD_clG"};
4263
4264 for(int i=0; i<4; i++) {
4265 fHFindCl[i] = new TH1D(nameFindCl[i], "", 30, -0.5, 29.5);
4266 }
69b55c55 4267}
7e448bcc 4268
4269//////////////////////////////////////////////////////////////////////////////////////////
4270
4271void AliTRDtracker::SaveLogHists() {
4272
4273 TDirectory *sav = gDirectory;
4274 TFile *logFile = 0;
4275
4276 TSeqCollection *col = gROOT->GetListOfFiles();
4277 int N = col->GetEntries();
4278 for(int i=0; i<N; i++) {
4279 logFile = (TFile*)col->At(i);
4280 if (strstr(logFile->GetName(), "AliESDs.root")) break;
4281 }
4282
4283 logFile->cd();
4284 fHBackfit->Write(fHBackfit->GetName(), TObject::kOverwrite);
4285 fHRefit->Write(fHRefit->GetName(), TObject::kOverwrite);
4286 fHClSearch->Write(fHClSearch->GetName(), TObject::kOverwrite);
4287 fHX->Write(fHX->GetName(), TObject::kOverwrite);
4288 fHNCl->Write(fHNCl->GetName(), TObject::kOverwrite);
4289 fHNClTrack->Write(fHNClTrack->GetName(), TObject::kOverwrite);
4290
4291 fHMinYPos->Write(fHMinYPos->GetName(), TObject::kOverwrite);
4292 fHMinYNeg->Write(fHMinYNeg->GetName(), TObject::kOverwrite);
4293 fHMinD->Write(fHMinD->GetName(), TObject::kOverwrite);
4294 fHMinZ->Write(fHMinZ->GetName(), TObject::kOverwrite);
4295
4296 fHDeltaX->Write(fHDeltaX->GetName(), TObject::kOverwrite);
4297 fHXCl->Write(fHXCl->GetName(), TObject::kOverwrite);
4298
4299
4300 for(int i=0; i<4; i++)
4301 fHFindCl[i]->Write(fHFindCl[i]->GetName(), TObject::kOverwrite);
4302
4303 logFile->Flush();
4304
4305 sav->cd();
4306}
4307
4308//////////////////////////////////////////////////////////////////////////////////////////