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