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