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