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