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