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