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