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