]> 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
kZDiElectron forcing Z->e+e- added.
[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) ||
6c23ffed 709 (track->Pt() > 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) &&
6c23ffed 729 (track->Pt() < 1.5)) {
4551ea7c 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;
306e3dc1 981 //AliTRDtrack &t = *pt;
43bfc8af 982 FollowProlongation(*pt);
983 if (pt->GetNumberOfClusters() >= foundMin) {
4551ea7c 984 //UseClusters(&t);
6c94f330 985 //CookLabel(pt, 1-fgkLabelFraction);
43bfc8af 986 pt->CookdEdx();
987 pt->CookdEdxTimBin();
306e3dc1 988 pt->CookPID(seed);
989 //pt->Calibrate(); // slot for calibration
1e9bb598 990 }
991 found++;
33744e5d 992
4551ea7c 993 Double_t xTPC = 250.0;
43bfc8af 994 if (PropagateToX(*pt,xTPC,fgkMaxStep)) {
7e448bcc 995
4551ea7c 996 seed->UpdateTrackParams(pt,AliESDtrack::kTRDrefit);
7e448bcc 997 fHRefit->Fill(5);
998
4551ea7c 999 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
1000 for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
6d45eaef 1001 seed->SetTRDsignals(pt->GetPIDsignals(i,j),i,j);
1002 }
7ad19338 1003 seed->SetTRDTimBin(pt->GetPIDTimBin(i),i);
1004 }
4551ea7c 1005 }
1006 else {
1007 // If not prolongation to TPC - propagate without update
7e448bcc 1008 fHRefit->Fill(5);
4551ea7c 1009 AliTRDtrack *seed2 = new AliTRDtrack(*seed);
1010 seed2->ResetCovariance(5.0);
1011 AliTRDtrack *pt2 = new AliTRDtrack(*seed2,seed2->GetAlpha());
7bed16a7 1012 delete seed2;
59393e34 1013 if (PropagateToX(*pt2,xTPC,fgkMaxStep)) {
4551ea7c 1014 pt2->CookdEdx( );
306e3dc1 1015 pt2->CookdEdxTimBin();
1016 seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit);
1017 fHRefit->Fill(6);
7e448bcc 1018
4551ea7c 1019 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
1020 for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
6d45eaef 1021 seed->SetTRDsignals(pt2->GetPIDsignals(i,j),i,j);
43bfc8af 1022 }
7ad19338 1023 seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
1024 }
eab5961e 1025 }
306e3dc1 1026
1027 // Add TRD track to ESDfriendTrack - maybe this tracks are not useful for post-processing - TODO make decission
1028 if (AliTRDReconstructor::StreamLevel() > 0) {
1029 seed->AddCalibObject(new AliTRDtrack(*pt2/*, kTRUE*/));
1030 }
7bed16a7 1031 delete pt2;
4551ea7c 1032 }
1033
306e3dc1 1034 // Add TRD track to ESDfriendTrack
1035 if (AliTRDReconstructor::StreamLevel() > 0) {
1036 seed->AddCalibObject(new AliTRDtrack(*pt/*, kTRUE*/));
1037 }
1e9bb598 1038 delete pt;
4551ea7c 1039
eab5961e 1040 }
1e9bb598 1041
33744e5d 1042 AliInfo(Form("Number of loaded seeds: %d",nseed));
1043 AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
1e9bb598 1044
7e448bcc 1045 SaveLogHists();
c630aafd 1046 return 0;
c630aafd 1047}
1048
75bd7f81 1049//_____________________________________________________________________________
4551ea7c 1050Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t)
8979685e 1051{
43bfc8af 1052 //
1053 // Starting from current position on track=t this function tries
1054 // to extrapolate the track up to timeBin=0 and to confirm prolongation
1055 // if a close cluster is found. Returns the number of clusters
1056 // expected to be found in sensitive layers
1057 // GeoManager used to estimate mean density
1058 //
4551ea7c 1059
43bfc8af 1060 Int_t sector;
1061 Int_t lastplane = GetLastPlane(&t);
97ccb02b 1062 Double_t xx0 = 0.0;
1063 Double_t xrho= 0.0;
43bfc8af 1064 Int_t expectedNumberOfClusters = 0;
4551ea7c 1065
43bfc8af 1066 for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
4551ea7c 1067
43bfc8af 1068 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
1069 Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
4551ea7c 1070
43bfc8af 1071 //
1072 // Propagate track close to the plane if neccessary
1073 //
1074 Double_t currentx = fTrSec[0]->GetLayer(rowlast)->GetX();
1075 if (currentx < (-fgkMaxStep + t.GetX())) {
1076 // Propagate closer to chamber - safety space fgkMaxStep
1077 if (!PropagateToX(t,currentx+fgkMaxStep,fgkMaxStep)) {
1078 break;
1079 }
1080 }
4551ea7c 1081
43bfc8af 1082 if (!AdjustSector(&t)) {
1083 break;
1084 }
4551ea7c 1085
43bfc8af 1086 //
1087 // Get material budget
1088 //
1089 Double_t xyz0[3];
1090 Double_t xyz1[3];
1091 Double_t param[7];
1092 Double_t x;
1093 Double_t y;
1094 Double_t z;
1095
1096 // Starting global position
1097 t.GetXYZ(xyz0);
1098 // End global position
1099 x = fTrSec[0]->GetLayer(row0)->GetX();
1100 if (!t.GetProlongation(x,y,z)) {
1101 break;
1102 }
1103 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1104 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1105 xyz1[2] = z;
97ccb02b 1106 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1107 xrho= param[0]*param[4];
1108 xx0 = param[1]; // Get mean propagation parameters
43bfc8af 1109 // A.Bercuci 25.07.07
1110 // Flags for marking the track momentum and direction. The track is
1111 // marked only if it has at least 1 cluster picked up in the current
1112 // chamber.
1113 Bool_t kUPDATED = kFALSE, kMARKED = kFALSE;
4551ea7c 1114
43bfc8af 1115 //
1116 // Propagate and update
1117 //
1118 sector = t.GetSector();
1119 //for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
1120 for (Int_t itime = 0 ; itime < GetTimeBinsPerPlane(); itime++) {
1121 // A.Bercuci 25.07.07
1122 // Mark track kinematics
1123 if(itime > 10 && kUPDATED && !kMARKED){
1124 t.SetTrackSegmentDirMom(iplane);
1125 kMARKED = kTRUE;
1126 }
1127
1128 Int_t ilayer = GetGlobalTimeBin(0,iplane,itime);
1129 expectedNumberOfClusters++;
1130 t.SetNExpected(t.GetNExpected() + 1);
1131 if (t.GetX() > 345.0) {
1132 t.SetNExpectedLast(t.GetNExpectedLast() + 1);
1133 }
1134 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
1135 AliTRDcluster *cl = 0;
1136 UInt_t index = 0;
1137 Double_t maxChi2 = fgkMaxChi2;
1138 x = timeBin.GetX();
1139
1140 if (timeBin) {
1141
1142 AliTRDcluster *cl0 = timeBin[0];
1143 if (!cl0) {
1144 // No clusters in given time bin
1145 continue;
1146 }
1147
1148 Int_t plane = fGeom->GetPlane(cl0->GetDetector());
1149 if (plane > lastplane) {
1150 continue;
1151 }
1152
1153 Int_t timebin = cl0->GetLocalTimeBin();
1154 AliTRDcluster *cl2 = GetCluster(&t,plane,timebin,index);
1155
1156 if (cl2) {
1157 cl = cl2;
1158 //Double_t h01 = GetTiltFactor(cl); //I.B's fix
1159 //maxChi2=t.GetPredictedChi2(cl,h01);
1160 }
1161 if (cl) {
1162 //if (cl->GetNPads()<5)
1163 Double_t dxsample = timeBin.GetdX();
1164 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
1165 Double_t h01 = GetTiltFactor(cl);
1166 Int_t det = cl->GetDetector();
1167 Int_t plane = fGeom->GetPlane(det);
1168 if (t.GetX() > 345.0) {
1169 t.SetNLast(t.GetNLast() + 1);
1170 t.SetChi2Last(t.GetChi2Last() + maxChi2);
1171 }
1172
1173 Double_t xcluster = cl->GetX();
97ccb02b 1174 t.PropagateTo(xcluster,xx0,xrho);
43bfc8af 1175
1176 if (!AdjustSector(&t)) {
1177 break; //I.B's fix
1178 }
1179 maxChi2 = t.GetPredictedChi2(cl,h01);
1180
1181 if (maxChi2<1e+10)
1182 if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1183 // ????
1184 } else {
1185 // A.Bercuci 25.07.07
1186 //SetCluster(cl, GetNumberOfClusters()-1);
1187 kUPDATED = kTRUE;
1188 }
1189 }
1190 }
1191 }
4551ea7c 1192 }
1193
43bfc8af 1194 return expectedNumberOfClusters;
75bd7f81 1195}
7b580082 1196
75bd7f81 1197//_____________________________________________________________________________
4551ea7c 1198Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
8979685e 1199{
75bd7f81 1200 //
8979685e 1201 // Starting from current radial position of track <t> this function
1202 // extrapolates the track up to outer timebin and in the sensitive
1203 // layers confirms prolongation if a close cluster is found.
1204 // Returns the number of clusters expected to be found in sensitive layers
1205 // Use GEO manager for material Description
75bd7f81 1206 //
7e448bcc 1207 // return number of assigned clusters ?
1208 //
1209
59393e34 1210
4551ea7c 1211 Int_t sector;
1212 Int_t clusters[1000];
97ccb02b 1213 Double_t xx0 = 0.0;
1214 Double_t xrho = 0.0;
4551ea7c 1215 Int_t expectedNumberOfClusters = 0;
1216 Float_t ratio0 = 0.0;
8979685e 1217 AliTRDtracklet tracklet;
33744e5d 1218
77566f2a 1219 // Calibration fill 2D
082a2283 1220 AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
77566f2a 1221 if (!calibra) {
1222 AliInfo("Could not get Calibra instance\n");
1223 }
8ec526a4 1224 if (calibra->GetMITracking()) {
1225 calibra->ResetTrack();
77566f2a 1226 }
1227
33744e5d 1228 for (Int_t i = 0; i < 1000; i++) {
1229 clusters[i] = -1;
1230 }
1231
4551ea7c 1232 for (Int_t iplane = 0; iplane < AliESDtrack::kNPlane; iplane++) {
1233
7e448bcc 1234 int hb = iplane * 10;
1235 fHClSearch->Fill(hb);
1236
4551ea7c 1237 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
1238 Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
1239 Double_t currentx = fTrSec[0]->GetLayer(row0)->GetX();
1240 if (currentx < t.GetX()) {
7e448bcc 1241 fHClSearch->Fill(hb+1);
4551ea7c 1242 continue;
1243 }
1244
6c94f330 1245 //
4551ea7c 1246 // Propagate closer to chamber if neccessary
6c94f330 1247 //
4551ea7c 1248 if (currentx > (fgkMaxStep + t.GetX())) {
1249 if (!PropagateToX(t,currentx-fgkMaxStep,fgkMaxStep)) {
7e448bcc 1250 fHClSearch->Fill(hb+2);
4551ea7c 1251 break;
1252 }
1253 }
1254 if (!AdjustSector(&t)) {
7e448bcc 1255 fHClSearch->Fill(hb+3);
4551ea7c 1256 break;
8979685e 1257 }
4551ea7c 1258 if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) {
7e448bcc 1259 fHClSearch->Fill(hb+4);
4551ea7c 1260 break;
1261 }
1262
8979685e 1263 //
4551ea7c 1264 // Get material budget inside of chamber
59393e34 1265 //
4551ea7c 1266 Double_t xyz0[3];
1267 Double_t xyz1[3];
1268 Double_t param[7];
1269 Double_t x;
1270 Double_t y;
1271 Double_t z;
1272 // Starting global position
1273 t.GetXYZ(xyz0);
1274 // End global position
7b580082 1275 x = fTrSec[0]->GetLayer(rowlast)->GetX();
4551ea7c 1276 if (!t.GetProlongation(x,y,z)) {
7e448bcc 1277 fHClSearch->Fill(hb+5);
4551ea7c 1278 break;
1279 }
1280 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1281 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1282 xyz1[2] = z;
97ccb02b 1283 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1284 xrho = param[0]*param[4];
1285 xx0 = param[1]; // Get mean propagation parameters
4551ea7c 1286
8979685e 1287 //
7b580082 1288 // Find clusters
8979685e 1289 //
6c94f330 1290 sector = t.GetSector();
4551ea7c 1291 Float_t ncl = FindClusters(sector,row0,rowlast,&t,clusters,tracklet);
7e448bcc 1292 fHNCl->Fill(tracklet.GetN());
1293
4551ea7c 1294 if (tracklet.GetN() < GetTimeBinsPerPlane()/3) {
7e448bcc 1295 fHClSearch->Fill(hb+6);
4551ea7c 1296 continue;
1297 }
1298
8979685e 1299 //
7b580082 1300 // Propagate and update track
8979685e 1301 //
4551ea7c 1302 for (Int_t itime = GetTimeBinsPerPlane()-1; itime >= 0; itime--) {
1303
1304 Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
8979685e 1305 expectedNumberOfClusters++;
27eaf44b 1306 t.SetNExpected(t.GetNExpected() + 1);
4551ea7c 1307 if (t.GetX() > 345.0) {
27eaf44b 1308 t.SetNExpectedLast(t.GetNExpectedLast() + 1);
4551ea7c 1309 }
1310 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
7e448bcc 1311 AliTRDcluster *cl = 0;
4551ea7c 1312 UInt_t index = 0;
1313 Double_t maxChi2 = fgkMaxChi2;
8979685e 1314 x = timeBin.GetX();
4551ea7c 1315
8979685e 1316 if (timeBin) {
4551ea7c 1317
1318 if (clusters[ilayer] > 0) {
8979685e 1319 index = clusters[ilayer];
4551ea7c 1320 cl = (AliTRDcluster *)GetCluster(index);
33744e5d 1321 //Double_t h01 = GetTiltFactor(cl); // I.B's fix
1322 //maxChi2=t.GetPredictedChi2(cl,h01); //
8979685e 1323 }
1324
1325 if (cl) {
4551ea7c 1326
1327 //if (cl->GetNPads() < 5)
59393e34 1328 Double_t dxsample = timeBin.GetdX();
1329 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
4551ea7c 1330 Double_t h01 = GetTiltFactor(cl);
1331 Int_t det = cl->GetDetector();
1332 Int_t plane = fGeom->GetPlane(det);
1333 if (t.GetX() > 345.0) {
27eaf44b 1334 t.SetNLast(t.GetNLast() + 1);
1335 t.SetChi2Last(t.GetChi2Last() + maxChi2);
8979685e 1336 }
59393e34 1337 Double_t xcluster = cl->GetX();
97ccb02b 1338 t.PropagateTo(xcluster,xx0,xrho);
4551ea7c 1339 maxChi2 = t.GetPredictedChi2(cl,h01);
6c94f330 1340
c2239ca8 1341 if (maxChi2<1e+10)
1342 if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1343 if (!t.Update(cl,maxChi2,index,h01)) {
1344 // ????
1345 }
43bfc8af 1346 } // else SetCluster(cl, GetNumberOfClusters()-1); // A.Bercuci 25.07.07
1347
4551ea7c 1348
8ec526a4 1349 if (calibra->GetMITracking()) {
77566f2a 1350 calibra->UpdateHistograms(cl,&t);
1351 }
1352
4551ea7c 1353 // Reset material budget if 2 consecutive gold
1354 if (plane > 0) {
27eaf44b 1355 if ((t.GetTracklets(plane).GetN() + t.GetTracklets(plane-1).GetN()) > 20) {
1356 t.SetBudget(2,0.0);
4551ea7c 1357 }
1358 }
1359
1360 }
1361
8979685e 1362 }
4551ea7c 1363
59393e34 1364 }
4551ea7c 1365
1366 ratio0 = ncl / Float_t(fTimeBinsPerPlane);
27eaf44b 1367 Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);
1368 if ((tracklet.GetChi2() < 18.0) &&
1369 (ratio0 > 0.8) &&
1370 (ratio1 > 0.6) &&
1371 (ratio0+ratio1 > 1.5) &&
1372 (t.GetNCross() == 0) &&
1373 (TMath::Abs(t.GetSnp()) < 0.85) &&
1374 (t.GetNumberOfClusters() > 20)){
7e448bcc 1375 //if (ratio0 > 0.8) {
4551ea7c 1376 t.MakeBackupTrack(); // Make backup of the track until is gold
59393e34 1377 }
7b580082 1378
8979685e 1379 }
5443e65e 1380
75bd7f81 1381 return expectedNumberOfClusters;
75bd7f81 1382}
1e9bb598 1383
75bd7f81 1384//_____________________________________________________________________________
4551ea7c 1385Int_t AliTRDtracker::PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep)
5443e65e 1386{
75bd7f81 1387 //
5443e65e 1388 // Starting from current radial position of track <t> this function
1389 // extrapolates the track up to radial position <xToGo>.
1390 // Returns 1 if track reaches the plane, and 0 otherwise
75bd7f81 1391 //
1392
59393e34 1393 const Double_t kEpsilon = 0.00001;
4551ea7c 1394 //Double_t tanmax = TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
1395 Double_t xpos = t.GetX();
1396 Double_t dir = (xpos<xToGo) ? 1.0 : -1.0;
1397
1398 while (((xToGo-xpos)*dir) > kEpsilon) {
1399
1400 Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos),maxStep);
1401
1402 Double_t xyz0[3];
1403 Double_t xyz1[3];
1404 Double_t param[7];
1405 Double_t x;
1406 Double_t y;
1407 Double_t z;
1408 // Starting global position
1409 t.GetXYZ(xyz0);
1410 x = xpos + step;
1411
1412 if (!t.GetProlongation(x,y,z)) {
1413 return 0; // No prolongation
1414 }
1415
1416 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1417 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1418 xyz1[2] = z;
1419
97ccb02b 1420 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1421 if (!t.PropagateTo(x,param[1],param[0]*param[4])) {
4551ea7c 1422 return 0;
1423 }
59393e34 1424 AdjustSector(&t);
1425 xpos = t.GetX();
4551ea7c 1426
5443e65e 1427 }
75bd7f81 1428
5443e65e 1429 return 1;
5443e65e 1430
59393e34 1431}
5443e65e 1432
c630aafd 1433//_____________________________________________________________________________
1434Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1435{
75bd7f81 1436 //
c630aafd 1437 // Fills clusters into TRD tracking_sectors
1438 // Note that the numbering scheme for the TRD tracking_sectors
1439 // differs from that of TRD sectors
75bd7f81 1440 //
1441
c630aafd 1442 if (ReadClusters(fClusters,cTree)) {
4551ea7c 1443 AliError("Problem with reading the clusters !");
c630aafd 1444 return 1;
1445 }
4551ea7c 1446 Int_t ncl = fClusters->GetEntriesFast();
1447 fNclusters = ncl;
1448 AliInfo(Form("Sorting %d clusters",ncl));
c630aafd 1449
1450 UInt_t index;
4551ea7c 1451 for (Int_t ichamber = 0; ichamber < 5; ichamber++) {
1452 for (Int_t isector = 0; isector < 18; isector++) {
1453 fHoles[ichamber][isector] = kTRUE;
3c625a9b 1454 }
4551ea7c 1455 }
ad670fba 1456
6c94f330 1457 while (ncl--) {
33744e5d 1458
4551ea7c 1459 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(ncl);
1460 Int_t detector = c->GetDetector();
1461 Int_t localTimeBin = c->GetLocalTimeBin();
1462 Int_t sector = fGeom->GetSector(detector);
1463 Int_t plane = fGeom->GetPlane(detector);
439c63c8 1464 Int_t trackingSector = sector;
4551ea7c 1465
7e448bcc 1466 //if (c->GetLabel(0) > 0) {
1467 if (c->GetQ() > 10) {
3c625a9b 1468 Int_t chamber = fGeom->GetChamber(detector);
4551ea7c 1469 fHoles[chamber][trackingSector] = kFALSE;
3c625a9b 1470 }
c630aafd 1471
6c94f330 1472 Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
4551ea7c 1473 if (gtb < 0) {
1474 continue;
1475 }
029cd327 1476 Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
c630aafd 1477
4551ea7c 1478 index = ncl;
33744e5d 1479
1480 // Apply pos correction
7e448bcc 1481 Transform(c);
1482 fHXCl->Fill(c->GetX());
33744e5d 1483
7e448bcc 1484 fTrSec[trackingSector]->GetLayer(layer)->SetX(c->GetX());
1485 fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
4551ea7c 1486 }
75bd7f81 1487
c630aafd 1488 return 0;
75bd7f81 1489
c630aafd 1490}
1491
5443e65e 1492//_____________________________________________________________________________
b7a0917f 1493void AliTRDtracker::UnloadClusters()
5443e65e 1494{
1495 //
1496 // Clears the arrays of clusters and tracks. Resets sectors and timebins
1497 //
1498
4551ea7c 1499 Int_t i;
1500 Int_t nentr;
5443e65e 1501
1502 nentr = fClusters->GetEntriesFast();
4551ea7c 1503 for (i = 0; i < nentr; i++) {
1504 delete fClusters->RemoveAt(i);
1505 }
b7a0917f 1506 fNclusters = 0;
5443e65e 1507
1508 nentr = fSeeds->GetEntriesFast();
4551ea7c 1509 for (i = 0; i < nentr; i++) {
1510 delete fSeeds->RemoveAt(i);
1511 }
5443e65e 1512
1513 nentr = fTracks->GetEntriesFast();
4551ea7c 1514 for (i = 0; i < nentr; i++) {
1515 delete fTracks->RemoveAt(i);
1516 }
5443e65e 1517
1518 Int_t nsec = AliTRDgeometry::kNsect;
5443e65e 1519 for (i = 0; i < nsec; i++) {
1520 for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1521 fTrSec[i]->GetLayer(pl)->Clear();
1522 }
1523 }
1524
1525}
1526
75bd7f81 1527//_____________________________________________________________________________
af885e0f 1528void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESDEvent *esd)
7ad19338 1529{
1530 //
1531 // Creates seeds using clusters between position inner plane and outer plane
1532 //
75bd7f81 1533
4551ea7c 1534 const Double_t kMaxTheta = 1.0;
1535 const Double_t kMaxPhi = 2.0;
1536
1537 const Double_t kRoad0y = 6.0; // Road for middle cluster
1538 const Double_t kRoad0z = 8.5; // Road for middle cluster
1539
1540 const Double_t kRoad1y = 2.0; // Road in y for seeded cluster
1541 const Double_t kRoad1z = 20.0; // Road in z for seeded cluster
1542
1543 const Double_t kRoad2y = 3.0; // Road in y for extrapolated cluster
1544 const Double_t kRoad2z = 20.0; // Road in z for extrapolated cluster
c6f438c0 1545 const Int_t kMaxSeed = 3000;
7ad19338 1546
4551ea7c 1547 Int_t maxSec = AliTRDgeometry::kNsect;
1548
1549 // Linear fitters in planes
1550 TLinearFitter fitterTC(2,"hyp2"); // Fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
1551 TLinearFitter fitterT2(4,"hyp4"); // Fitting with tilting pads - kz not fixed
69b55c55 1552 fitterTC.StoreData(kTRUE);
1553 fitterT2.StoreData(kTRUE);
4551ea7c 1554 AliRieman rieman(1000); // Rieman fitter
1555 AliRieman rieman2(1000); // Rieman fitter
1556
1557 // Find the maximal and minimal layer for the planes
7ad19338 1558 Int_t layers[6][2];
4551ea7c 1559 AliTRDpropagationLayer *reflayers[6];
1560 for (Int_t i = 0; i < 6; i++) {
1561 layers[i][0] = 10000;
1562 layers[i][1] = 0;
1563 }
1564 for (Int_t ns = 0; ns < maxSec; ns++) {
1565 for (Int_t ilayer = 0; ilayer < fTrSec[ns]->GetNumberOfLayers(); ilayer++) {
1566 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(ilayer));
1567 if (layer == 0) {
1568 continue;
1569 }
7ad19338 1570 Int_t det = layer[0]->GetDetector();
1571 Int_t plane = fGeom->GetPlane(det);
4551ea7c 1572 if (ilayer < layers[plane][0]) {
1573 layers[plane][0] = ilayer;
1574 }
1575 if (ilayer > layers[plane][1]) {
1576 layers[plane][1] = ilayer;
1577 }
7ad19338 1578 }
1579 }
4551ea7c 1580
f162af62 1581 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0);
6c94f330 1582 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
4551ea7c 1583 Double_t hL[6]; // Tilting angle
1584 Double_t xcl[6]; // X - position of reference cluster
1585 Double_t ycl[6]; // Y - position of reference cluster
1586 Double_t zcl[6]; // Z - position of reference cluster
1587
1588 AliTRDcluster *cl[6] = { 0, 0, 0, 0, 0, 0 }; // Seeding clusters
1589 Float_t padlength[6] = { 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 }; // Current pad-length
1590
1591 Double_t chi2R = 0.0;
1592 Double_t chi2Z = 0.0;
1593 Double_t chi2RF = 0.0;
1594 Double_t chi2ZF = 0.0;
1595
1596 Int_t nclusters; // Total number of clusters
1597 for (Int_t i = 0; i < 6; i++) {
1598 hL[i] = h01;
1599 if (i%2==1) {
1600 hL[i]*=-1.0;
1601 }
1602 }
1603
1604 // Registered seed
c6f438c0 1605 AliTRDseed *pseed = new AliTRDseed[kMaxSeed*6];
1606 AliTRDseed *seed[kMaxSeed];
4551ea7c 1607 for (Int_t iseed = 0; iseed < kMaxSeed; iseed++) {
1608 seed[iseed]= &pseed[iseed*6];
1609 }
69b55c55 1610 AliTRDseed *cseed = seed[0];
4551ea7c 1611
1612 Double_t seedquality[kMaxSeed];
1613 Double_t seedquality2[kMaxSeed];
1614 Double_t seedparams[kMaxSeed][7];
1615 Int_t seedlayer[kMaxSeed];
1616 Int_t registered = 0;
1617 Int_t sort[kMaxSeed];
1618
1619 //
1620 // Seeding part
1621 //
1622 for (Int_t ns = 0; ns < maxSec; ns++) { // Loop over sectors
1623 //for (Int_t ns = 0; ns < 5; ns++) { // Loop over sectors
1624
1625 registered = 0; // Reset registerd seed counter
1626 cseed = seed[registered];
1627 Float_t iter = 0.0;
1628
1629 for (Int_t sLayer = 2; sLayer >= 0; sLayer--) {
1630 //for (Int_t dseed = 5; dseed < 15; dseed += 3) {
1631
1632 iter += 1.0;
1633 Int_t dseed = 5 + Int_t(iter) * 3;
1634
69b55c55 1635 // Initialize seeding layers
4551ea7c 1636 for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
69b55c55 1637 reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
1638 xcl[ilayer] = reflayers[ilayer]->GetX();
4551ea7c 1639 }
1640
1641 Double_t xref = (xcl[sLayer+1] + xcl[sLayer+2]) * 0.5;
1642 AliTRDpropagationLayer &layer0 = *reflayers[sLayer+0];
1643 AliTRDpropagationLayer &layer1 = *reflayers[sLayer+1];
1644 AliTRDpropagationLayer &layer2 = *reflayers[sLayer+2];
1645 AliTRDpropagationLayer &layer3 = *reflayers[sLayer+3];
1646
1647 Int_t maxn3 = layer3;
1648 for (Int_t icl3 = 0; icl3 < maxn3; icl3++) {
33744e5d 1649
69b55c55 1650 AliTRDcluster *cl3 = layer3[icl3];
4551ea7c 1651 if (!cl3) {
1652 continue;
1653 }
1654 padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2() * 12.0);
69b55c55 1655 ycl[sLayer+3] = cl3->GetY();
1656 zcl[sLayer+3] = cl3->GetZ();
4551ea7c 1657 Float_t yymin0 = ycl[sLayer+3] - 1.0 - kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1658 Float_t yymax0 = ycl[sLayer+3] + 1.0 + kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1659 Int_t maxn0 = layer0;
1660
1661 for (Int_t icl0 = layer0.Find(yymin0); icl0 < maxn0; icl0++) {
1662
69b55c55 1663 AliTRDcluster *cl0 = layer0[icl0];
4551ea7c 1664 if (!cl0) {
1665 continue;
1666 }
1667 if (cl3->IsUsed() && cl0->IsUsed()) {
1668 continue;
1669 }
69b55c55 1670 ycl[sLayer+0] = cl0->GetY();
1671 zcl[sLayer+0] = cl0->GetZ();
4551ea7c 1672 if (ycl[sLayer+0] > yymax0) {
1673 break;
1674 }
1675 Double_t tanphi = (ycl[sLayer+3]-ycl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1676 if (TMath::Abs(tanphi) > kMaxPhi) {
1677 continue;
1678 }
1679 Double_t tantheta = (zcl[sLayer+3]-zcl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1680 if (TMath::Abs(tantheta) > kMaxTheta) {
1681 continue;
1682 }
1683 padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2() * 12.0);
1684
1685 // Expected position in 1 layer
1686 Double_t y1exp = ycl[sLayer+0] + (tanphi) * (xcl[sLayer+1]-xcl[sLayer+0]);
1687 Double_t z1exp = zcl[sLayer+0] + (tantheta) * (xcl[sLayer+1]-xcl[sLayer+0]);
1688 Float_t yymin1 = y1exp - kRoad0y - tanphi;
1689 Float_t yymax1 = y1exp + kRoad0y + tanphi;
1690 Int_t maxn1 = layer1;
1691
1692 for (Int_t icl1 = layer1.Find(yymin1); icl1 < maxn1; icl1++) {
1693
69b55c55 1694 AliTRDcluster *cl1 = layer1[icl1];
4551ea7c 1695 if (!cl1) {
1696 continue;
1697 }
69b55c55 1698 Int_t nusedCl = 0;
1699 if (cl3->IsUsed()) nusedCl++;
1700 if (cl0->IsUsed()) nusedCl++;
1701 if (cl1->IsUsed()) nusedCl++;
4551ea7c 1702 if (nusedCl > 1) {
1703 continue;
1704 }
69b55c55 1705 ycl[sLayer+1] = cl1->GetY();
1706 zcl[sLayer+1] = cl1->GetZ();
4551ea7c 1707 if (ycl[sLayer+1] > yymax1) {
1708 break;
1709 }
1710 if (TMath::Abs(ycl[sLayer+1]-y1exp) > kRoad0y+tanphi) {
1711 continue;
1712 }
1713 if (TMath::Abs(zcl[sLayer+1]-z1exp) > kRoad0z) {
1714 continue;
1715 }
1716 padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2() * 12.0);
1717
1718 Double_t y2exp = ycl[sLayer+0]+(tanphi) * (xcl[sLayer+2]-xcl[sLayer+0]) + (ycl[sLayer+1]-y1exp);
1719 Double_t z2exp = zcl[sLayer+0]+(tantheta) * (xcl[sLayer+2]-xcl[sLayer+0]);
1720 Int_t index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y,kRoad1z);
1721 if (index2 <= 0) {
1722 continue;
1723 }
1724 AliTRDcluster *cl2 = (AliTRDcluster *) GetCluster(index2);
1725 padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2() * 12.0);
1726 ycl[sLayer+2] = cl2->GetY();
1727 zcl[sLayer+2] = cl2->GetZ();
1728 if (TMath::Abs(cl2->GetZ()-z2exp) > kRoad0z) {
1729 continue;
1730 }
1731
69b55c55 1732 rieman.Reset();
1733 rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
1734 rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
1735 rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);
1736 rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
1737 rieman.Update();
4551ea7c 1738
1739 // Reset fitter
1740 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
69b55c55 1741 cseed[iLayer].Reset();
1742 }
4551ea7c 1743 chi2Z = 0.0;
1744 chi2R = 0.0;
1745
1746 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
ca1e1e5b 1747 cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1748 chi2Z += (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer])
1749 * (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer]);
1750 cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
1751 cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1752 chi2R += (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer])
1753 * (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer]);
1754 cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
69b55c55 1755 }
4551ea7c 1756 if (TMath::Sqrt(chi2R) > 1.0/iter) {
1757 continue;
1758 }
1759 if (TMath::Sqrt(chi2Z) > 7.0/iter) {
1760 continue;
69b55c55 1761 }
4551ea7c 1762
1763 Float_t minmax[2] = { -100.0, 100.0 };
1764 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1765 Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer] * 0.5
ca1e1e5b 1766 + 1.0 - cseed[sLayer+iLayer].GetZref(0);
4551ea7c 1767 if (max < minmax[1]) {
1768 minmax[1] = max;
1769 }
1770 Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer] * 0.5
ca1e1e5b 1771 - 1.0 - cseed[sLayer+iLayer].GetZref(0);
4551ea7c 1772 if (min > minmax[0]) {
1773 minmax[0] = min;
1774 }
1775 }
1776
69b55c55 1777 Bool_t isFake = kFALSE;
4551ea7c 1778 if (cl0->GetLabel(0) != cl3->GetLabel(0)) {
1779 isFake = kTRUE;
1780 }
1781 if (cl1->GetLabel(0) != cl3->GetLabel(0)) {
1782 isFake = kTRUE;
1783 }
1784 if (cl2->GetLabel(0) != cl3->GetLabel(0)) {
1785 isFake = kTRUE;
1786 }
1787
1788 if (AliTRDReconstructor::StreamLevel() > 0) {
1789 if ((!isFake) || ((icl3%10) == 0)) { // Debugging print
1790 TTreeSRedirector &cstream = *fDebugStreamer;
1791 cstream << "Seeds0"
1792 << "isFake=" << isFake
1793 << "Cl0.=" << cl0
1794 << "Cl1.=" << cl1
1795 << "Cl2.=" << cl2
1796 << "Cl3.=" << cl3
1797 << "Xref=" << xref
1798 << "X0=" << xcl[sLayer+0]
1799 << "X1=" << xcl[sLayer+1]
1800 << "X2=" << xcl[sLayer+2]
1801 << "X3=" << xcl[sLayer+3]
1802 << "Y2exp=" << y2exp
1803 << "Z2exp=" << z2exp
1804 << "Chi2R=" << chi2R
1805 << "Chi2Z=" << chi2Z
1806 << "Seed0.=" << &cseed[sLayer+0]
1807 << "Seed1.=" << &cseed[sLayer+1]
1808 << "Seed2.=" << &cseed[sLayer+2]
1809 << "Seed3.=" << &cseed[sLayer+3]
1810 << "Zmin=" << minmax[0]
1811 << "Zmax=" << minmax[1]
1812 << "\n";
d337ef8d 1813 }
69b55c55 1814 }
33744e5d 1815
4551ea7c 1816 ////////////////////////////////////////////////////////////////////////////////////
33744e5d 1817 //
4551ea7c 1818 // Fit seeding part
33744e5d 1819 //
4551ea7c 1820 ////////////////////////////////////////////////////////////////////////////////////
33744e5d 1821
69b55c55 1822 cl[sLayer+0] = cl0;
1823 cl[sLayer+1] = cl1;
1824 cl[sLayer+2] = cl2;
1825 cl[sLayer+3] = cl3;
4551ea7c 1826 Bool_t isOK = kTRUE;
1827
1828 for (Int_t jLayer = 0; jLayer < 4; jLayer++) {
1829
ca1e1e5b 1830 cseed[sLayer+jLayer].SetTilt(hL[sLayer+jLayer]);
1831 cseed[sLayer+jLayer].SetPadLength(padlength[sLayer+jLayer]);
1832 cseed[sLayer+jLayer].SetX0(xcl[sLayer+jLayer]);
4551ea7c 1833
1834 for (Int_t iter = 0; iter < 2; iter++) {
1835
69b55c55 1836 //
4551ea7c 1837 // In iteration 0 we try only one pad-row
1838 // If quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
69b55c55 1839 //
1840 AliTRDseed tseed = cseed[sLayer+jLayer];
4551ea7c 1841 Float_t roadz = padlength[sLayer+jLayer] * 0.5;
1842 if (iter > 0) {
1843 roadz = padlength[sLayer+jLayer];
1844 }
1845
1846 Float_t quality = 10000.0;
1847
1848 for (Int_t iTime = 2; iTime < 20; iTime++) {
1849
1850 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
1851 Double_t dxlayer = layer.GetX() - xcl[sLayer+jLayer];
1852 Double_t zexp = cl[sLayer+jLayer]->GetZ();
1853
1854 if (iter > 0) {
1855 // Try 2 pad-rows in second iteration
ca1e1e5b 1856 zexp = tseed.GetZref(0) + tseed.GetZref(1) * dxlayer;
4551ea7c 1857 if (zexp > cl[sLayer+jLayer]->GetZ()) {
1858 zexp = cl[sLayer+jLayer]->GetZ() + padlength[sLayer+jLayer]*0.5;
1859 }
1860 if (zexp < cl[sLayer+jLayer]->GetZ()) {
1861 zexp = cl[sLayer+jLayer]->GetZ() - padlength[sLayer+jLayer]*0.5;
1862 }
69b55c55 1863 }
4551ea7c 1864
ca1e1e5b 1865 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
4551ea7c 1866 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
1867 if (index <= 0) {
1868 continue;
1869 }
1870 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
1871
ca1e1e5b 1872 tseed.SetIndexes(iTime,index);
1873 tseed.SetClusters(iTime,cl); // Register cluster
1874 tseed.SetX(iTime,dxlayer); // Register cluster
1875 tseed.SetY(iTime,cl->GetY()); // Register cluster
1876 tseed.SetZ(iTime,cl->GetZ()); // Register cluster
4551ea7c 1877
1878 }
1879
69b55c55 1880 tseed.Update();
4551ea7c 1881
1882 // Count the number of clusters and distortions into quality
ca1e1e5b 1883 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1884 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
1885 + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
1886 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 1887 if ((iter == 0) && tseed.IsOK()) {
1888 cseed[sLayer+jLayer] = tseed;
1889 quality = tquality;
1890 if (tquality < 5) {
1891 break;
1892 }
1893 }
1894 if (tseed.IsOK() && (tquality < quality)) {
69b55c55 1895 cseed[sLayer+jLayer] = tseed;
69b55c55 1896 }
4551ea7c 1897
1898 } // Loop: iter
1899
1900 if (!cseed[sLayer+jLayer].IsOK()) {
69b55c55 1901 isOK = kFALSE;
1902 break;
4551ea7c 1903 }
1904
69b55c55 1905 cseed[sLayer+jLayer].CookLabels();
1906 cseed[sLayer+jLayer].UpdateUsed();
ca1e1e5b 1907 nusedCl += cseed[sLayer+jLayer].GetNUsed();
4551ea7c 1908 if (nusedCl > 25) {
69b55c55 1909 isOK = kFALSE;
1910 break;
4551ea7c 1911 }
1912
1913 } // Loop: jLayer
1914
1915 if (!isOK) {
1916 continue;
69b55c55 1917 }
4551ea7c 1918 nclusters = 0;
1919 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1920 if (cseed[sLayer+iLayer].IsOK()) {
ca1e1e5b 1921 nclusters += cseed[sLayer+iLayer].GetN2();
69b55c55 1922 }
1923 }
4551ea7c 1924
1925 // Iteration 0
69b55c55 1926 rieman.Reset();
4551ea7c 1927 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1928 rieman.AddPoint(xcl[sLayer+iLayer]
ca1e1e5b 1929 ,cseed[sLayer+iLayer].GetYfitR(0)
1930 ,cseed[sLayer+iLayer].GetZProb()
4551ea7c 1931 ,1
1932 ,10);
69b55c55 1933 }
1934 rieman.Update();
4551ea7c 1935
1936 chi2R = 0.0;
1937 chi2Z = 0.0;
1938
1939 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
ca1e1e5b 1940 cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1941 chi2R += (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0))
1942 * (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0));
1943 cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
1944 cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1945 chi2Z += (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz())
1946 * (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz());
1947 cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
69b55c55 1948 }
1949 Double_t curv = rieman.GetC();
4551ea7c 1950
69b55c55 1951 //
4551ea7c 1952 // Likelihoods
6c94f330 1953 //
ca1e1e5b 1954 Double_t sumda = TMath::Abs(cseed[sLayer+0].GetYfitR(1) - cseed[sLayer+0].GetYref(1))
1955 + TMath::Abs(cseed[sLayer+1].GetYfitR(1) - cseed[sLayer+1].GetYref(1))
1956 + TMath::Abs(cseed[sLayer+2].GetYfitR(1) - cseed[sLayer+2].GetYref(1))
1957 + TMath::Abs(cseed[sLayer+3].GetYfitR(1) - cseed[sLayer+3].GetYref(1));
4551ea7c 1958 Double_t likea = TMath::Exp(-sumda*10.6);
1959 Double_t likechi2 = 0.0000000001;
1960 if (chi2R < 0.5) {
1961 likechi2 += TMath::Exp(-TMath::Sqrt(chi2R) * 7.73);
1962 }
1963 Double_t likechi2z = TMath::Exp(-chi2Z * 0.088) / TMath::Exp(-chi2Z * 0.019);
1964 Double_t likeN = TMath::Exp(-(72 - nclusters) * 0.19);
1965 Double_t like = likea * likechi2 * likechi2z * likeN;
ca1e1e5b 1966 Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetYref(1) - 130.0*curv) * 1.9);
1967 Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetZref(1)
1968 - cseed[sLayer+0].GetZref(0) / xcl[sLayer+0]) * 5.9);
6c94f330 1969 Double_t likePrim = TMath::Max(likePrimY*likePrimZ,0.0005);
69b55c55 1970
4551ea7c 1971 seedquality[registered] = like;
1972 seedlayer[registered] = sLayer;
1973 if (TMath::Log(0.000000000000001 + like) < -15) {
1974 continue;
1975 }
69b55c55 1976 AliTRDseed seedb[6];
4551ea7c 1977 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
69b55c55 1978 seedb[iLayer] = cseed[iLayer];
1979 }
4551ea7c 1980
1981 ////////////////////////////////////////////////////////////////////////////////////
ad670fba 1982 //
4551ea7c 1983 // Full track fit part
33744e5d 1984 //
4551ea7c 1985 ////////////////////////////////////////////////////////////////////////////////////
1986
1987 Int_t nlayers = 0;
1988 Int_t nusedf = 0;
1989 Int_t findable = 0;
1990
69b55c55 1991 //
4551ea7c 1992 // Add new layers - avoid long extrapolation
69b55c55 1993 //
4551ea7c 1994 Int_t tLayer[2] = { 0, 0 };
1995 if (sLayer == 2) {
1996 tLayer[0] = 1;
1997 tLayer[1] = 0;
1998 }
1999 if (sLayer == 1) {
2000 tLayer[0] = 5;
2001 tLayer[1] = 0;
2002 }
2003 if (sLayer == 0) {
2004 tLayer[0] = 4;
2005 tLayer[1] = 5;
2006 }
2007
2008 for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
2009 Int_t jLayer = tLayer[iLayer]; // Set tracking layer
69b55c55 2010 cseed[jLayer].Reset();
ca1e1e5b 2011 cseed[jLayer].SetTilt(hL[jLayer]);
2012 cseed[jLayer].SetPadLength(padlength[jLayer]);
2013 cseed[jLayer].SetX0(xcl[jLayer]);
4551ea7c 2014 // Get pad length and rough cluster
ca1e1e5b 2015 Int_t indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].GetYref(0)
2016 ,cseed[jLayer].GetZref(0)
4551ea7c 2017 ,kRoad2y
2018 ,kRoad2z);
2019 if (indexdummy <= 0) {
2020 continue;
2021 }
2022 AliTRDcluster *cldummy = (AliTRDcluster *) GetCluster(indexdummy);
2023 padlength[jLayer] = TMath::Sqrt(cldummy->GetSigmaZ2() * 12.0);
69b55c55 2024 }
4551ea7c 2025 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
2026
2027 for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
2028
2029 Int_t jLayer = tLayer[iLayer]; // set tracking layer
2030 if ((jLayer == 0) && !(cseed[1].IsOK())) {
2031 continue; // break not allowed
2032 }
2033 if ((jLayer == 5) && !(cseed[4].IsOK())) {
2034 continue; // break not allowed
2035 }
ca1e1e5b 2036 Float_t zexp = cseed[jLayer].GetZref(0);
4551ea7c 2037 Double_t zroad = padlength[jLayer] * 0.5 + 1.0;
2038
2039 for (Int_t iter = 0; iter < 2; iter++) {
2040
69b55c55 2041 AliTRDseed tseed = cseed[jLayer];
4551ea7c 2042 Float_t quality = 10000.0;
2043
2044 for (Int_t iTime = 2; iTime < 20; iTime++) {
2045 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
2046 Double_t dxlayer = layer.GetX()-xcl[jLayer];
ca1e1e5b 2047 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
4551ea7c 2048 Float_t yroad = kRoad1y;
2049 Int_t index = layer.FindNearestCluster(yexp,zexp,yroad,zroad);
2050 if (index <= 0) {
2051 continue;
2052 }
2053 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
ca1e1e5b 2054 tseed.SetIndexes(iTime,index);
2055 tseed.SetClusters(iTime,cl); // Register cluster
2056 tseed.SetX(iTime,dxlayer); // Register cluster
2057 tseed.SetY(iTime,cl->GetY()); // Register cluster
2058 tseed.SetZ(iTime,cl->GetZ()); // Register cluster
4551ea7c 2059 }
2060
69b55c55 2061 tseed.Update();
4551ea7c 2062 if (tseed.IsOK()) {
ca1e1e5b 2063 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
2064 Float_t tquality = (18.0 - tseed.GetN2())/2.0 + TMath::Abs(dangle) / 0.1
2065 + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
2066 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 2067 if (tquality < quality) {
2068 cseed[jLayer] = tseed;
2069 quality = tquality;
69b55c55 2070 }
2071 }
4551ea7c 2072
2073 zroad *= 2.0;
2074
2075 } // Loop: iter
2076
2077 if ( cseed[jLayer].IsOK()) {
69b55c55 2078 cseed[jLayer].CookLabels();
2079 cseed[jLayer].UpdateUsed();
ca1e1e5b 2080 nusedf += cseed[jLayer].GetNUsed();
4551ea7c 2081 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
69b55c55 2082 }
4551ea7c 2083
2084 } // Loop: iLayer
2085
2086 // Make copy
69b55c55 2087 AliTRDseed bseed[6];
4551ea7c 2088 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
69b55c55 2089 bseed[jLayer] = cseed[jLayer];
4551ea7c 2090 }
2091 Float_t lastquality = 10000.0;
2092 Float_t lastchi2 = 10000.0;
2093 Float_t chi2 = 1000.0;
2094
2095 for (Int_t iter = 0; iter < 4; iter++) {
2096
2097 // Sort tracklets according "quality", try to "improve" 4 worst
2098 Float_t sumquality = 0.0;
69b55c55 2099 Float_t squality[6];
2100 Int_t sortindexes[6];
4551ea7c 2101
2102 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2103 if (bseed[jLayer].IsOK()) {
69b55c55 2104 AliTRDseed &tseed = bseed[jLayer];
ca1e1e5b 2105 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
2106 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
2107 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
2108 + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
2109 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 2110 squality[jLayer] = tquality;
2111 }
2112 else {
2113 squality[jLayer] = -1.0;
ad670fba 2114 }
6c94f330 2115 sumquality +=squality[jLayer];
69b55c55 2116 }
2117
4551ea7c 2118 if ((sumquality >= lastquality) ||
2119 (chi2 > lastchi2)) {
2120 break;
2121 }
69b55c55 2122 lastquality = sumquality;
2123 lastchi2 = chi2;
4551ea7c 2124 if (iter > 0) {
2125 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
69b55c55 2126 cseed[jLayer] = bseed[jLayer];
2127 }
2128 }
2129 TMath::Sort(6,squality,sortindexes,kFALSE);
4551ea7c 2130
2131 for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
2132
6c94f330 2133 Int_t bLayer = sortindexes[jLayer];
2134 AliTRDseed tseed = bseed[bLayer];
4551ea7c 2135
2136 for (Int_t iTime = 2; iTime < 20; iTime++) {
2137
2138 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
2139 Double_t dxlayer = layer.GetX() - xcl[bLayer];
ca1e1e5b 2140 Double_t zexp = tseed.GetZref(0);
2141 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
4551ea7c 2142 Float_t roadz = padlength[bLayer] + 1;
ca1e1e5b 2143 if (TMath::Abs(tseed.GetZProb() - zexp) > 0.5*padlength[bLayer]) {
4551ea7c 2144 roadz = padlength[bLayer] * 0.5;
2145 }
ca1e1e5b 2146 if (tseed.GetZfit(1)*tseed.GetZref(1) < 0.0) {
4551ea7c 2147 roadz = padlength[bLayer] * 0.5;
2148 }
ca1e1e5b 2149 if (TMath::Abs(tseed.GetZProb() - zexp) < 0.1*padlength[bLayer]) {
2150 zexp = tseed.GetZProb();
4551ea7c 2151 roadz = padlength[bLayer] * 0.5;
2152 }
2153
ca1e1e5b 2154 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer - zcor;
4551ea7c 2155 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
2156 if (index <= 0) {
2157 continue;
69b55c55 2158 }
4551ea7c 2159 AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
2160
ca1e1e5b 2161 tseed.SetIndexes(iTime,index);
2162 tseed.SetClusters(iTime,cl); // Register cluster
2163 tseed.SetX(iTime,dxlayer); // Register cluster
2164 tseed.SetY(iTime,cl->GetY()); // Register cluster
2165 tseed.SetZ(iTime,cl->GetZ()); // Register cluster
4551ea7c 2166
2167 }
2168
69b55c55 2169 tseed.Update();
c6f438c0 2170 if (tseed.IsOK()) {
ca1e1e5b 2171 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
2172 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
2173 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0
2174 + TMath::Abs(dangle) / 0.1
2175 + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
2176 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
4551ea7c 2177 if (tquality<squality[bLayer]) {
69b55c55 2178 bseed[bLayer] = tseed;
4551ea7c 2179 }
69b55c55 2180 }
4551ea7c 2181
2182 } // Loop: jLayer
2183
2184 chi2 = AliTRDseed::FitRiemanTilt(bseed,kTRUE);
2185
2186 } // Loop: iter
2187
2188 nclusters = 0;
2189 nlayers = 0;
2190 findable = 0;
2191 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ca1e1e5b 2192 if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) {
69b55c55 2193 findable++;
4551ea7c 2194 }
2195 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2196 nclusters += cseed[iLayer].GetN2();
69b55c55 2197 nlayers++;
2198 }
2199 }
4551ea7c 2200 if (nlayers < 3) {
2201 continue;
2202 }
69b55c55 2203 rieman.Reset();
4551ea7c 2204 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2205 if (cseed[iLayer].IsOK()) {
2206 rieman.AddPoint(xcl[iLayer]
ca1e1e5b 2207 ,cseed[iLayer].GetYfitR(0)
2208 ,cseed[iLayer].GetZProb()
4551ea7c 2209 ,1
2210 ,10);
2211 }
69b55c55 2212 }
2213 rieman.Update();
4551ea7c 2214
2215 chi2RF = 0.0;
2216 chi2ZF = 0.0;
2217 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2218 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2219 cseed[iLayer].SetYref(0,rieman.GetYat(xcl[iLayer]));
2220 chi2RF += (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0))
2221 * (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0));
2222 cseed[iLayer].SetYref(1,rieman.GetDYat(xcl[iLayer]));
2223 cseed[iLayer].SetZref(0,rieman.GetZat(xcl[iLayer]));
2224 chi2ZF += (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz())
2225 * (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz());
2226 cseed[iLayer].SetZref(1,rieman.GetDZat(xcl[iLayer]));
69b55c55 2227 }
2228 }
4551ea7c 2229 chi2RF /= TMath::Max((nlayers - 3.0),1.0);
2230 chi2ZF /= TMath::Max((nlayers - 3.0),1.0);
2231 curv = rieman.GetC();
2232
2233 Double_t xref2 = (xcl[2] + xcl[3]) * 0.5; // Middle of the chamber
2234 Double_t dzmf = rieman.GetDZat(xref2);
2235 Double_t zmf = rieman.GetZat(xref2);
69b55c55 2236
69b55c55 2237 //
4551ea7c 2238 // Fit hyperplane
69b55c55 2239 //
4551ea7c 2240 Int_t npointsT = 0;
69b55c55 2241 fitterTC.ClearPoints();
2242 fitterT2.ClearPoints();
2243 rieman2.Reset();
4551ea7c 2244
2245 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2246
2247 if (!cseed[iLayer].IsOK()) {
2248 continue;
2249 }
2250
2251 for (Int_t itime = 0; itime < 25; itime++) {
2252
ca1e1e5b 2253 if (!cseed[iLayer].IsUsable(itime)) {
4551ea7c 2254 continue;
2255 }
2256 // X relative to the middle chamber
ca1e1e5b 2257 Double_t x = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;
2258 Double_t y = cseed[iLayer].GetY(itime);
2259 Double_t z = cseed[iLayer].GetZ(itime);
69b55c55 2260 // ExB correction to the correction
4551ea7c 2261 // Tilted rieman
69b55c55 2262 Double_t uvt[6];
4551ea7c 2263 // Global x
ca1e1e5b 2264 Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();
4551ea7c 2265 Double_t t = 1.0 / (x2*x2 + y*y);
2266 uvt[1] = t; // t
2267 uvt[0] = 2.0 * x2 * uvt[1]; // u
2268 uvt[2] = 2.0 * hL[iLayer] * uvt[1];
2269 uvt[3] = 2.0 * hL[iLayer] * x * uvt[1];
2270 uvt[4] = 2.0 * (y + hL[iLayer]*z) * uvt[1];
2271 Double_t error = 2.0 * 0.2 * uvt[1];
69b55c55 2272 fitterT2.AddPoint(uvt,uvt[4],error);
4551ea7c 2273
69b55c55 2274 //
4551ea7c 2275 // Constrained rieman
69b55c55 2276 //
ca1e1e5b 2277 z = cseed[iLayer].GetZ(itime);
4551ea7c 2278 uvt[0] = 2.0 * x2 * t; // u
2279 uvt[1] = 2.0 * hL[iLayer] * x2 * uvt[1];
2280 uvt[2] = 2.0 * (y + hL[iLayer] * (z - GetZ())) * t;
69b55c55 2281 fitterTC.AddPoint(uvt,uvt[2],error);
69b55c55 2282 rieman2.AddPoint(x2,y,z,1,10);
2283 npointsT++;
4551ea7c 2284
69b55c55 2285 }
4551ea7c 2286
2287 } // Loop: iLayer
2288
69b55c55 2289 rieman2.Update();
2290 fitterTC.Eval();
2291 fitterT2.Eval();
2292 Double_t rpolz0 = fitterT2.GetParameter(3);
2293 Double_t rpolz1 = fitterT2.GetParameter(4);
4551ea7c 2294
69b55c55 2295 //
4551ea7c 2296 // Linear fitter - not possible to make boundaries
2297 // Do not accept non possible z and dzdx combinations
69b55c55 2298 //
4551ea7c 2299 Bool_t acceptablez = kTRUE;
2300 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2301 if (cseed[iLayer].IsOK()) {
2302 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
ca1e1e5b 2303 if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > padlength[iLayer] * 0.5 + 1.0) {
69b55c55 2304 acceptablez = kFALSE;
4551ea7c 2305 }
69b55c55 2306 }
2307 }
4551ea7c 2308 if (!acceptablez) {
69b55c55 2309 fitterT2.FixParameter(3,zmf);
2310 fitterT2.FixParameter(4,dzmf);
2311 fitterT2.Eval();
2312 fitterT2.ReleaseParameter(3);
2313 fitterT2.ReleaseParameter(4);
2314 rpolz0 = fitterT2.GetParameter(3);
2315 rpolz1 = fitterT2.GetParameter(4);
2316 }
4551ea7c 2317
2318 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
2319 Double_t chi2TC = fitterTC.GetChisquare() / Float_t(npointsT);
69b55c55 2320 Double_t polz1c = fitterTC.GetParameter(2);
4551ea7c 2321 Double_t polz0c = polz1c * xref2;
6c94f330 2322 Double_t aC = fitterTC.GetParameter(0);
2323 Double_t bC = fitterTC.GetParameter(1);
4551ea7c 2324 Double_t cC = aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
6c94f330 2325 Double_t aR = fitterT2.GetParameter(0);
2326 Double_t bR = fitterT2.GetParameter(1);
2327 Double_t dR = fitterT2.GetParameter(2);
4551ea7c 2328 Double_t cR = 1.0 + bR*bR - dR*aR;
2329 Double_t dca = 0.0;
2330 if (cR > 0.0) {
2331 dca = -dR / (TMath::Sqrt(1.0 + bR*bR - dR*aR) + TMath::Sqrt(1.0 + bR*bR));
2332 cR = aR / TMath::Sqrt(cR);
69b55c55 2333 }
4551ea7c 2334
2335 Double_t chi2ZT2 = 0.0;
2336 Double_t chi2ZTC = 0.0;
2337 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2338 if (cseed[iLayer].IsOK()) {
2339 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
2340 Double_t zTC = polz0c + polz1c * (xcl[iLayer] - xref2);
ca1e1e5b 2341 chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
2342 chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
69b55c55 2343 }
2344 }
4551ea7c 2345 chi2ZT2 /= TMath::Max((nlayers - 3.0),1.0);
2346 chi2ZTC /= TMath::Max((nlayers - 3.0),1.0);
2347
2348 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
2349 Float_t sumdaf = 0.0;
2350 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2351 if (cseed[iLayer].IsOK()) {
ca1e1e5b 2352 sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))
2353 / cseed[iLayer].GetSigmaY2());
4551ea7c 2354 }
2355 }
2356 sumdaf /= Float_t (nlayers - 2.0);
2357
69b55c55 2358 //
4551ea7c 2359 // Likelihoods for full track
69b55c55 2360 //
4551ea7c 2361 Double_t likezf = TMath::Exp(-chi2ZF * 0.14);
2362 Double_t likechi2C = TMath::Exp(-chi2TC * 0.677);
2363 Double_t likechi2TR = TMath::Exp(-chi2TR * 0.78);
2364 Double_t likeaf = TMath::Exp(-sumdaf * 3.23);
2365 seedquality2[registered] = likezf * likechi2TR * likeaf;
33744e5d 2366
4551ea7c 2367 // Still needed ????
69b55c55 2368// Bool_t isGold = kFALSE;
2369//
6c94f330 2370// if (nlayers == 6 && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE; // gold
2371// if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE; // gold
69b55c55 2372// if (isGold &&nusedf<10){
2373// for (Int_t jLayer=0;jLayer<6;jLayer++){
6c94f330 2374// if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
69b55c55 2375// seed[index][jLayer].UseClusters(); //sign gold
2376// }
2377// }
4551ea7c 2378
2379 Int_t index0 = 0;
2380 if (!cseed[0].IsOK()) {
69b55c55 2381 index0 = 1;
4551ea7c 2382 if (!cseed[1].IsOK()) {
2383 index0 = 2;
2384 }
69b55c55 2385 }
ca1e1e5b 2386 seedparams[registered][0] = cseed[index0].GetX0();
2387 seedparams[registered][1] = cseed[index0].GetYref(0);
2388 seedparams[registered][2] = cseed[index0].GetZref(0);
c6f438c0 2389 seedparams[registered][5] = cR;
ca1e1e5b 2390 seedparams[registered][3] = cseed[index0].GetX0() * cR - TMath::Sin(TMath::ATan(cseed[0].GetYref(1)));
2391 seedparams[registered][4] = cseed[index0].GetZref(1)
2392 / TMath::Sqrt(1.0 + cseed[index0].GetYref(1) * cseed[index0].GetYref(1));
69b55c55 2393 seedparams[registered][6] = ns;
4551ea7c 2394
2395 Int_t labels[12];
2396 Int_t outlab[24];
2397 Int_t nlab = 0;
2398 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2399 if (!cseed[iLayer].IsOK()) {
2400 continue;
2401 }
ca1e1e5b 2402 if (cseed[iLayer].GetLabels(0) >= 0) {
2403 labels[nlab] = cseed[iLayer].GetLabels(0);
69b55c55 2404 nlab++;
2405 }
ca1e1e5b 2406 if (cseed[iLayer].GetLabels(1) >= 0) {
2407 labels[nlab] = cseed[iLayer].GetLabels(1);
69b55c55 2408 nlab++;
4551ea7c 2409 }
69b55c55 2410 }
2411 Freq(nlab,labels,outlab,kFALSE);
4551ea7c 2412 Int_t label = outlab[0];
2413 Int_t frequency = outlab[1];
2414 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ca1e1e5b 2415 cseed[iLayer].SetFreq(frequency);
2416 cseed[iLayer].SetC(cR);
2417 cseed[iLayer].SetCC(cC);
2418 cseed[iLayer].SetChi2(chi2TR);
2419 cseed[iLayer].SetChi2Z(chi2ZF);
69b55c55 2420 }
33744e5d 2421
2422 // Debugging print
4551ea7c 2423 if (1 || (!isFake)) {
69b55c55 2424 Float_t zvertex = GetZ();
4551ea7c 2425 TTreeSRedirector &cstream = *fDebugStreamer;
2426 if (AliTRDReconstructor::StreamLevel() > 0) {
2427 cstream << "Seeds1"
2428 << "isFake=" << isFake
2429 << "Vertex=" << zvertex
2430 << "Rieman2.=" << &rieman2
2431 << "Rieman.=" << &rieman
2432 << "Xref=" << xref
2433 << "X0=" << xcl[0]
2434 << "X1=" << xcl[1]
2435 << "X2=" << xcl[2]
2436 << "X3=" << xcl[3]
2437 << "X4=" << xcl[4]
2438 << "X5=" << xcl[5]
2439 << "Chi2R=" << chi2R
2440 << "Chi2Z=" << chi2Z
2441 << "Chi2RF=" << chi2RF // Chi2 of trackletes on full track
2442 << "Chi2ZF=" << chi2ZF // Chi2 z on tracklets on full track
2443 << "Chi2ZT2=" << chi2ZT2 // Chi2 z on tracklets on full track - rieman tilt
2444 << "Chi2ZTC=" << chi2ZTC // Chi2 z on tracklets on full track - rieman tilt const
2445 << "Chi2TR=" << chi2TR // Chi2 without vertex constrain
2446 << "Chi2TC=" << chi2TC // Chi2 with vertex constrain
2447 << "C=" << curv // Non constrained - no tilt correction
2448 << "DR=" << dR // DR parameter - tilt correction
2449 << "DCA=" << dca // DCA - tilt correction
2450 << "CR=" << cR // Non constrained curvature - tilt correction
2451 << "CC=" << cC // Constrained curvature
2452 << "Polz0=" << polz0c
2453 << "Polz1=" << polz1c
2454 << "RPolz0=" << rpolz0
2455 << "RPolz1=" << rpolz1
2456 << "Ncl=" << nclusters
2457 << "Nlayers=" << nlayers
2458 << "NUsedS=" << nusedCl
2459 << "NUsed=" << nusedf
2460 << "Findable=" << findable
2461 << "Like=" << like
2462 << "LikePrim=" << likePrim
2463 << "Likechi2C=" << likechi2C
2464 << "Likechi2TR=" << likechi2TR
2465 << "Likezf=" << likezf
2466 << "LikeF=" << seedquality2[registered]
2467 << "S0.=" << &cseed[0]
2468 << "S1.=" << &cseed[1]
2469 << "S2.=" << &cseed[2]
2470 << "S3.=" << &cseed[3]
2471 << "S4.=" << &cseed[4]
2472 << "S5.=" << &cseed[5]
2473 << "SB0.=" << &seedb[0]
2474 << "SB1.=" << &seedb[1]
2475 << "SB2.=" << &seedb[2]
2476 << "SB3.=" << &seedb[3]
2477 << "SB4.=" << &seedb[4]
2478 << "SB5.=" << &seedb[5]
2479 << "Label=" << label
2480 << "Freq=" << frequency
2481 << "sLayer=" << sLayer
2482 << "\n";
2483 }
69b55c55 2484 }
4551ea7c 2485
2486 if (registered<kMaxSeed - 1) {
69b55c55 2487 registered++;
2488 cseed = seed[registered];
2489 }
4551ea7c 2490
2491 } // End of loop over layer 1
2492
2493 } // End of loop over layer 0
2494
2495 } // End of loop over layer 3
2496
2497 } // End of loop over seeding time bins
2498
69b55c55 2499 //
4551ea7c 2500 // Choose best
69b55c55 2501 //
4551ea7c 2502
69b55c55 2503 TMath::Sort(registered,seedquality2,sort,kTRUE);
c6f438c0 2504 Bool_t signedseed[kMaxSeed];
4551ea7c 2505 for (Int_t i = 0; i < registered; i++) {
2506 signedseed[i] = kFALSE;
69b55c55 2507 }
4551ea7c 2508
2509 for (Int_t iter = 0; iter < 5; iter++) {
2510
2511 for (Int_t iseed = 0; iseed < registered; iseed++) {
2512
69b55c55 2513 Int_t index = sort[iseed];
4551ea7c 2514 if (signedseed[index]) {
2515 continue;
2516 }
69b55c55 2517 Int_t labelsall[1000];
4551ea7c 2518 Int_t nlabelsall = 0;
2519 Int_t naccepted = 0;;
2520 Int_t sLayer = seedlayer[index];
2521 Int_t ncl = 0;
2522 Int_t nused = 0;
2523 Int_t nlayers = 0;
69b55c55 2524 Int_t findable = 0;
4551ea7c 2525
2526 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2527
ca1e1e5b 2528 if (TMath::Abs(seed[index][jLayer].GetYref(0) / xcl[jLayer]) < 0.15) {
69b55c55 2529 findable++;
4551ea7c 2530 }
2531 if (seed[index][jLayer].IsOK()) {
69b55c55 2532 seed[index][jLayer].UpdateUsed();
ca1e1e5b 2533 ncl +=seed[index][jLayer].GetN2();
2534 nused +=seed[index][jLayer].GetNUsed();
69b55c55 2535 nlayers++;
4551ea7c 2536 // Cooking label
2537 for (Int_t itime = 0; itime < 25; itime++) {
ca1e1e5b 2538 if (seed[index][jLayer].IsUsable(itime)) {
69b55c55 2539 naccepted++;
4551ea7c 2540 for (Int_t ilab = 0; ilab < 3; ilab++) {
ca1e1e5b 2541 Int_t tindex = seed[index][jLayer].GetClusters(itime)->GetLabel(ilab);
4551ea7c 2542 if (tindex >= 0) {
69b55c55 2543 labelsall[nlabelsall] = tindex;
2544 nlabelsall++;
2545 }
2546 }
2547 }
2548 }
2549 }
4551ea7c 2550
69b55c55 2551 }
4551ea7c 2552
2553 if (nused > 30) {
2554 continue;
69b55c55 2555 }
4551ea7c 2556
2557 if (iter == 0) {
2558 if (nlayers < 6) {
2559 continue;
2560 }
2561 if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2562 continue; // Gold
2563 }
7ad19338 2564 }
4551ea7c 2565
2566 if (iter == 1) {
2567 if (nlayers < findable) {
2568 continue;
2569 }
2570 if (TMath::Log(0.000000001+seedquality2[index]) < -4.0) {
2571 continue;
2572 }
2573 }
2574
2575 if (iter == 2) {
2576 if ((nlayers == findable) ||
2577 (nlayers == 6)) {
2578 continue;
2579 }
2580 if (TMath::Log(0.000000001+seedquality2[index]) < -6.0) {
2581 continue;
2582 }
69b55c55 2583 }
4551ea7c 2584
2585 if (iter == 3) {
2586 if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2587 continue;
2588 }
69b55c55 2589 }
4551ea7c 2590
2591 if (iter == 4) {
2592 if (TMath::Log(0.000000001+seedquality2[index]) - nused/(nlayers-3.0) < -15.0) {
2593 continue;
2594 }
69b55c55 2595 }
4551ea7c 2596
69b55c55 2597 signedseed[index] = kTRUE;
4551ea7c 2598
2599 Int_t labels[1000];
2600 Int_t outlab[1000];
2601 Int_t nlab = 0;
2602 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2603 if (seed[index][iLayer].IsOK()) {
ca1e1e5b 2604 if (seed[index][iLayer].GetLabels(0) >= 0) {
2605 labels[nlab] = seed[index][iLayer].GetLabels(0);
69b55c55 2606 nlab++;
2607 }
ca1e1e5b 2608 if (seed[index][iLayer].GetLabels(1) >= 0) {
2609 labels[nlab] = seed[index][iLayer].GetLabels(1);
69b55c55 2610 nlab++;
4551ea7c 2611 }
2612 }
7ad19338 2613 }
69b55c55 2614 Freq(nlab,labels,outlab,kFALSE);
4551ea7c 2615 Int_t label = outlab[0];
2616 Int_t frequency = outlab[1];
69b55c55 2617 Freq(nlabelsall,labelsall,outlab,kFALSE);
4551ea7c 2618 Int_t label1 = outlab[0];
2619 Int_t label2 = outlab[2];
2620 Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
2621 Float_t ratio = Float_t(nused) / Float_t(ncl);
2622 if (ratio < 0.25) {
2623 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2624 if ((seed[index][jLayer].IsOK()) &&
ca1e1e5b 2625 (TMath::Abs(seed[index][jLayer].GetYfit(1) - seed[index][jLayer].GetYfit(1)) < 0.2)) {
4551ea7c 2626 seed[index][jLayer].UseClusters(); // Sign gold
2627 }
69b55c55 2628 }
7ad19338 2629 }
4551ea7c 2630
31fd97b2 2631 Int_t eventNrInFile = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
4551ea7c 2632 TTreeSRedirector &cstream = *fDebugStreamer;
2633
69b55c55 2634 //
4551ea7c 2635 // Register seed
69b55c55 2636 //
4551ea7c 2637 AliTRDtrack *track = RegisterSeed(seed[index],seedparams[index]);
2638 AliTRDtrack dummy;
2639 if (!track) {
2640 track = &dummy;
2641 }
2642 else {
69b55c55 2643 AliESDtrack esdtrack;
4551ea7c 2644 esdtrack.UpdateTrackParams(track,AliESDtrack::kTRDout);
69b55c55 2645 esdtrack.SetLabel(label);
2646 esd->AddTrack(&esdtrack);
4551ea7c 2647 TTreeSRedirector &cstream = *fDebugStreamer;
2648 if (AliTRDReconstructor::StreamLevel() > 0) {
2649 cstream << "Tracks"
31fd97b2 2650 << "EventNrInFile=" << eventNrInFile
4551ea7c 2651 << "ESD.=" << &esdtrack
2652 << "trd.=" << track
2653 << "trdback.=" << track
2654 << "\n";
2655 }
7ad19338 2656 }
4551ea7c 2657
2658 if (AliTRDReconstructor::StreamLevel() > 0) {
2659 cstream << "Seeds2"
2660 << "Iter=" << iter
2661 << "Track.=" << track
2662 << "Like=" << seedquality[index]
2663 << "LikeF=" << seedquality2[index]
2664 << "S0.=" << &seed[index][0]
2665 << "S1.=" << &seed[index][1]
2666 << "S2.=" << &seed[index][2]
2667 << "S3.=" << &seed[index][3]
2668 << "S4.=" << &seed[index][4]
2669 << "S5.=" << &seed[index][5]
2670 << "Label=" << label
2671 << "Label1=" << label1
2672 << "Label2=" << label2
2673 << "FakeRatio=" << fakeratio
2674 << "Freq=" << frequency
2675 << "Ncl=" << ncl
2676 << "Nlayers=" << nlayers
2677 << "Findable=" << findable
2678 << "NUsed=" << nused
2679 << "sLayer=" << sLayer
31fd97b2 2680 << "EventNrInFile=" << eventNrInFile
4551ea7c 2681 << "\n";
2682 }
2683
2684 } // Loop: iseed
2685
2686 } // Loop: iter
2687
2688 } // End of loop over sectors
75bd7f81 2689
69b55c55 2690 delete [] pseed;
75bd7f81 2691
69b55c55 2692}
2693
5443e65e 2694//_____________________________________________________________________________
4551ea7c 2695Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *clusterTree) const
5443e65e 2696{
2697 //
a819a5f7 2698 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
2699 // from the file. The names of the cluster tree and branches
2700 // should match the ones used in AliTRDclusterizer::WriteClusters()
2701 //
75bd7f81 2702
4551ea7c 2703 Int_t nsize = Int_t(clusterTree->GetTotBytes() / (sizeof(AliTRDcluster)));
6c94f330 2704 TObjArray *clusterArray = new TObjArray(nsize+1000);
5443e65e 2705
4551ea7c 2706 TBranch *branch = clusterTree->GetBranch("TRDcluster");
c630aafd 2707 if (!branch) {
4551ea7c 2708 AliError("Can't get the branch !");
c630aafd 2709 return 1;
2710 }
029cd327 2711 branch->SetAddress(&clusterArray);
5443e65e 2712
a819a5f7 2713 // Loop through all entries in the tree
4551ea7c 2714 Int_t nEntries = (Int_t) clusterTree->GetEntries();
2715 Int_t nbytes = 0;
a819a5f7 2716 AliTRDcluster *c = 0;
a819a5f7 2717 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2718
2719 // Import the tree
4551ea7c 2720 nbytes += clusterTree->GetEvent(iEntry);
5443e65e 2721
a819a5f7 2722 // Get the number of points in the detector
029cd327 2723 Int_t nCluster = clusterArray->GetEntriesFast();
5443e65e 2724
a819a5f7 2725 // Loop through all TRD digits
2726 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
4551ea7c 2727 c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
4f1c04d3 2728 AliTRDcluster *co = c;
a819a5f7 2729 array->AddLast(co);
4f1c04d3 2730 clusterArray->RemoveAt(iCluster);
a819a5f7 2731 }
4551ea7c 2732
a819a5f7 2733 }
2734
029cd327 2735 delete clusterArray;
5443e65e 2736
c630aafd 2737 return 0;
75bd7f81 2738
a819a5f7 2739}
2740
75bd7f81 2741//_____________________________________________________________________________
4551ea7c 2742Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint &p) const
3551db50 2743{
2744 //
2745 // Get track space point with index i
2746 // Origin: C.Cheshkov
2747 //
2748
4551ea7c 2749 AliTRDcluster *cl = (AliTRDcluster *) fClusters->UncheckedAt(index);
2750 Int_t idet = cl->GetDetector();
2751 Int_t isector = fGeom->GetSector(idet);
2752 Int_t ichamber = fGeom->GetChamber(idet);
2753 Int_t iplan = fGeom->GetPlane(idet);
3551db50 2754 Double_t local[3];
4551ea7c 2755 local[0] = GetX(isector,iplan,cl->GetLocalTimeBin());
2756 local[1] = cl->GetY();
2757 local[2] = cl->GetZ();
3551db50 2758 Double_t global[3];
2759 fGeom->RotateBack(idet,local,global);
2760 p.SetXYZ(global[0],global[1],global[2]);
ae079791 2761 AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
3551db50 2762 switch (iplan) {
2763 case 0:
ae079791 2764 iLayer = AliGeomManager::kTRD1;
3551db50 2765 break;
2766 case 1:
ae079791 2767 iLayer = AliGeomManager::kTRD2;
3551db50 2768 break;
2769 case 2:
ae079791 2770 iLayer = AliGeomManager::kTRD3;
3551db50 2771 break;
2772 case 3:
ae079791 2773 iLayer = AliGeomManager::kTRD4;
3551db50 2774 break;
2775 case 4:
ae079791 2776 iLayer = AliGeomManager::kTRD5;
3551db50 2777 break;
2778 case 5:
ae079791 2779 iLayer = AliGeomManager::kTRD6;
3551db50 2780 break;
2781 };
4551ea7c 2782 Int_t modId = isector * fGeom->Ncham() + ichamber;
ae079791 2783 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,modId);
3551db50 2784 p.SetVolumeID(volid);
2785
2786 return kTRUE;
2787
2788}
2789
75bd7f81 2790//_____________________________________________________________________________
4551ea7c 2791void AliTRDtracker::CookLabel(AliKalmanTrack *pt, Float_t wrong) const
029cd327 2792{
2793 //
2794 // This cooks a label. Mmmmh, smells good...
2795 //
46d29e70 2796
4551ea7c 2797 Int_t label = 123456789;
2798 Int_t index;
2799 Int_t i;
2800 Int_t j;
2801 Int_t ncl = pt->GetNumberOfClusters();
2802
2803 const Int_t kRange = fTrSec[0]->GetOuterTimeBin() + 1;
5443e65e 2804
029cd327 2805 Bool_t labelAdded;
46d29e70 2806
6c94f330 2807 Int_t **s = new Int_t* [kRange];
4551ea7c 2808 for (i = 0; i < kRange; i++) {
d1b06c24 2809 s[i] = new Int_t[2];
2810 }
4551ea7c 2811 for (i = 0; i < kRange; i++) {
2812 s[i][0] = -1;
2813 s[i][1] = 0;
46d29e70 2814 }
2815
4551ea7c 2816 Int_t t0;
2817 Int_t t1;
2818 Int_t t2;
2819
2820 for (i = 0; i < ncl; i++) {
2821 index = pt->GetClusterIndex(i);
2822 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
6c94f330 2823 t0=c->GetLabel(0);
2824 t1=c->GetLabel(1);
2825 t2=c->GetLabel(2);
46d29e70 2826 }
2827
4551ea7c 2828 for (i = 0; i < ncl; i++) {
2829 index = pt->GetClusterIndex(i);
2830 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2831 for (Int_t k = 0; k < 3; k++) {
2832 label = c->GetLabel(k);
2833 labelAdded = kFALSE;
2834 j = 0;
46d29e70 2835 if (label >= 0) {
4551ea7c 2836 while ((!labelAdded) && (j < kRange)) {
2837 if ((s[j][0] == label) ||
2838 (s[j][1] == 0)) {
2839 s[j][0] = label;
2840 s[j][1] = s[j][1] + 1;
2841 labelAdded = kTRUE;
a9814c09 2842 }
2843 j++;
2844 }
46d29e70 2845 }
2846 }
2847 }
2848
4551ea7c 2849 Int_t max = 0;
6c94f330 2850 label = -123456789;
46d29e70 2851
4551ea7c 2852 for (i = 0; i < kRange; i++) {
2853 if (s[i][1] > max) {
2854 max = s[i][1];
2855 label = s[i][0];
46d29e70 2856 }
2857 }
5443e65e 2858
4551ea7c 2859 for (i = 0; i < kRange; i++) {
6c94f330 2860 delete []s[i];
5443e65e 2861 }
2862
6c94f330 2863 delete []s;
5443e65e 2864
4551ea7c 2865 if ((1.0 - Float_t(max)/ncl) > wrong) {
2866 label = -label;
2867 }
5443e65e 2868
2869 pt->SetLabel(label);
2870
46d29e70 2871}
2872
75bd7f81 2873//_____________________________________________________________________________
4551ea7c 2874void AliTRDtracker::UseClusters(const AliKalmanTrack *t, Int_t from) const
029cd327 2875{
2876 //
2877 // Use clusters, but don't abuse them!
2878 //
75bd7f81 2879
4551ea7c 2880 const Float_t kmaxchi2 = 18;
2881 const Float_t kmincl = 10;
2882
2883 AliTRDtrack *track = (AliTRDtrack *) t;
2884
2885 Int_t ncl = t->GetNumberOfClusters();
2886 for (Int_t i = from; i < ncl; i++) {
2887 Int_t index = t->GetClusterIndex(i);
2888 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
69b55c55 2889 Int_t iplane = fGeom->GetPlane(c->GetDetector());
27eaf44b 2890 if (track->GetTracklets(iplane).GetChi2() > kmaxchi2) {
4551ea7c 2891 continue;
2892 }
27eaf44b 2893 if (track->GetTracklets(iplane).GetN() < kmincl) {
4551ea7c 2894 continue;
2895 }
2896 if (!(c->IsUsed())) {
2897 c->Use();
2898 }
5443e65e 2899 }
5443e65e 2900
75bd7f81 2901}
5443e65e 2902
75bd7f81 2903//_____________________________________________________________________________
029cd327 2904Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
5443e65e 2905{
75bd7f81 2906 //
5443e65e 2907 // Parametrised "expected" error of the cluster reconstruction in Y
75bd7f81 2908 //
5443e65e 2909
2910 Double_t s = 0.08 * 0.08;
2911 return s;
75bd7f81 2912
5443e65e 2913}
2914
75bd7f81 2915//_____________________________________________________________________________
029cd327 2916Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
0a29d0f1 2917{
75bd7f81 2918 //
5443e65e 2919 // Parametrised "expected" error of the cluster reconstruction in Z
75bd7f81 2920 //
5443e65e 2921
4551ea7c 2922 Double_t s = 9.0 * 9.0 / 12.0;
5443e65e 2923 return s;
75bd7f81 2924
5443e65e 2925}
2926
75bd7f81 2927//_____________________________________________________________________________
029cd327 2928Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
5443e65e 2929{
2930 //
029cd327 2931 // Returns radial position which corresponds to time bin <localTB>
5443e65e 2932 // in tracking sector <sector> and plane <plane>
2933 //
2934
6c94f330 2935 Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
4551ea7c 2936 Int_t pl = fTrSec[sector]->GetLayerNumber(index);
2937
5443e65e 2938 return fTrSec[sector]->GetLayer(pl)->GetX();
2939
2940}
2941
75bd7f81 2942//_____________________________________________________________________________
2943AliTRDtracker::AliTRDpropagationLayer
2944 ::AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
2945 , Double_t radLength, Int_t tbIndex, Int_t plane)
ad670fba 2946 :fN(0)
2947 ,fSec(0)
2948 ,fClusters(NULL)
2949 ,fIndex(NULL)
2950 ,fX(x)
2951 ,fdX(dx)
2952 ,fRho(rho)
2953 ,fX0(radLength)
2954 ,fTimeBinIndex(tbIndex)
2955 ,fPlane(plane)
2956 ,fYmax(0)
2957 ,fYmaxSensitive(0)
2958 ,fHole(kFALSE)
2959 ,fHoleZc(0)
2960 ,fHoleZmax(0)
2961 ,fHoleYc(0)
2962 ,fHoleYmax(0)
2963 ,fHoleRho(0)
2964 ,fHoleX0(0)
5443e65e 2965{
0a29d0f1 2966 //
5443e65e 2967 // AliTRDpropagationLayer constructor
0a29d0f1 2968 //
46d29e70 2969
ad670fba 2970 for (Int_t i = 0; i < (Int_t) kZones; i++) {
2971 fZc[i] = 0;
2972 fZmax[i] = 0;
a819a5f7 2973 }
5443e65e 2974
ad670fba 2975 if (fTimeBinIndex >= 0) {
029cd327 2976 fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
ad670fba 2977 fIndex = new UInt_t[kMaxClusterPerTimeBin];
a819a5f7 2978 }
46d29e70 2979
ad670fba 2980 for (Int_t i = 0; i < 5; i++) {
2981 fIsHole[i] = kFALSE;
2982 }
5443e65e 2983
2984}
2985
75bd7f81 2986//_____________________________________________________________________________
2987void AliTRDtracker::AliTRDpropagationLayer
2988 ::SetHole(Double_t Zmax, Double_t Ymax, Double_t rho
2989 , Double_t radLength, Double_t Yc, Double_t Zc)
5443e65e 2990{
2991 //
2992 // Sets hole in the layer
2993 //
75bd7f81 2994
4551ea7c 2995 fHole = kTRUE;
2996 fHoleZc = Zc;
5443e65e 2997 fHoleZmax = Zmax;
4551ea7c 2998 fHoleYc = Yc;
5443e65e 2999 fHoleYmax = Ymax;
4551ea7c 3000 fHoleRho = rho;
3001 fHoleX0 = radLength;
75bd7f81 3002
5443e65e 3003}
46d29e70 3004
75bd7f81 3005//_____________________________________________________________________________
3006AliTRDtracker::AliTRDtrackingSector
ad670fba 3007 ::AliTRDtrackingSector(AliTRDgeometry *geo, Int_t gs)
3008 :fN(0)
3009 ,fGeom(geo)
3010 ,fGeomSector(gs)
5443e65e 3011{
3012 //
3013 // AliTRDtrackingSector Constructor
3014 //
75bd7f81 3015
7e448bcc 3016 AliTRDpadPlane *padPlane = 0;
3017 AliTRDpropagationLayer *ppl = 0;
a5cadd36 3018
ad670fba 3019 // Get holes description from geometry
3c625a9b 3020 Bool_t holes[AliTRDgeometry::kNcham];
ad670fba 3021 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3c625a9b 3022 holes[icham] = fGeom->IsHole(0,icham,gs);
3c625a9b 3023 }
3c625a9b 3024
ad670fba 3025 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
3026 fTimeBinIndex[i] = -1;
3027 }
5443e65e 3028
ad670fba 3029 Double_t x;
3030 Double_t dx;
3031 Double_t rho;
3032 Double_t radLength;
5443e65e 3033
ad670fba 3034 // Add layers for each of the planes
3035 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
a305677e 3036 //Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
5443e65e 3037
ad670fba 3038 const Int_t kNchambers = AliTRDgeometry::Ncham();
a305677e 3039 Int_t tbIndex;
ad670fba 3040 Double_t ymax = 0;
3041 Double_t ymaxsensitive = 0;
3042 Double_t *zc = new Double_t[kNchambers];
3043 Double_t *zmax = new Double_t[kNchambers];
3c625a9b 3044 Double_t *zmaxsensitive = new Double_t[kNchambers];
3551db50 3045
ad670fba 3046 for (Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
3047
3048 ymax = fGeom->GetChamberWidth(plane) / 2.0;
f162af62 3049 padPlane = fGeom->GetPadPlane(plane,0);
ad670fba 3050 ymaxsensitive = (padPlane->GetColSize(1) * padPlane->GetNcols() - 4.0) / 2.0;
3051
3052 for (Int_t ch = 0; ch < kNchambers; ch++) {
3053 zmax[ch] = fGeom->GetChamberLength(plane,ch) / 2.0;
3054 Float_t pad = padPlane->GetRowSize(1);
f162af62 3055 Float_t row0 = fGeom->GetRow0(plane,ch,0);
3056 Int_t nPads = fGeom->GetRowMax(plane,ch,0);
ad670fba 3057 zmaxsensitive[ch] = Float_t(nPads) * pad / 2.0;
3058 zc[ch] = -(pad * nPads) / 2.0 + row0;
5443e65e 3059 }
3060
ad670fba 3061 dx = AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
f162af62 3062 / AliTRDCommonParam::Instance()->GetSamplingFrequency();
ad670fba 3063 rho = 0.00295 * 0.85; //????
3064 radLength = 11.0;
5443e65e 3065
3551db50 3066 Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(plane);
a305677e 3067 //Double_t xbottom = x0 - dxDrift;
ad670fba 3068 //Double_t xtop = x0 + dxAmp;
3069
59393e34 3070 Int_t nTimeBins = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
ad670fba 3071 for (Int_t iTime = 0; iTime < nTimeBins; iTime++) {
3072
3073 Double_t xlayer = iTime * dx - dxAmp;
3074 //if (xlayer<0) xlayer = dxAmp / 2.0;
59393e34 3075 x = x0 - xlayer;
ad670fba 3076
3077 tbIndex = CookTimeBinIndex(plane,iTime);
3078 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex,plane);
3c625a9b 3079 ppl->SetYmax(ymax,ymaxsensitive);
ad670fba 3080 ppl->SetZ(zc,zmax,zmaxsensitive);
3c625a9b 3081 ppl->SetHoles(holes);
59393e34 3082 InsertLayer(ppl);
ad670fba 3083
5443e65e 3084 }
ad670fba 3085
5443e65e 3086 }
3087
5443e65e 3088 MapTimeBinLayers();
ad670fba 3089
029cd327 3090 delete [] zc;
3091 delete [] zmax;
4f1c04d3 3092 delete [] zmaxsensitive;
5443e65e 3093
3094}
3095
ad670fba 3096//_____________________________________________________________________________
3097AliTRDtracker::AliTRDtrackingSector
3098 ::AliTRDtrackingSector(const AliTRDtrackingSector &/*t*/)
3099 :fN(0)
3100 ,fGeom(0)
3101 ,fGeomSector(0)
3102{
3103 //
3104 // Copy constructor
3105 //
3106
3107}
3108
75bd7f81 3109//_____________________________________________________________________________
3110Int_t AliTRDtracker::AliTRDtrackingSector
3111 ::CookTimeBinIndex(Int_t plane, Int_t localTB) const
5443e65e 3112{
3113 //
6c94f330 3114 // depending on the digitization parameters calculates "global"
029cd327 3115 // time bin index for timebin <localTB> in plane <plane>
5443e65e 3116 //
59393e34 3117 //
75bd7f81 3118
59393e34 3119 Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
4551ea7c 3120 Int_t gtb = (plane+1) * tbPerPlane - localTB - 1;
3121 if (localTB < 0) {
3122 return -1;
3123 }
3124 if (gtb < 0) {
3125 return -1;
3126 }
75bd7f81 3127
5443e65e 3128 return gtb;
5443e65e 3129
75bd7f81 3130}
5443e65e 3131
75bd7f81 3132//_____________________________________________________________________________
3133void AliTRDtracker::AliTRDtrackingSector
3134 ::MapTimeBinLayers()
5443e65e 3135{
3136 //
3137 // For all sensitive time bins sets corresponding layer index
3138 // in the array fTimeBins
3139 //
3140
3141 Int_t index;
3142
4551ea7c 3143 for (Int_t i = 0; i < fN; i++) {
3144
5443e65e 3145 index = fLayers[i]->GetTimeBinIndex();
6c94f330 3146
4551ea7c 3147 if (index < 0) {
3148 continue;
3149 }
3150 if (index >= (Int_t) kMaxTimeBinIndex) {
3151 //AliWarning(Form("Index %d exceeds allowed maximum of %d!\n"
3152 // ,index,kMaxTimeBinIndex-1));
5443e65e 3153 continue;
3154 }
4551ea7c 3155
5443e65e 3156 fTimeBinIndex[index] = i;
4551ea7c 3157
5443e65e 3158 }
5443e65e 3159
75bd7f81 3160}
5443e65e 3161
75bd7f81 3162//_____________________________________________________________________________
3163Int_t AliTRDtracker::AliTRDtrackingSector
3164 ::GetLayerNumber(Double_t x) const
5443e65e 3165{
3166 //
3167 // Returns the number of time bin which in radial position is closest to <x>
3168 //
3169
4551ea7c 3170 if (x >= fLayers[fN-1]->GetX()) {
3171 return fN - 1;
3172 }
3173 if (x <= fLayers[ 0]->GetX()) {
3174 return 0;
3175 }
3176
3177 Int_t b = 0;
3178 Int_t e = fN - 1;
3179 Int_t m = (b + e) / 2;
ad670fba 3180
4551ea7c 3181 for ( ; b < e; m = (b + e) / 2) {
3182 if (x > fLayers[m]->GetX()) {
3183 b = m + 1;
3184 }
3185 else {
3186 e = m;
3187 }
5443e65e 3188 }
75bd7f81 3189
4551ea7c 3190 if (TMath::Abs(x - fLayers[m]->GetX()) > TMath::Abs(x - fLayers[m+1]->GetX())) {
3191 return m + 1;
3192 }
3193 else {
3194 return m;
3195 }
5443e65e 3196
3197}
3198
75bd7f81 3199//_____________________________________________________________________________
3200Int_t AliTRDtracker::AliTRDtrackingSector
3201 ::GetInnerTimeBin() const
5443e65e 3202{
3203 //
3204 // Returns number of the innermost SENSITIVE propagation layer
3205 //
3206
3207 return GetLayerNumber(0);
5443e65e 3208
75bd7f81 3209}
5443e65e 3210
75bd7f81 3211//_____________________________________________________________________________
3212Int_t AliTRDtracker::AliTRDtrackingSector
3213 ::GetOuterTimeBin() const
5443e65e 3214{
3215 //
3216 // Returns number of the outermost SENSITIVE time bin
3217 //
3218
3219 return GetLayerNumber(GetNumberOfTimeBins() - 1);
46d29e70 3220
75bd7f81 3221}
5443e65e 3222
75bd7f81 3223//_____________________________________________________________________________
3224Int_t AliTRDtracker::AliTRDtrackingSector
3225 ::GetNumberOfTimeBins() const
5443e65e 3226{
3227 //
3228 // Returns number of SENSITIVE time bins
3229 //
3230
4551ea7c 3231 Int_t tb;
3232 Int_t layer;
3233
3234 for (tb = kMaxTimeBinIndex - 1; tb >= 0; tb--) {
5443e65e 3235 layer = GetLayerNumber(tb);
4551ea7c 3236 if (layer >= 0) {
3237 break;
3238 }
5443e65e 3239 }
75bd7f81 3240
4551ea7c 3241 return tb + 1;
5443e65e 3242
75bd7f81 3243}
5443e65e 3244
75bd7f81 3245//_____________________________________________________________________________
3246void AliTRDtracker::AliTRDtrackingSector
4551ea7c 3247 ::InsertLayer(AliTRDpropagationLayer *pl)
5443e65e 3248{
3249 //
3250 // Insert layer <pl> in fLayers array.
3251 // Layers are sorted according to X coordinate.
75bd7f81 3252 //
5443e65e 3253
4551ea7c 3254 if (fN == ((Int_t) kMaxLayersPerSector)) {
3255 //AliWarning("Too many layers !\n");
3256 return;
3257 }
3258
3259 if (fN == 0) {
3260 fLayers[fN++] = pl;
ad670fba 3261 return;
3262 }
3263
4551ea7c 3264 Int_t i = Find(pl->GetX());
3265
3266 memmove(fLayers+i+1,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
3267
3268 fLayers[i] = pl;
3269 fN++;
5443e65e 3270
3271}
3272
75bd7f81 3273//_____________________________________________________________________________
3274Int_t AliTRDtracker::AliTRDtrackingSector
3275 ::Find(Double_t x) const
5443e65e 3276{
3277 //
3278 // Returns index of the propagation layer nearest to X
3279 //
3280
4551ea7c 3281 if (x <= fLayers[0]->GetX()) {
3282 return 0;
3283 }
3284
3285 if (x > fLayers[fN-1]->GetX()) {
3286 return fN;
3287 }
3288
3289 Int_t b = 0;
3290 Int_t e = fN-1;
3291 Int_t m = (b + e) / 2;
3292
3293 for (; b < e; m = (b + e) / 2) {
3294 if (x > fLayers[m]->GetX()) {
3295 b = m + 1;
3296 }
3297 else {
3298 e = m;
3299 }
5443e65e 3300 }
7ad19338 3301
75bd7f81 3302 return m;
7ad19338 3303
75bd7f81 3304}
7ad19338 3305
75bd7f81 3306//_____________________________________________________________________________
3307void AliTRDtracker::AliTRDpropagationLayer
4551ea7c 3308 ::SetZ(Double_t *center, Double_t *w, Double_t *wsensitive )
3c625a9b 3309{
3310 //
6c94f330 3311 // set centers and the width of sectors
75bd7f81 3312 //
3313
4551ea7c 3314 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3315 fZc[icham] = center[icham];
3316 fZmax[icham] = w[icham];
3c625a9b 3317 fZmaxSensitive[icham] = wsensitive[icham];
3c625a9b 3318 }
75bd7f81 3319
3c625a9b 3320}
5443e65e 3321
75bd7f81 3322//_____________________________________________________________________________
3c625a9b 3323void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
3324{
3325 //
6c94f330 3326 // set centers and the width of sectors
75bd7f81 3327 //
3328
3c625a9b 3329 fHole = kFALSE;
4551ea7c 3330
3331 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3c625a9b 3332 fIsHole[icham] = holes[icham];
4551ea7c 3333 if (holes[icham]) {
3334 fHole = kTRUE;
3335 }
3c625a9b 3336 }
5443e65e 3337
75bd7f81 3338}
5443e65e 3339
75bd7f81 3340//_____________________________________________________________________________
3341void AliTRDtracker::AliTRDpropagationLayer
4551ea7c 3342 ::InsertCluster(AliTRDcluster *c, UInt_t index)
75bd7f81 3343{
3344 //
3345 // Insert cluster in cluster array.
3346 // Clusters are sorted according to Y coordinate.
3347 //
5443e65e 3348
4551ea7c 3349 if (fTimeBinIndex < 0) {
3350 //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");
3351 return;
3352 }
3353
3354 if (fN == (Int_t) kMaxClusterPerTimeBin) {
3355 //AliWarning("Too many clusters !\n");
5443e65e 3356 return;
3357 }
ad670fba 3358
4551ea7c 3359 if (fN == 0) {
3360 fIndex[0] = index;
3361 fClusters[fN++] = c;
ad670fba 3362 return;
3363 }
4551ea7c 3364
3365 Int_t i = Find(c->GetY());
3366 memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
3367 memmove(fIndex +i+1,fIndex +i,(fN-i)*sizeof(UInt_t));
3368 fIndex[i] = index;
3369 fClusters[i] = c;
3370 fN++;
5443e65e 3371
75bd7f81 3372}
5443e65e 3373
75bd7f81 3374//_____________________________________________________________________________
3375Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const
3376{
3377 //
3378 // Returns index of the cluster nearest in Y
3379 //
5443e65e 3380
4551ea7c 3381 if (fN <= 0) {
3382 return 0;
3383 }
3384 if (y <= fClusters[0]->GetY()) {
3385 return 0;
3386 }
3387 if (y > fClusters[fN-1]->GetY()) {
3388 return fN;
3389 }
3390
3391 Int_t b = 0;
3392 Int_t e = fN - 1;
3393 Int_t m = (b + e) / 2;
3394
3395 for ( ; b < e; m = (b + e) / 2) {
3396 if (y > fClusters[m]->GetY()) {
3397 b = m + 1;
3398 }
3399 else {
3400 e = m;
3401 }
5443e65e 3402 }
75bd7f81 3403
5443e65e 3404 return m;
75bd7f81 3405
5443e65e 3406}
3407
75bd7f81 3408//_____________________________________________________________________________
3409Int_t AliTRDtracker::AliTRDpropagationLayer
3410 ::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad
3411 , Float_t maxroadz) const
7ad19338 3412{
3413 //
3414 // Returns index of the cluster nearest to the given y,z
3415 //
75bd7f81 3416
4551ea7c 3417 Int_t index = -1;
3418 Int_t maxn = fN;
69b55c55 3419 Float_t mindist = maxroad;
4551ea7c 3420
3421 for (Int_t i = Find(y-maxroad); i < maxn; i++) {
3422 AliTRDcluster *c = (AliTRDcluster *) (fClusters[i]);
69b55c55 3423 Float_t ycl = c->GetY();
4551ea7c 3424 if (ycl > (y + maxroad)) {
3425 break;
3426 }
3427 if (TMath::Abs(c->GetZ() - z) > maxroadz) {
3428 continue;
3429 }
3430 if (TMath::Abs(ycl - y) < mindist) {
3431 mindist = TMath::Abs(ycl - y);
3432 index = fIndex[i];
3433 }
6c94f330 3434 }
75bd7f81 3435
7ad19338 3436 return index;
7ad19338 3437
75bd7f81 3438}
7ad19338 3439
75bd7f81 3440//_____________________________________________________________________________
4551ea7c 3441Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c)
75bd7f81 3442{
3443 //
3444 // Returns correction factor for tilted pads geometry
3445 //
5443e65e 3446
4551ea7c 3447 Int_t det = c->GetDetector();
3448 Int_t plane = fGeom->GetPlane(det);
f162af62 3449 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(plane,0);
4551ea7c 3450 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
b8dc2353 3451
4551ea7c 3452 if (fNoTilt) {
3453 h01 = 0;
3454 }
75bd7f81 3455
fd621f36 3456 return h01;
5443e65e 3457
75bd7f81 3458}
c630aafd 3459
3460
75bd7f81 3461//_____________________________________________________________________________
3462Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
4551ea7c 3463 , AliTRDtrack *track
3464 , Int_t *clusters, AliTRDtracklet &tracklet)
4f1c04d3 3465{
6c94f330 3466 //
4f1c04d3 3467 //
75bd7f81 3468 // Try to find nearest clusters to the track in timebins from t0 to t1
4f1c04d3 3469 //
4551ea7c 3470 // Correction coeficients - depend on TRD parameters - to be changed accordingly
3471 //
3472
3473 Double_t x[100];
3474 Double_t yt[100];
3475 Double_t zt[100];
3476 Double_t xmean = 0.0; // Reference x
3477 Double_t dz[10][100];
3478 Double_t dy[10][100];
3479 Float_t zmean[100];
3480 Float_t nmean[100];
3481 Int_t clfound = 0;
3482 Int_t indexes[10][100]; // Indexes of the clusters in the road
3483 Int_t best[10][100]; // Index of best matching cluster
3484 AliTRDcluster *cl[10][100]; // Pointers to the clusters in the road
3485
3486 for (Int_t it = 0; it < 100; it++) {
3487 x[it] = 0.0;
3488 yt[it] = 0.0;
3489 zt[it] = 0.0;
3490 clusters[it] = -2;
3491 zmean[it] = 0.0;
3492 nmean[it] = 0.0;
3493 for (Int_t ih = 0; ih < 10;ih++) {
3494 indexes[ih][it] = -2; // Reset indexes1
3495 cl[ih][it] = 0;
3496 dz[ih][it] = -100.0;
3497 dy[ih][it] = -100.0;
3498 best[ih][it] = 0;
3499 }
3500 }
ad670fba 3501
4551ea7c 3502 Double_t x0 = track->GetX();
3503 Double_t sigmaz = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
3504 Int_t nall = 0;
3505 Int_t nfound = 0;
3506 Double_t h01 = 0.0;
3507 Int_t plane = -1;
3508 Int_t detector = -1;
3509 Float_t padlength = 0.0;
3510 AliTRDtrack track2(* track);
3511 Float_t snpy = track->GetSnp();
3512 Float_t tany = TMath::Sqrt(snpy*snpy / (1.0 - snpy*snpy));
3513 if (snpy < 0.0) {
3514 tany *= -1.0;
3515 }
ad670fba 3516
6c23ffed 3517 Double_t sy2 = ExpectedSigmaY2(x0,track->GetTgl(),track->GetSignedPt());
4551ea7c 3518 Double_t sz2 = ExpectedSigmaZ2(x0,track->GetTgl());
3519 Double_t road = 15.0 * TMath::Sqrt(track->GetSigmaY2() + sy2);
3520 if (road > 6.0) {
3521 road = 6.0;
3522 }
7e448bcc 3523 //road = 20.0;
4551ea7c 3524
3525 for (Int_t it = 0; it < t1-t0; it++) {
3526
3527 Double_t maxChi2[2] = { fgkMaxChi2, fgkMaxChi2 };
3528 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(it+t0));
3529 if (timeBin == 0) {
3530 continue; // No indexes1
6c94f330 3531 }
4551ea7c 3532
3533 Int_t maxn = timeBin;
3534 x[it] = timeBin.GetX();
7ad19338 3535 track2.PropagateTo(x[it]);
6c94f330 3536 yt[it] = track2.GetY();
3537 zt[it] = track2.GetZ();
7ad19338 3538
4551ea7c 3539 Double_t y = yt[it];
3540 Double_t z = zt[it];
3541 Double_t chi2 = 1000000.0;
4f1c04d3 3542 nall++;
4551ea7c 3543
4f1c04d3 3544 //
4551ea7c 3545 // Find 2 nearest cluster at given time bin
6c94f330 3546 //
7e448bcc 3547 int checkPoint[4] = {0,0,0,0};
3548 double minY = 123456789;
3549 double minD[2] = {1,1};
3550
4551ea7c 3551 for (Int_t i = timeBin.Find(y - road); i < maxn; i++) {
7e448bcc 3552 //for (Int_t i = 0; i < maxn; i++) {
4551ea7c 3553
3554 AliTRDcluster *c = (AliTRDcluster *) (timeBin[i]);
7ad19338 3555 h01 = GetTiltFactor(c);
4551ea7c 3556 if (plane < 0) {
c6f438c0 3557 Int_t det = c->GetDetector();
4551ea7c 3558 plane = fGeom->GetPlane(det);
3559 padlength = TMath::Sqrt(c->GetSigmaZ2() * 12.0);
3560 }
7e448bcc 3561
4551ea7c 3562 //if (c->GetLocalTimeBin()==0) continue;
3563 if (c->GetY() > (y + road)) {
3564 break;
3565 }
7e448bcc 3566
3567 fHDeltaX->Fill(c->GetX() - x[it]);
3568 //printf("%f\t%f\t%f \n", c->GetX(), x[it], c->GetX()-x[it]);
3569
3570 if (TMath::Abs(c->GetY()-y) < TMath::Abs(minY)) {
3571 minY = c->GetY()-y;
3572 minD[0] = c->GetY()-y;
3573 minD[1] = c->GetZ()-z;
3574 }
3575
3576 checkPoint[0]++;
3577
3578 fHMinZ->Fill(c->GetZ() - z);
3579 if ((c->GetZ() - z) * (c->GetZ() - z) > 2 * (12.0 * sz2)) {
4551ea7c 3580 continue;
7ad19338 3581 }
7e448bcc 3582 checkPoint[1]++;
7ad19338 3583
4551ea7c 3584 Double_t dist = TMath::Abs(c->GetZ() - z);
7e448bcc 3585 if (dist > (0.5 * padlength + 6.0 * sigmaz)) { // 0.5
4551ea7c 3586 continue; // 6 sigma boundary cut
3587 }
7e448bcc 3588 checkPoint[2]++;
3589
4551ea7c 3590 Double_t cost = 0.0;
3591 // Sigma boundary cost function
3592 if (dist> (0.5 * padlength - sigmaz)){
3593 cost = (dist - 0.5*padlength) / (2.0 * sigmaz);
3594 if (cost > -1) {
3595 cost = (cost + 1.0) * (cost + 1.0);
3596 }
3597 else {
3598 cost = 0.0;
3599 }
3600 }
3601 //Int_t label = TMath::Abs(track->GetLabel());
3602 //if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue;
3603 chi2 = track2.GetPredictedChi2(c,h01) + cost;
7ad19338 3604 clfound++;
6c94f330 3605
4551ea7c 3606 if (chi2 > maxChi2[1]) {
3607 continue;
3608 }
7e448bcc 3609 checkPoint[3]++;
3610
4551ea7c 3611 detector = c->GetDetector();
3612 // Store the clusters in the road
3613 for (Int_t ih = 2; ih < 9; ih++) {
3614 if (cl[ih][it] == 0) {
3615 cl[ih][it] = c;
3616 indexes[ih][it] = timeBin.GetIndex(i); // Index - 9 - reserved for outliers
7ad19338 3617 break;
3618 }
4f1c04d3 3619 }
4551ea7c 3620
3621 if (chi2 < maxChi2[0]) {
7ad19338 3622 maxChi2[1] = maxChi2[0];
3623 maxChi2[0] = chi2;
3624 indexes[1][it] = indexes[0][it];
3625 cl[1][it] = cl[0][it];
3626 indexes[0][it] = timeBin.GetIndex(i);
3627 cl[0][it] = c;
3628 continue;
3629 }
4551ea7c 3630 maxChi2[1] = chi2;
3631 cl[1][it] = c;
3632 indexes[1][it] = timeBin.GetIndex(i);
3633
3634 }
7e448bcc 3635
3636 for(int iCheckPoint = 0; iCheckPoint<4; iCheckPoint++)
3637 fHFindCl[iCheckPoint]->Fill(checkPoint[iCheckPoint]);
3638
3639 if (checkPoint[3]) {
6c23ffed 3640 if (track->GetSignedPt() > 0) fHMinYPos->Fill(minY);
7e448bcc 3641 else fHMinYNeg->Fill(minY);
3642
3643 fHMinD->Fill(minD[0], minD[1]);
3644 }
4551ea7c 3645
3646 if (cl[0][it]) {
7ad19338 3647 nfound++;
3648 xmean += x[it];
3649 }
4551ea7c 3650
4f1c04d3 3651 }
4551ea7c 3652
3653 if (nfound < 4) {
3654 return 0;
3655 }
3656 xmean /= Float_t(nfound); // Middle x
3657 track2.PropagateTo(xmean); // Propagate track to the center
3658
4f1c04d3 3659 //
4551ea7c 3660 // Choose one of the variants
7ad19338 3661 //
4551ea7c 3662 Int_t changes[10];
3663 Float_t sumz = 0.0;
3664 Float_t sum = 0.0;
3665 Double_t sumdy = 0.0;
3666 Double_t sumdy2 = 0.0;
3667 Double_t sumx = 0.0;
3668 Double_t sumxy = 0.0;
3669 Double_t sumx2 = 0.0;
3670 Double_t mpads = 0.0;
3671
3672 Int_t ngood[10];
3673 Int_t nbad[10];
3674
7ad19338 3675 Double_t meanz[10];
4551ea7c 3676 Double_t moffset[10]; // Mean offset
3677 Double_t mean[10]; // Mean value
3678 Double_t angle[10]; // Angle
3679
3680 Double_t smoffset[10]; // Sigma of mean offset
3681 Double_t smean[10]; // Sigma of mean value
3682 Double_t sangle[10]; // Sigma of angle
3683 Double_t smeanangle[10]; // Correlation
3684
7ad19338 3685 Double_t sigmas[10];
4551ea7c 3686 Double_t tchi2s[10]; // Chi2s for tracklet
ad670fba 3687
4551ea7c 3688 for (Int_t it = 0; it < 10; it++) {
ad670fba 3689
4551ea7c 3690 ngood[it] = 0;
3691 nbad[it] = 0;
3692
3693 meanz[it] = 0.0;
3694 moffset[it] = 0.0; // Mean offset
3695 mean[it] = 0.0; // Mean value
3696 angle[it] = 0.0; // Angle
3697
7e448bcc 3698 smoffset[it] = 1.0e5; // Sigma of mean offset
3699 smean[it] = 1.0e5; // Sigma of mean value
3700 sangle[it] = 1.0e5; // Sigma of angle
4551ea7c 3701 smeanangle[it] = 0.0; // Correlation
3702
7e448bcc 3703 sigmas[it] = 1.0e5;
3704 tchi2s[it] = 1.0e5; // Chi2s for tracklet
75bd7f81 3705
3706 }
3707
7ad19338 3708 //
4551ea7c 3709 // Calculate zmean
7ad19338 3710 //
4551ea7c 3711 for (Int_t it = 0; it < t1 - t0; it++) {
3712 if (!cl[0][it]) {
3713 continue;
3714 }
3715 for (Int_t dt = -3; dt <= 3; dt++) {
3716 if (it+dt < 0) {
3717 continue;
3718 }
3719 if (it+dt > t1-t0) {
3720 continue;
3721 }
3722 if (!cl[0][it+dt]) {
3723 continue;
3724 }
3725 zmean[it] += cl[0][it+dt]->GetZ();
3726 nmean[it] += 1.0;
7ad19338 3727 }
4551ea7c 3728 zmean[it] /= nmean[it];
7ad19338 3729 }
4551ea7c 3730
3731 for (Int_t it = 0; it < t1 - t0; it++) {
3732
3733 best[0][it] = 0;
3734
3735 for (Int_t ih = 0; ih < 10; ih++) {
3736 dz[ih][it] = -100.0;
3737 dy[ih][it] = -100.0;
3738 if (!cl[ih][it]) {
3739 continue;
3740 }
3741 Double_t xcluster = cl[ih][it]->GetX();
3742 Double_t ytrack;
3743 Double_t ztrack;
3744 track2.GetProlongation(xcluster,ytrack,ztrack );
3745 dz[ih][it] = cl[ih][it]->GetZ()- ztrack; // Calculate distance from track in z
3746 dy[ih][it] = cl[ih][it]->GetY() + dz[ih][it]*h01 - ytrack; // and in y
7ad19338 3747 }
4551ea7c 3748
3749 // Minimize changes
3750 if (!cl[0][it]) {
3751 continue;
3752 }
3753 if ((TMath::Abs(cl[0][it]->GetZ()-zmean[it]) > padlength * 0.8) &&
3754 (cl[1][it])) {
3755 if (TMath::Abs(cl[1][it]->GetZ()-zmean[it]) < padlength * 0.5) {
3756 best[0][it] = 1;
4f1c04d3 3757 }
4551ea7c 3758 }
3759
7ad19338 3760 }
4551ea7c 3761
7ad19338 3762 //
4551ea7c 3763 // Iterative choice of "best path"
6c94f330 3764 //
4551ea7c 3765 Int_t label = TMath::Abs(track->GetLabel());
3766 Int_t bestiter = 0;
3767
3768 for (Int_t iter = 0; iter < 9; iter++) {
3769
3770 changes[iter] = 0;
3771 sumz = 0;
3772 sum = 0;
3773 sumdy = 0;
3774 sumdy2 = 0;
3775 sumx = 0;
3776 sumx2 = 0;
3777 sumxy = 0;
3778 mpads = 0;
3779 ngood[iter] = 0;
3780 nbad[iter] = 0;
3781
3782 // Linear fit
3783 for (Int_t it = 0; it < t1 - t0; it++) {
3784
3785 if (!cl[best[iter][it]][it]) {
3786 continue;
3787 }
3788
3789 // Calculates pad-row changes
3790 Double_t zbefore = cl[best[iter][it]][it]->GetZ();
3791 Double_t zafter = cl[best[iter][it]][it]->GetZ();
3792 for (Int_t itd = it - 1; itd >= 0; itd--) {
7ad19338 3793 if (cl[best[iter][itd]][itd]) {
4551ea7c 3794 zbefore = cl[best[iter][itd]][itd]->GetZ();
7ad19338 3795 break;
3796 }
3797 }
4551ea7c 3798 for (Int_t itd = it + 1; itd < t1 - t0; itd++) {
7ad19338 3799 if (cl[best[iter][itd]][itd]) {
4551ea7c 3800 zafter = cl[best[iter][itd]][itd]->GetZ();
7ad19338 3801 break;
3802 }
3803 }
4551ea7c 3804 if ((TMath::Abs(cl[best[iter][it]][it]->GetZ()-zbefore) > 0.1) &&
3805 (TMath::Abs(cl[best[iter][it]][it]->GetZ()- zafter) > 0.1)) {
3806 changes[iter]++;
3807 }
3808
3809 Double_t dx = x[it]-xmean; // Distance to reference x
3810 sumz += cl[best[iter][it]][it]->GetZ();
4f1c04d3 3811 sum++;
4551ea7c 3812 sumdy += dy[best[iter][it]][it];
3813 sumdy2 += dy[best[iter][it]][it]*dy[best[iter][it]][it];
3814 sumx += dx;
3815 sumx2 += dx*dx;
7ad19338 3816 sumxy += dx*dy[best[iter][it]][it];
4551ea7c 3817 mpads += cl[best[iter][it]][it]->GetNPads();
3818 if ((cl[best[iter][it]][it]->GetLabel(0) == label) ||
3819 (cl[best[iter][it]][it]->GetLabel(1) == label) ||
3820 (cl[best[iter][it]][it]->GetLabel(2) == label)) {
7ad19338 3821 ngood[iter]++;
4f1c04d3 3822 }
4551ea7c 3823 else {
7ad19338 3824 nbad[iter]++;
4f1c04d3 3825 }
4551ea7c 3826
4f1c04d3 3827 }
4551ea7c 3828
7ad19338 3829 //
6c94f330 3830 // calculates line parameters
7ad19338 3831 //
4551ea7c 3832 Double_t det = sum*sumx2 - sumx*sumx;
3833 angle[iter] = (sum*sumxy - sumx*sumdy) / det;
3834 mean[iter] = (sumx2*sumdy - sumx*sumxy) / det;
3835 meanz[iter] = sumz / sum;
3836 moffset[iter] = sumdy / sum;
3837 mpads /= sum; // Mean number of pads
3838
3839 Double_t sigma2 = 0.0; // Normalized residuals - for line fit
3840 Double_t sigma1 = 0.0; // Normalized residuals - constant fit
3841
3842 for (Int_t it = 0; it < t1 - t0; it++) {
3843 if (!cl[best[iter][it]][it]) {
3844 continue;
3845 }
3846 Double_t dx = x[it] - xmean;
3847 Double_t ytr = mean[iter] + angle[iter] * dx;
3848 sigma2 += (dy[best[iter][it]][it] - ytr)
3849 * (dy[best[iter][it]][it] - ytr);
3850 sigma1 += (dy[best[iter][it]][it] - moffset[iter])
3851 * (dy[best[iter][it]][it] - moffset[iter]);
7ad19338 3852 sum++;
4f1c04d3 3853 }
4551ea7c 3854 sigma2 /= (sum - 2); // Normalized residuals
3855 sigma1 /= (sum - 1); // Normalized residuals
3856 smean[iter] = sigma2 * (sumx2 / det); // Estimated error2 of mean
3857 sangle[iter] = sigma2 * ( sum / det); // Estimated error2 of angle
3858 smeanangle[iter] = sigma2 * (-sumx / det); // Correlation
3859 sigmas[iter] = TMath::Sqrt(sigma1);
3860 smoffset[iter] = (sigma1 / sum) + 0.01*0.01; // Sigma of mean offset + unisochronity sigma
3861
6c94f330 3862 //
4551ea7c 3863 // Iterative choice of "better path"
6c94f330 3864 //
4551ea7c 3865 for (Int_t it = 0; it < t1 - t0; it++) {
3866
3867 if (!cl[best[iter][it]][it]) {
3868 continue;
3869 }
3870
3871 // Add unisochronity + angular effect contribution
3872 Double_t sigmatr2 = smoffset[iter] + 0.5*tany*tany;
3873 Double_t sweight = 1.0/sigmatr2 + 1.0/track->GetSigmaY2();
3874 Double_t weighty = (moffset[iter] / sigmatr2) / sweight; // Weighted mean
3875 Double_t sigmacl = TMath::Sqrt(sigma1*sigma1 + track->GetSigmaY2());
3876 Double_t mindist = 100000.0;
3877 Int_t ihbest = 0;
3878
3879 for (Int_t ih = 0; ih < 10; ih++) {
3880 if (!cl[ih][it]) {
3881 break;
3882 }
3883 Double_t dist2 = (dy[ih][it] - weighty) / sigmacl;
3884 dist2 *= dist2; // Chi2 distance
3885 if (dist2 < mindist) {
7ad19338 3886 mindist = dist2;
4551ea7c 3887 ihbest = ih;
7ad19338 3888 }
3889 }
4551ea7c 3890 best[iter+1][it] = ihbest;
4f1c04d3 3891 }
4551ea7c 3892
4f1c04d3 3893 //
4551ea7c 3894 // Update best hypothesy if better chi2 according tracklet position and angle
7ad19338 3895 //
69b55c55 3896 Double_t sy2 = smean[iter] + track->GetSigmaY2();
4551ea7c 3897 Double_t sa2 = sangle[iter] + track->GetSigmaSnp2(); // track->fCee;
3898 Double_t say = track->GetSigmaSnpY(); // track->fCey;
3899 //Double_t chi20 = mean[bestiter]*mean[bestiter ] / sy2+angle[bestiter]*angle[bestiter]/sa2;
3900 //Double_t chi21 = mean[iter]*mean[iter] / sy2+angle[iter]*angle[iter]/sa2;
6c94f330 3901
4551ea7c 3902 Double_t detchi = sy2*sa2 - say*say;
3903 Double_t invers[3] = {sa2/detchi,sy2/detchi,-say/detchi}; // Inverse value of covariance matrix
4f1c04d3 3904
4551ea7c 3905 Double_t chi20 = mean[bestiter] * mean[bestiter] * invers[0]
3906 + angle[bestiter] * angle[bestiter] * invers[1]
3907 + 2.0 * mean[bestiter] * angle[bestiter] * invers[2];
3908 Double_t chi21 = mean[iter] * mean[iter] * invers[0]
3909 + angle[iter] * angle[iter] * invers[1]
3910 + 2.0 * mean[iter] * angle[iter] * invers[2];
3911 tchi2s[iter] = chi21;
3912
3913 if ((changes[iter] <= changes[bestiter]) &&
3914 (chi21 < chi20)) {
3915 bestiter = iter;
7ad19338 3916 }
4551ea7c 3917
7ad19338 3918 }
4551ea7c 3919
7ad19338 3920 //
4551ea7c 3921 // Set clusters
7ad19338 3922 //
4551ea7c 3923 Double_t sigma2 = sigmas[0]; // Choose as sigma from 0 iteration
3924 Short_t maxpos = -1;
3925 Float_t maxcharge = 0.0;
3926 Short_t maxpos4 = -1;
3927 Float_t maxcharge4 = 0.0;
3928 Short_t maxpos5 = -1;
3929 Float_t maxcharge5 = 0.0;
8979685e 3930
7ad19338 3931 //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
3932 //if (tchi2s[bestiter]>25.) sigma2=1000.; // dont'accept
3933
4e009ce4 3934 Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
3935 ,-AliTracker::GetBz()*0.1);
4551ea7c 3936 Double_t expectederr = sigma2*sigma2 + 0.01*0.01;
3937 if (mpads > 3.5) {
3938 expectederr += (mpads - 3.5) * 0.04;
3939 }
3940 if (changes[bestiter] > 1) {
3941 expectederr += changes[bestiter] * 0.01;
3942 }
3943 expectederr += (0.03 * (tany-exB)*(tany-exB)) * 15.0;
3944 //if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.;
7ad19338 3945 //expectederr+=10000;
4551ea7c 3946
3947 for (Int_t it = 0; it < t1 - t0; it++) {
3948
3949 if (!cl[best[bestiter][it]][it]) {
3950 continue;
69b55c55 3951 }
4551ea7c 3952
3953 cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); // Set cluster error
3954 if (!cl[best[bestiter][it]][it]->IsUsed()) {
3955 cl[best[bestiter][it]][it]->SetY(cl[best[bestiter][it]][it]->GetY());
3956 //cl[best[bestiter][it]][it]->Use();
3957 }
3958
3959 // Time bins with maximal charge
3960 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge) {
69b55c55 3961 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3962 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 3963 }
6c94f330 3964
4551ea7c 3965 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
3966 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
69b55c55 3967 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3968 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 3969 }
3970 }
4551ea7c 3971
3972 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
3973 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
69b55c55 3974 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3975 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 3976 }
7ad19338 3977 }
4551ea7c 3978
3979 // Time bins with maximal charge
3980 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge) {
8979685e 3981 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3982 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 3983 }
6c94f330 3984
4551ea7c 3985 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
3986 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
8979685e 3987 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3988 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 3989 }
3990 }
4551ea7c 3991
3992 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
3993 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
8979685e 3994 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 3995 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 3996 }
3997 }
4551ea7c 3998
7ad19338 3999 clusters[it+t0] = indexes[best[bestiter][it]][it];
4551ea7c 4000
4001 // Still needed ????
4002 //if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>4 &&
4003 //cl[best[bestiter][it]][it]->GetLocalTimeBin()<18) clusters[it+t0]
4004 // = indexes[best[bestiter][it]][it]; //Test
4005
7ad19338 4006 }
4551ea7c 4007
7ad19338 4008 //
4551ea7c 4009 // Set tracklet parameters
6c94f330 4010 //
4551ea7c 4011 Double_t trackleterr2 = smoffset[bestiter] + 0.01*0.01;
4012 if (mpads > 3.5) {
4013 trackleterr2 += (mpads - 3.5) * 0.04;
4014 }
4015 trackleterr2 += changes[bestiter] * 0.01;
4016 trackleterr2 *= TMath::Max(14.0 - nfound,1.0);
4017 trackleterr2 += 0.2 * (tany-exB)*(tany-exB);
4018
4019 // Set tracklet parameters
4020 tracklet.Set(xmean
4021 ,track2.GetY() + moffset[bestiter]
4022 ,meanz[bestiter]
4023 ,track2.GetAlpha()
4024 ,trackleterr2);
7ad19338 4025 tracklet.SetTilt(h01);
4026 tracklet.SetP0(mean[bestiter]);
4027 tracklet.SetP1(angle[bestiter]);
4028 tracklet.SetN(nfound);
4029 tracklet.SetNCross(changes[bestiter]);
4030 tracklet.SetPlane(plane);
4031 tracklet.SetSigma2(expectederr);
4032 tracklet.SetChi2(tchi2s[bestiter]);
8979685e 4033 tracklet.SetMaxPos(maxpos,maxpos4,maxpos5);
27eaf44b 4034 track->SetTracklets(plane,tracklet);
4035 track->SetNWrong(track->GetNWrong() + nbad[0]);
4551ea7c 4036
7ad19338 4037 //
4038 // Debuging part
4039 //
69b55c55 4040 TClonesArray array0("AliTRDcluster");
4041 TClonesArray array1("AliTRDcluster");
4551ea7c 4042 array0.ExpandCreateFast(t1 - t0 + 1);
4043 array1.ExpandCreateFast(t1 - t0 + 1);
4044 TTreeSRedirector &cstream = *fDebugStreamer;
7ad19338 4045 AliTRDcluster dummy;
4046 Double_t dy0[100];
8979685e 4047 Double_t dyb[100];
4048
4551ea7c 4049 for (Int_t it = 0; it < t1 - t0; it++) {
7ad19338 4050 dy0[it] = dy[0][it];
4051 dyb[it] = dy[best[bestiter][it]][it];
4551ea7c 4052 if (cl[0][it]) {
7ad19338 4053 new(array0[it]) AliTRDcluster(*cl[0][it]);
4054 }
4551ea7c 4055 else {
7ad19338 4056 new(array0[it]) AliTRDcluster(dummy);
4057 }
4058 if(cl[best[bestiter][it]][it]) {
4059 new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
4060 }
4061 else{
4062 new(array1[it]) AliTRDcluster(dummy);
4063 }
4f1c04d3 4064 }
7ad19338 4065 TGraph graph0(t1-t0,x,dy0);
4066 TGraph graph1(t1-t0,x,dyb);
4067 TGraph graphy(t1-t0,x,yt);
4068 TGraph graphz(t1-t0,x,zt);
4551ea7c 4069
4070 if (AliTRDReconstructor::StreamLevel() > 0) {
4071 cstream << "tracklet"
4072 << "track.=" << track // Track parameters
4073 << "tany=" << tany // Tangent of the local track angle
4074 << "xmean=" << xmean // Xmean - reference x of tracklet
4075 << "tilt=" << h01 // Tilt angle
4076 << "nall=" << nall // Number of foundable clusters
4077 << "nfound=" << nfound // Number of found clusters
4078 << "clfound=" << clfound // Total number of found clusters in road
4079 << "mpads=" << mpads // Mean number of pads per cluster
4080 << "plane=" << plane // Plane number
4081 << "detector=" << detector // Detector number
4082 << "road=" << road // The width of the used road
4083 << "graph0.=" << &graph0 // x - y = dy for closest cluster
4084 << "graph1.=" << &graph1 // x - y = dy for second closest cluster
4085 << "graphy.=" << &graphy // y position of the track
4086 << "graphz.=" << &graphz // z position of the track
4087 //<< "fCl.=" << &array0 // closest cluster
4088 //<< "fCl2.=" << &array1 // second closest cluster
4089 << "maxpos=" << maxpos // Maximal charge postion
4090 << "maxcharge=" << maxcharge // Maximal charge
4091 << "maxpos4=" << maxpos4 // Maximal charge postion - after bin 4
4092 << "maxcharge4=" << maxcharge4 // Maximal charge - after bin 4
4093 << "maxpos5=" << maxpos5 // Maximal charge postion - after bin 5
4094 << "maxcharge5=" << maxcharge5 // Maximal charge - after bin 5
4095 << "bestiter=" << bestiter // Best iteration number
4096 << "tracklet.=" << &tracklet // Corrspond to the best iteration
4097 << "tchi20=" << tchi2s[0] // Chi2 of cluster in the 0 iteration
4098 << "tchi2b=" << tchi2s[bestiter] // Chi2 of cluster in the best iteration
4099 << "sigmas0=" << sigmas[0] // Residuals sigma
4100 << "sigmasb=" << sigmas[bestiter] // Residulas sigma
4101 << "ngood0=" << ngood[0] // Number of good clusters in 0 iteration
4102 << "nbad0=" << nbad[0] // Number of bad clusters in 0 iteration
4103 << "ngoodb=" << ngood[bestiter] // in best iteration
4104 << "nbadb=" << nbad[bestiter] // in best iteration
4105 << "changes0=" << changes[0] // Changes of pardrows in iteration number 0
4106 << "changesb=" << changes[bestiter] // Changes of pardrows in best iteration
4107 << "moffset0=" << moffset[0] // Offset fixing angle in iter=0
4108 << "smoffset0=" << smoffset[0] // Sigma of offset fixing angle in iter=0
4109 << "moffsetb=" << moffset[bestiter] // Offset fixing angle in iter=best
4110 << "smoffsetb=" << smoffset[bestiter] // Sigma of offset fixing angle in iter=best
4111 << "mean0=" << mean[0] // Mean dy in iter=0;
4112 << "smean0=" << smean[0] // Sigma of mean dy in iter=0
4113 << "meanb=" << mean[bestiter] // Mean dy in iter=best
4114 << "smeanb=" << smean[bestiter] // Sigma of mean dy in iter=best
4115 << "angle0=" << angle[0] // Angle deviation in the iteration number 0
4116 << "sangle0=" << sangle[0] // Sigma of angular deviation in iteration number 0
4117 << "angleb=" << angle[bestiter] // Angle deviation in the best iteration
4118 << "sangleb=" << sangle[bestiter] // Sigma of angle deviation in the best iteration
4119 << "expectederr=" << expectederr // Expected error of cluster position
4120 << "\n";
4121 }
4122
4f1c04d3 4123 return nfound;
4f1c04d3 4124
75bd7f81 4125}
4f1c04d3 4126
75bd7f81 4127//_____________________________________________________________________________
4128Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist
4129 , Int_t *outlist, Bool_t down)
69b55c55 4130{
4131 //
4551ea7c 4132 // Sort eleements according occurancy
4133 // The size of output array has is 2*n
69b55c55 4134 //
75bd7f81 4135
4551ea7c 4136 Int_t *sindexS = new Int_t[n]; // Temporary array for sorting
4137 Int_t *sindexF = new Int_t[2*n];
4138 for (Int_t i = 0; i < n; i++) {
4139 sindexF[i] = 0;
4140 }
4141
4142 TMath::Sort(n,inlist,sindexS,down);
4143
4144 Int_t last = inlist[sindexS[0]];
4145 Int_t val = last;
4146 sindexF[0] = 1;
4147 sindexF[0+n] = last;
4148 Int_t countPos = 0;
4149
4150 // Find frequency
4151 for (Int_t i = 1; i < n; i++) {
69b55c55 4152 val = inlist[sindexS[i]];
4551ea7c 4153 if (last == val) {
4154 sindexF[countPos]++;
4155 }
4156 else {
69b55c55 4157 countPos++;
4158 sindexF[countPos+n] = val;
4159 sindexF[countPos]++;
4551ea7c 4160 last = val;
69b55c55 4161 }
4162 }
4551ea7c 4163 if (last == val) {
4164 countPos++;
4165 }
4166
4167 // Sort according frequency
4168 TMath::Sort(countPos,sindexF,sindexS,kTRUE);
4169
4170 for (Int_t i = 0; i < countPos; i++) {
69b55c55 4171 outlist[2*i ] = sindexF[sindexS[i]+n];
4172 outlist[2*i+1] = sindexF[sindexS[i]];
4173 }
4551ea7c 4174
69b55c55 4175 delete [] sindexS;
4176 delete [] sindexF;
4177
4178 return countPos;
75bd7f81 4179
69b55c55 4180}
4181
75bd7f81 4182//_____________________________________________________________________________
4551ea7c 4183AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params)
69b55c55 4184{
4185 //
75bd7f81 4186 // Register a seed
69b55c55 4187 //
4551ea7c 4188
4189 Double_t alpha = AliTRDgeometry::GetAlpha();
4190 Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
69b55c55 4191 Double_t c[15];
4551ea7c 4192
4193 c[ 0] = 0.2;
4194 c[ 1] = 0.0; c[ 2] = 2.0;
4195 c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
4196 c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0; c[ 9] = 0.1;
4197 c[10] = 0.0; c[11] = 0.0; c[12] = 0.0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
4198
4199 Int_t index = 0;
4200 AliTRDcluster *cl = 0;
4201
4202 for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
4203 if (seeds[ilayer].IsOK()) {
4204 for (Int_t itime = 22; itime > 0; itime--) {
ca1e1e5b 4205 if (seeds[ilayer].GetIndexes(itime) > 0) {
4206 index = seeds[ilayer].GetIndexes(itime);
4207 cl = seeds[ilayer].GetClusters(itime);
69b55c55 4208 break;
4209 }
4210 }
4211 }
4551ea7c 4212 if (index > 0) {
4213 break;
4214 }
69b55c55 4215 }
4551ea7c 4216 if (cl == 0) {
4217 return 0;
4218 }
4219
4220 AliTRDtrack *track = new AliTRDtrack(cl
4221 ,index
4222 ,&params[1]
4223 ,c
4224 ,params[0]
4225 ,params[6]*alpha+shift);
43bfc8af 4226 // SetCluster(cl, 0); // A. Bercuci
4227 track->PropagateTo(params[0]-5.0);
69b55c55 4228 track->ResetCovariance(1);
4551ea7c 4229
4230 Int_t rc = FollowBackProlongation(*track);
4231 if (rc < 30) {
69b55c55 4232 delete track;
4551ea7c 4233 track = 0;
4234 }
4235 else {
69b55c55 4236 track->CookdEdx();
43bfc8af 4237 track->CookdEdxTimBin();
4551ea7c 4238 CookLabel(track,0.9);
69b55c55 4239 }
69b55c55 4240
75bd7f81 4241 return track;
7e448bcc 4242}
4243
4244//////////////////////////////////////////////////////////////////////////////////////////
69b55c55 4245
7e448bcc 4246void AliTRDtracker::InitLogHists() {
4247
4248 fHBackfit = new TH1D("logTRD_backfit", "", 40, -0.5, 39.5);
4249 fHRefit = new TH1D("logTRD_refit", "", 40, -0.5, 39.5);
4250 fHClSearch = new TH1D("logTRD_clSearch", "", 60, -0.5, 59.5);
4251
4252 fHX = new TH1D("logTRD_X", ";x (cm)", 200, 50, 400);
4253 fHNCl = new TH1D("logTRD_ncl", "", 40, -0.5, 39.5);
4254 fHNClTrack = new TH1D("logTRD_nclTrack", "", 180, -0.5, 179.5);
4255
4256 fHMinYPos = new TH1D("logTRD_minYPos", ";#delta Y (cm)", 400, -6, 6);
4257 fHMinYNeg = new TH1D("logTRD_minYNeg", ";#delta Y (cm)", 400, -6, 6);
4258 fHMinZ = new TH1D("logTRD_minZ", ";#delta Z (cm)", 400, -20, 20);
4259 fHMinD = new TH2D("logTRD_minD", ";#delta Y (cm);#delta Z (cm)", 100, -6, 6, 100, -50, 50);
4260
4261 fHDeltaX = new TH1D("logTRD_deltaX", ";#delta X (cm)", 100, -5, 5);
4262 fHXCl = new TH1D("logTRD_xCl", ";cluster x position (cm)", 1000, 280, 380);
4263
4264 const char *nameFindCl[4] = {"logTRD_clY", "logTRD_clZ", "logTRD_clB", "logTRD_clG"};
4265
4266 for(int i=0; i<4; i++) {
4267 fHFindCl[i] = new TH1D(nameFindCl[i], "", 30, -0.5, 29.5);
4268 }
69b55c55 4269}
7e448bcc 4270
4271//////////////////////////////////////////////////////////////////////////////////////////
4272
4273void AliTRDtracker::SaveLogHists() {
4274
4275 TDirectory *sav = gDirectory;
4276 TFile *logFile = 0;
4277
4278 TSeqCollection *col = gROOT->GetListOfFiles();
4279 int N = col->GetEntries();
4280 for(int i=0; i<N; i++) {
4281 logFile = (TFile*)col->At(i);
4282 if (strstr(logFile->GetName(), "AliESDs.root")) break;
4283 }
4284
4285 logFile->cd();
4286 fHBackfit->Write(fHBackfit->GetName(), TObject::kOverwrite);
4287 fHRefit->Write(fHRefit->GetName(), TObject::kOverwrite);
4288 fHClSearch->Write(fHClSearch->GetName(), TObject::kOverwrite);
4289 fHX->Write(fHX->GetName(), TObject::kOverwrite);
4290 fHNCl->Write(fHNCl->GetName(), TObject::kOverwrite);
4291 fHNClTrack->Write(fHNClTrack->GetName(), TObject::kOverwrite);
4292
4293 fHMinYPos->Write(fHMinYPos->GetName(), TObject::kOverwrite);
4294 fHMinYNeg->Write(fHMinYNeg->GetName(), TObject::kOverwrite);
4295 fHMinD->Write(fHMinD->GetName(), TObject::kOverwrite);
4296 fHMinZ->Write(fHMinZ->GetName(), TObject::kOverwrite);
4297
4298 fHDeltaX->Write(fHDeltaX->GetName(), TObject::kOverwrite);
4299 fHXCl->Write(fHXCl->GetName(), TObject::kOverwrite);
4300
4301
4302 for(int i=0; i<4; i++)
4303 fHFindCl[i]->Write(fHFindCl[i]->GetName(), TObject::kOverwrite);
4304
4305 logFile->Flush();
4306
4307 sav->cd();
4308}
4309
4310//////////////////////////////////////////////////////////////////////////////////////////