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