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