]> 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
Move pad planes from AliTRDCommomParam to AliTRDgeometry
[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"
082a2283 57#include "AliTRDCalibraFillHisto.h"
58
7e448bcc 59ClassImp(AliTRDtracker)
46d29e70 60
b08a2e7e 61const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5; //
62const Float_t AliTRDtracker::fgkLabelFraction = 0.8; //
63const Double_t AliTRDtracker::fgkMaxChi2 = 12.0; //
7e448bcc 64const Double_t AliTRDtracker::fgkMaxSnp = 0.95; // Corresponds to tan = 3
65const Double_t AliTRDtracker::fgkMaxStep = 2.0; // Maximal step size in propagation
7ad19338 66
75bd7f81 67//_____________________________________________________________________________
4551ea7c 68AliTRDtracker::AliTRDtracker()
4e009ce4 69 :AliTracker()
70 ,fHBackfit(0x0)
71 ,fHClSearch(0x0)
72 ,fHRefit(0x0)
73 ,fHX(0x0)
74 ,fHNCl(0x0)
75 ,fHNClTrack(0x0)
76 ,fHMinYPos(0x0)
77 ,fHMinYNeg(0x0)
78 ,fHMinZ(0x0)
79 ,fHMinD(0x0)
80 ,fHDeltaX(0x0)
81 ,fHXCl(0x0)
82 ,fGeom(0)
83 ,fNclusters(0)
84 ,fClusters(0)
85 ,fNseeds(0)
86 ,fSeeds(0)
87 ,fNtracks(0)
88 ,fTracks(0)
89 ,fTimeBinsPerPlane(0)
90 ,fAddTRDseeds(kFALSE)
91 ,fNoTilt(kFALSE)
92 ,fDebugStreamer(0)
89f05372 93{
75bd7f81 94 //
b7a0917f 95 // Default constructor
75bd7f81 96 //
b7a0917f 97
4551ea7c 98 for (Int_t i = 0; i < kTrackingSectors; i++) {
99 fTrSec[i] = 0;
100 }
101 for (Int_t j = 0; j < 5; j++) {
102 for (Int_t k = 0; k < 18; k++) {
103 fHoles[j][k] = kFALSE;
104 }
105 }
75bd7f81 106
7e448bcc 107 InitLogHists();
4e009ce4 108
89f05372 109}
75bd7f81 110
111//_____________________________________________________________________________
6c94f330 112AliTRDtracker::AliTRDtracker(const AliTRDtracker &t)
113 :AliTracker(t)
4e009ce4 114 ,fHBackfit(0x0)
115 ,fHClSearch(0x0)
116 ,fHRefit(0x0)
117 ,fHX(0x0)
118 ,fHNCl(0x0)
119 ,fHNClTrack(0x0)
120 ,fHMinYPos(0x0)
121 ,fHMinYNeg(0x0)
122 ,fHMinZ(0x0)
123 ,fHMinD(0x0)
124 ,fHDeltaX(0x0)
125 ,fHXCl(0x0)
ad670fba 126 ,fGeom(0)
127 ,fNclusters(0)
128 ,fClusters(0)
129 ,fNseeds(0)
130 ,fSeeds(0)
131 ,fNtracks(0)
132 ,fTracks(0)
133 ,fTimeBinsPerPlane(0)
134 ,fAddTRDseeds(kFALSE)
135 ,fNoTilt(kFALSE)
136 ,fDebugStreamer(0)
6c94f330 137{
ad670fba 138 //
139 // Copy constructor
140 //
4e009ce4 141
ad670fba 142}
143
144//_____________________________________________________________________________
145AliTRDtracker::AliTRDtracker(const TFile *geomfile)
146 :AliTracker()
4e009ce4 147 ,fHBackfit(0x0)
148 ,fHClSearch(0x0)
149 ,fHRefit(0x0)
150 ,fHX(0x0)
151 ,fHNCl(0x0)
152 ,fHNClTrack(0x0)
153 ,fHMinYPos(0x0)
154 ,fHMinYNeg(0x0)
155 ,fHMinZ(0x0)
156 ,fHMinD(0x0)
157 ,fHDeltaX(0x0)
158 ,fHXCl(0x0)
ad670fba 159 ,fGeom(0)
160 ,fNclusters(0)
161 ,fClusters(new TObjArray(2000))
162 ,fNseeds(0)
163 ,fSeeds(new TObjArray(2000))
164 ,fNtracks(0)
165 ,fTracks(new TObjArray(1000))
166 ,fTimeBinsPerPlane(0)
167 ,fAddTRDseeds(kFALSE)
168 ,fNoTilt(kFALSE)
169 ,fDebugStreamer(0)
46d29e70 170{
5443e65e 171 //
172 // Main constructor
173 //
b8dc2353 174
4551ea7c 175 TDirectory *savedir = gDirectory;
46d29e70 176
f162af62 177 //TFile *in = (TFile *) geomfile;
178 //
179 //if (!in->IsOpen()) {
180 // AliWarning("geometry file is not open!\n");
181 // AliWarning("FULL TRD geometry and DEFAULT TRD parameter will be used\n");
182 //}
183 //else {
184 // in->cd();
185 // fGeom = (AliTRDgeometry *) in->Get("TRDgeometry");
186 //}
187 //
188 //if (!fGeom) {
189 // AliWarning("Cannot find TRD geometry!\n");
190 // fGeom = new AliTRDgeometry();
191 //}
192 fGeom = new AliTRDgeometry();
c6f438c0 193 fGeom->ReadGeoMatrices();
c630aafd 194
f162af62 195 //savedir->cd();
46d29e70 196
4551ea7c 197 for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
439c63c8 198 Int_t trS = geomS;
4551ea7c 199 fTrSec[trS] = new AliTRDtrackingSector(fGeom,geomS);
200 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
201 fHoles[icham][trS] = fGeom->IsHole(0,icham,geomS);
3c625a9b 202 }
5443e65e 203 }
4551ea7c 204
f162af62 205 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0);
7ad19338 206 Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
4551ea7c 207 if (tiltAngle < 0.1) {
b8dc2353 208 fNoTilt = kTRUE;
209 }
210
fc88348b 211 if (!AliTRDcalibDB::Instance()) {
212 AliFatal("Could not get calibration object");
213 }
214 fTimeBinsPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
46d29e70 215
4551ea7c 216 fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
0a29d0f1 217
9c9d2487 218 savedir->cd();
7e448bcc 219
220 InitLogHists();
75bd7f81 221
5443e65e 222}
46d29e70 223
75bd7f81 224//_____________________________________________________________________________
5443e65e 225AliTRDtracker::~AliTRDtracker()
46d29e70 226{
029cd327 227 //
228 // Destructor of AliTRDtracker
229 //
230
89f05372 231 if (fClusters) {
232 fClusters->Delete();
233 delete fClusters;
234 }
4551ea7c 235
89f05372 236 if (fTracks) {
237 fTracks->Delete();
238 delete fTracks;
239 }
4551ea7c 240
89f05372 241 if (fSeeds) {
242 fSeeds->Delete();
243 delete fSeeds;
244 }
4551ea7c 245
f162af62 246 if (fGeom) {
247 delete fGeom;
248 }
0a29d0f1 249
4551ea7c 250 for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
029cd327 251 delete fTrSec[geomS];
5443e65e 252 }
4551ea7c 253
7ad19338 254 if (fDebugStreamer) {
7ad19338 255 delete fDebugStreamer;
256 }
9c9d2487 257
75bd7f81 258}
59393e34 259
75bd7f81 260//_____________________________________________________________________________
261Int_t AliTRDtracker::LocalToGlobalID(Int_t lid)
262{
59393e34 263 //
75bd7f81 264 // Transform internal TRD ID to global detector ID
59393e34 265 //
75bd7f81 266
4551ea7c 267 Int_t isector = fGeom->GetSector(lid);
268 Int_t ichamber = fGeom->GetChamber(lid);
269 Int_t iplan = fGeom->GetPlane(lid);
270
ae079791 271 AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
59393e34 272 switch (iplan) {
273 case 0:
ae079791 274 iLayer = AliGeomManager::kTRD1;
59393e34 275 break;
276 case 1:
ae079791 277 iLayer = AliGeomManager::kTRD2;
59393e34 278 break;
279 case 2:
ae079791 280 iLayer = AliGeomManager::kTRD3;
59393e34 281 break;
282 case 3:
ae079791 283 iLayer = AliGeomManager::kTRD4;
59393e34 284 break;
285 case 4:
ae079791 286 iLayer = AliGeomManager::kTRD5;
59393e34 287 break;
288 case 5:
ae079791 289 iLayer = AliGeomManager::kTRD6;
59393e34 290 break;
291 };
4551ea7c 292
293 Int_t modId = isector * fGeom->Ncham() + ichamber;
ae079791 294 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,modId);
75bd7f81 295
59393e34 296 return volid;
75bd7f81 297
59393e34 298}
299
75bd7f81 300//_____________________________________________________________________________
301Int_t AliTRDtracker::GlobalToLocalID(Int_t gid)
302{
59393e34 303 //
75bd7f81 304 // Transform global detector ID to local detector ID
59393e34 305 //
75bd7f81 306
4551ea7c 307 Int_t modId = 0;
ae079791 308 AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(gid,modId);
4551ea7c 309
310 Int_t isector = modId / fGeom->Ncham();
311 Int_t ichamber = modId % fGeom->Ncham();
312 Int_t iLayer = -1;
313
59393e34 314 switch (layerId) {
ae079791 315 case AliGeomManager::kTRD1:
6c94f330 316 iLayer = 0;
59393e34 317 break;
ae079791 318 case AliGeomManager::kTRD2:
6c94f330 319 iLayer = 1;
59393e34 320 break;
ae079791 321 case AliGeomManager::kTRD3:
6c94f330 322 iLayer = 2;
59393e34 323 break;
ae079791 324 case AliGeomManager::kTRD4:
6c94f330 325 iLayer = 3;
59393e34 326 break;
ae079791 327 case AliGeomManager::kTRD5:
6c94f330 328 iLayer = 4;
59393e34 329 break;
ae079791 330 case AliGeomManager::kTRD6:
6c94f330 331 iLayer = 5;
59393e34 332 break;
333 default:
6c94f330 334 iLayer =-1;
59393e34 335 }
4551ea7c 336
337 if (iLayer < 0) {
338 return -1;
339 }
340
59393e34 341 Int_t lid = fGeom->GetDetector(iLayer,ichamber,isector);
75bd7f81 342
59393e34 343 return lid;
59393e34 344
75bd7f81 345}
59393e34 346
75bd7f81 347//_____________________________________________________________________________
4551ea7c 348Bool_t AliTRDtracker::Transform(AliTRDcluster *cluster)
75bd7f81 349{
59393e34 350 //
eb37e0af 351 // Transform from cluster system to tracking system
c6f438c0 352 //
75bd7f81 353
33744e5d 354 // Magic constants for geo manager transformation
4551ea7c 355 const Double_t kX0shift = 2.52;
33744e5d 356
6c94f330 357 //
33744e5d 358 // Apply alignment and calibration to transform cluster
6c94f330 359 //
360 Int_t detector = cluster->GetDetector();
4551ea7c 361 Int_t plane = fGeom->GetPlane(cluster->GetDetector());
362 Int_t chamber = fGeom->GetChamber(cluster->GetDetector());
363 Int_t sector = fGeom->GetSector(cluster->GetDetector());
364
365 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
366 Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.0); // Drift distance
c6f438c0 367
59393e34 368 //
59393e34 369 // ExB correction
370 //
371 Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
4e009ce4 372 Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift,-AliTracker::GetBz()*0.1);
4551ea7c 373
f162af62 374 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(plane,chamber);
4551ea7c 375 Double_t zshiftIdeal = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
376 Double_t localPos[3];
377 Double_t localPosTracker[3];
c6f438c0 378 localPos[0] = -cluster->GetX();
4551ea7c 379 localPos[1] = cluster->GetY() - driftX * exB;
380 localPos[2] = cluster->GetZ() - zshiftIdeal;
381
c6f438c0 382 cluster->SetY(cluster->GetY() - driftX*exB);
383 Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
2fa01832 384 cluster->SetX(xplane - cluster->GetX());
4551ea7c 385
2fa01832 386 TGeoHMatrix *matrix = fGeom->GetCorrectionMatrix(cluster->GetDetector());
4551ea7c 387 if (!matrix) {
388 // No matrix found - if somebody used geometry with holes
c6f438c0 389 AliError("Invalid Geometry - Default Geometry used\n");
390 return kTRUE;
391 }
4551ea7c 392 matrix->LocalToMaster(localPos,localPosTracker);
393
394 if (AliTRDReconstructor::StreamLevel() > 1) {
395 (* fDebugStreamer) << "Transform"
396 << "Cl.=" << cluster
397 << "matrix.=" << matrix
398 << "Detector=" << detector
399 << "Sector=" << sector
400 << "Plane=" << plane
401 << "Chamber=" << chamber
402 << "lx0=" << localPosTracker[0]
403 << "ly0=" << localPosTracker[1]
404 << "lz0=" << localPosTracker[2]
405 << "\n";
c6f438c0 406 }
4551ea7c 407
eb37e0af 408 cluster->SetX(localPosTracker[0]+kX0shift);
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;
31fd97b2 794 Int_t eventNrInFile = event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
4551ea7c 795 if (AliTRDReconstructor::StreamLevel() > 0) {
796 if (track->GetBackupTrack()) {
797 cstream << "Tracks"
31fd97b2 798 << "EventNrInFile=" << eventNrInFile
4551ea7c 799 << "ESD.=" << seed
800 << "trd.=" << track
801 << "trdback.=" << track->GetBackupTrack()
802 << "\n";
803 }
804 else {
805 cstream << "Tracks"
31fd97b2 806 << "EventNrInFile=" << eventNrInFile
4551ea7c 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
082a2283 1217 AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
77566f2a 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
f162af62 1577 AliTRDpadPlane *padPlane = fGeom->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
31fd97b2 2627 Int_t eventNrInFile = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
4551ea7c 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"
31fd97b2 2646 << "EventNrInFile=" << eventNrInFile
4551ea7c 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
31fd97b2 2676 << "EventNrInFile=" << eventNrInFile
4551ea7c 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]);
ae079791 2757 AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
3551db50 2758 switch (iplan) {
2759 case 0:
ae079791 2760 iLayer = AliGeomManager::kTRD1;
3551db50 2761 break;
2762 case 1:
ae079791 2763 iLayer = AliGeomManager::kTRD2;
3551db50 2764 break;
2765 case 2:
ae079791 2766 iLayer = AliGeomManager::kTRD3;
3551db50 2767 break;
2768 case 3:
ae079791 2769 iLayer = AliGeomManager::kTRD4;
3551db50 2770 break;
2771 case 4:
ae079791 2772 iLayer = AliGeomManager::kTRD5;
3551db50 2773 break;
2774 case 5:
ae079791 2775 iLayer = AliGeomManager::kTRD6;
3551db50 2776 break;
2777 };
4551ea7c 2778 Int_t modId = isector * fGeom->Ncham() + ichamber;
ae079791 2779 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,modId);
3551db50 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];
3551db50 3041
ad670fba 3042 for (Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
3043
3044 ymax = fGeom->GetChamberWidth(plane) / 2.0;
f162af62 3045 padPlane = fGeom->GetPadPlane(plane,0);
ad670fba 3046 ymaxsensitive = (padPlane->GetColSize(1) * padPlane->GetNcols() - 4.0) / 2.0;
3047
3048 for (Int_t ch = 0; ch < kNchambers; ch++) {
3049 zmax[ch] = fGeom->GetChamberLength(plane,ch) / 2.0;
3050 Float_t pad = padPlane->GetRowSize(1);
f162af62 3051 Float_t row0 = fGeom->GetRow0(plane,ch,0);
3052 Int_t nPads = fGeom->GetRowMax(plane,ch,0);
ad670fba 3053 zmaxsensitive[ch] = Float_t(nPads) * pad / 2.0;
3054 zc[ch] = -(pad * nPads) / 2.0 + row0;
5443e65e 3055 }
3056
ad670fba 3057 dx = AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
f162af62 3058 / AliTRDCommonParam::Instance()->GetSamplingFrequency();
ad670fba 3059 rho = 0.00295 * 0.85; //????
3060 radLength = 11.0;
5443e65e 3061
3551db50 3062 Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(plane);
a305677e 3063 //Double_t xbottom = x0 - dxDrift;
ad670fba 3064 //Double_t xtop = x0 + dxAmp;
3065
59393e34 3066 Int_t nTimeBins = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
ad670fba 3067 for (Int_t iTime = 0; iTime < nTimeBins; iTime++) {
3068
3069 Double_t xlayer = iTime * dx - dxAmp;
3070 //if (xlayer<0) xlayer = dxAmp / 2.0;
59393e34 3071 x = x0 - xlayer;
ad670fba 3072
3073 tbIndex = CookTimeBinIndex(plane,iTime);
3074 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex,plane);
3c625a9b 3075 ppl->SetYmax(ymax,ymaxsensitive);
ad670fba 3076 ppl->SetZ(zc,zmax,zmaxsensitive);
3c625a9b 3077 ppl->SetHoles(holes);
59393e34 3078 InsertLayer(ppl);
ad670fba 3079
5443e65e 3080 }
ad670fba 3081
5443e65e 3082 }
3083
5443e65e 3084 MapTimeBinLayers();
ad670fba 3085
029cd327 3086 delete [] zc;
3087 delete [] zmax;
4f1c04d3 3088 delete [] zmaxsensitive;
5443e65e 3089
3090}
3091
ad670fba 3092//_____________________________________________________________________________
3093AliTRDtracker::AliTRDtrackingSector
3094 ::AliTRDtrackingSector(const AliTRDtrackingSector &/*t*/)
3095 :fN(0)
3096 ,fGeom(0)
3097 ,fGeomSector(0)
3098{
3099 //
3100 // Copy constructor
3101 //
3102
3103}
3104
75bd7f81 3105//_____________________________________________________________________________
3106Int_t AliTRDtracker::AliTRDtrackingSector
3107 ::CookTimeBinIndex(Int_t plane, Int_t localTB) const
5443e65e 3108{
3109 //
6c94f330 3110 // depending on the digitization parameters calculates "global"
029cd327 3111 // time bin index for timebin <localTB> in plane <plane>
5443e65e 3112 //
59393e34 3113 //
75bd7f81 3114
59393e34 3115 Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
4551ea7c 3116 Int_t gtb = (plane+1) * tbPerPlane - localTB - 1;
3117 if (localTB < 0) {
3118 return -1;
3119 }
3120 if (gtb < 0) {
3121 return -1;
3122 }
75bd7f81 3123
5443e65e 3124 return gtb;
5443e65e 3125
75bd7f81 3126}
5443e65e 3127
75bd7f81 3128//_____________________________________________________________________________
3129void AliTRDtracker::AliTRDtrackingSector
3130 ::MapTimeBinLayers()
5443e65e 3131{
3132 //
3133 // For all sensitive time bins sets corresponding layer index
3134 // in the array fTimeBins
3135 //
3136
3137 Int_t index;
3138
4551ea7c 3139 for (Int_t i = 0; i < fN; i++) {
3140
5443e65e 3141 index = fLayers[i]->GetTimeBinIndex();
6c94f330 3142
4551ea7c 3143 if (index < 0) {
3144 continue;
3145 }
3146 if (index >= (Int_t) kMaxTimeBinIndex) {
3147 //AliWarning(Form("Index %d exceeds allowed maximum of %d!\n"
3148 // ,index,kMaxTimeBinIndex-1));
5443e65e 3149 continue;
3150 }
4551ea7c 3151
5443e65e 3152 fTimeBinIndex[index] = i;
4551ea7c 3153
5443e65e 3154 }
5443e65e 3155
75bd7f81 3156}
5443e65e 3157
75bd7f81 3158//_____________________________________________________________________________
3159Int_t AliTRDtracker::AliTRDtrackingSector
3160 ::GetLayerNumber(Double_t x) const
5443e65e 3161{
3162 //
3163 // Returns the number of time bin which in radial position is closest to <x>
3164 //
3165
4551ea7c 3166 if (x >= fLayers[fN-1]->GetX()) {
3167 return fN - 1;
3168 }
3169 if (x <= fLayers[ 0]->GetX()) {
3170 return 0;
3171 }
3172
3173 Int_t b = 0;
3174 Int_t e = fN - 1;
3175 Int_t m = (b + e) / 2;
ad670fba 3176
4551ea7c 3177 for ( ; b < e; m = (b + e) / 2) {
3178 if (x > fLayers[m]->GetX()) {
3179 b = m + 1;
3180 }
3181 else {
3182 e = m;
3183 }
5443e65e 3184 }
75bd7f81 3185
4551ea7c 3186 if (TMath::Abs(x - fLayers[m]->GetX()) > TMath::Abs(x - fLayers[m+1]->GetX())) {
3187 return m + 1;
3188 }
3189 else {
3190 return m;
3191 }
5443e65e 3192
3193}
3194
75bd7f81 3195//_____________________________________________________________________________
3196Int_t AliTRDtracker::AliTRDtrackingSector
3197 ::GetInnerTimeBin() const
5443e65e 3198{
3199 //
3200 // Returns number of the innermost SENSITIVE propagation layer
3201 //
3202
3203 return GetLayerNumber(0);
5443e65e 3204
75bd7f81 3205}
5443e65e 3206
75bd7f81 3207//_____________________________________________________________________________
3208Int_t AliTRDtracker::AliTRDtrackingSector
3209 ::GetOuterTimeBin() const
5443e65e 3210{
3211 //
3212 // Returns number of the outermost SENSITIVE time bin
3213 //
3214
3215 return GetLayerNumber(GetNumberOfTimeBins() - 1);
46d29e70 3216
75bd7f81 3217}
5443e65e 3218
75bd7f81 3219//_____________________________________________________________________________
3220Int_t AliTRDtracker::AliTRDtrackingSector
3221 ::GetNumberOfTimeBins() const
5443e65e 3222{
3223 //
3224 // Returns number of SENSITIVE time bins
3225 //
3226
4551ea7c 3227 Int_t tb;
3228 Int_t layer;
3229
3230 for (tb = kMaxTimeBinIndex - 1; tb >= 0; tb--) {
5443e65e 3231 layer = GetLayerNumber(tb);
4551ea7c 3232 if (layer >= 0) {
3233 break;
3234 }
5443e65e 3235 }
75bd7f81 3236
4551ea7c 3237 return tb + 1;
5443e65e 3238
75bd7f81 3239}
5443e65e 3240
75bd7f81 3241//_____________________________________________________________________________
3242void AliTRDtracker::AliTRDtrackingSector
4551ea7c 3243 ::InsertLayer(AliTRDpropagationLayer *pl)
5443e65e 3244{
3245 //
3246 // Insert layer <pl> in fLayers array.
3247 // Layers are sorted according to X coordinate.
75bd7f81 3248 //
5443e65e 3249
4551ea7c 3250 if (fN == ((Int_t) kMaxLayersPerSector)) {
3251 //AliWarning("Too many layers !\n");
3252 return;
3253 }
3254
3255 if (fN == 0) {
3256 fLayers[fN++] = pl;
ad670fba 3257 return;
3258 }
3259
4551ea7c 3260 Int_t i = Find(pl->GetX());
3261
3262 memmove(fLayers+i+1,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
3263
3264 fLayers[i] = pl;
3265 fN++;
5443e65e 3266
3267}
3268
75bd7f81 3269//_____________________________________________________________________________
3270Int_t AliTRDtracker::AliTRDtrackingSector
3271 ::Find(Double_t x) const
5443e65e 3272{
3273 //
3274 // Returns index of the propagation layer nearest to X
3275 //
3276
4551ea7c 3277 if (x <= fLayers[0]->GetX()) {
3278 return 0;
3279 }
3280
3281 if (x > fLayers[fN-1]->GetX()) {
3282 return fN;
3283 }
3284
3285 Int_t b = 0;
3286 Int_t e = fN-1;
3287 Int_t m = (b + e) / 2;
3288
3289 for (; b < e; m = (b + e) / 2) {
3290 if (x > fLayers[m]->GetX()) {
3291 b = m + 1;
3292 }
3293 else {
3294 e = m;
3295 }
5443e65e 3296 }
7ad19338 3297
75bd7f81 3298 return m;
7ad19338 3299
75bd7f81 3300}
7ad19338 3301
75bd7f81 3302//_____________________________________________________________________________
3303void AliTRDtracker::AliTRDpropagationLayer
4551ea7c 3304 ::SetZ(Double_t *center, Double_t *w, Double_t *wsensitive )
3c625a9b 3305{
3306 //
6c94f330 3307 // set centers and the width of sectors
75bd7f81 3308 //
3309
4551ea7c 3310 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3311 fZc[icham] = center[icham];
3312 fZmax[icham] = w[icham];
3c625a9b 3313 fZmaxSensitive[icham] = wsensitive[icham];
3c625a9b 3314 }
75bd7f81 3315
3c625a9b 3316}
5443e65e 3317
75bd7f81 3318//_____________________________________________________________________________
3c625a9b 3319void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
3320{
3321 //
6c94f330 3322 // set centers and the width of sectors
75bd7f81 3323 //
3324
3c625a9b 3325 fHole = kFALSE;
4551ea7c 3326
3327 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3c625a9b 3328 fIsHole[icham] = holes[icham];
4551ea7c 3329 if (holes[icham]) {
3330 fHole = kTRUE;
3331 }
3c625a9b 3332 }
5443e65e 3333
75bd7f81 3334}
5443e65e 3335
75bd7f81 3336//_____________________________________________________________________________
3337void AliTRDtracker::AliTRDpropagationLayer
4551ea7c 3338 ::InsertCluster(AliTRDcluster *c, UInt_t index)
75bd7f81 3339{
3340 //
3341 // Insert cluster in cluster array.
3342 // Clusters are sorted according to Y coordinate.
3343 //
5443e65e 3344
4551ea7c 3345 if (fTimeBinIndex < 0) {
3346 //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");
3347 return;
3348 }
3349
3350 if (fN == (Int_t) kMaxClusterPerTimeBin) {
3351 //AliWarning("Too many clusters !\n");
5443e65e 3352 return;
3353 }
ad670fba 3354
4551ea7c 3355 if (fN == 0) {
3356 fIndex[0] = index;
3357 fClusters[fN++] = c;
ad670fba 3358 return;
3359 }
4551ea7c 3360
3361 Int_t i = Find(c->GetY());
3362 memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
3363 memmove(fIndex +i+1,fIndex +i,(fN-i)*sizeof(UInt_t));
3364 fIndex[i] = index;
3365 fClusters[i] = c;
3366 fN++;
5443e65e 3367
75bd7f81 3368}
5443e65e 3369
75bd7f81 3370//_____________________________________________________________________________
3371Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const
3372{
3373 //
3374 // Returns index of the cluster nearest in Y
3375 //
5443e65e 3376
4551ea7c 3377 if (fN <= 0) {
3378 return 0;
3379 }
3380 if (y <= fClusters[0]->GetY()) {
3381 return 0;
3382 }
3383 if (y > fClusters[fN-1]->GetY()) {
3384 return fN;
3385 }
3386
3387 Int_t b = 0;
3388 Int_t e = fN - 1;
3389 Int_t m = (b + e) / 2;
3390
3391 for ( ; b < e; m = (b + e) / 2) {
3392 if (y > fClusters[m]->GetY()) {
3393 b = m + 1;
3394 }
3395 else {
3396 e = m;
3397 }
5443e65e 3398 }
75bd7f81 3399
5443e65e 3400 return m;
75bd7f81 3401
5443e65e 3402}
3403
75bd7f81 3404//_____________________________________________________________________________
3405Int_t AliTRDtracker::AliTRDpropagationLayer
3406 ::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad
3407 , Float_t maxroadz) const
7ad19338 3408{
3409 //
3410 // Returns index of the cluster nearest to the given y,z
3411 //
75bd7f81 3412
4551ea7c 3413 Int_t index = -1;
3414 Int_t maxn = fN;
69b55c55 3415 Float_t mindist = maxroad;
4551ea7c 3416
3417 for (Int_t i = Find(y-maxroad); i < maxn; i++) {
3418 AliTRDcluster *c = (AliTRDcluster *) (fClusters[i]);
69b55c55 3419 Float_t ycl = c->GetY();
4551ea7c 3420 if (ycl > (y + maxroad)) {
3421 break;
3422 }
3423 if (TMath::Abs(c->GetZ() - z) > maxroadz) {
3424 continue;
3425 }
3426 if (TMath::Abs(ycl - y) < mindist) {
3427 mindist = TMath::Abs(ycl - y);
3428 index = fIndex[i];
3429 }
6c94f330 3430 }
75bd7f81 3431
7ad19338 3432 return index;
7ad19338 3433
75bd7f81 3434}
7ad19338 3435
75bd7f81 3436//_____________________________________________________________________________
4551ea7c 3437Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c)
75bd7f81 3438{
3439 //
3440 // Returns correction factor for tilted pads geometry
3441 //
5443e65e 3442
4551ea7c 3443 Int_t det = c->GetDetector();
3444 Int_t plane = fGeom->GetPlane(det);
f162af62 3445 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(plane,0);
4551ea7c 3446 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
b8dc2353 3447
4551ea7c 3448 if (fNoTilt) {
3449 h01 = 0;
3450 }
75bd7f81 3451
fd621f36 3452 return h01;
5443e65e 3453
75bd7f81 3454}
c630aafd 3455
75bd7f81 3456//_____________________________________________________________________________
4551ea7c 3457void AliTRDtracker::CookdEdxTimBin(AliTRDtrack &TRDtrack)
eab5961e 3458{
75bd7f81 3459 //
b08a2e7e 3460 // This is setting fdEdxPlane and fTimBinPlane
3461 // Sums up the charge in each plane for track TRDtrack and also get the
3462 // Time bin for Max. Cluster
eab5961e 3463 // Prashant Shukla (shukla@physi.uni-heidelberg.de)
b08a2e7e 3464 //
eab5961e 3465
6d45eaef 3466 Double_t clscharge[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
3467 Double_t maxclscharge[AliESDtrack::kNPlane];
3468 Int_t nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
3469 Int_t timebin[AliESDtrack::kNPlane];
3470
b08a2e7e 3471 // Initialization of cluster charge per plane.
6d45eaef 3472 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
3473 for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
b08a2e7e 3474 clscharge[iPlane][iSlice] = 0.0;
4551ea7c 3475 nCluster[iPlane][iSlice] = 0;
6d45eaef 3476 }
b08a2e7e 3477 }
3478
3479 // Initialization of cluster charge per plane.
3480 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
4551ea7c 3481 timebin[iPlane] = -1;
eab5961e 3482 maxclscharge[iPlane] = 0.0;
3483 }
3484
3485 // Loop through all clusters associated to track TRDtrack
b08a2e7e 3486 Int_t nClus = TRDtrack.GetNumberOfClusters(); // from Kalmantrack
3487 for (Int_t iClus = 0; iClus < nClus; iClus++) {
eab5961e 3488 Double_t charge = TRDtrack.GetClusterdQdl(iClus);
4551ea7c 3489 Int_t index = TRDtrack.GetClusterIndex(iClus);
b08a2e7e 3490 AliTRDcluster *pTRDcluster = (AliTRDcluster *) GetCluster(index);
3491 if (!pTRDcluster) {
3492 continue;
3493 }
3494 Int_t tb = pTRDcluster->GetLocalTimeBin();
3495 if (!tb) {
3496 continue;
3497 }
4551ea7c 3498 Int_t detector = pTRDcluster->GetDetector();
3499 Int_t iPlane = fGeom->GetPlane(detector);
b08a2e7e 3500 if (iPlane >= AliESDtrack::kNPlane) {
3501 AliError(Form("Wrong plane %d",iPlane));
3502 continue;
3503 }
3504 Int_t iSlice = tb * AliESDtrack::kNSlice
3505 / AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
3506 if (iSlice >= AliESDtrack::kNSlice) {
3507 AliError(Form("Wrong slice %d",iSlice));
3508 continue;
3509 }
3510 clscharge[iPlane][iSlice] = clscharge[iPlane][iSlice] + charge;
4551ea7c 3511 if (charge > maxclscharge[iPlane]) {
eab5961e 3512 maxclscharge[iPlane] = charge;
4551ea7c 3513 timebin[iPlane] = tb;
eab5961e 3514 }
6d45eaef 3515 nCluster[iPlane][iSlice]++;
4551ea7c 3516 } // End of loop over cluster
eab5961e 3517
b08a2e7e 3518 // Setting the fdEdxPlane and fTimBinPlane variabales
3519 Double_t totalCharge = 0.0;
3520
f122c485 3521 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
b08a2e7e 3522 for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
3523 if (nCluster[iPlane][iSlice]) {
3524 clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice];
3525 }
4551ea7c 3526 TRDtrack.SetPIDsignals(clscharge[iPlane][iSlice],iPlane,iSlice);
b08a2e7e 3527 totalCharge = totalCharge+clscharge[iPlane][iSlice];
bd50219c 3528 }
b08a2e7e 3529 TRDtrack.SetPIDTimBin(timebin[iPlane],iPlane);
eab5961e 3530 }
b08a2e7e 3531
3532 // Still needed ????
3533 // Int_t i;
3534 // Int_t nc=TRDtrack.GetNumberOfClusters();
3535 // Float_t dedx=0;
3536 // for (i=0; i<nc; i++) dedx += TRDtrack.GetClusterdQdl(i);
3537 // dedx /= nc;
3538 // for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
3539 // TRDtrack.SetPIDsignals(dedx, iPlane);
3540 // TRDtrack.SetPIDTimBin(timbin[iPlane], iPlane);
3541 // }
3542
75bd7f81 3543}
c630aafd 3544
75bd7f81 3545//_____________________________________________________________________________
3546Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
4551ea7c 3547 , AliTRDtrack *track
3548 , Int_t *clusters, AliTRDtracklet &tracklet)
4f1c04d3 3549{
6c94f330 3550 //
4f1c04d3 3551 //
75bd7f81 3552 // Try to find nearest clusters to the track in timebins from t0 to t1
4f1c04d3 3553 //
4551ea7c 3554 // Correction coeficients - depend on TRD parameters - to be changed accordingly
3555 //
3556
3557 Double_t x[100];
3558 Double_t yt[100];
3559 Double_t zt[100];
3560 Double_t xmean = 0.0; // Reference x
3561 Double_t dz[10][100];
3562 Double_t dy[10][100];
3563 Float_t zmean[100];
3564 Float_t nmean[100];
3565 Int_t clfound = 0;
3566 Int_t indexes[10][100]; // Indexes of the clusters in the road
3567 Int_t best[10][100]; // Index of best matching cluster
3568 AliTRDcluster *cl[10][100]; // Pointers to the clusters in the road
3569
3570 for (Int_t it = 0; it < 100; it++) {
3571 x[it] = 0.0;
3572 yt[it] = 0.0;
3573 zt[it] = 0.0;
3574 clusters[it] = -2;
3575 zmean[it] = 0.0;
3576 nmean[it] = 0.0;
3577 for (Int_t ih = 0; ih < 10;ih++) {
3578 indexes[ih][it] = -2; // Reset indexes1
3579 cl[ih][it] = 0;
3580 dz[ih][it] = -100.0;
3581 dy[ih][it] = -100.0;
3582 best[ih][it] = 0;
3583 }
3584 }
ad670fba 3585
4551ea7c 3586 Double_t x0 = track->GetX();
3587 Double_t sigmaz = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
3588 Int_t nall = 0;
3589 Int_t nfound = 0;
3590 Double_t h01 = 0.0;
3591 Int_t plane = -1;
3592 Int_t detector = -1;
3593 Float_t padlength = 0.0;
3594 AliTRDtrack track2(* track);
3595 Float_t snpy = track->GetSnp();
3596 Float_t tany = TMath::Sqrt(snpy*snpy / (1.0 - snpy*snpy));
3597 if (snpy < 0.0) {
3598 tany *= -1.0;
3599 }
ad670fba 3600
4551ea7c 3601 Double_t sy2 = ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt());
3602 Double_t sz2 = ExpectedSigmaZ2(x0,track->GetTgl());
3603 Double_t road = 15.0 * TMath::Sqrt(track->GetSigmaY2() + sy2);
3604 if (road > 6.0) {
3605 road = 6.0;
3606 }
7e448bcc 3607 //road = 20.0;
4551ea7c 3608
3609 for (Int_t it = 0; it < t1-t0; it++) {
3610
3611 Double_t maxChi2[2] = { fgkMaxChi2, fgkMaxChi2 };
3612 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(it+t0));
3613 if (timeBin == 0) {
3614 continue; // No indexes1
6c94f330 3615 }
4551ea7c 3616
3617 Int_t maxn = timeBin;
3618 x[it] = timeBin.GetX();
7ad19338 3619 track2.PropagateTo(x[it]);
6c94f330 3620 yt[it] = track2.GetY();
3621 zt[it] = track2.GetZ();
7ad19338 3622
4551ea7c 3623 Double_t y = yt[it];
3624 Double_t z = zt[it];
3625 Double_t chi2 = 1000000.0;
4f1c04d3 3626 nall++;
4551ea7c 3627
4f1c04d3 3628 //
4551ea7c 3629 // Find 2 nearest cluster at given time bin
6c94f330 3630 //
7e448bcc 3631 int checkPoint[4] = {0,0,0,0};
3632 double minY = 123456789;
3633 double minD[2] = {1,1};
3634
4551ea7c 3635 for (Int_t i = timeBin.Find(y - road); i < maxn; i++) {
7e448bcc 3636 //for (Int_t i = 0; i < maxn; i++) {
4551ea7c 3637
3638 AliTRDcluster *c = (AliTRDcluster *) (timeBin[i]);
7ad19338 3639 h01 = GetTiltFactor(c);
4551ea7c 3640 if (plane < 0) {
c6f438c0 3641 Int_t det = c->GetDetector();
4551ea7c 3642 plane = fGeom->GetPlane(det);
3643 padlength = TMath::Sqrt(c->GetSigmaZ2() * 12.0);
3644 }
7e448bcc 3645
4551ea7c 3646 //if (c->GetLocalTimeBin()==0) continue;
3647 if (c->GetY() > (y + road)) {
3648 break;
3649 }
7e448bcc 3650
3651 fHDeltaX->Fill(c->GetX() - x[it]);
3652 //printf("%f\t%f\t%f \n", c->GetX(), x[it], c->GetX()-x[it]);
3653
3654 if (TMath::Abs(c->GetY()-y) < TMath::Abs(minY)) {
3655 minY = c->GetY()-y;
3656 minD[0] = c->GetY()-y;
3657 minD[1] = c->GetZ()-z;
3658 }
3659
3660 checkPoint[0]++;
3661
3662 fHMinZ->Fill(c->GetZ() - z);
3663 if ((c->GetZ() - z) * (c->GetZ() - z) > 2 * (12.0 * sz2)) {
4551ea7c 3664 continue;
7ad19338 3665 }
7e448bcc 3666 checkPoint[1]++;
7ad19338 3667
4551ea7c 3668 Double_t dist = TMath::Abs(c->GetZ() - z);
7e448bcc 3669 if (dist > (0.5 * padlength + 6.0 * sigmaz)) { // 0.5
4551ea7c 3670 continue; // 6 sigma boundary cut
3671 }
7e448bcc 3672 checkPoint[2]++;
3673
4551ea7c 3674 Double_t cost = 0.0;
3675 // Sigma boundary cost function
3676 if (dist> (0.5 * padlength - sigmaz)){
3677 cost = (dist - 0.5*padlength) / (2.0 * sigmaz);
3678 if (cost > -1) {
3679 cost = (cost + 1.0) * (cost + 1.0);
3680 }
3681 else {
3682 cost = 0.0;
3683 }
3684 }
3685 //Int_t label = TMath::Abs(track->GetLabel());
3686 //if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue;
3687 chi2 = track2.GetPredictedChi2(c,h01) + cost;
7ad19338 3688 clfound++;
6c94f330 3689
4551ea7c 3690 if (chi2 > maxChi2[1]) {
3691 continue;
3692 }
7e448bcc 3693 checkPoint[3]++;
3694
4551ea7c 3695 detector = c->GetDetector();
3696 // Store the clusters in the road
3697 for (Int_t ih = 2; ih < 9; ih++) {
3698 if (cl[ih][it] == 0) {
3699 cl[ih][it] = c;
3700 indexes[ih][it] = timeBin.GetIndex(i); // Index - 9 - reserved for outliers
7ad19338 3701 break;
3702 }
4f1c04d3 3703 }
4551ea7c 3704
3705 if (chi2 < maxChi2[0]) {
7ad19338 3706 maxChi2[1] = maxChi2[0];
3707 maxChi2[0] = chi2;
3708 indexes[1][it] = indexes[0][it];
3709 cl[1][it] = cl[0][it];
3710 indexes[0][it] = timeBin.GetIndex(i);
3711 cl[0][it] = c;
3712 continue;
3713 }
4551ea7c 3714 maxChi2[1] = chi2;
3715 cl[1][it] = c;
3716 indexes[1][it] = timeBin.GetIndex(i);
3717
3718 }
7e448bcc 3719
3720 for(int iCheckPoint = 0; iCheckPoint<4; iCheckPoint++)
3721 fHFindCl[iCheckPoint]->Fill(checkPoint[iCheckPoint]);
3722
3723 if (checkPoint[3]) {
3724 if (track->GetPt() > 0) fHMinYPos->Fill(minY);
3725 else fHMinYNeg->Fill(minY);
3726
3727 fHMinD->Fill(minD[0], minD[1]);
3728 }
4551ea7c 3729
3730 if (cl[0][it]) {
7ad19338 3731 nfound++;
3732 xmean += x[it];
3733 }
4551ea7c 3734
4f1c04d3 3735 }
4551ea7c 3736
3737 if (nfound < 4) {
3738 return 0;
3739 }
3740 xmean /= Float_t(nfound); // Middle x
3741 track2.PropagateTo(xmean); // Propagate track to the center
3742
4f1c04d3 3743 //
4551ea7c 3744 // Choose one of the variants
7ad19338 3745 //
4551ea7c 3746 Int_t changes[10];
3747 Float_t sumz = 0.0;
3748 Float_t sum = 0.0;
3749 Double_t sumdy = 0.0;
3750 Double_t sumdy2 = 0.0;
3751 Double_t sumx = 0.0;
3752 Double_t sumxy = 0.0;
3753 Double_t sumx2 = 0.0;
3754 Double_t mpads = 0.0;
3755
3756 Int_t ngood[10];
3757 Int_t nbad[10];
3758
7ad19338 3759 Double_t meanz[10];
4551ea7c 3760 Double_t moffset[10]; // Mean offset
3761 Double_t mean[10]; // Mean value
3762 Double_t angle[10]; // Angle
3763
3764 Double_t smoffset[10]; // Sigma of mean offset
3765 Double_t smean[10]; // Sigma of mean value
3766 Double_t sangle[10]; // Sigma of angle
3767 Double_t smeanangle[10]; // Correlation
3768
7ad19338 3769 Double_t sigmas[10];
4551ea7c 3770 Double_t tchi2s[10]; // Chi2s for tracklet
ad670fba 3771
4551ea7c 3772 for (Int_t it = 0; it < 10; it++) {
ad670fba 3773
4551ea7c 3774 ngood[it] = 0;
3775 nbad[it] = 0;
3776
3777 meanz[it] = 0.0;
3778 moffset[it] = 0.0; // Mean offset
3779 mean[it] = 0.0; // Mean value
3780 angle[it] = 0.0; // Angle
3781
7e448bcc 3782 smoffset[it] = 1.0e5; // Sigma of mean offset
3783 smean[it] = 1.0e5; // Sigma of mean value
3784 sangle[it] = 1.0e5; // Sigma of angle
4551ea7c 3785 smeanangle[it] = 0.0; // Correlation
3786
7e448bcc 3787 sigmas[it] = 1.0e5;
3788 tchi2s[it] = 1.0e5; // Chi2s for tracklet
75bd7f81 3789
3790 }
3791
7ad19338 3792 //
4551ea7c 3793 // Calculate zmean
7ad19338 3794 //
4551ea7c 3795 for (Int_t it = 0; it < t1 - t0; it++) {
3796 if (!cl[0][it]) {
3797 continue;
3798 }
3799 for (Int_t dt = -3; dt <= 3; dt++) {
3800 if (it+dt < 0) {
3801 continue;
3802 }
3803 if (it+dt > t1-t0) {
3804 continue;
3805 }
3806 if (!cl[0][it+dt]) {
3807 continue;
3808 }
3809 zmean[it] += cl[0][it+dt]->GetZ();
3810 nmean[it] += 1.0;
7ad19338 3811 }
4551ea7c 3812 zmean[it] /= nmean[it];
7ad19338 3813 }
4551ea7c 3814
3815 for (Int_t it = 0; it < t1 - t0; it++) {
3816
3817 best[0][it] = 0;
3818
3819 for (Int_t ih = 0; ih < 10; ih++) {
3820 dz[ih][it] = -100.0;
3821 dy[ih][it] = -100.0;
3822 if (!cl[ih][it]) {
3823 continue;
3824 }
3825 Double_t xcluster = cl[ih][it]->GetX();
3826 Double_t ytrack;
3827 Double_t ztrack;
3828 track2.GetProlongation(xcluster,ytrack,ztrack );
3829 dz[ih][it] = cl[ih][it]->GetZ()- ztrack; // Calculate distance from track in z
3830 dy[ih][it] = cl[ih][it]->GetY() + dz[ih][it]*h01 - ytrack; // and in y
7ad19338 3831 }
4551ea7c 3832
3833 // Minimize changes
3834 if (!cl[0][it]) {
3835 continue;
3836 }
3837 if ((TMath::Abs(cl[0][it]->GetZ()-zmean[it]) > padlength * 0.8) &&
3838 (cl[1][it])) {
3839 if (TMath::Abs(cl[1][it]->GetZ()-zmean[it]) < padlength * 0.5) {
3840 best[0][it] = 1;
4f1c04d3 3841 }
4551ea7c 3842 }
3843
7ad19338 3844 }
4551ea7c 3845
7ad19338 3846 //
4551ea7c 3847 // Iterative choice of "best path"
6c94f330 3848 //
4551ea7c 3849 Int_t label = TMath::Abs(track->GetLabel());
3850 Int_t bestiter = 0;
3851
3852 for (Int_t iter = 0; iter < 9; iter++) {
3853
3854 changes[iter] = 0;
3855 sumz = 0;
3856 sum = 0;
3857 sumdy = 0;
3858 sumdy2 = 0;
3859 sumx = 0;
3860 sumx2 = 0;
3861 sumxy = 0;
3862 mpads = 0;
3863 ngood[iter] = 0;
3864 nbad[iter] = 0;
3865
3866 // Linear fit
3867 for (Int_t it = 0; it < t1 - t0; it++) {
3868
3869 if (!cl[best[iter][it]][it]) {
3870 continue;
3871 }
3872
3873 // Calculates pad-row changes
3874 Double_t zbefore = cl[best[iter][it]][it]->GetZ();
3875 Double_t zafter = cl[best[iter][it]][it]->GetZ();
3876 for (Int_t itd = it - 1; itd >= 0; itd--) {
7ad19338 3877 if (cl[best[iter][itd]][itd]) {
4551ea7c 3878 zbefore = cl[best[iter][itd]][itd]->GetZ();
7ad19338 3879 break;
3880 }
3881 }
4551ea7c 3882 for (Int_t itd = it + 1; itd < t1 - t0; itd++) {
7ad19338 3883 if (cl[best[iter][itd]][itd]) {
4551ea7c 3884 zafter = cl[best[iter][itd]][itd]->GetZ();
7ad19338 3885 break;
3886 }
3887 }
4551ea7c 3888 if ((TMath::Abs(cl[best[iter][it]][it]->GetZ()-zbefore) > 0.1) &&
3889 (TMath::Abs(cl[best[iter][it]][it]->GetZ()- zafter) > 0.1)) {
3890 changes[iter]++;
3891 }
3892
3893 Double_t dx = x[it]-xmean; // Distance to reference x
3894 sumz += cl[best[iter][it]][it]->GetZ();
4f1c04d3 3895 sum++;
4551ea7c 3896 sumdy += dy[best[iter][it]][it];
3897 sumdy2 += dy[best[iter][it]][it]*dy[best[iter][it]][it];
3898 sumx += dx;
3899 sumx2 += dx*dx;
7ad19338 3900 sumxy += dx*dy[best[iter][it]][it];
4551ea7c 3901 mpads += cl[best[iter][it]][it]->GetNPads();
3902 if ((cl[best[iter][it]][it]->GetLabel(0) == label) ||
3903 (cl[best[iter][it]][it]->GetLabel(1) == label) ||
3904 (cl[best[iter][it]][it]->GetLabel(2) == label)) {
7ad19338 3905 ngood[iter]++;
4f1c04d3 3906 }
4551ea7c 3907 else {
7ad19338 3908 nbad[iter]++;
4f1c04d3 3909 }
4551ea7c 3910
4f1c04d3 3911 }
4551ea7c 3912
7ad19338 3913 //
6c94f330 3914 // calculates line parameters
7ad19338 3915 //
4551ea7c 3916 Double_t det = sum*sumx2 - sumx*sumx;
3917 angle[iter] = (sum*sumxy - sumx*sumdy) / det;
3918 mean[iter] = (sumx2*sumdy - sumx*sumxy) / det;
3919 meanz[iter] = sumz / sum;
3920 moffset[iter] = sumdy / sum;
3921 mpads /= sum; // Mean number of pads
3922
3923 Double_t sigma2 = 0.0; // Normalized residuals - for line fit
3924 Double_t sigma1 = 0.0; // Normalized residuals - constant fit
3925
3926 for (Int_t it = 0; it < t1 - t0; it++) {
3927 if (!cl[best[iter][it]][it]) {
3928 continue;
3929 }
3930 Double_t dx = x[it] - xmean;
3931 Double_t ytr = mean[iter] + angle[iter] * dx;
3932 sigma2 += (dy[best[iter][it]][it] - ytr)
3933 * (dy[best[iter][it]][it] - ytr);
3934 sigma1 += (dy[best[iter][it]][it] - moffset[iter])
3935 * (dy[best[iter][it]][it] - moffset[iter]);
7ad19338 3936 sum++;
4f1c04d3 3937 }
4551ea7c 3938 sigma2 /= (sum - 2); // Normalized residuals
3939 sigma1 /= (sum - 1); // Normalized residuals
3940 smean[iter] = sigma2 * (sumx2 / det); // Estimated error2 of mean
3941 sangle[iter] = sigma2 * ( sum / det); // Estimated error2 of angle
3942 smeanangle[iter] = sigma2 * (-sumx / det); // Correlation
3943 sigmas[iter] = TMath::Sqrt(sigma1);
3944 smoffset[iter] = (sigma1 / sum) + 0.01*0.01; // Sigma of mean offset + unisochronity sigma
3945
6c94f330 3946 //
4551ea7c 3947 // Iterative choice of "better path"
6c94f330 3948 //
4551ea7c 3949 for (Int_t it = 0; it < t1 - t0; it++) {
3950
3951 if (!cl[best[iter][it]][it]) {
3952 continue;
3953 }
3954
3955 // Add unisochronity + angular effect contribution
3956 Double_t sigmatr2 = smoffset[iter] + 0.5*tany*tany;
3957 Double_t sweight = 1.0/sigmatr2 + 1.0/track->GetSigmaY2();
3958 Double_t weighty = (moffset[iter] / sigmatr2) / sweight; // Weighted mean
3959 Double_t sigmacl = TMath::Sqrt(sigma1*sigma1 + track->GetSigmaY2());
3960 Double_t mindist = 100000.0;
3961 Int_t ihbest = 0;
3962
3963 for (Int_t ih = 0; ih < 10; ih++) {
3964 if (!cl[ih][it]) {
3965 break;
3966 }
3967 Double_t dist2 = (dy[ih][it] - weighty) / sigmacl;
3968 dist2 *= dist2; // Chi2 distance
3969 if (dist2 < mindist) {
7ad19338 3970 mindist = dist2;
4551ea7c 3971 ihbest = ih;
7ad19338 3972 }
3973 }
4551ea7c 3974 best[iter+1][it] = ihbest;
4f1c04d3 3975 }
4551ea7c 3976
4f1c04d3 3977 //
4551ea7c 3978 // Update best hypothesy if better chi2 according tracklet position and angle
7ad19338 3979 //
69b55c55 3980 Double_t sy2 = smean[iter] + track->GetSigmaY2();
4551ea7c 3981 Double_t sa2 = sangle[iter] + track->GetSigmaSnp2(); // track->fCee;
3982 Double_t say = track->GetSigmaSnpY(); // track->fCey;
3983 //Double_t chi20 = mean[bestiter]*mean[bestiter ] / sy2+angle[bestiter]*angle[bestiter]/sa2;
3984 //Double_t chi21 = mean[iter]*mean[iter] / sy2+angle[iter]*angle[iter]/sa2;
6c94f330 3985
4551ea7c 3986 Double_t detchi = sy2*sa2 - say*say;
3987 Double_t invers[3] = {sa2/detchi,sy2/detchi,-say/detchi}; // Inverse value of covariance matrix
4f1c04d3 3988
4551ea7c 3989 Double_t chi20 = mean[bestiter] * mean[bestiter] * invers[0]
3990 + angle[bestiter] * angle[bestiter] * invers[1]
3991 + 2.0 * mean[bestiter] * angle[bestiter] * invers[2];
3992 Double_t chi21 = mean[iter] * mean[iter] * invers[0]
3993 + angle[iter] * angle[iter] * invers[1]
3994 + 2.0 * mean[iter] * angle[iter] * invers[2];
3995 tchi2s[iter] = chi21;
3996
3997 if ((changes[iter] <= changes[bestiter]) &&
3998 (chi21 < chi20)) {
3999 bestiter = iter;
7ad19338 4000 }
4551ea7c 4001
7ad19338 4002 }
4551ea7c 4003
7ad19338 4004 //
4551ea7c 4005 // Set clusters
7ad19338 4006 //
4551ea7c 4007 Double_t sigma2 = sigmas[0]; // Choose as sigma from 0 iteration
4008 Short_t maxpos = -1;
4009 Float_t maxcharge = 0.0;
4010 Short_t maxpos4 = -1;
4011 Float_t maxcharge4 = 0.0;
4012 Short_t maxpos5 = -1;
4013 Float_t maxcharge5 = 0.0;
8979685e 4014
7ad19338 4015 //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
4016 //if (tchi2s[bestiter]>25.) sigma2=1000.; // dont'accept
4017
4e009ce4 4018 Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
4019 ,-AliTracker::GetBz()*0.1);
4551ea7c 4020 Double_t expectederr = sigma2*sigma2 + 0.01*0.01;
4021 if (mpads > 3.5) {
4022 expectederr += (mpads - 3.5) * 0.04;
4023 }
4024 if (changes[bestiter] > 1) {
4025 expectederr += changes[bestiter] * 0.01;
4026 }
4027 expectederr += (0.03 * (tany-exB)*(tany-exB)) * 15.0;
4028 //if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.;
7ad19338 4029 //expectederr+=10000;
4551ea7c 4030
4031 for (Int_t it = 0; it < t1 - t0; it++) {
4032
4033 if (!cl[best[bestiter][it]][it]) {
4034 continue;
69b55c55 4035 }
4551ea7c 4036
4037 cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); // Set cluster error
4038 if (!cl[best[bestiter][it]][it]->IsUsed()) {
4039 cl[best[bestiter][it]][it]->SetY(cl[best[bestiter][it]][it]->GetY());
4040 //cl[best[bestiter][it]][it]->Use();
4041 }
4042
4043 // Time bins with maximal charge
4044 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge) {
69b55c55 4045 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 4046 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 4047 }
6c94f330 4048
4551ea7c 4049 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
4050 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
69b55c55 4051 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 4052 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 4053 }
4054 }
4551ea7c 4055
4056 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
4057 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
69b55c55 4058 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 4059 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 4060 }
7ad19338 4061 }
4551ea7c 4062
4063 // Time bins with maximal charge
4064 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge) {
8979685e 4065 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 4066 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 4067 }
6c94f330 4068
4551ea7c 4069 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
4070 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
8979685e 4071 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 4072 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 4073 }
4074 }
4551ea7c 4075
4076 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
4077 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
8979685e 4078 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4551ea7c 4079 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 4080 }
4081 }
4551ea7c 4082
7ad19338 4083 clusters[it+t0] = indexes[best[bestiter][it]][it];
4551ea7c 4084
4085 // Still needed ????
4086 //if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>4 &&
4087 //cl[best[bestiter][it]][it]->GetLocalTimeBin()<18) clusters[it+t0]
4088 // = indexes[best[bestiter][it]][it]; //Test
4089
7ad19338 4090 }
4551ea7c 4091
7ad19338 4092 //
4551ea7c 4093 // Set tracklet parameters
6c94f330 4094 //
4551ea7c 4095 Double_t trackleterr2 = smoffset[bestiter] + 0.01*0.01;
4096 if (mpads > 3.5) {
4097 trackleterr2 += (mpads - 3.5) * 0.04;
4098 }
4099 trackleterr2 += changes[bestiter] * 0.01;
4100 trackleterr2 *= TMath::Max(14.0 - nfound,1.0);
4101 trackleterr2 += 0.2 * (tany-exB)*(tany-exB);
4102
4103 // Set tracklet parameters
4104 tracklet.Set(xmean
4105 ,track2.GetY() + moffset[bestiter]
4106 ,meanz[bestiter]
4107 ,track2.GetAlpha()
4108 ,trackleterr2);
7ad19338 4109 tracklet.SetTilt(h01);
4110 tracklet.SetP0(mean[bestiter]);
4111 tracklet.SetP1(angle[bestiter]);
4112 tracklet.SetN(nfound);
4113 tracklet.SetNCross(changes[bestiter]);
4114 tracklet.SetPlane(plane);
4115 tracklet.SetSigma2(expectederr);
4116 tracklet.SetChi2(tchi2s[bestiter]);
8979685e 4117 tracklet.SetMaxPos(maxpos,maxpos4,maxpos5);
27eaf44b 4118 track->SetTracklets(plane,tracklet);
4119 track->SetNWrong(track->GetNWrong() + nbad[0]);
4551ea7c 4120
7ad19338 4121 //
4122 // Debuging part
4123 //
69b55c55 4124 TClonesArray array0("AliTRDcluster");
4125 TClonesArray array1("AliTRDcluster");
4551ea7c 4126 array0.ExpandCreateFast(t1 - t0 + 1);
4127 array1.ExpandCreateFast(t1 - t0 + 1);
4128 TTreeSRedirector &cstream = *fDebugStreamer;
7ad19338 4129 AliTRDcluster dummy;
4130 Double_t dy0[100];
8979685e 4131 Double_t dyb[100];
4132
4551ea7c 4133 for (Int_t it = 0; it < t1 - t0; it++) {
7ad19338 4134 dy0[it] = dy[0][it];
4135 dyb[it] = dy[best[bestiter][it]][it];
4551ea7c 4136 if (cl[0][it]) {
7ad19338 4137 new(array0[it]) AliTRDcluster(*cl[0][it]);
4138 }
4551ea7c 4139 else {
7ad19338 4140 new(array0[it]) AliTRDcluster(dummy);
4141 }
4142 if(cl[best[bestiter][it]][it]) {
4143 new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
4144 }
4145 else{
4146 new(array1[it]) AliTRDcluster(dummy);
4147 }
4f1c04d3 4148 }
7ad19338 4149 TGraph graph0(t1-t0,x,dy0);
4150 TGraph graph1(t1-t0,x,dyb);
4151 TGraph graphy(t1-t0,x,yt);
4152 TGraph graphz(t1-t0,x,zt);
4551ea7c 4153
4154 if (AliTRDReconstructor::StreamLevel() > 0) {
4155 cstream << "tracklet"
4156 << "track.=" << track // Track parameters
4157 << "tany=" << tany // Tangent of the local track angle
4158 << "xmean=" << xmean // Xmean - reference x of tracklet
4159 << "tilt=" << h01 // Tilt angle
4160 << "nall=" << nall // Number of foundable clusters
4161 << "nfound=" << nfound // Number of found clusters
4162 << "clfound=" << clfound // Total number of found clusters in road
4163 << "mpads=" << mpads // Mean number of pads per cluster
4164 << "plane=" << plane // Plane number
4165 << "detector=" << detector // Detector number
4166 << "road=" << road // The width of the used road
4167 << "graph0.=" << &graph0 // x - y = dy for closest cluster
4168 << "graph1.=" << &graph1 // x - y = dy for second closest cluster
4169 << "graphy.=" << &graphy // y position of the track
4170 << "graphz.=" << &graphz // z position of the track
4171 //<< "fCl.=" << &array0 // closest cluster
4172 //<< "fCl2.=" << &array1 // second closest cluster
4173 << "maxpos=" << maxpos // Maximal charge postion
4174 << "maxcharge=" << maxcharge // Maximal charge
4175 << "maxpos4=" << maxpos4 // Maximal charge postion - after bin 4
4176 << "maxcharge4=" << maxcharge4 // Maximal charge - after bin 4
4177 << "maxpos5=" << maxpos5 // Maximal charge postion - after bin 5
4178 << "maxcharge5=" << maxcharge5 // Maximal charge - after bin 5
4179 << "bestiter=" << bestiter // Best iteration number
4180 << "tracklet.=" << &tracklet // Corrspond to the best iteration
4181 << "tchi20=" << tchi2s[0] // Chi2 of cluster in the 0 iteration
4182 << "tchi2b=" << tchi2s[bestiter] // Chi2 of cluster in the best iteration
4183 << "sigmas0=" << sigmas[0] // Residuals sigma
4184 << "sigmasb=" << sigmas[bestiter] // Residulas sigma
4185 << "ngood0=" << ngood[0] // Number of good clusters in 0 iteration
4186 << "nbad0=" << nbad[0] // Number of bad clusters in 0 iteration
4187 << "ngoodb=" << ngood[bestiter] // in best iteration
4188 << "nbadb=" << nbad[bestiter] // in best iteration
4189 << "changes0=" << changes[0] // Changes of pardrows in iteration number 0
4190 << "changesb=" << changes[bestiter] // Changes of pardrows in best iteration
4191 << "moffset0=" << moffset[0] // Offset fixing angle in iter=0
4192 << "smoffset0=" << smoffset[0] // Sigma of offset fixing angle in iter=0
4193 << "moffsetb=" << moffset[bestiter] // Offset fixing angle in iter=best
4194 << "smoffsetb=" << smoffset[bestiter] // Sigma of offset fixing angle in iter=best
4195 << "mean0=" << mean[0] // Mean dy in iter=0;
4196 << "smean0=" << smean[0] // Sigma of mean dy in iter=0
4197 << "meanb=" << mean[bestiter] // Mean dy in iter=best
4198 << "smeanb=" << smean[bestiter] // Sigma of mean dy in iter=best
4199 << "angle0=" << angle[0] // Angle deviation in the iteration number 0
4200 << "sangle0=" << sangle[0] // Sigma of angular deviation in iteration number 0
4201 << "angleb=" << angle[bestiter] // Angle deviation in the best iteration
4202 << "sangleb=" << sangle[bestiter] // Sigma of angle deviation in the best iteration
4203 << "expectederr=" << expectederr // Expected error of cluster position
4204 << "\n";
4205 }
4206
4f1c04d3 4207 return nfound;
4f1c04d3 4208
75bd7f81 4209}
4f1c04d3 4210
75bd7f81 4211//_____________________________________________________________________________
4212Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist
4213 , Int_t *outlist, Bool_t down)
69b55c55 4214{
4215 //
4551ea7c 4216 // Sort eleements according occurancy
4217 // The size of output array has is 2*n
69b55c55 4218 //
75bd7f81 4219
4551ea7c 4220 Int_t *sindexS = new Int_t[n]; // Temporary array for sorting
4221 Int_t *sindexF = new Int_t[2*n];
4222 for (Int_t i = 0; i < n; i++) {
4223 sindexF[i] = 0;
4224 }
4225
4226 TMath::Sort(n,inlist,sindexS,down);
4227
4228 Int_t last = inlist[sindexS[0]];
4229 Int_t val = last;
4230 sindexF[0] = 1;
4231 sindexF[0+n] = last;
4232 Int_t countPos = 0;
4233
4234 // Find frequency
4235 for (Int_t i = 1; i < n; i++) {
69b55c55 4236 val = inlist[sindexS[i]];
4551ea7c 4237 if (last == val) {
4238 sindexF[countPos]++;
4239 }
4240 else {
69b55c55 4241 countPos++;
4242 sindexF[countPos+n] = val;
4243 sindexF[countPos]++;
4551ea7c 4244 last = val;
69b55c55 4245 }
4246 }
4551ea7c 4247 if (last == val) {
4248 countPos++;
4249 }
4250
4251 // Sort according frequency
4252 TMath::Sort(countPos,sindexF,sindexS,kTRUE);
4253
4254 for (Int_t i = 0; i < countPos; i++) {
69b55c55 4255 outlist[2*i ] = sindexF[sindexS[i]+n];
4256 outlist[2*i+1] = sindexF[sindexS[i]];
4257 }
4551ea7c 4258
69b55c55 4259 delete [] sindexS;
4260 delete [] sindexF;
4261
4262 return countPos;
75bd7f81 4263
69b55c55 4264}
4265
75bd7f81 4266//_____________________________________________________________________________
4551ea7c 4267AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params)
69b55c55 4268{
4269 //
75bd7f81 4270 // Register a seed
69b55c55 4271 //
4551ea7c 4272
4273 Double_t alpha = AliTRDgeometry::GetAlpha();
4274 Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
69b55c55 4275 Double_t c[15];
4551ea7c 4276
4277 c[ 0] = 0.2;
4278 c[ 1] = 0.0; c[ 2] = 2.0;
4279 c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
4280 c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0; c[ 9] = 0.1;
4281 c[10] = 0.0; c[11] = 0.0; c[12] = 0.0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
4282
4283 Int_t index = 0;
4284 AliTRDcluster *cl = 0;
4285
4286 for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
4287 if (seeds[ilayer].IsOK()) {
4288 for (Int_t itime = 22; itime > 0; itime--) {
ca1e1e5b 4289 if (seeds[ilayer].GetIndexes(itime) > 0) {
4290 index = seeds[ilayer].GetIndexes(itime);
4291 cl = seeds[ilayer].GetClusters(itime);
69b55c55 4292 break;
4293 }
4294 }
4295 }
4551ea7c 4296 if (index > 0) {
4297 break;
4298 }
69b55c55 4299 }
4551ea7c 4300 if (cl == 0) {
4301 return 0;
4302 }
4303
4304 AliTRDtrack *track = new AliTRDtrack(cl
4305 ,index
4306 ,&params[1]
4307 ,c
4308 ,params[0]
4309 ,params[6]*alpha+shift);
4310 track->PropagateTo(params[0]-5.0);
69b55c55 4311 track->ResetCovariance(1);
4551ea7c 4312
4313 Int_t rc = FollowBackProlongation(*track);
4314 if (rc < 30) {
69b55c55 4315 delete track;
4551ea7c 4316 track = 0;
4317 }
4318 else {
69b55c55 4319 track->CookdEdx();
4320 CookdEdxTimBin(*track);
4551ea7c 4321 CookLabel(track,0.9);
69b55c55 4322 }
69b55c55 4323
75bd7f81 4324 return track;
7e448bcc 4325}
4326
4327//////////////////////////////////////////////////////////////////////////////////////////
69b55c55 4328
7e448bcc 4329void AliTRDtracker::InitLogHists() {
4330
4331 fHBackfit = new TH1D("logTRD_backfit", "", 40, -0.5, 39.5);
4332 fHRefit = new TH1D("logTRD_refit", "", 40, -0.5, 39.5);
4333 fHClSearch = new TH1D("logTRD_clSearch", "", 60, -0.5, 59.5);
4334
4335 fHX = new TH1D("logTRD_X", ";x (cm)", 200, 50, 400);
4336 fHNCl = new TH1D("logTRD_ncl", "", 40, -0.5, 39.5);
4337 fHNClTrack = new TH1D("logTRD_nclTrack", "", 180, -0.5, 179.5);
4338
4339 fHMinYPos = new TH1D("logTRD_minYPos", ";#delta Y (cm)", 400, -6, 6);
4340 fHMinYNeg = new TH1D("logTRD_minYNeg", ";#delta Y (cm)", 400, -6, 6);
4341 fHMinZ = new TH1D("logTRD_minZ", ";#delta Z (cm)", 400, -20, 20);
4342 fHMinD = new TH2D("logTRD_minD", ";#delta Y (cm);#delta Z (cm)", 100, -6, 6, 100, -50, 50);
4343
4344 fHDeltaX = new TH1D("logTRD_deltaX", ";#delta X (cm)", 100, -5, 5);
4345 fHXCl = new TH1D("logTRD_xCl", ";cluster x position (cm)", 1000, 280, 380);
4346
4347 const char *nameFindCl[4] = {"logTRD_clY", "logTRD_clZ", "logTRD_clB", "logTRD_clG"};
4348
4349 for(int i=0; i<4; i++) {
4350 fHFindCl[i] = new TH1D(nameFindCl[i], "", 30, -0.5, 29.5);
4351 }
69b55c55 4352}
7e448bcc 4353
4354//////////////////////////////////////////////////////////////////////////////////////////
4355
4356void AliTRDtracker::SaveLogHists() {
4357
4358 TDirectory *sav = gDirectory;
4359 TFile *logFile = 0;
4360
4361 TSeqCollection *col = gROOT->GetListOfFiles();
4362 int N = col->GetEntries();
4363 for(int i=0; i<N; i++) {
4364 logFile = (TFile*)col->At(i);
4365 if (strstr(logFile->GetName(), "AliESDs.root")) break;
4366 }
4367
4368 logFile->cd();
4369 fHBackfit->Write(fHBackfit->GetName(), TObject::kOverwrite);
4370 fHRefit->Write(fHRefit->GetName(), TObject::kOverwrite);
4371 fHClSearch->Write(fHClSearch->GetName(), TObject::kOverwrite);
4372 fHX->Write(fHX->GetName(), TObject::kOverwrite);
4373 fHNCl->Write(fHNCl->GetName(), TObject::kOverwrite);
4374 fHNClTrack->Write(fHNClTrack->GetName(), TObject::kOverwrite);
4375
4376 fHMinYPos->Write(fHMinYPos->GetName(), TObject::kOverwrite);
4377 fHMinYNeg->Write(fHMinYNeg->GetName(), TObject::kOverwrite);
4378 fHMinD->Write(fHMinD->GetName(), TObject::kOverwrite);
4379 fHMinZ->Write(fHMinZ->GetName(), TObject::kOverwrite);
4380
4381 fHDeltaX->Write(fHDeltaX->GetName(), TObject::kOverwrite);
4382 fHXCl->Write(fHXCl->GetName(), TObject::kOverwrite);
4383
4384
4385 for(int i=0; i<4; i++)
4386 fHFindCl[i]->Write(fHFindCl[i]->GetName(), TObject::kOverwrite);
4387
4388 logFile->Flush();
4389
4390 sav->cd();
4391}
4392
4393//////////////////////////////////////////////////////////////////////////////////////////