Fix Coverity defects
[u/mrichter/AliRoot.git] / TRD / AliTRDtracker.cxx
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 //                                                                           //
26 ///////////////////////////////////////////////////////////////////////////////
27
28 #include <TBranch.h>
29 #include <TFile.h>
30 #include <TGraph.h>
31 #include <TLinearFitter.h>
32 #include <TObjArray.h> 
33 #include <TROOT.h>
34 #include <TTree.h>  
35 #include <TTreeStream.h>
36
37 #include "AliESDEvent.h"
38 #include "AliESDtrack.h"
39 //#include "AliAlignObj.h"
40 #include "AliGeomManager.h"
41 #include "AliRieman.h"
42 #include "AliTrackPointArray.h"
43
44 #include "AliTRDgeometry.h"
45 #include "AliTRDpadPlane.h"
46 #include "AliTRDgeometry.h"
47 #include "AliTRDcluster.h" 
48 #include "AliTRDtrack.h"
49 #include "AliTRDseed.h"
50 #include "AliTRDcalibDB.h"
51 #include "AliTRDCommonParam.h"
52 #include "AliTRDtracker.h"
53 #include "AliTRDReconstructor.h"
54 #include "AliTRDrecoParam.h"
55 #include "AliTRDCalibraFillHisto.h"
56
57 ClassImp(AliTRDtracker)
58
59 const  Float_t  AliTRDtracker::fgkMinClustersInTrack =  0.5;  //
60 const  Float_t  AliTRDtracker::fgkLabelFraction      =  0.8;  //
61 const  Double_t AliTRDtracker::fgkMaxChi2            = 12.0;  //
62 const  Double_t AliTRDtracker::fgkMaxSnp             =  0.95; // Maximum local sine of the azimuthal angle
63 const  Double_t AliTRDtracker::fgkMaxStep            =  2.0;  // Maximal step size in propagation 
64
65 //_____________________________________________________________________________
66 AliTRDtracker::AliTRDtracker(AliTRDReconstructor *rec)
67   :AliTracker()
68   ,fReconstructor(rec)
69   ,fGeom(0)
70   ,fNclusters(0)
71   ,fClusters(0)
72   ,fNseeds(0)
73   ,fSeeds(0)
74   ,fNtracks(0)
75   ,fTracks(0)
76   ,fTimeBinsPerPlane(0)
77   ,fAddTRDseeds(kFALSE)
78   ,fNoTilt(kFALSE)
79   ,fHBackfit(NULL)
80   ,fHClSearch(NULL)
81   ,fHRefit(NULL)
82   ,fHX(NULL)
83   ,fHNCl(NULL)
84   ,fHNClTrack(NULL)
85   ,fHMinYPos(NULL)
86   ,fHMinYNeg(NULL)
87   ,fHMinZ(NULL)
88   ,fHMinD(NULL)
89   ,fHDeltaX(NULL)
90   ,fHXCl(NULL)
91   ,fDebugStreamer(0)
92 {
93   //
94   // Default constructor
95   //
96
97   for (Int_t i = 0; i < kTrackingSectors; i++) {
98     fTrSec[i] = 0;
99   }
100
101   InitLogHists();
102
103
104
105 //_____________________________________________________________________________
106 AliTRDtracker::AliTRDtracker(const AliTRDtracker &t)
107   :AliTracker(t)
108   ,fReconstructor(t.fReconstructor)
109   ,fGeom(0)
110   ,fNclusters(0)
111   ,fClusters(0)
112   ,fNseeds(0)
113   ,fSeeds(0)
114   ,fNtracks(0)
115   ,fTracks(0)
116   ,fTimeBinsPerPlane(0)
117   ,fAddTRDseeds(kFALSE)
118   ,fNoTilt(kFALSE)
119   ,fHBackfit(NULL)
120   ,fHClSearch(NULL)
121   ,fHRefit(NULL)
122   ,fHX(NULL)
123   ,fHNCl(NULL)
124   ,fHNClTrack(NULL)
125   ,fHMinYPos(NULL)
126   ,fHMinYNeg(NULL)
127   ,fHMinZ(NULL)
128   ,fHMinD(NULL)
129   ,fHDeltaX(NULL)
130   ,fHXCl(NULL)
131   ,fDebugStreamer(0)
132 {
133   //
134   // Copy constructor
135   //
136
137 }
138
139 //_____________________________________________________________________________
140 AliTRDtracker::AliTRDtracker(const TFile */*geomfile*/, AliTRDReconstructor *rec)
141   :AliTracker()
142   ,fReconstructor(rec)
143   ,fGeom(0)
144   ,fNclusters(0)
145   ,fClusters(new TObjArray(2000))
146   ,fNseeds(0)
147   ,fSeeds(new TObjArray(2000))
148   ,fNtracks(0)
149   ,fTracks(new TObjArray(1000))
150   ,fTimeBinsPerPlane(0)
151   ,fAddTRDseeds(kFALSE)
152   ,fNoTilt(kFALSE)
153   ,fHBackfit(NULL)
154   ,fHClSearch(NULL)
155   ,fHRefit(NULL)
156   ,fHX(NULL)
157   ,fHNCl(NULL)
158   ,fHNClTrack(NULL)
159   ,fHMinYPos(NULL)
160   ,fHMinYNeg(NULL)
161   ,fHMinZ(NULL)
162   ,fHMinD(NULL)
163   ,fHDeltaX(NULL)
164   ,fHXCl(NULL)
165   ,fDebugStreamer(0)
166 {
167   // 
168   //  Main constructor
169   //  
170    
171   TDirectory *savedir = gDirectory; 
172
173   fGeom = new AliTRDgeometry();
174
175   for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
176     Int_t trS   = geomS;
177     fTrSec[trS] = new AliTRDtrackingSector(fGeom,geomS);
178   }
179
180   AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0);
181   Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
182   if (tiltAngle < 0.1) {
183     fNoTilt = kTRUE;
184   }
185
186   if (!AliTRDcalibDB::Instance()) {
187     AliFatal("Could not get calibration object");
188   }
189   fTimeBinsPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS();
190
191   fDebugStreamer    = new TTreeSRedirector("TRDdebug.root");
192
193   savedir->cd();
194   
195   InitLogHists();
196
197 }   
198
199 //_____________________________________________________________________________
200 AliTRDtracker::~AliTRDtracker()
201 {
202   //
203   // Destructor of AliTRDtracker 
204   //
205
206   if (fClusters) {
207     fClusters->Delete();
208     delete fClusters;
209   }
210
211   if (fTracks) {
212     fTracks->Delete();
213     delete fTracks;
214   }
215
216   if (fSeeds) {
217     fSeeds->Delete();
218     delete fSeeds;
219   }
220
221   if (fGeom) {
222     delete fGeom;  
223   }
224
225   for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
226     delete fTrSec[geomS];
227   }
228
229   if (fDebugStreamer) {    
230     delete fDebugStreamer;
231   }
232
233 }   
234
235 //_____________________________________________________________________________
236 Int_t  AliTRDtracker::LocalToGlobalID(Int_t lid)
237 {
238   //
239   // Transform internal TRD ID to global detector ID
240   //
241
242   Int_t  isector = fGeom->GetSector(lid);
243   Int_t  istack  = fGeom->GetStack(lid);
244   Int_t  ilayer  = fGeom->GetLayer(lid);
245
246   AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
247   switch (ilayer) {
248   case 0:
249     iLayer = AliGeomManager::kTRD1;
250     break;
251   case 1:
252     iLayer = AliGeomManager::kTRD2;
253     break;
254   case 2:
255     iLayer = AliGeomManager::kTRD3;
256     break;
257   case 3:
258     iLayer = AliGeomManager::kTRD4;
259     break;
260   case 4:
261     iLayer = AliGeomManager::kTRD5;
262     break;
263   case 5:
264     iLayer = AliGeomManager::kTRD6;
265     break;
266   };
267
268   Int_t    modId = isector * fGeom->Nstack() + istack;
269   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,modId);
270
271   return volid;
272
273 }
274
275 //_____________________________________________________________________________
276 Int_t  AliTRDtracker::GlobalToLocalID(Int_t gid)
277 {
278   //
279   // Transform global detector ID to local detector ID
280   // 
281
282   Int_t modId    = 0;
283   AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(gid,modId);
284
285   Int_t isector = modId / fGeom->Nstack();
286   Int_t istack  = modId % fGeom->Nstack();
287   Int_t iLayer  = -1;
288
289   switch (layerId) {
290   case AliGeomManager::kTRD1:
291     iLayer = 0;
292     break;
293   case AliGeomManager::kTRD2:
294     iLayer = 1;
295     break;
296   case AliGeomManager::kTRD3:
297     iLayer = 2;
298     break;
299   case AliGeomManager::kTRD4:
300     iLayer = 3;
301     break;
302   case AliGeomManager::kTRD5:
303     iLayer = 4;
304     break;
305   case AliGeomManager::kTRD6:
306     iLayer = 5;
307     break;
308   default:
309     iLayer =-1;
310   }
311
312   if (iLayer < 0) {
313     return -1;
314   }
315
316   Int_t lid = fGeom->GetDetector(iLayer,istack,isector);
317
318   return lid;
319
320 }
321
322 //_____________________________________________________________________________
323 Bool_t AliTRDtracker::AdjustSector(AliTRDtrack* const track) const 
324 {
325   //
326   // Rotates the track when necessary
327   //
328
329   Double_t alpha = AliTRDgeometry::GetAlpha(); 
330   Double_t y     = track->GetY();
331   Double_t ymax  = track->GetX()*TMath::Tan(0.5*alpha);
332
333   if      (y >  ymax) {
334     if (!track->Rotate( alpha)) {
335       return kFALSE;
336     }
337   } 
338   else if (y < -ymax) {
339     if (!track->Rotate(-alpha)) {
340       return kFALSE;   
341     }
342   } 
343
344   return kTRUE;
345
346 }
347
348 //_____________________________________________________________________________
349 AliTRDcluster *AliTRDtracker::GetCluster(AliTRDtrack * const track, Int_t plane
350                                        , Int_t timebin, UInt_t &index)
351 {
352   //
353   // Try to find cluster in the backup list
354   //
355
356   AliTRDcluster *cl =0;
357   const Int_t *indexes = track->GetBackupIndexes();
358
359   for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
360     if (indexes[i] == 0) {
361       break;  
362     }
363     AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
364     if (!cli) {
365       break;
366     }
367     if (cli->GetLocalTimeBin() != timebin) {
368       continue;
369     }
370     Int_t ilayer = fGeom->GetLayer(cli->GetDetector());
371     if (ilayer == plane) {
372       cl    = cli;
373       index = indexes[i];
374       break;
375     }
376   }
377
378   return cl;
379
380 }
381
382 //_____________________________________________________________________________
383 Int_t  AliTRDtracker::GetLastPlane(AliTRDtrack * const track)
384 {
385   //
386   // Return last updated plane
387   //
388
389   Int_t  lastplane = 0;
390   const Int_t *indexes   = track->GetBackupIndexes();
391
392   for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
393     AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
394     if (!cli) {
395       break;
396     }
397     Int_t ilayer = fGeom->GetLayer(cli->GetDetector());
398     if (ilayer > lastplane) {
399       lastplane = ilayer;
400     }
401   }
402
403   return lastplane;
404
405 }
406      
407 //_____________________________________________________________________________
408 Int_t AliTRDtracker::PropagateBack(AliESDEvent *event) 
409 {
410   //
411   // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
412   // backpropagated by the TPC tracker. Each seed is first propagated 
413   // to the TRD, and then its prolongation is searched in the TRD.
414   // If sufficiently long continuation of the track is found in the TRD
415   // the track is updated, otherwise it's stored as originaly defined 
416   // by the TPC tracker.   
417   //  
418
419         Int_t   found    = 0;     // number of tracks found
420         Float_t foundMin = 20.0;
421         
422         Int_t    nSeed   = event->GetNumberOfTracks();
423         if(!nSeed){
424                 // run stand alone tracking
425                 if (fReconstructor->IsSeeding()) Clusters2Tracks(event);
426                 return 0;
427         }
428         
429         Float_t *quality = new Float_t[nSeed];
430         Int_t   *index   = new Int_t[nSeed];
431         for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
432                 AliESDtrack *seed = event->GetTrack(iSeed);
433                 Double_t covariance[15];
434                 seed->GetExternalCovariance(covariance);
435                 quality[iSeed] = covariance[0] + covariance[2];
436         }
437         // Sort tracks according to covariance of local Y and Z
438         TMath::Sort(nSeed,quality,index,kFALSE);
439         
440         // Backpropagate all seeds
441         for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
442         
443                 // Get the seeds in sorted sequence
444                 AliESDtrack *seed = event->GetTrack(index[iSeed]);
445                 fHBackfit->Fill(0);   // All seeds
446         
447                 // Check the seed status
448                 ULong_t status = seed->GetStatus();
449                 if ((status & AliESDtrack::kTPCout) == 0) {
450                         fHBackfit->Fill(1); // TPC outer edge reached
451                         continue;
452                 }
453                 if ((status & AliESDtrack::kTRDout) != 0) {
454                         fHBackfit->Fill(2); // TRD outer edge reached (does this happen ?)
455                         continue;
456                 }
457         
458                 // Do the back prolongation
459                 Int_t   lbl         = seed->GetLabel();
460                 AliTRDtrack *track  = new AliTRDtrack(*seed);
461                 track->SetSeedLabel(lbl);
462                 seed->UpdateTrackParams(track,AliESDtrack::kTRDbackup); // Make backup
463                 fNseeds++;
464                 Float_t p4          = track->GetC();
465                 Int_t   expectedClr = FollowBackProlongation(*track);
466                 fHBackfit->Fill(3);   // Back prolongation done
467                 fHX->Fill(track->GetX());
468         
469                 if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
470                                 (track->Pt()                                     > 0.8)) {
471                         
472                         fHBackfit->Fill(4);
473                         
474                         //
475                         // Make backup for back propagation
476                         //
477                         
478                         Int_t foundClr = track->GetNumberOfClusters();
479                         if (foundClr >= foundMin) {
480                                 track->CookdEdx();
481                                 track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
482                                 CookLabel(track,1 - fgkLabelFraction);
483                                 if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
484                                 
485                                 
486                                 // Sign only gold tracks
487                                 if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
488                                         if ((seed->GetKinkIndex(0)      ==   0) &&
489                                                         (track->Pt()                <  1.5)) {
490                                                                         UseClusters(track);
491                                         }
492                                 }
493                                 Bool_t isGold = kFALSE;
494         
495                                 // Full gold track
496                                 if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
497                                         //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
498                                         if (track->GetBackupTrack()) {
499                                                 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
500                                         }
501                                         isGold = kTRUE;
502                                         //fHBackfit->Fill()
503                                 }
504         
505                                 // Almost gold track
506                                 if ((!isGold)  && (track->GetNCross() == 0) &&
507                                                 (track->GetChi2() / track->GetNumberOfClusters()  < 7)) {
508                                         //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
509                                         if (track->GetBackupTrack()) {
510                                                 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
511                                         }
512                                         isGold = kTRUE;
513                                 }
514                                 
515                                 if ((!isGold) && (track->GetBackupTrack())) {
516                                         if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) && ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
517                                                 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
518                                                 isGold = kTRUE;
519                                         }
520                                 }
521         
522                                 if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected())  > 0.4)) {
523                                         //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
524                                 }
525                         }
526                 }
527                 /**/
528         
529                 /**/
530                 // Debug part of tracking
531                 TTreeSRedirector &cstream = *fDebugStreamer;
532                 Int_t eventNrInFile = event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
533                 if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
534                         if (track->GetBackupTrack()) {
535                                 cstream << "Tracks"
536                                 << "EventNrInFile="  << eventNrInFile
537                                 << "ESD.="     << seed
538                                 << "trd.="     << track
539                                 << "trdback.=" << track->GetBackupTrack()
540                                 << "\n";
541                         }
542                         else {
543                                 cstream << "Tracks"
544                                 << "EventNrInFile="  << eventNrInFile
545                                 << "ESD.="     << seed
546                                 << "trd.="     << track
547                                 << "trdback.=" << track
548                                 << "\n";
549                         }
550                 }
551                 /**/
552         
553                 // Propagation to the TOF (I.Belikov)
554                 if (track->GetStop() == kFALSE) {
555                         fHBackfit->Fill(5);
556         
557                         Double_t xtof  = 371.0;
558                         Double_t xTOF0 = 370.0;
559                 
560                         Double_t c2    = track->GetSnp() + track->GetC() * (xtof - track->GetX());
561                         if (TMath::Abs(c2) >= 0.99) {
562                                 fHBackfit->Fill(6);
563                                 delete track;
564                                 continue;
565                         }
566                         
567                         PropagateToX(*track,xTOF0,fgkMaxStep);
568         
569                         // Energy losses taken to the account - check one more time
570                         c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
571                         if (TMath::Abs(c2) >= 0.99) {
572                                 fHBackfit->Fill(7);
573                                 delete track;
574                                 continue;
575                         }
576                         
577                         //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
578                         //      fHBackfit->Fill(7);
579                         //delete track;
580                         //      continue;
581                         //}
582         
583                         Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
584                         Double_t y;
585                         track->GetYAt(xtof,GetBz(),y);
586                         if (y >  ymax) {
587                                 if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
588                                         fHBackfit->Fill(8);
589                                         delete track;
590                                         continue;
591                                 }
592                         }
593                         else if (y < -ymax) {
594                                 if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
595                                         fHBackfit->Fill(9);
596                                         delete track;
597                                         continue;
598                                 }
599                         }
600                                         
601                         if (track->PropagateTo(xtof)) {
602                                 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
603                                 fHBackfit->Fill(10);
604         
605         seed->SetNumberOfTRDslices(AliTRDCalPID::kNSlicesLQ);
606                                 for (Int_t i = 0; i < AliTRDtrack::kNplane; i++) {
607                                         for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ; j++) {
608                                                 seed->SetTRDslice(track->GetPIDsignals(i,j),i,j);
609                                         }
610                                         seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
611                                 }
612                                 //seed->SetTRDtrack(new AliTRDtrack(*track));
613                                 if (track->GetNumberOfClusters() > foundMin) {
614                                         fHBackfit->Fill(11);
615                                                                 found++;
616                                 }
617                         }
618                 }
619                 else {                  
620                         fHBackfit->Fill(12);
621                         
622                         if ((track->GetNumberOfClusters() >              15) &&
623                                         (track->GetNumberOfClusters() > 0.5*expectedClr)) {
624                                 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
625                                 fHBackfit->Fill(13);
626         
627                                 //seed->SetStatus(AliESDtrack::kTRDStop);
628
629         seed->SetNumberOfTRDslices(AliTRDCalPID::kNSlicesLQ);
630                                 for (Int_t i = 0; i < AliTRDtrack::kNplane; i++) {
631                                         for (Int_t j = 0; j <AliTRDCalPID::kNSlicesLQ; j++) {
632                                                 seed->SetTRDslice(track->GetPIDsignals(i,j),i,j);
633                                         }
634                                         seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
635                                 }
636                                 //seed->SetTRDtrack(new AliTRDtrack(*track));
637                                 found++;
638                         }
639                 }
640         
641                 seed->SetTRDQuality(track->StatusForTOF());
642                 seed->SetTRDBudget(track->GetBudget(0));
643         
644                 fHBackfit->Fill(14);
645                 delete track;
646         }
647         
648         AliInfo(Form("Number of seeds: %d",fNseeds));
649         AliInfo(Form("Number of back propagated TRD tracks: %d",found));
650                 
651         fSeeds->Clear();
652         fNseeds = 0;
653         
654         delete [] index;
655         delete [] quality;
656         
657         SaveLogHists();
658   
659   return 0;
660 }
661
662 //_____________________________________________________________________________
663 Int_t AliTRDtracker::RefitInward(AliESDEvent *event)
664 {
665   //
666   // Refits tracks within the TRD. The ESD event is expected to contain seeds 
667   // at the outer part of the TRD. 
668   // The tracks are propagated to the innermost time bin 
669   // of the TRD and the ESD event is updated
670   // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
671   //
672
673   //Int_t   timeBins = fTrSec[0]->GetNumberOfTimeBins();
674   //Float_t foundMin = fgkMinClustersInTrack * timeBins; 
675   Int_t   nseed    = 0;
676   Int_t   found    = 0;
677   Int_t   pidQ     = 0;
678   //Int_t innerTB    = fTrSec[0]->GetInnerTimeBin();
679   AliTRDtrack seed2;
680   
681   // Calibration fill 2D
682   AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
683   if (!calibra) {
684     AliInfo("Could not get Calibra instance\n");
685   }
686   
687   Int_t n = event->GetNumberOfTracks();
688   for (Int_t i = 0; i < n; i++) {
689
690     AliESDtrack *seed = event->GetTrack(i);
691     new (&seed2) AliTRDtrack(*seed);
692     fHRefit->Fill(0);
693
694     if (seed2.GetX() < 270.0) {
695       seed->UpdateTrackParams(&seed2,AliESDtrack::kTRDbackup); // Backup TPC track - only update
696       fHRefit->Fill(1);
697       continue;
698     }
699
700     ULong_t status = seed->GetStatus();
701     if ((status & AliESDtrack::kTRDout) == 0) {
702       fHRefit->Fill(2);
703       continue;
704     }
705     if ((status & AliESDtrack::kTRDin)  != 0) {
706       fHRefit->Fill(3);
707       continue;
708     }
709     
710     nseed++; 
711     fHRefit->Fill(4);
712
713     seed2.ResetCovariance(50.0); 
714
715     AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
716     const Int_t *indexes2 = seed2.GetIndexes();
717     for (Int_t l = 0; l < AliTRDtrack::kNplane;++l) {
718       for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ;j++) {
719         pt->SetPIDsignals(seed2.GetPIDsignals(l,j),l,j);
720       }
721       pt->SetPIDTimBin(seed2.GetPIDTimBin(l),l);
722     }
723
724     Int_t *indexes3 = const_cast<Int_t *>(pt->GetBackupIndexes());
725     for (Int_t l = 0; l < 200;++l) {
726       if (indexes2[l] == 0) {
727         break;
728       }
729       indexes3[l] = indexes2[l];
730     }  
731         
732     FollowProlongation(*pt); 
733     pt->CookdEdx();
734     pt->CookdEdxTimBin(seed->GetID());
735
736     //calculate PID methods    
737     pt->SetPIDMethod(AliTRDtrack::kLQ);
738     pt->CookPID(pidQ);
739     seed->SetTRDpid(pt->GetPID());
740     seed->SetTRDntracklets(pidQ<<3);
741
742     // update calibration
743     if(calibra->GetHisto2d()) calibra->UpdateHistograms(pt);
744     found++;
745
746     Double_t xTPC = 250.0;
747     if (PropagateToX(*pt,xTPC,fgkMaxStep)) {
748       seed->UpdateTrackParams(pt,AliESDtrack::kTRDrefit);
749       fHRefit->Fill(5);
750
751       for (Int_t l = 0; l < AliTRDtrack::kNplane; ++l) {
752         for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ; j++) {
753           seed->SetTRDslice(pt->GetPIDsignals(l,j),l,j);
754         }
755         seed->SetTRDTimBin(pt->GetPIDTimBin(l),l);
756       }
757     } else {
758       // If not prolongation to TPC - propagate without update
759       fHRefit->Fill(5);
760       AliTRDtrack *seed2t = new AliTRDtrack(*seed);
761       seed2t->ResetCovariance(5.0); 
762       AliTRDtrack *pt2   = new AliTRDtrack(*seed2t,seed2t->GetAlpha());
763       delete seed2t;
764
765       if (PropagateToX(*pt2,xTPC,fgkMaxStep)) { 
766         pt2->CookdEdx(); 
767         pt2->CookdEdxTimBin(seed->GetID()); 
768         seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit);
769         fHRefit->Fill(6);
770
771         for (Int_t l = 0; l < AliTRDtrack::kNplane; ++l) {
772           for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ; j++) {
773             seed->SetTRDslice(pt2->GetPIDsignals(l,j),l,j);
774           }
775           seed->SetTRDTimBin(pt2->GetPIDTimBin(l),l);
776         }
777       }
778
779       // Add TRD track to ESDfriendTrack - maybe this tracks are
780       // not useful for post-processing - TODO make decision
781       if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0)  {
782         seed->AddCalibObject(new AliTRDtrack(*pt2/*, kTRUE*/));
783       }
784       delete pt2;
785
786     }
787
788     // Add TRD track to ESDfriendTrack
789     if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0)  {
790       seed->AddCalibObject(new AliTRDtrack(*pt/*, kTRUE*/));
791     }
792     delete pt;
793
794   }   
795
796   AliInfo(Form("Number of loaded seeds: %d",nseed));
797   AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
798
799   SaveLogHists();
800
801   return 0;
802
803 }
804
805 //_____________________________________________________________________________
806 Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t)
807 {
808   //
809   // Starting from current position on track=t this function tries
810   // to extrapolate the track up to timeBin=0 and to confirm prolongation
811   // if a close cluster is found. Returns the number of clusters
812   // expected to be found in sensitive layers
813   // GeoManager used to estimate mean density
814   //
815
816   Int_t    sector;
817   Int_t    lastplane = GetLastPlane(&t);
818   Double_t xx0 = 0.0;
819   Double_t xrho= 0.0;
820   Int_t    expectedNumberOfClusters = 0;
821
822   for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
823
824     Int_t row0    = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
825     Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
826
827     // Propagate track close to the plane if neccessary
828     Double_t currentx  = fTrSec[0]->GetLayer(rowlast)->GetX();
829     if (currentx < (-fgkMaxStep + t.GetX())) {
830       // Propagate closer to chamber - safety space fgkMaxStep
831       if (!PropagateToX(t,currentx+fgkMaxStep,fgkMaxStep)) {
832         break;
833       }
834     }
835
836     if (!AdjustSector(&t)) {
837       break;
838     }
839
840     // Get material budget
841     Double_t xyz0[3];
842     Double_t xyz1[3];
843     Double_t param[7]; 
844     Double_t x;
845     Double_t y;
846     Double_t z;
847
848     // Starting global position
849     t.GetXYZ(xyz0);
850     // End global position
851     x = fTrSec[0]->GetLayer(row0)->GetX();
852     if (!t.GetProlongation(x,y,z)) {
853       break;
854     }
855     xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
856     xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
857     xyz1[2] =  z;
858     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
859     xrho= param[0]*param[4];
860     xx0 = param[1]; // Get mean propagation parameters
861
862     // Flags for marking the track momentum and direction. The track is
863     // marked only if it has at least 1 cluster picked up in the current
864     // chamber.
865     Bool_t kUPDATED = kFALSE;
866     Bool_t kMARKED  = kFALSE;
867
868     // Propagate and update
869     sector = t.GetSector();
870     //for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
871     for (Int_t itime = 0 ; itime < GetTimeBinsPerPlane(); itime++) {
872
873       // Mark track kinematics
874       if (itime > 10 && kUPDATED && !kMARKED) {
875         t.SetTrackSegmentDirMom(iplane);
876         kMARKED = kTRUE;
877       }
878
879       Int_t ilayer = GetGlobalTimeBin(0,iplane,itime);
880       expectedNumberOfClusters++;
881       t.SetNExpected(t.GetNExpected() + 1);
882       if (t.GetX() > 345.0) {
883         t.SetNExpectedLast(t.GetNExpectedLast() + 1);
884       }
885       AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
886       AliTRDcluster *cl = 0;
887       UInt_t   index   = 0;
888       Double_t maxChi2 = fgkMaxChi2;
889       x = timeBin.GetX();
890
891       if (timeBin) {
892
893         AliTRDcluster *cl0 = timeBin[0];
894         if (!cl0) {
895           // No clusters in given time bin
896           continue;
897         }
898
899         Int_t layer = fGeom->GetLayer(cl0->GetDetector());
900         if (layer > lastplane) {
901           continue;
902         }
903
904         Int_t timebin = cl0->GetLocalTimeBin();
905         AliTRDcluster *cl2 = GetCluster(&t,layer,timebin,index);
906
907         if (cl2) {
908           cl = cl2;
909           //Double_t h01 = GetTiltFactor(cl);    //I.B's fix
910           //maxChi2=t.GetPredictedChi2(cl,h01);
911         }
912
913         if (cl) {
914
915           //if (cl->GetNPads()<5)
916           Double_t dxsample = timeBin.GetdX();
917           t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
918           Double_t h01      = GetTiltFactor(cl);
919           Int_t    det      = cl->GetDetector();
920           Int_t    llayer   = fGeom->GetLayer(det);
921
922           if (t.GetX() > 345.0) {
923             t.SetNLast(t.GetNLast() + 1);
924             t.SetChi2Last(t.GetChi2Last() + maxChi2);
925           }
926
927           Double_t xcluster = cl->GetX();
928           t.PropagateTo(xcluster,xx0,xrho);                     
929           if (!AdjustSector(&t)) {
930             break;           //I.B's fix
931           }
932
933           maxChi2 = t.GetPredictedChi2(cl,h01);                                 
934           if (maxChi2 < 1e+10) {
935             if (!t.UpdateMI(cl,maxChi2,index,h01,llayer)) {
936               // ????
937             } 
938             else {
939               //SetCluster(cl, GetNumberOfClusters()-1);
940               kUPDATED = kTRUE;
941             }
942           }
943
944         }
945
946       }
947
948     }
949
950   }
951
952   return expectedNumberOfClusters;
953
954 }                
955
956 //_____________________________________________________________________________
957 Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
958 {
959   //  
960   // Starting from current radial position of track <t> this function
961   // extrapolates the track up to the outer timebin and in the sensitive
962   // layers confirms prolongation if a close cluster is found. 
963   // Returns the number of clusters expected to be found in sensitive layers
964   // Uses the geomanager for material description
965   //
966   // return number of assigned clusters ? 
967   //
968
969   Int_t    sector;
970
971   Double_t xx0    = 0.0;
972   Double_t xrho   = 0.0;
973
974   Float_t  ratio0 = 0.0;
975
976   Int_t    expectedNumberOfClusters = 0;
977
978   AliTRDtracklet tracklet;
979
980   const Int_t kNclusters = 1000;
981   Int_t       clusters[kNclusters];
982   for (Int_t i = 0; i < kNclusters; i++) {
983     clusters[i] = -1;
984   }
985
986 //   // Calibration fill 2D
987 //   AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
988 //   if (!calibra) {
989 //     AliInfo("Could not get Calibra instance\n");
990 //   }
991 //   if (calibra->GetMITracking()) {
992 //     calibra->ResetTrack();
993 //   }
994
995   // Loop through the TRD planes
996   for (Int_t iplane = 0; iplane < AliTRDtrack::kNplane; iplane++) {
997
998     Int_t hb = iplane * 10;
999     fHClSearch->Fill(hb);
1000
1001     // Get the global time bin numbers for the first an last 
1002     // local time bin of the given plane
1003     Int_t    row0     = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
1004     Int_t    rowlast  = GetGlobalTimeBin(0,iplane,0);
1005
1006     // Get the X coordinates of the propagation layer for the first time bin
1007     Double_t currentx = fTrSec[0]->GetLayer(row0)->GetX();
1008     if (currentx < t.GetX()) {
1009       fHClSearch->Fill(hb+1);
1010       continue;
1011     }
1012
1013     // Propagate closer to the current chamber if neccessary 
1014     if (currentx > (fgkMaxStep + t.GetX())) {
1015       if (!PropagateToX(t,currentx-fgkMaxStep,fgkMaxStep)) {
1016         fHClSearch->Fill(hb+2);
1017         break;
1018       }
1019     }
1020
1021     // Rotate track to adjacent sector if neccessary
1022     if (!AdjustSector(&t)) {
1023       fHClSearch->Fill(hb+3);
1024       break;
1025     }
1026
1027     // Check whether azimuthal angle is getting too large
1028     if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) {
1029       fHClSearch->Fill(hb+4);
1030       break;
1031     }
1032
1033     Double_t xyz0[3];
1034     Double_t xyz1[3];
1035     Double_t param[7];
1036     Double_t x;
1037     Double_t y;
1038     Double_t z;
1039     // Global start position (beginning of chamber)
1040     t.GetXYZ(xyz0);   
1041     // X-position of the end of the chamber
1042     x = fTrSec[0]->GetLayer(rowlast)->GetX();
1043     // Get local Y and Z at the X-position of the end of the chamber
1044     if (!t.GetProlongation(x,y,z)) {
1045       fHClSearch->Fill(hb+5);
1046       break;
1047     }
1048     // Global end position (end of chamber)
1049     xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha()); 
1050     xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1051     xyz1[2] =  z;
1052
1053     // Calculate the mean material budget along the path inside the chamber
1054     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);    
1055     // The mean propagation parameters (density*length and radiation length)
1056     xrho = param[0]*param[4];
1057     xx0  = param[1]; 
1058
1059     // Find the clusters and tracklet along the path inside the chamber
1060     sector      = t.GetSector();
1061     Float_t ncl = FindClusters(sector,row0,rowlast,&t,clusters,tracklet);
1062     fHNCl->Fill(tracklet.GetN());
1063
1064     // Discard if in less than 1/3 of the available timebins 
1065     // clusters are found 
1066     if (tracklet.GetN() < GetTimeBinsPerPlane()/3) {
1067       fHClSearch->Fill(hb+6);
1068       continue;
1069     }
1070
1071     //
1072     // Propagate and update track
1073     //
1074     for (Int_t itime = GetTimeBinsPerPlane()-1; itime >= 0; itime--) {
1075       Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
1076       expectedNumberOfClusters++;       
1077       t.SetNExpected(t.GetNExpected() + 1);
1078       if (t.GetX() > 345.0) {
1079         t.SetNExpectedLast(t.GetNExpectedLast() + 1);
1080       }
1081       AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
1082       AliTRDcluster *cl = 0;
1083       UInt_t   index   = 0;
1084       Double_t maxChi2 = fgkMaxChi2;
1085       x = timeBin.GetX();
1086
1087       if (timeBin) {    
1088                                 if (clusters[ilayer] > 0) {
1089                                         index = clusters[ilayer];
1090                                         cl    = (AliTRDcluster *)GetCluster(index);
1091                                         //Double_t h01 = GetTiltFactor(cl);   // I.B's fix
1092                                                                 //maxChi2=t.GetPredictedChi2(cl,h01); //
1093                                 }
1094         
1095         if (cl) {
1096                                         //if (cl->GetNPads() < 5)
1097                                         Double_t dxsample = timeBin.GetdX();
1098                                         t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
1099                                                                 Double_t h01      = GetTiltFactor(cl);
1100                                         Int_t    det      = cl->GetDetector();
1101                                         Int_t    layer    = fGeom->GetLayer(det);
1102                                         if (t.GetX() > 345.0) {
1103                                                 t.SetNLast(t.GetNLast() + 1);
1104                                                 t.SetChi2Last(t.GetChi2Last() + maxChi2);
1105                                         }
1106                                         Double_t xcluster = cl->GetX();
1107                                         t.PropagateTo(xcluster,xx0,xrho);
1108                                         maxChi2 = t.GetPredictedChi2(cl,h01);
1109
1110                                         if (maxChi2<1e+10)
1111                                                 if (!t.UpdateMI(cl,maxChi2,index,h01,layer)) {
1112                                                         if (!t.Update(cl,maxChi2,index,h01)) {
1113                                         // ????
1114                                                         }
1115                                                 } // else SetCluster(cl, GetNumberOfClusters()-1); // A.Bercuci 25.07.07
1116   
1117
1118 //           if (calibra->GetMITracking()) {
1119 //             calibra->UpdateHistograms(cl,&t);
1120 //           }
1121
1122                                         // Reset material budget if 2 consecutive gold
1123                                         if (layer > 0) {
1124                                                 if ((t.GetTracklets(layer).GetN() + t.GetTracklets(layer-1).GetN()) > 20) {
1125                                                         t.SetBudget(2,0.0);
1126                                                 }
1127                                         }
1128                         
1129                                 }
1130
1131       }
1132
1133     }
1134
1135     ratio0 = ncl / Float_t(fTimeBinsPerPlane);
1136     Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);  
1137     if ((tracklet.GetChi2()      <  18.0) && 
1138         (ratio0                  >   0.8) && 
1139         (ratio1                  >   0.6) && 
1140         (ratio0+ratio1           >   1.5) && 
1141         (t.GetNCross()           ==    0) && 
1142         (TMath::Abs(t.GetSnp())  <  0.85) &&
1143         (t.GetNumberOfClusters() >    20)){
1144       //if (ratio0 > 0.8) {
1145       t.MakeBackupTrack(); // Make backup of the track until is gold
1146     }
1147     
1148   }
1149
1150   return expectedNumberOfClusters;  
1151
1152 }         
1153
1154 //_____________________________________________________________________________
1155 Int_t AliTRDtracker::PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep)
1156 {
1157   //
1158   // Starting from current X-position of track <t> this function
1159   // extrapolates the track up to radial position <xToGo>. 
1160   // Returns 1 if track reaches the plane, and 0 otherwise 
1161   //
1162
1163   const Double_t kEpsilon = 0.00001;
1164
1165   // Current track X-position
1166   Double_t xpos = t.GetX();
1167
1168   // Direction: inward or outward
1169   Double_t dir  = (xpos < xToGo) ? 1.0 : -1.0;
1170
1171   while (((xToGo - xpos) * dir) > kEpsilon) {
1172
1173     Double_t xyz0[3];
1174     Double_t xyz1[3];
1175     Double_t param[7];
1176     Double_t x;
1177     Double_t y;
1178     Double_t z;
1179
1180     // The next step size
1181     Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos),maxStep);
1182
1183     // Get the global position of the starting point
1184     t.GetXYZ(xyz0);
1185
1186     // X-position after next step
1187     x = xpos + step;
1188
1189     // Get local Y and Z at the X-position of the next step
1190     if (!t.GetProlongation(x,y,z)) {
1191       return 0; // No prolongation possible
1192     }
1193
1194     // The global position of the end point of this prolongation step
1195     xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha()); 
1196     xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1197     xyz1[2] =  z;
1198
1199     // Calculate the mean material budget between start and
1200     // end point of this prolongation step
1201     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1202
1203     // Propagate the track to the X-position after the next step
1204     if (!t.PropagateTo(x,param[1],param[0]*param[4])) {
1205       return 0;
1206     }
1207
1208     // Rotate the track if necessary
1209     AdjustSector(&t);
1210
1211     // New track X-position
1212     xpos = t.GetX();
1213
1214   }
1215
1216   return 1;
1217
1218 }
1219
1220 //_____________________________________________________________________________
1221 Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1222 {
1223   //
1224   // Fills clusters into TRD tracking_sectors 
1225   // Note that the numbering scheme for the TRD tracking_sectors 
1226   // differs from that of TRD sectors
1227   //
1228
1229         
1230   if (ReadClusters(fClusters, cTree)) {
1231     AliError("Problem with reading the clusters !");
1232     return 1;
1233   }
1234   Int_t ncl  = fClusters->GetEntriesFast();
1235   fNclusters = ncl;
1236   AliInfo(Form("Sorting %d clusters",ncl));
1237               
1238   UInt_t index;
1239   while (ncl--) {
1240
1241     AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(ncl);
1242     Int_t detector       = c->GetDetector();
1243     Int_t localTimeBin   = c->GetLocalTimeBin();
1244     Int_t sector         = fGeom->GetSector(detector);
1245     Int_t layer          = fGeom->GetLayer(detector);
1246     Int_t trackingSector = sector;
1247
1248     //if (c->GetQ() > 10) {
1249     //  Int_t stack = fGeom->GetStack(detector);
1250     //}
1251
1252     Int_t gtb   = fTrSec[trackingSector]->CookTimeBinIndex(layer,localTimeBin);
1253     if (gtb < 0) {
1254       continue; 
1255     }
1256     Int_t trLayer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1257
1258     index = ncl;
1259
1260     fHXCl->Fill(c->GetX());
1261
1262     fTrSec[trackingSector]->GetLayer(trLayer)->SetX(c->GetX());
1263     fTrSec[trackingSector]->GetLayer(trLayer)->InsertCluster(c,index);
1264
1265   }
1266
1267   return 0;
1268
1269 }
1270
1271 //_____________________________________________________________________________
1272 void AliTRDtracker::UnloadClusters() 
1273
1274   //
1275   // Clears the arrays of clusters and tracks. Resets sectors and timebins 
1276   //
1277
1278   Int_t i;
1279   Int_t nentr;
1280
1281   nentr = fClusters->GetEntriesFast();
1282   for (i = 0; i < nentr; i++) {
1283     delete fClusters->RemoveAt(i);
1284   }
1285   fNclusters = 0;
1286
1287   nentr = fSeeds->GetEntriesFast();
1288   for (i = 0; i < nentr; i++) {
1289     delete fSeeds->RemoveAt(i);
1290   }
1291
1292   nentr = fTracks->GetEntriesFast();
1293   for (i = 0; i < nentr; i++) {
1294     delete fTracks->RemoveAt(i);
1295   }
1296
1297   Int_t nsec = AliTRDgeometry::kNsector;
1298   for (i = 0; i < nsec; i++) {    
1299     for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1300       fTrSec[i]->GetLayer(pl)->Clear();
1301     }
1302   }
1303
1304 }
1305
1306 //_____________________________________________________________________________
1307 Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *esd)
1308 {
1309   //
1310   // Creates  seeds using clusters between  position inner plane  and outer plane 
1311   //
1312
1313   const Double_t kMaxTheta =  1.0;
1314   const Double_t kMaxPhi   =  2.0;
1315
1316   const Double_t kRoad0y   =  6.0; // Road for middle cluster 
1317   const Double_t kRoad0z   =  8.5; // Road for middle cluster 
1318
1319   const Double_t kRoad1y   =  2.0; // Road in y for seeded cluster
1320   const Double_t kRoad1z   = 20.0; // Road in z for seeded cluster
1321
1322   const Double_t kRoad2y   =  3.0; // Road in y for extrapolated cluster
1323   const Double_t kRoad2z   = 20.0; // Road in z for extrapolated cluster
1324   const Int_t    kMaxSeed  = 3000;
1325
1326   Int_t maxSec = AliTRDgeometry::kNsector;  
1327
1328   // Linear fitters in planes
1329   TLinearFitter fitterTC(2,"hyp2");  // Fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
1330   TLinearFitter fitterT2(4,"hyp4");  // Fitting with tilting pads - kz not fixed
1331   fitterTC.StoreData(kTRUE);
1332   fitterT2.StoreData(kTRUE);
1333   AliRieman rieman(1000);            // Rieman fitter
1334   AliRieman rieman2(1000);           // Rieman fitter
1335
1336   // Find the maximal and minimal layer for the planes
1337   Int_t layers[6][2];
1338   AliTRDpropagationLayer *reflayers[6];
1339   for (Int_t i = 0; i < 6; i++) {
1340     layers[i][0] = 10000; 
1341     layers[i][1] = 0;
1342   }
1343   for (Int_t ns = 0; ns < maxSec; ns++) {
1344     for (Int_t ilayer = 0; ilayer < fTrSec[ns]->GetNumberOfLayers(); ilayer++) {
1345       AliTRDpropagationLayer &propLayer = *(fTrSec[ns]->GetLayer(ilayer));
1346       if (propLayer == 0) {
1347         continue;
1348       }
1349       Int_t det   = propLayer[0]->GetDetector();    
1350       Int_t layer = fGeom->GetLayer(det);
1351       if (ilayer < layers[layer][0]) {
1352         layers[layer][0] = ilayer;
1353       }
1354       if (ilayer > layers[layer][1]) {
1355         layers[layer][1] = ilayer;
1356       }
1357     }
1358   }
1359
1360   AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0);
1361   Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
1362   Double_t hL[6];                                                // Tilting angle
1363   Double_t xcl[6];                                               // X - position of reference cluster
1364   Double_t ycl[6];                                               // Y - position of reference cluster
1365   Double_t zcl[6];                                               // Z - position of reference cluster
1366
1367   AliTRDcluster *cl[6] = {    0,    0,    0,    0,    0,    0 }; // Seeding clusters
1368   Float_t padlength[6] = { 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 }; // Current pad-length 
1369
1370   Double_t chi2R  = 0.0;
1371   Double_t chi2Z  = 0.0;
1372   Double_t chi2RF = 0.0;
1373   Double_t chi2ZF = 0.0;
1374
1375   Int_t nclusters;                                               // Total number of clusters
1376   for (Int_t i = 0; i < 6; i++) {
1377     hL[i] = h01; 
1378     if (i%2==1) {
1379       hL[i]*=-1.0;
1380     }
1381   }
1382
1383   // Registered seed
1384   AliTRDseed *pseed = new AliTRDseed[kMaxSeed*6];
1385   AliTRDseed *seed[kMaxSeed];
1386   for (Int_t iseed = 0; iseed < kMaxSeed; iseed++) {
1387     seed[iseed]= &pseed[iseed*6];
1388   }
1389   AliTRDseed *cseed = seed[0];
1390
1391   Double_t seedquality[kMaxSeed];  
1392   Double_t seedquality2[kMaxSeed];  
1393   Double_t seedparams[kMaxSeed][7];
1394   Int_t    seedlayer[kMaxSeed];
1395   Int_t    registered = 0;
1396   Int_t    sort[kMaxSeed];
1397
1398   //
1399   // Seeding part
1400   //
1401   for (Int_t ns = 0; ns < maxSec; ns++) { // Loop over sectors
1402   //for (Int_t ns = 0; ns < 5; ns++) { // Loop over sectors
1403
1404     registered   = 0; // Reset registerd seed counter
1405     cseed        = seed[registered];
1406     Float_t iter = 0.0;
1407
1408     for (Int_t sLayer = 2; sLayer >= 0; sLayer--) {
1409       //for (Int_t dseed = 5; dseed < 15; dseed += 3) {
1410
1411       iter += 1.0;
1412       Int_t dseed = 5 + Int_t(iter) * 3;
1413
1414       // Initialize seeding layers
1415       for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
1416         reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
1417         xcl[ilayer]       = reflayers[ilayer]->GetX();
1418       }
1419
1420       Double_t xref = (xcl[sLayer+1] + xcl[sLayer+2]) * 0.5;      
1421       AliTRDpropagationLayer &layer0 = *reflayers[sLayer+0];
1422       AliTRDpropagationLayer &layer1 = *reflayers[sLayer+1];
1423       AliTRDpropagationLayer &layer2 = *reflayers[sLayer+2];
1424       AliTRDpropagationLayer &layer3 = *reflayers[sLayer+3];
1425
1426       Int_t maxn3 = layer3;
1427       for (Int_t icl3 = 0; icl3 < maxn3; icl3++) {
1428
1429         AliTRDcluster *cl3 = layer3[icl3];
1430         if (!cl3) {
1431           continue;
1432         }
1433         padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2() * 12.0);
1434         ycl[sLayer+3] = cl3->GetY();
1435         zcl[sLayer+3] = cl3->GetZ();
1436         Float_t yymin0 = ycl[sLayer+3] - 1.0 - kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1437         Float_t yymax0 = ycl[sLayer+3] + 1.0 + kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1438         Int_t   maxn0  = layer0;
1439
1440         for (Int_t icl0 = layer0.Find(yymin0); icl0 < maxn0; icl0++) {
1441
1442           AliTRDcluster *cl0 = layer0[icl0];
1443           if (!cl0) {
1444             continue;
1445           }
1446           if (cl3->IsUsed() && cl0->IsUsed()) {
1447             continue;
1448           }
1449           ycl[sLayer+0] = cl0->GetY();
1450           zcl[sLayer+0] = cl0->GetZ();
1451           if (ycl[sLayer+0] > yymax0) {
1452             break;
1453           }
1454           Double_t tanphi     = (ycl[sLayer+3]-ycl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]); 
1455           if (TMath::Abs(tanphi)   > kMaxPhi) {
1456             continue;
1457           }
1458           Double_t tantheta   = (zcl[sLayer+3]-zcl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]); 
1459           if (TMath::Abs(tantheta) > kMaxTheta) {
1460             continue; 
1461           }
1462           padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2() * 12.0);
1463
1464           // Expected position in 1 layer
1465           Double_t y1exp = ycl[sLayer+0] + (tanphi)   * (xcl[sLayer+1]-xcl[sLayer+0]);    
1466           Double_t z1exp = zcl[sLayer+0] + (tantheta) * (xcl[sLayer+1]-xcl[sLayer+0]);    
1467           Float_t yymin1 = y1exp - kRoad0y - tanphi;
1468           Float_t yymax1 = y1exp + kRoad0y + tanphi;
1469           Int_t   maxn1  = layer1;
1470
1471           for (Int_t icl1 = layer1.Find(yymin1); icl1 < maxn1; icl1++) {
1472
1473             AliTRDcluster *cl1 = layer1[icl1];
1474             if (!cl1) {
1475               continue;
1476             }
1477             Int_t nusedCl = 0;
1478             if (cl3->IsUsed()) nusedCl++;
1479             if (cl0->IsUsed()) nusedCl++;
1480             if (cl1->IsUsed()) nusedCl++;
1481             if (nusedCl > 1) {
1482               continue;
1483             }
1484             ycl[sLayer+1] = cl1->GetY();
1485             zcl[sLayer+1] = cl1->GetZ();
1486             if (ycl[sLayer+1] > yymax1) {
1487               break;
1488             }
1489             if (TMath::Abs(ycl[sLayer+1]-y1exp) > kRoad0y+tanphi) {
1490               continue;
1491             }
1492             if (TMath::Abs(zcl[sLayer+1]-z1exp) > kRoad0z) {
1493               continue;
1494             }
1495             padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2() * 12.0);
1496
1497             Double_t y2exp  = ycl[sLayer+0]+(tanphi)   * (xcl[sLayer+2]-xcl[sLayer+0]) + (ycl[sLayer+1]-y1exp);   
1498             Double_t z2exp  = zcl[sLayer+0]+(tantheta) * (xcl[sLayer+2]-xcl[sLayer+0]);
1499             Int_t    index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y,kRoad1z);
1500             if (index2 <= 0) {
1501               continue; 
1502             }
1503             AliTRDcluster *cl2 = (AliTRDcluster *) GetCluster(index2);
1504             padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2() * 12.0);
1505             ycl[sLayer+2]       = cl2->GetY();
1506             zcl[sLayer+2]       = cl2->GetZ();
1507             if (TMath::Abs(cl2->GetZ()-z2exp) > kRoad0z) {
1508               continue;
1509             }
1510
1511             rieman.Reset();
1512             rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
1513             rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
1514             rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);        
1515             rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
1516             rieman.Update();
1517
1518             // Reset fitter
1519             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1520               cseed[iLayer].Reset();
1521             }     
1522             chi2Z = 0.0; 
1523             chi2R = 0.0;
1524
1525             for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1526               cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1527               chi2Z += (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer])
1528                      * (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer]);
1529               cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));             
1530               cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1531               chi2R += (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer])
1532                      * (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer]);
1533               cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
1534             }
1535             if (TMath::Sqrt(chi2R) > 1.0/iter) {
1536               continue;
1537             }
1538             if (TMath::Sqrt(chi2Z) > 7.0/iter) {
1539               continue;
1540             }
1541
1542             Float_t minmax[2] = { -100.0,  100.0 };
1543             for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1544               Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer] * 0.5
1545                           + 1.0 - cseed[sLayer+iLayer].GetZref(0);
1546               if (max < minmax[1]) {
1547                 minmax[1] = max; 
1548               }
1549               Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer] * 0.5
1550                           - 1.0 - cseed[sLayer+iLayer].GetZref(0);
1551               if (min > minmax[0]) {
1552                 minmax[0] = min; 
1553               }
1554             }
1555
1556             Bool_t isFake = kFALSE; 
1557             if (cl0->GetLabel(0) != cl3->GetLabel(0)) {
1558               isFake = kTRUE;
1559             }
1560             if (cl1->GetLabel(0) != cl3->GetLabel(0)) {
1561               isFake = kTRUE;
1562             }
1563             if (cl2->GetLabel(0) != cl3->GetLabel(0)) {
1564               isFake = kTRUE;
1565             }
1566
1567             if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
1568               if ((!isFake) || ((icl3%10) == 0)) {  // Debugging print
1569                 TTreeSRedirector &cstream = *fDebugStreamer;
1570                 cstream << "Seeds0"
1571                         << "isFake=" << isFake
1572                         << "Cl0.="   << cl0
1573                         << "Cl1.="   << cl1
1574                         << "Cl2.="   << cl2
1575                         << "Cl3.="   << cl3
1576                         << "Xref="   << xref
1577                         << "X0="     << xcl[sLayer+0]
1578                         << "X1="     << xcl[sLayer+1]
1579                         << "X2="     << xcl[sLayer+2]
1580                         << "X3="     << xcl[sLayer+3]
1581                         << "Y2exp="  << y2exp
1582                         << "Z2exp="  << z2exp
1583                         << "Chi2R="  << chi2R
1584                         << "Chi2Z="  << chi2Z
1585                         << "Seed0.=" << &cseed[sLayer+0]
1586                         << "Seed1.=" << &cseed[sLayer+1]
1587                         << "Seed2.=" << &cseed[sLayer+2]
1588                         << "Seed3.=" << &cseed[sLayer+3]
1589                         << "Zmin="   << minmax[0]
1590                         << "Zmax="   << minmax[1]
1591                         << "\n";
1592               }
1593             }
1594
1595             ////////////////////////////////////////////////////////////////////////////////////
1596             //
1597             // Fit seeding part
1598             //
1599             ////////////////////////////////////////////////////////////////////////////////////
1600
1601             cl[sLayer+0] = cl0;
1602             cl[sLayer+1] = cl1;
1603             cl[sLayer+2] = cl2;
1604             cl[sLayer+3] = cl3;
1605             Bool_t isOK = kTRUE;
1606
1607             for (Int_t jLayer = 0; jLayer < 4; jLayer++) {
1608
1609               cseed[sLayer+jLayer].SetTilt(hL[sLayer+jLayer]);
1610               cseed[sLayer+jLayer].SetPadLength(padlength[sLayer+jLayer]);
1611               cseed[sLayer+jLayer].SetX0(xcl[sLayer+jLayer]);
1612
1613               for (Int_t jter = 0; jter < 2; jter++) {
1614
1615                 //
1616                 // In iteration 0 we try only one pad-row
1617                 // If quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
1618                 //
1619                 AliTRDseed tseed = cseed[sLayer+jLayer];
1620                 Float_t    roadz = padlength[sLayer+jLayer] * 0.5;
1621                 if (jter > 0) {
1622                   roadz = padlength[sLayer+jLayer];
1623                 }
1624
1625                 Float_t quality = 10000.0;
1626
1627                 for (Int_t iTime = 2; iTime < 20; iTime++) { 
1628
1629                   AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
1630                   Double_t dxlayer = layer.GetX() - xcl[sLayer+jLayer];          
1631                   Double_t zexp    = cl[sLayer+jLayer]->GetZ();
1632
1633                   if (jter > 0) {
1634                     // Try 2 pad-rows in second iteration
1635                     zexp = tseed.GetZref(0) + tseed.GetZref(1) * dxlayer;
1636                     if (zexp > cl[sLayer+jLayer]->GetZ()) {
1637                       zexp = cl[sLayer+jLayer]->GetZ() + padlength[sLayer+jLayer]*0.5;
1638                     }
1639                     if (zexp < cl[sLayer+jLayer]->GetZ()) {
1640                       zexp = cl[sLayer+jLayer]->GetZ() - padlength[sLayer+jLayer]*0.5;
1641                     }
1642                   }
1643
1644                   Double_t yexp  = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
1645                   Int_t    index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
1646                   if (index <= 0) {
1647                     continue; 
1648                   }
1649                   AliTRDcluster *clu = (AliTRDcluster *) GetCluster(index);           
1650
1651                   tseed.SetIndexes(iTime,index);
1652                   tseed.SetClusters(iTime,clu);  // Register cluster
1653                   tseed.SetX(iTime,dxlayer);     // Register cluster
1654                   tseed.SetY(iTime,clu->GetY()); // Register cluster
1655                   tseed.SetZ(iTime,clu->GetZ()); // Register cluster
1656
1657                 }
1658
1659                 tseed.Update();
1660
1661                 // Count the number of clusters and distortions into quality
1662                 Float_t dangle   = tseed.GetYfit(1) - tseed.GetYref(1);
1663                 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
1664                                  + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0))   / 0.2
1665                                  + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
1666                 if ((jter == 0) && tseed.IsOK()) {
1667                   cseed[sLayer+jLayer] = tseed;
1668                   quality              = tquality;
1669                   if (tquality < 5) {
1670                     break;  
1671                   }
1672                 }
1673                 if (tseed.IsOK() && (tquality < quality)) {
1674                   cseed[sLayer+jLayer] = tseed;
1675                 }
1676         
1677               } // Loop: jter
1678
1679               if (!cseed[sLayer+jLayer].IsOK()) {
1680                 isOK = kFALSE;
1681                 break;
1682               }
1683
1684               cseed[sLayer+jLayer].CookLabels();
1685               cseed[sLayer+jLayer].UpdateUsed();
1686               nusedCl += cseed[sLayer+jLayer].GetNUsed();
1687               if (nusedCl > 25) {
1688                 isOK = kFALSE;
1689                 break;
1690               }
1691
1692             } // Loop: jLayer
1693
1694             if (!isOK) {
1695               continue;
1696             }
1697             nclusters = 0;
1698             for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1699               if (cseed[sLayer+iLayer].IsOK()) {
1700                 nclusters += cseed[sLayer+iLayer].GetN2();          
1701               }
1702             }
1703
1704             // Iteration 0
1705             rieman.Reset();
1706             for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1707               rieman.AddPoint(xcl[sLayer+iLayer]
1708                              ,cseed[sLayer+iLayer].GetYfitR(0)
1709                              ,cseed[sLayer+iLayer].GetZProb()
1710                              ,1
1711                              ,10);
1712             }
1713             rieman.Update();
1714
1715             chi2R = 0.0; 
1716             chi2Z = 0.0;
1717
1718             for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1719               cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1720               chi2R += (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0))
1721                      * (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0));
1722               cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
1723               cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1724               chi2Z += (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz())
1725                      * (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz());
1726               cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
1727             }
1728             Double_t curv = rieman.GetC();
1729
1730             //
1731             // Likelihoods
1732             //
1733             Double_t sumda     = TMath::Abs(cseed[sLayer+0].GetYfitR(1) - cseed[sLayer+0].GetYref(1))
1734                                + TMath::Abs(cseed[sLayer+1].GetYfitR(1) - cseed[sLayer+1].GetYref(1))
1735                                + TMath::Abs(cseed[sLayer+2].GetYfitR(1) - cseed[sLayer+2].GetYref(1))
1736                                + TMath::Abs(cseed[sLayer+3].GetYfitR(1) - cseed[sLayer+3].GetYref(1));
1737             Double_t likea     = TMath::Exp(-sumda*10.6);
1738             Double_t likechi2  = 0.0000000001;
1739             if (chi2R < 0.5) {
1740               likechi2 += TMath::Exp(-TMath::Sqrt(chi2R) * 7.73);
1741             }
1742             Double_t likechi2z = TMath::Exp(-chi2Z * 0.088) / TMath::Exp(-chi2Z * 0.019);
1743             Double_t likeN     = TMath::Exp(-(72 - nclusters) * 0.19);
1744             Double_t like      = likea * likechi2 * likechi2z * likeN;
1745             Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetYref(1) - 130.0*curv) * 1.9);
1746             Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetZref(1)
1747                                                       - cseed[sLayer+0].GetZref(0) / xcl[sLayer+0]) * 5.9);
1748             Double_t likePrim  = TMath::Max(likePrimY*likePrimZ,0.0005);
1749                                             
1750             seedquality[registered] = like; 
1751             seedlayer[registered]   = sLayer;
1752             if (TMath::Log(0.000000000000001 + like) < -15) {
1753               continue;
1754             }
1755             AliTRDseed seedb[6];
1756             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1757               seedb[iLayer] = cseed[iLayer]; 
1758             }
1759
1760             ////////////////////////////////////////////////////////////////////////////////////
1761             //
1762             // Full track fit part
1763             //
1764             ////////////////////////////////////////////////////////////////////////////////////
1765
1766             Int_t nlayers  = 0;
1767             Int_t nusedf   = 0;
1768             Int_t findable = 0;
1769
1770             //
1771             // Add new layers  - avoid long extrapolation
1772             //
1773             Int_t tLayer[2] = { 0, 0 };
1774             if (sLayer == 2) {
1775               tLayer[0] = 1; 
1776               tLayer[1] = 0;
1777             }
1778             if (sLayer == 1) {
1779               tLayer[0] = 5; 
1780               tLayer[1] = 0;
1781             }
1782             if (sLayer == 0) {
1783               tLayer[0] = 4; 
1784               tLayer[1] = 5;
1785             }
1786
1787             for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
1788               Int_t jLayer = tLayer[iLayer]; // Set tracking layer            
1789               cseed[jLayer].Reset();
1790               cseed[jLayer].SetTilt(hL[jLayer]);
1791               cseed[jLayer].SetPadLength(padlength[jLayer]);
1792               cseed[jLayer].SetX0(xcl[jLayer]);
1793               // Get pad length and rough cluster
1794               Int_t indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].GetYref(0)
1795                                                                       ,cseed[jLayer].GetZref(0)
1796                                                                       ,kRoad2y
1797                                                                       ,kRoad2z);
1798               if (indexdummy <= 0) {
1799                 continue; 
1800               }
1801               AliTRDcluster *cldummy = (AliTRDcluster *) GetCluster(indexdummy);
1802               padlength[jLayer]      = TMath::Sqrt(cldummy->GetSigmaZ2() * 12.0);
1803             }
1804             AliTRDseed::FitRiemanTilt(cseed,kTRUE);
1805
1806             for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
1807
1808               Int_t    jLayer = tLayer[iLayer];      // set tracking layer      
1809               if ((jLayer == 0) && !(cseed[1].IsOK())) {
1810                 continue;  // break not allowed
1811               }
1812               if ((jLayer == 5) && !(cseed[4].IsOK())) {
1813                 continue;  // break not allowed
1814               }
1815               Float_t  zexp   = cseed[jLayer].GetZref(0);
1816               Double_t zroad  = padlength[jLayer] * 0.5 + 1.0;
1817
1818               for (Int_t jter = 0; jter < 2; jter++) {
1819
1820                 AliTRDseed tseed = cseed[jLayer];
1821                 Float_t quality = 10000.0;
1822
1823                 for (Int_t iTime = 2; iTime < 20; iTime++) { 
1824                   AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
1825                   Double_t dxlayer = layer.GetX()-xcl[jLayer];
1826                   Double_t yexp    = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
1827                   Float_t  yroad   = kRoad1y;
1828                   Int_t    index   = layer.FindNearestCluster(yexp,zexp,yroad,zroad);
1829                   if (index <= 0) {
1830                     continue; 
1831                   }
1832                   AliTRDcluster *clu = (AliTRDcluster *) GetCluster(index);           
1833                   tseed.SetIndexes(iTime,index);
1834                   tseed.SetClusters(iTime,clu);  // Register cluster
1835                   tseed.SetX(iTime,dxlayer);     // Register cluster
1836                   tseed.SetY(iTime,clu->GetY()); // Register cluster
1837                   tseed.SetZ(iTime,clu->GetZ()); // Register cluster
1838                 }
1839
1840                 tseed.Update();
1841                 if (tseed.IsOK()) {
1842                   Float_t dangle   = tseed.GetYfit(1) - tseed.GetYref(1);
1843                   Float_t tquality = (18.0 - tseed.GetN2())/2.0 + TMath::Abs(dangle) / 0.1
1844                                    + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
1845                                    + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
1846                   if (tquality < quality) {
1847                     cseed[jLayer] = tseed;
1848                     quality       = tquality;
1849                   }
1850                 }
1851
1852                 zroad *= 2.0;
1853
1854               } // Loop: jter
1855
1856               if ( cseed[jLayer].IsOK()) {
1857                 cseed[jLayer].CookLabels();
1858                 cseed[jLayer].UpdateUsed();
1859                 nusedf += cseed[jLayer].GetNUsed();
1860                 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
1861               }
1862
1863             } // Loop: iLayer
1864
1865             // Make copy
1866             AliTRDseed bseed[6];
1867             for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1868               bseed[jLayer] = cseed[jLayer];
1869             }
1870             Float_t lastquality = 10000.0;
1871             Float_t lastchi2    = 10000.0;
1872             Float_t chi2        =  1000.0;
1873
1874             for (Int_t jter = 0; jter < 4; jter++) {
1875
1876               // Sort tracklets according "quality", try to "improve" 4 worst 
1877               Float_t sumquality = 0.0;
1878               Float_t squality[6];
1879               Int_t   sortindexes[6];
1880
1881               for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1882                 if (bseed[jLayer].IsOK()) { 
1883                   AliTRDseed &tseed = bseed[jLayer];
1884                   Double_t zcor     = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
1885                   Float_t  dangle   = tseed.GetYfit(1) - tseed.GetYref(1);
1886                   Float_t  tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle)        / 0.1
1887                                     + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
1888                                     + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0))    / padlength[jLayer];
1889                   squality[jLayer]  = tquality;
1890                 }
1891                 else {
1892                   squality[jLayer] = -1.0;
1893                 }
1894                 sumquality +=squality[jLayer];
1895               }
1896
1897               if ((sumquality >= lastquality) ||  
1898                   (chi2       >     lastchi2)) {
1899                 break;
1900               }
1901               lastquality = sumquality;  
1902               lastchi2    = chi2;
1903               if (jter > 0) {
1904                 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1905                   cseed[jLayer] = bseed[jLayer];
1906                 }
1907               }
1908               TMath::Sort(6,squality,sortindexes,kFALSE);
1909
1910               for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
1911
1912                 Int_t bLayer = sortindexes[jLayer];
1913                 AliTRDseed tseed = bseed[bLayer];
1914
1915                 for (Int_t iTime = 2; iTime < 20; iTime++) { 
1916
1917                   AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
1918                   Double_t dxlayer = layer.GetX() - xcl[bLayer];
1919                   Double_t zexp    = tseed.GetZref(0);
1920                   Double_t zcor    = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
1921                   Float_t  roadz   = padlength[bLayer] + 1;
1922                   if (TMath::Abs(tseed.GetZProb() - zexp) > 0.5*padlength[bLayer]) {
1923                     roadz = padlength[bLayer] * 0.5;
1924                   }
1925                   if (tseed.GetZfit(1)*tseed.GetZref(1)   < 0.0) {
1926                     roadz = padlength[bLayer] * 0.5;
1927                   }
1928                   if (TMath::Abs(tseed.GetZProb() - zexp) < 0.1*padlength[bLayer]) {
1929                     zexp  = tseed.GetZProb(); 
1930                     roadz = padlength[bLayer] * 0.5;
1931                   }
1932
1933                   Double_t yexp  = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer - zcor;
1934                   Int_t    index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
1935                   if (index <= 0) {
1936                     continue; 
1937                   }
1938                   AliTRDcluster *clu = (AliTRDcluster *) GetCluster(index);           
1939
1940                   tseed.SetIndexes(iTime,index);
1941                   tseed.SetClusters(iTime,clu);  // Register cluster
1942                   tseed.SetX(iTime,dxlayer);     // Register cluster
1943                   tseed.SetY(iTime,clu->GetY()); // Register cluster
1944                   tseed.SetZ(iTime,clu->GetZ()); // Register cluster
1945
1946                 }
1947
1948                 tseed.Update();
1949                 if (tseed.IsOK()) {
1950                   Float_t  dangle   = tseed.GetYfit(1) - tseed.GetYref(1);
1951                   Double_t zcor     = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
1952                   Float_t  tquality = (18.0 - tseed.GetN2()) / 2.0 
1953                                     + TMath::Abs(dangle)     / 0.1
1954                                     + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
1955                                     + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0))    / padlength[jLayer];
1956                   if (tquality<squality[bLayer]) {
1957                     bseed[bLayer] = tseed;
1958                   }
1959                 }
1960
1961               } // Loop: jLayer
1962
1963               chi2 = AliTRDseed::FitRiemanTilt(bseed,kTRUE);
1964
1965             } // Loop: jter
1966
1967             nclusters = 0;
1968             nlayers   = 0;
1969             findable  = 0;
1970             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1971               if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) {
1972                 findable++;
1973               }
1974               if (cseed[iLayer].IsOK()) {
1975                 nclusters += cseed[iLayer].GetN2();         
1976                 nlayers++;
1977               }
1978             }
1979             if (nlayers < 3) {
1980               continue;
1981             }
1982             rieman.Reset();
1983             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1984               if (cseed[iLayer].IsOK()) {
1985                 rieman.AddPoint(xcl[iLayer]
1986                                ,cseed[iLayer].GetYfitR(0)
1987                                ,cseed[iLayer].GetZProb()
1988                                ,1
1989                                ,10);
1990               }
1991             }
1992             rieman.Update();
1993
1994             chi2RF = 0.0;
1995             chi2ZF = 0.0;
1996             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1997               if (cseed[iLayer].IsOK()) {
1998                 cseed[iLayer].SetYref(0,rieman.GetYat(xcl[iLayer]));
1999                 chi2RF += (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0))
2000                         * (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0));
2001                 cseed[iLayer].SetYref(1,rieman.GetDYat(xcl[iLayer]));
2002                 cseed[iLayer].SetZref(0,rieman.GetZat(xcl[iLayer]));
2003                 chi2ZF += (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz())
2004                         * (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz());
2005                 cseed[iLayer].SetZref(1,rieman.GetDZat(xcl[iLayer]));
2006               }
2007             }
2008             chi2RF /= TMath::Max((nlayers - 3.0),1.0);
2009             chi2ZF /= TMath::Max((nlayers - 3.0),1.0);
2010             curv    = rieman.GetC();
2011
2012             Double_t xref2 = (xcl[2] + xcl[3]) * 0.5; // Middle of the chamber
2013             Double_t dzmf  = rieman.GetDZat(xref2);
2014             Double_t zmf   = rieman.GetZat(xref2);
2015
2016             //
2017             // Fit hyperplane
2018             //
2019             Int_t npointsT = 0;
2020             fitterTC.ClearPoints();
2021             fitterT2.ClearPoints();
2022             rieman2.Reset();
2023
2024             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2025
2026               if (!cseed[iLayer].IsOK()) {
2027                 continue;
2028               }
2029
2030               for (Int_t itime = 0; itime < 25; itime++) {
2031
2032                 if (!cseed[iLayer].IsUsable(itime)) {
2033                   continue;
2034                 }
2035                 // X relative to the middle chamber
2036                 Double_t x  = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;  
2037                 Double_t y  = cseed[iLayer].GetY(itime);
2038                 Double_t z  = cseed[iLayer].GetZ(itime);
2039                 // ExB correction to the correction
2040                 // Tilted rieman
2041                 Double_t uvt[6];
2042                 // Global x
2043                 Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();      
2044                 Double_t t  = 1.0 / (x2*x2 + y*y);
2045                 uvt[1] = t;                 // t
2046                 uvt[0] = 2.0 * x2 * uvt[1]; // u 
2047                 uvt[2] = 2.0 * hL[iLayer] * uvt[1];
2048                 uvt[3] = 2.0 * hL[iLayer] * x * uvt[1];       
2049                 uvt[4] = 2.0 * (y + hL[iLayer]*z) * uvt[1];
2050                 Double_t error = 2.0 * 0.2 * uvt[1];
2051                 fitterT2.AddPoint(uvt,uvt[4],error);
2052
2053                 //
2054                 // Constrained rieman
2055                 // 
2056                 z = cseed[iLayer].GetZ(itime);
2057                 uvt[0] = 2.0 * x2 * t; // u 
2058                 uvt[1] = 2.0 * hL[iLayer] * x2 * uvt[1];              
2059                 uvt[2] = 2.0 * (y + hL[iLayer] * (z - GetZ())) * t;
2060                 fitterTC.AddPoint(uvt,uvt[2],error);
2061                 rieman2.AddPoint(x2,y,z,1,10);
2062                 npointsT++;
2063
2064               }
2065
2066             } // Loop: iLayer
2067
2068             rieman2.Update();
2069             fitterTC.Eval();
2070             fitterT2.Eval();
2071             Double_t rpolz0 = fitterT2.GetParameter(3);
2072             Double_t rpolz1 = fitterT2.GetParameter(4);     
2073
2074             //
2075             // Linear fitter  - not possible to make boundaries
2076             // Do not accept non possible z and dzdx combinations
2077             //      
2078             Bool_t acceptablez = kTRUE;
2079             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2080               if (cseed[iLayer].IsOK()) {
2081                 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
2082                 if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > padlength[iLayer] * 0.5 + 1.0) {
2083                   acceptablez = kFALSE;
2084                 }
2085               }
2086             }
2087             if (!acceptablez) {
2088               fitterT2.FixParameter(3,zmf);
2089               fitterT2.FixParameter(4,dzmf);
2090               fitterT2.Eval();
2091               fitterT2.ReleaseParameter(3);
2092               fitterT2.ReleaseParameter(4);
2093               rpolz0 = fitterT2.GetParameter(3);
2094               rpolz1 = fitterT2.GetParameter(4);
2095             }
2096
2097             Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
2098             Double_t chi2TC = fitterTC.GetChisquare() / Float_t(npointsT);
2099             Double_t polz1c = fitterTC.GetParameter(2);
2100             Double_t polz0c = polz1c * xref2;
2101             Double_t aC     =  fitterTC.GetParameter(0);
2102             Double_t bC     =  fitterTC.GetParameter(1);
2103             Double_t cC     =  aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
2104             Double_t aR     =  fitterT2.GetParameter(0);
2105             Double_t bR     =  fitterT2.GetParameter(1);
2106             Double_t dR     =  fitterT2.GetParameter(2);            
2107             Double_t cR     =  1.0 + bR*bR - dR*aR;
2108             Double_t dca    =  0.0;         
2109             if (cR > 0.0) {
2110               dca = -dR / (TMath::Sqrt(1.0 + bR*bR - dR*aR) + TMath::Sqrt(1.0 + bR*bR)); 
2111               cR  =  aR / TMath::Sqrt(cR);
2112             }
2113
2114             Double_t chi2ZT2 = 0.0;
2115             Double_t chi2ZTC = 0.0;
2116             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2117               if (cseed[iLayer].IsOK()) {
2118                 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
2119                 Double_t zTC = polz0c + polz1c * (xcl[iLayer] - xref2);
2120                 chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
2121                 chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
2122               }
2123             }
2124             chi2ZT2 /= TMath::Max((nlayers - 3.0),1.0);
2125             chi2ZTC /= TMath::Max((nlayers - 3.0),1.0);     
2126
2127             AliTRDseed::FitRiemanTilt(cseed,kTRUE);
2128             Float_t sumdaf = 0.0;
2129             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2130               if (cseed[iLayer].IsOK()) {
2131                 sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))
2132                                     / cseed[iLayer].GetSigmaY2());
2133               }
2134             }
2135             sumdaf /= Float_t (nlayers - 2.0);
2136
2137             //
2138             // Likelihoods for full track
2139             //
2140             Double_t likezf     = TMath::Exp(-chi2ZF * 0.14);
2141             Double_t likechi2C  = TMath::Exp(-chi2TC * 0.677);
2142             Double_t likechi2TR = TMath::Exp(-chi2TR * 0.78);
2143             Double_t likeaf     = TMath::Exp(-sumdaf * 3.23);
2144             seedquality2[registered] = likezf * likechi2TR * likeaf; 
2145
2146             // Still needed ????
2147 //          Bool_t isGold = kFALSE;
2148 //          
2149 //          if (nlayers == 6        && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE;   // gold
2150 //          if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE;   // gold
2151 //          if (isGold &&nusedf<10){
2152 //            for (Int_t jLayer=0;jLayer<6;jLayer++){
2153 //              if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
2154 //                seed[index][jLayer].UseClusters();  //sign gold
2155 //            }
2156 //          }
2157
2158             Int_t index0 = 0;
2159             if (!cseed[0].IsOK()) {
2160               index0 = 1;
2161               if (!cseed[1].IsOK()) {
2162                 index0 = 2;
2163               }
2164             }
2165             seedparams[registered][0] = cseed[index0].GetX0();
2166             seedparams[registered][1] = cseed[index0].GetYref(0);
2167             seedparams[registered][2] = cseed[index0].GetZref(0);
2168             seedparams[registered][5] = cR;
2169             seedparams[registered][3] = cseed[index0].GetX0() * cR - TMath::Sin(TMath::ATan(cseed[0].GetYref(1)));
2170             seedparams[registered][4] = cseed[index0].GetZref(1)
2171                                       / TMath::Sqrt(1.0 + cseed[index0].GetYref(1) * cseed[index0].GetYref(1));
2172             seedparams[registered][6] = ns;
2173
2174             Int_t labels[12];
2175             Int_t outlab[24];
2176             Int_t nlab = 0;
2177             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2178               if (!cseed[iLayer].IsOK()) {
2179                 continue;
2180               }
2181               if (cseed[iLayer].GetLabels(0) >= 0) {
2182                 labels[nlab] = cseed[iLayer].GetLabels(0);
2183                 nlab++;
2184               }
2185               if (cseed[iLayer].GetLabels(1) >= 0) {
2186                 labels[nlab] = cseed[iLayer].GetLabels(1);
2187                 nlab++;
2188               }
2189             }
2190             Freq(nlab,labels,outlab,kFALSE);
2191             Int_t label     = outlab[0];
2192             Int_t frequency = outlab[1];
2193             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2194               cseed[iLayer].SetFreq(frequency);
2195               cseed[iLayer].SetC(cR);
2196               cseed[iLayer].SetCC(cC);
2197               cseed[iLayer].SetChi2(chi2TR);
2198               cseed[iLayer].SetChi2Z(chi2ZF);
2199             }
2200
2201             // Debugging print
2202             if (1 || (!isFake)) {
2203               Float_t zvertex = GetZ();
2204               TTreeSRedirector &cstream = *fDebugStreamer;
2205               if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
2206                 cstream << "Seeds1"
2207                         << "isFake="     << isFake
2208                         << "Vertex="     << zvertex
2209                         << "Rieman2.="   << &rieman2
2210                         << "Rieman.="    << &rieman
2211                         << "Xref="       << xref
2212                         << "X0="         << xcl[0]
2213                         << "X1="         << xcl[1]
2214                         << "X2="         << xcl[2]
2215                         << "X3="         << xcl[3]
2216                         << "X4="         << xcl[4]
2217                         << "X5="         << xcl[5]
2218                         << "Chi2R="      << chi2R
2219                         << "Chi2Z="      << chi2Z
2220                         << "Chi2RF="     << chi2RF  // Chi2 of trackletes on full track
2221                         << "Chi2ZF="     << chi2ZF  // Chi2 z on tracklets on full track
2222                         << "Chi2ZT2="    << chi2ZT2 // Chi2 z on tracklets on full track  - rieman tilt
2223                         << "Chi2ZTC="    << chi2ZTC // Chi2 z on tracklets on full track  - rieman tilt const
2224                         << "Chi2TR="     << chi2TR  // Chi2 without vertex constrain
2225                         << "Chi2TC="     << chi2TC  // Chi2 with    vertex constrain
2226                         << "C="          << curv    // Non constrained - no tilt correction
2227                         << "DR="         << dR      // DR parameter    - tilt correction
2228                         << "DCA="        << dca     // DCA             - tilt correction
2229                         << "CR="         << cR      // Non constrained curvature - tilt correction
2230                         << "CC="         << cC      // Constrained curvature
2231                         << "Polz0="      << polz0c
2232                         << "Polz1="      << polz1c
2233                         << "RPolz0="     << rpolz0
2234                         << "RPolz1="     << rpolz1
2235                         << "Ncl="        << nclusters
2236                         << "Nlayers="    << nlayers
2237                         << "NUsedS="     << nusedCl
2238                         << "NUsed="      << nusedf
2239                         << "Findable="   << findable
2240                         << "Like="       << like
2241                         << "LikePrim="   << likePrim
2242                         << "Likechi2C="  << likechi2C
2243                         << "Likechi2TR=" << likechi2TR
2244                         << "Likezf="     << likezf
2245                         << "LikeF="      << seedquality2[registered]
2246                         << "S0.="        << &cseed[0]
2247                         << "S1.="        << &cseed[1]
2248                         << "S2.="        << &cseed[2]
2249                         << "S3.="        << &cseed[3]
2250                         << "S4.="        << &cseed[4]
2251                         << "S5.="        << &cseed[5]
2252                         << "SB0.="       << &seedb[0]
2253                         << "SB1.="       << &seedb[1]
2254                         << "SB2.="       << &seedb[2]
2255                         << "SB3.="       << &seedb[3]
2256                         << "SB4.="       << &seedb[4]
2257                         << "SB5.="       << &seedb[5]
2258                         << "Label="      << label
2259                         << "Freq="       << frequency
2260                         << "sLayer="     << sLayer
2261                         << "\n";
2262               }
2263             }
2264
2265             if (registered<kMaxSeed - 1) {
2266               registered++;
2267               cseed = seed[registered];
2268             }
2269
2270           } // End of loop over layer 1
2271
2272         } // End of loop over layer 0 
2273
2274       } // End of loop over layer 3     
2275
2276     } // End of loop over seeding time bins 
2277
2278     //
2279     // Choose best
2280     //
2281
2282     TMath::Sort(registered,seedquality2,sort,kTRUE);
2283     Bool_t signedseed[kMaxSeed];
2284     for (Int_t i = 0; i < registered; i++) {
2285       signedseed[i] = kFALSE;
2286     }
2287
2288     for (Int_t jter = 0; jter < 5; jter++) {
2289
2290       for (Int_t iseed = 0; iseed < registered; iseed++) {
2291       
2292         Int_t index = sort[iseed];
2293         if (signedseed[index]) {
2294           continue;
2295         }
2296         Int_t labelsall[1000];
2297         Int_t nlabelsall = 0;
2298         Int_t naccepted  = 0;;
2299         Int_t sLayer     = seedlayer[index];
2300         Int_t ncl        = 0;
2301         Int_t nused      = 0;
2302         Int_t nlayers    = 0;
2303         Int_t findable   = 0;
2304
2305         for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2306
2307           if (TMath::Abs(seed[index][jLayer].GetYref(0) / xcl[jLayer]) < 0.15) {
2308             findable++;
2309           }
2310           if (seed[index][jLayer].IsOK()) {
2311             seed[index][jLayer].UpdateUsed();
2312             ncl   +=seed[index][jLayer].GetN2();
2313             nused +=seed[index][jLayer].GetNUsed();
2314             nlayers++;
2315             // Cooking label
2316             for (Int_t itime = 0; itime < 25; itime++) {
2317               if (seed[index][jLayer].IsUsable(itime)) {
2318                 naccepted++;
2319                 for (Int_t ilab = 0; ilab < 3; ilab++) {
2320                   Int_t tindex = seed[index][jLayer].GetClusters(itime)->GetLabel(ilab);
2321                   if (tindex >= 0) {
2322                     labelsall[nlabelsall] = tindex;
2323                     nlabelsall++;
2324                   }
2325                 }
2326               }
2327             }
2328           }
2329
2330         }
2331
2332         if (nused > 30) {
2333           continue;
2334         }
2335
2336         if (jter == 0) {
2337           if (nlayers < 6) {
2338             continue;
2339           }
2340           if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2341             continue; // Gold
2342           }
2343         }
2344
2345         if (jter == 1) {
2346           if (nlayers < findable) {
2347             continue;
2348           }
2349           if (TMath::Log(0.000000001+seedquality2[index]) < -4.0) {
2350             continue;
2351           }
2352         }
2353
2354         if (jter == 2) {
2355           if ((nlayers == findable) || 
2356               (nlayers ==        6)) {
2357             continue;
2358           }
2359           if (TMath::Log(0.000000001+seedquality2[index]) < -6.0) {
2360             continue;
2361           }
2362         }
2363
2364         if (jter == 3) {
2365           if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2366             continue;
2367           }
2368         }
2369
2370         if (jter == 4) {
2371           if (TMath::Log(0.000000001+seedquality2[index]) - nused/(nlayers-3.0) < -15.0) {
2372             continue;
2373           }
2374         }
2375
2376         signedseed[index] = kTRUE;
2377
2378         Int_t labels[1000];
2379         Int_t outlab[1000];
2380         Int_t nlab = 0;
2381         for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2382           if (seed[index][iLayer].IsOK()) {
2383             if (seed[index][iLayer].GetLabels(0) >= 0) {
2384               labels[nlab] = seed[index][iLayer].GetLabels(0);
2385               nlab++;
2386             }
2387             if (seed[index][iLayer].GetLabels(1) >= 0) {
2388               labels[nlab] = seed[index][iLayer].GetLabels(1);
2389               nlab++;
2390             }
2391           }   
2392         }
2393         Freq(nlab,labels,outlab,kFALSE);
2394         Int_t   label     = outlab[0];
2395         Int_t   frequency = outlab[1];
2396         Freq(nlabelsall,labelsall,outlab,kFALSE);
2397         Int_t   label1    = outlab[0];
2398         Int_t   label2    = outlab[2];
2399         Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
2400         Float_t ratio     = Float_t(nused) / Float_t(ncl);
2401         if (ratio < 0.25) {
2402           for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2403             if ((seed[index][jLayer].IsOK()) && 
2404                 (TMath::Abs(seed[index][jLayer].GetYfit(1) - seed[index][jLayer].GetYfit(1)) < 0.2)) {
2405               seed[index][jLayer].UseClusters(); // Sign gold
2406             }
2407           }
2408         }
2409
2410         Int_t eventNrInFile = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
2411         TTreeSRedirector &cstream = *fDebugStreamer;
2412
2413         //
2414         // Register seed
2415         //
2416         AliTRDtrack *track = RegisterSeed(seed[index],seedparams[index]);
2417         AliTRDtrack  dummy;
2418         if (!track) {
2419           track = &dummy;
2420         }
2421         else {
2422           AliESDtrack esdtrack;
2423           esdtrack.UpdateTrackParams(track,AliESDtrack::kTRDout);
2424           esdtrack.SetLabel(label);
2425           esd->AddTrack(&esdtrack);     
2426           if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
2427             cstream << "Tracks"
2428                     << "EventNrInFile="  << eventNrInFile
2429                     << "ESD.="     << &esdtrack
2430                     << "trd.="     << track
2431                     << "trdback.=" << track
2432                     << "\n";
2433           }
2434         }
2435
2436         if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
2437           cstream << "Seeds2"
2438                   << "Iter="      << jter
2439                   << "Track.="    << track
2440                   << "Like="      << seedquality[index]
2441                   << "LikeF="     << seedquality2[index]
2442                   << "S0.="       << &seed[index][0]
2443                   << "S1.="       << &seed[index][1]
2444                   << "S2.="       << &seed[index][2]
2445                   << "S3.="       << &seed[index][3]
2446                   << "S4.="       << &seed[index][4]
2447                   << "S5.="       << &seed[index][5]
2448                   << "Label="     << label
2449                   << "Label1="    << label1
2450                   << "Label2="    << label2
2451                   << "FakeRatio=" << fakeratio
2452                   << "Freq="      << frequency
2453                   << "Ncl="       << ncl
2454                   << "Nlayers="   << nlayers
2455                   << "Findable="  << findable
2456                   << "NUsed="     << nused
2457                   << "sLayer="    << sLayer
2458                   << "EventNrInFile="   << eventNrInFile
2459                   << "\n";
2460         }
2461
2462       } // Loop: iseed
2463
2464     } // Loop: jter
2465
2466   } // End of loop over sectors
2467
2468   delete [] pseed;
2469
2470         return 0;
2471 }
2472           
2473 //_____________________________________________________________________________
2474 Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *clusterTree) const
2475 {
2476   //
2477   // Reads AliTRDclusters from the file. 
2478   // The names of the cluster tree and branches 
2479   // should match the ones used in AliTRDclusterizer::WriteClusters()
2480   //
2481
2482   Int_t nsize = Int_t(clusterTree->GetTotBytes() / (sizeof(AliTRDcluster))); 
2483   TObjArray *clusterArray = new TObjArray(nsize+1000); 
2484   
2485   TBranch *branch = clusterTree->GetBranch("TRDcluster");
2486   if (!branch) {
2487     AliError("Can't get the branch !");
2488     return 1;
2489   }
2490   branch->SetAddress(&clusterArray); 
2491   
2492   // Loop through all entries in the tree
2493   Int_t nEntries   = (Int_t) clusterTree->GetEntries();
2494   Int_t nbytes     = 0;
2495   AliTRDcluster *c = NULL;
2496   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {    
2497     
2498     // Import the tree
2499     nbytes += clusterTree->GetEvent(iEntry);  
2500     
2501     // Get the number of points in the detector
2502     Int_t nCluster = clusterArray->GetEntriesFast();  
2503     // Loop through all TRD digits
2504     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
2505       if(!(c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster))) continue;
2506       array->AddLast(c);
2507                         //printf("Add cluster 0x%x.\n", c);
2508       clusterArray->RemoveAt(iCluster); 
2509     }
2510
2511   }
2512   delete clusterArray;
2513
2514   return 0;
2515
2516 }
2517
2518 //_____________________________________________________________________________
2519 Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint &p) const
2520 {
2521   //
2522   // Get track space point with index i
2523   // Origin: C.Cheshkov
2524   //
2525
2526   AliTRDcluster *cl = (AliTRDcluster *) fClusters->UncheckedAt(index);
2527   Int_t  idet    = cl->GetDetector();
2528   Int_t  isector = fGeom->GetSector(idet);
2529   Int_t  istack  = fGeom->GetStack(idet);
2530   Int_t  ilayer  = fGeom->GetLayer(idet);
2531   Double_t local[3];
2532   local[0] = GetX(isector,ilayer,cl->GetLocalTimeBin());
2533   local[1] = cl->GetY();
2534   local[2] = cl->GetZ();
2535   Double_t global[3];
2536   fGeom->RotateBack(idet,local,global);
2537   p.SetXYZ(global[0],global[1],global[2]);
2538   AliGeomManager::ELayerID iGeoLayer = AliGeomManager::kTRD1;
2539   switch (ilayer) {
2540   case 0:
2541     iGeoLayer = AliGeomManager::kTRD1;
2542     break;
2543   case 1:
2544     iGeoLayer = AliGeomManager::kTRD2;
2545     break;
2546   case 2:
2547     iGeoLayer = AliGeomManager::kTRD3;
2548     break;
2549   case 3:
2550     iGeoLayer = AliGeomManager::kTRD4;
2551     break;
2552   case 4:
2553     iGeoLayer = AliGeomManager::kTRD5;
2554     break;
2555   case 5:
2556     iGeoLayer = AliGeomManager::kTRD6;
2557     break;
2558   };
2559   Int_t    modId = isector * fGeom->Nstack() + istack;
2560   UShort_t volid = AliGeomManager::LayerToVolUID(iGeoLayer,modId);
2561   p.SetVolumeID(volid);
2562
2563   return kTRUE;
2564
2565 }
2566
2567 //_____________________________________________________________________________
2568 void AliTRDtracker::CookLabel(AliKalmanTrack *pt, Float_t wrong) const 
2569 {
2570   //
2571   // This cooks a label. Mmmmh, smells good...
2572   //
2573
2574   Int_t label = 123456789;
2575   Int_t index;
2576   Int_t i;
2577   Int_t j;
2578   Int_t ncl   = pt->GetNumberOfClusters();
2579
2580   const Int_t kRange = fTrSec[0]->GetOuterTimeBin() + 1;
2581
2582   Bool_t labelAdded;
2583
2584   Int_t **s = new Int_t* [kRange];
2585   for (i = 0; i < kRange; i++) {
2586     s[i] = new Int_t[2];
2587   }
2588   for (i = 0; i < kRange; i++) {
2589     s[i][0] = -1;
2590     s[i][1] =  0;
2591   }
2592
2593   Int_t t0;
2594   Int_t t1;
2595   Int_t t2;
2596
2597   for (i = 0; i < ncl; i++) {
2598     index = pt->GetClusterIndex(i);
2599     AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2600     t0=c->GetLabel(0);
2601     t1=c->GetLabel(1);
2602     t2=c->GetLabel(2);
2603   }
2604
2605   for (i = 0; i < ncl; i++) {
2606     index = pt->GetClusterIndex(i);
2607     AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2608     for (Int_t k = 0; k < 3; k++) { 
2609       label      = c->GetLabel(k);
2610       labelAdded = kFALSE; 
2611       j = 0;
2612       if (label >= 0) {
2613         while ((!labelAdded) && (j < kRange)) {
2614           if ((s[j][0] == label) || 
2615               (s[j][1] ==     0)) {
2616             s[j][0] = label; 
2617             s[j][1] = s[j][1] + 1; 
2618             labelAdded = kTRUE;
2619           }
2620           j++;
2621         }
2622       }
2623     }
2624   }
2625
2626   Int_t max = 0;
2627   label = -123456789;
2628
2629   for (i = 0; i < kRange; i++) {
2630     if (s[i][1] > max) {
2631       max   = s[i][1]; 
2632       label = s[i][0];
2633     }
2634   }
2635
2636   for (i = 0; i < kRange; i++) {
2637     delete []s[i];
2638   }        
2639
2640   delete []s;
2641
2642   if ((1.0 - Float_t(max)/ncl) > wrong) {
2643     label = -label;   
2644   }
2645
2646   pt->SetLabel(label); 
2647
2648 }
2649
2650 //_____________________________________________________________________________
2651 void AliTRDtracker::UseClusters(const AliKalmanTrack *t, Int_t from) const 
2652 {
2653   //
2654   // Use clusters, but don't abuse them!
2655   //
2656
2657   const Float_t kmaxchi2 = 18;
2658   const Float_t kmincl   = 10;
2659
2660   AliTRDtrack *track = (AliTRDtrack *) t;
2661
2662   Int_t ncl = t->GetNumberOfClusters();
2663   for (Int_t i = from; i < ncl; i++) {
2664     Int_t index  = t->GetClusterIndex(i);
2665     AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2666     Int_t ilayer = fGeom->GetLayer(c->GetDetector());
2667     if (track->GetTracklets(ilayer).GetChi2() > kmaxchi2) {
2668       continue; 
2669     }
2670     if (track->GetTracklets(ilayer).GetN()    <   kmincl) {
2671       continue; 
2672     }
2673     if (!(c->IsUsed())) {
2674       c->Use();
2675     }
2676   }
2677
2678 }
2679
2680 //_____________________________________________________________________________
2681 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
2682 {
2683   //
2684   // Parametrised "expected" error of the cluster reconstruction in Y 
2685   //
2686
2687   Double_t s = 0.08 * 0.08;    
2688   return s;
2689
2690 }
2691
2692 //_____________________________________________________________________________
2693 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
2694 {
2695   //
2696   // Parametrised "expected" error of the cluster reconstruction in Z 
2697   //
2698
2699   Double_t s = 9.0 * 9.0 / 12.0;  
2700   return s;
2701
2702 }                  
2703
2704 //_____________________________________________________________________________
2705 Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const 
2706 {
2707   //
2708   // Returns radial position which corresponds to time bin <localTB>
2709   // in tracking sector <sector> and plane <plane>
2710   //
2711
2712   Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB); 
2713   Int_t pl    = fTrSec[sector]->GetLayerNumber(index);
2714
2715   return fTrSec[sector]->GetLayer(pl)->GetX();
2716
2717 }
2718
2719
2720 //_____________________________________________________________________________
2721 AliTRDtracker::AliTRDtrackingSector
2722              ::AliTRDtrackingSector(AliTRDgeometry *geo, Int_t gs)
2723   :fN(0)
2724   ,fGeom(geo)
2725   ,fGeomSector(gs)
2726 {
2727   //
2728   // AliTRDtrackingSector Constructor
2729   //
2730
2731   AliTRDpadPlane *padPlane = 0;
2732   AliTRDpropagationLayer *ppl = 0;
2733
2734   // Get holes description from geometry
2735   Bool_t holes[AliTRDgeometry::kNstack];
2736   for (Int_t istack = 0; istack < AliTRDgeometry::kNstack; istack++) {
2737     holes[istack] = fGeom->IsHole(0,istack,gs);
2738   } 
2739   
2740   for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
2741     fTimeBinIndex[i] = -1;
2742   }
2743
2744   Double_t x;
2745   Double_t dx;
2746   Double_t rho;
2747   Double_t radLength;
2748
2749   // Add layers for each of the planes
2750   Double_t dxAmp   = (Double_t) fGeom->CamHght(); // Amplification region
2751   //Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region  
2752
2753   const Int_t kNstacks  = AliTRDgeometry::Nstack();
2754   Int_t    tbIndex;
2755   Double_t ymax           = 0;
2756   Double_t ymaxsensitive  = 0;
2757   Double_t *zc            = new Double_t[kNstacks];
2758   Double_t *zmax          = new Double_t[kNstacks];
2759   Double_t *zmaxsensitive = new Double_t[kNstacks];  
2760     
2761   for (Int_t layer = 0; layer < AliTRDgeometry::Nlayer(); layer++) {
2762
2763     ymax          = fGeom->GetChamberWidth(layer) / 2.0;
2764     padPlane      = fGeom->GetPadPlane(layer,0);
2765     ymaxsensitive = (padPlane->GetColSize(1) * padPlane->GetNcols() - 4.0) / 2.0;    
2766
2767     for (Int_t st = 0; st < kNstacks; st++) {
2768       zmax[st]          = fGeom->GetChamberLength(layer,st) / 2.0;
2769       Float_t pad       = padPlane->GetRowSize(1);
2770       Float_t row0      = fGeom->GetRow0(layer,st,0);
2771       Int_t   nPads     = fGeom->GetRowMax(layer,st,0);
2772       zmaxsensitive[st] = Float_t(nPads) * pad / 2.0;      
2773       zc[st]            = -(pad * nPads) / 2.0 + row0;
2774     }
2775
2776                 AliTRDcalibDB *fCalibration = AliTRDcalibDB::Instance();
2777     dx        = fCalibration->GetVdrift(0,0,0)
2778               / AliTRDCommonParam::Instance()->GetSamplingFrequency();
2779     rho       = 0.00295 * 0.85; //????
2780     radLength = 11.0;  
2781
2782     Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(layer);
2783     //Double_t xbottom = x0 - dxDrift;
2784     //Double_t xtop    = x0 + dxAmp;
2785                 
2786                 //temporary !! (A.Bercuci)
2787     Int_t t0 = (Int_t)fCalibration->GetT0Average(AliTRDgeometry::GetDetector(layer, 2, gs));
2788
2789     Int_t nTimeBins =  AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS();    
2790     for (Int_t iTime = 0; iTime < nTimeBins; iTime++) {
2791
2792       Double_t xlayer = iTime * dx - dxAmp;
2793       //if (xlayer<0) xlayer = dxAmp / 2.0;
2794       x = x0 - xlayer;
2795
2796       tbIndex = CookTimeBinIndex(layer,iTime);
2797       ppl     = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex,layer);
2798       ppl->SetYmax(ymax,ymaxsensitive);
2799       ppl->SetZ(zc,zmax,zmaxsensitive);
2800       ppl->SetHoles(holes);
2801       if(iTime == t0) ppl->SetT0();
2802                         
2803                         InsertLayer(ppl);      
2804
2805     }
2806
2807   }    
2808
2809   MapTimeBinLayers();
2810
2811   delete [] zc;
2812   delete [] zmax;
2813   delete [] zmaxsensitive;
2814
2815 }
2816
2817 //_____________________________________________________________________________
2818 AliTRDtracker::AliTRDtrackingSector
2819              ::AliTRDtrackingSector(const AliTRDtrackingSector &/*t*/)
2820   :fN(0)
2821   ,fGeom(0)
2822   ,fGeomSector(0)
2823 {
2824   //
2825   // Copy constructor
2826   //
2827
2828 }
2829
2830 //_____________________________________________________________________________
2831 AliTRDtracker::AliTRDtrackingSector
2832              ::~AliTRDtrackingSector()
2833 {
2834   //
2835   // Destructor
2836   //
2837
2838   for (Int_t i = 0; i < fN; i++) {
2839     delete fLayers[i]; 
2840   }
2841
2842 }
2843
2844 //_____________________________________________________________________________
2845 Int_t  AliTRDtracker::AliTRDtrackingSector
2846                     ::CookTimeBinIndex(Int_t plane, Int_t localTB) const
2847 {
2848   //
2849   // Depending on the digitization parameters calculates global
2850   // (i.e. for a whole TRD stack of 6 planes) time bin index for 
2851   // timebin <localTB> in plane <plane>
2852   //
2853
2854   Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS();
2855   Int_t gtb        = (plane+1) * tbPerPlane - localTB - 1;
2856   if (localTB < 0) {
2857     return -1;
2858   }
2859   if (gtb     < 0) {
2860     return -1;
2861   }
2862
2863   return gtb;
2864
2865 }
2866
2867 //_____________________________________________________________________________
2868 void AliTRDtracker::AliTRDtrackingSector
2869                   ::MapTimeBinLayers() 
2870 {
2871   //
2872   // For all sensitive time bins sets corresponding layer index
2873   // in the array fTimeBins 
2874   //
2875
2876   Int_t index;
2877
2878   for (Int_t i = 0; i < fN; i++) {
2879
2880     index = fLayers[i]->GetTimeBinIndex();
2881     
2882     if (index < 0) {
2883       continue;
2884     }
2885     if (index >= (Int_t) kMaxTimeBinIndex) {
2886       //AliWarning(Form("Index %d exceeds allowed maximum of %d!\n"
2887       //               ,index,kMaxTimeBinIndex-1));
2888       continue;
2889     }
2890
2891     fTimeBinIndex[index] = i;
2892
2893   }
2894
2895 }
2896
2897 //_____________________________________________________________________________
2898 Int_t AliTRDtracker::AliTRDtrackingSector
2899                    ::GetLayerNumber(Double_t x) const
2900 {
2901   // 
2902   // Returns the number of time bin which in radial position is closest to <x>
2903   //
2904
2905   if (x >= fLayers[fN-1]->GetX()) {
2906     return fN - 1; 
2907   }
2908   if (x <= fLayers[   0]->GetX()) {
2909     return 0; 
2910   }
2911
2912   Int_t b = 0;
2913   Int_t e = fN - 1;
2914   Int_t m = (b + e) / 2;
2915
2916   for ( ; b < e; m = (b + e) / 2) {
2917     if (x > fLayers[m]->GetX()) {
2918       b = m + 1;
2919     }
2920     else {
2921       e = m;
2922     }
2923   }
2924
2925   if (TMath::Abs(x - fLayers[m]->GetX()) > TMath::Abs(x - fLayers[m+1]->GetX())) {
2926     return m + 1;
2927   }
2928   else {
2929     return m;
2930   }
2931
2932 }
2933
2934 //_____________________________________________________________________________
2935 Int_t AliTRDtracker::AliTRDtrackingSector
2936                    ::GetInnerTimeBin() const 
2937 {
2938   // 
2939   // Returns number of the innermost SENSITIVE propagation layer
2940   //
2941
2942   return GetLayerNumber(0);
2943
2944 }
2945
2946 //_____________________________________________________________________________
2947 Int_t AliTRDtracker::AliTRDtrackingSector
2948                    ::GetOuterTimeBin() const 
2949 {
2950   // 
2951   // Returns number of the outermost SENSITIVE time bin
2952   //
2953
2954   return GetLayerNumber(GetNumberOfTimeBins() - 1);
2955
2956 }
2957
2958 //_____________________________________________________________________________
2959 Int_t AliTRDtracker::AliTRDtrackingSector
2960                    ::GetNumberOfTimeBins() const 
2961 {
2962   // 
2963   // Returns number of SENSITIVE time bins
2964   //
2965
2966   Int_t tb;
2967   Int_t layer;
2968
2969   for (tb = kMaxTimeBinIndex - 1; tb >= 0; tb--) {
2970     layer = GetLayerNumber(tb);
2971     if (layer >= 0) {
2972       break;
2973     }
2974   }
2975
2976   return tb + 1;
2977
2978 }
2979
2980 //_____________________________________________________________________________
2981 void AliTRDtracker::AliTRDtrackingSector
2982                   ::InsertLayer(AliTRDpropagationLayer *pl)
2983
2984   //
2985   // Insert layer <pl> in fLayers array.
2986   // Layers are sorted according to X coordinate.
2987   //
2988
2989   if (fN == ((Int_t) kMaxLayersPerSector)) {
2990     //AliWarning("Too many layers !\n");
2991     return;
2992   }
2993
2994   if (fN == 0) {
2995     fLayers[fN++] = pl; 
2996     return;
2997   }
2998
2999   Int_t i = Find(pl->GetX());
3000
3001   memmove(fLayers+i+1,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
3002
3003   fLayers[i] = pl; 
3004   fN++;
3005
3006 }              
3007
3008 //_____________________________________________________________________________
3009 Int_t AliTRDtracker::AliTRDtrackingSector
3010                    ::Find(Double_t x) const 
3011 {
3012   //
3013   // Returns index of the propagation layer nearest to X 
3014   //
3015
3016   if (x <= fLayers[0]->GetX()) {
3017     return 0;
3018   }
3019
3020   if (x > fLayers[fN-1]->GetX()) {
3021     return fN;
3022   }
3023
3024   Int_t b = 0;
3025   Int_t e = fN-1;
3026   Int_t m = (b + e) / 2;
3027
3028   for (; b < e; m = (b + e) / 2) {
3029     if (x > fLayers[m]->GetX()) {
3030       b = m + 1;
3031     }
3032     else {
3033       e = m;
3034     }
3035   }
3036
3037   return m;
3038
3039 }             
3040
3041 //_____________________________________________________________________________
3042 Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c) 
3043 {
3044   //
3045   // Returns correction factor for tilted pads geometry 
3046   //
3047
3048   Int_t    det   = c->GetDetector();    
3049   Int_t    layer = fGeom->GetLayer(det);
3050   AliTRDpadPlane *padPlane = fGeom->GetPadPlane(layer,0);
3051   Double_t h01   = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
3052
3053   if (fNoTilt) {
3054     h01 = 0;
3055   }
3056
3057   return h01;
3058
3059 }
3060
3061 //_____________________________________________________________________________
3062 Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
3063                                 , AliTRDtrack *track
3064                                 , Int_t *clusters
3065                                 , AliTRDtracklet &tracklet)
3066 {
3067   //
3068   // Try to find the nearest clusters to the track in the time bins 
3069   // between <t0> and <t1>. 
3070   // Also the corresponding tracklet is calculated
3071   // Correction coeficients - depend on TRD parameters - to be changed accordingly
3072   //
3073
3074   const Int_t    kN1 = 100;
3075   const Int_t    kN2 = 10;
3076
3077   Double_t       x[kN1];
3078   Double_t       yt[kN1];
3079   Double_t       zt[kN1];
3080   Float_t        zmean[kN1];
3081   Float_t        nmean[kN1];
3082
3083   Double_t       dz[kN2][kN1];
3084   Double_t       dy[kN2][kN1];
3085   Int_t          indexes[kN2][kN1];  // Indexes of the clusters in the road
3086   Int_t          best[kN2][kN1];     // Index of best matching cluster 
3087   AliTRDcluster *cl[kN2][kN1];       // Pointers to the clusters in the road
3088
3089   Double_t xmean   = 0.0;            // Reference x
3090   Int_t    clfound = 0;
3091
3092   // Initialize the arrays
3093   for (Int_t it = 0; it < kN1; it++) {
3094
3095     x[it]        =  0.0;
3096     yt[it]       =  0.0;
3097     zt[it]       =  0.0;
3098     clusters[it] = -2;
3099     zmean[it]    =  0.0;
3100     nmean[it]    =  0.0;
3101
3102     for (Int_t ih = 0; ih < kN2; ih++) {
3103       indexes[ih][it] = -2;
3104       cl[ih][it]      =  0;
3105       dz[ih][it]      = -100.0;
3106       dy[ih][it]      = -100.0;
3107       best[ih][it]    =  0;
3108     }
3109
3110   }  
3111
3112   Double_t x0        = track->GetX();
3113   Double_t sigmaz    = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
3114   Int_t    nall      = 0;
3115   Int_t    nfound    = 0;
3116   Double_t h01       = 0.0;
3117   Int_t    layer     = -1;
3118   Int_t    detector  = -1;
3119   Float_t  padlength = 0.0;
3120
3121   AliTRDtrack track2(* track);
3122   Float_t  snpy      = track->GetSnp();
3123   Float_t  tany      = TMath::Sqrt(snpy*snpy / ((1.-snpy)*(1.+snpy))); 
3124   if (snpy < 0.0) {
3125     tany *= -1.0;
3126   }
3127
3128   Double_t sy2       = ExpectedSigmaY2(x0,track->GetTgl(),track->GetSignedPt());
3129   Double_t sz2       = ExpectedSigmaZ2(x0,track->GetTgl());
3130   Double_t road      = 15.0 * TMath::Sqrt(track->GetSigmaY2() + sy2);
3131   if (road > 6.0) {
3132     road = 6.0;
3133   }
3134
3135   for (Int_t it = 0; it < t1-t0; it++) {
3136
3137     Double_t maxChi2[2] = { fgkMaxChi2, fgkMaxChi2 };      
3138     AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(it+t0));
3139     if (timeBin == 0) {
3140       continue; // No indexes1
3141     }
3142
3143     Int_t    maxn = timeBin;
3144     x[it]  = timeBin.GetX();
3145     track2.PropagateTo(x[it]);
3146     yt[it] = track2.GetY();
3147     zt[it] = track2.GetZ();
3148     
3149     Double_t y    = yt[it];
3150     Double_t z    = zt[it];
3151     Double_t chi2 = 1000000.0;
3152     nall++;
3153
3154     //
3155     // Find 2 nearest cluster at given time bin
3156     // 
3157     Int_t    checkPoint[4] = { 0, 0, 0, 0 };
3158     Double_t minY          = 123456789.0;
3159     Double_t minD[2]       = { 1.0, 1.0 };
3160
3161     for (Int_t i = timeBin.Find(y - road); i < maxn; i++) {
3162
3163       AliTRDcluster *c = (AliTRDcluster *) (timeBin[i]);
3164       h01 = GetTiltFactor(c);
3165       if (layer < 0) {
3166         Int_t det = c->GetDetector();
3167         layer     = fGeom->GetLayer(det);
3168         padlength = TMath::Sqrt(c->GetSigmaZ2() * 12.0);
3169       }
3170
3171       if (c->GetY() > (y + road)) {
3172         break;
3173       }
3174
3175       fHDeltaX->Fill(c->GetX() - x[it]);
3176
3177       if (TMath::Abs(c->GetY()-y) < TMath::Abs(minY)) {
3178         minY    = c->GetY() - y;
3179         minD[0] = c->GetY() - y;
3180         minD[1] = c->GetZ() - z;
3181       }
3182
3183       checkPoint[0]++;
3184
3185       fHMinZ->Fill(c->GetZ() - z);
3186       if ((c->GetZ() - z) * (c->GetZ() - z) > 2 * (12.0 * sz2)) {
3187         continue;
3188       }
3189       checkPoint[1]++;
3190
3191       Double_t dist = TMath::Abs(c->GetZ() - z);
3192       if (dist > (0.5 * padlength + 6.0 * sigmaz)) {
3193         continue;   // 6 sigma boundary cut
3194       }
3195       checkPoint[2]++;
3196
3197       // Sigma boundary cost function
3198       Double_t cost = 0.0;
3199       if (dist> (0.5 * padlength - sigmaz)){ 
3200         cost =  (dist - 0.5*padlength) / (2.0 * sigmaz);
3201         if (cost > -1) {
3202           cost = (cost + 1.0) * (cost + 1.0);
3203         }
3204         else {
3205           cost = 0.0;
3206         }
3207       }
3208       chi2 = track2.GetPredictedChi2(c,h01) + cost;
3209       clfound++;      
3210
3211       if (chi2 > maxChi2[1]) {
3212         continue;
3213       }
3214       checkPoint[3]++;
3215
3216       // Store the clusters in the road
3217       detector = c->GetDetector();
3218       for (Int_t ih = 2; ih < 9; ih++) {  
3219         if (cl[ih][it] == 0) {
3220           cl[ih][it]      = c;
3221           indexes[ih][it] = timeBin.GetIndex(i); // Index - 9 - reserved for outliers
3222           break;
3223         }
3224       }
3225
3226       if (chi2 < maxChi2[0]) {
3227         maxChi2[1]     = maxChi2[0];
3228         maxChi2[0]     = chi2;
3229         indexes[1][it] = indexes[0][it];
3230         cl[1][it]      = cl[0][it];
3231         indexes[0][it] = timeBin.GetIndex(i);
3232         cl[0][it]      = c;
3233         continue;
3234       }
3235       maxChi2[1]     = chi2;
3236       cl[1][it]      = c;
3237       indexes[1][it] = timeBin.GetIndex(i); 
3238
3239     }
3240     
3241     for(int iCheckPoint = 0; iCheckPoint<4; iCheckPoint++) {
3242       fHFindCl[iCheckPoint]->Fill(checkPoint[iCheckPoint]);
3243     }
3244
3245     if (checkPoint[3]) {
3246       if (track->GetSignedPt() > 0) {
3247         fHMinYPos->Fill(minY);
3248       }
3249       else {
3250         fHMinYNeg->Fill(minY);
3251       }
3252       fHMinD->Fill(minD[0],minD[1]);
3253     }
3254
3255     if (cl[0][it]) {
3256       nfound++;
3257       xmean += x[it];
3258     }
3259
3260   }
3261
3262   if (nfound < 4) {
3263     return 0;  
3264   }
3265   xmean /= Float_t(nfound);  // Middle x
3266   track2.PropagateTo(xmean); // Propagate track to the center
3267
3268   //
3269   // Choose one of the variants
3270   //
3271   Float_t  sumz   = 0.0;
3272   Float_t  sum    = 0.0;
3273   Double_t sumdy  = 0.0;
3274   Double_t sumdy2 = 0.0;
3275   Double_t sumx   = 0.0;
3276   Double_t sumxy  = 0.0;
3277   Double_t sumx2  = 0.0;
3278   Double_t mpads  = 0.0;
3279
3280   Int_t    changes[kN2];
3281
3282   Int_t    ngood[kN2];
3283   Int_t    nbad[kN2];
3284
3285   Double_t meanz[kN2];
3286   Double_t moffset[kN2];    // Mean offset
3287   Double_t mean[kN2];       // Mean value
3288   Double_t angle[kN2];      // Angle
3289
3290   Double_t smoffset[kN2];   // Sigma of mean offset
3291   Double_t smean[kN2];      // Sigma of mean value
3292   Double_t sangle[kN2];     // Sigma of angle
3293   Double_t smeanangle[kN2]; // Correlation
3294
3295   Double_t sigmas[kN2];     
3296   Double_t tchi2s[kN2];     // Chi2s for tracklet
3297
3298   for (Int_t it = 0; it < kN2; it++) {
3299
3300     ngood[it]      = 0;
3301     nbad[it]       = 0;
3302
3303     meanz[it]      = 0.0;
3304     moffset[it]    = 0.0;   // Mean offset
3305     mean[it]       = 0.0;   // Mean value
3306     angle[it]      = 0.0;   // Angle
3307
3308     smoffset[it]   = 1.0e5; // Sigma of mean offset
3309     smean[it]      = 1.0e5; // Sigma of mean value
3310     sangle[it]     = 1.0e5; // Sigma of angle
3311     smeanangle[it] = 0.0;   // Correlation
3312
3313     sigmas[it]     = 1.0e5;     
3314     tchi2s[it]     = 1.0e5; // Chi2s for tracklet
3315
3316   }
3317
3318   //
3319   // Calculate zmean
3320