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