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