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