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