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