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