bd1ab8385355a944f055c43c54fc3cf129dd1a4e
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackerV1.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 /* $Id$ */
18
19 ///////////////////////////////////////////////////////////////////////////////
20 //                                                                           //
21 //  Track finder                                                             //
22 //                                                                           //
23 //  Authors:                                                                 //
24 //    Alex Bercuci <A.Bercuci@gsi.de>                                        //
25 //    Markus Fasel <M.Fasel@gsi.de>                                          //
26 //                                                                           //
27 ///////////////////////////////////////////////////////////////////////////////
28
29 #include <Riostream.h>
30 #include <stdio.h>
31 #include <string.h>
32
33 #include <TBranch.h>
34 #include <TFile.h>
35 #include <TGraph.h>
36 #include <TH1D.h>
37 #include <TH2D.h>
38 #include <TLinearFitter.h>
39 #include <TObjArray.h> 
40 #include <TROOT.h>
41 #include <TTree.h>  
42 #include <TClonesArray.h>
43 #include <TRandom.h>
44 #include <TTreeStream.h>
45
46 #include "AliLog.h"
47 #include "AliESDEvent.h"
48 #include "AliAlignObj.h"
49 #include "AliRieman.h"
50 #include "AliTrackPointArray.h"
51
52 #include "AliTRDtracker.h"
53 #include "AliTRDtrackerV1.h"
54 #include "AliTRDgeometry.h"
55 #include "AliTRDpadPlane.h"
56 #include "AliTRDgeometry.h"
57 #include "AliTRDcluster.h" 
58 #include "AliTRDtrack.h"
59 #include "AliTRDseed.h"
60 #include "AliTRDcalibDB.h"
61 #include "AliTRDCommonParam.h"
62 #include "AliTRDReconstructor.h"
63 #include "AliTRDCalibraFillHisto.h"
64 #include "AliTRDtrackerFitter.h"
65 #include "AliTRDstackLayer.h"
66 #include "AliTRDrecoParam.h"
67 #include "AliTRDseedV1.h"
68 #include "AliTRDtrackV1.h"
69 #include "Cal/AliTRDCalDet.h"
70
71 #define DEBUG
72
73 ClassImp(AliTRDtrackerV1)
74 Double_t AliTRDtrackerV1::fgTopologicQA[kNConfigs] = {
75                 0.1112, 0.1112, 0.1112, 0.0786, 0.0786,
76                 0.0786, 0.0786, 0.0579, 0.0579, 0.0474,
77                 0.0474, 0.0408, 0.0335, 0.0335, 0.0335
78 };
79
80 //____________________________________________________________________
81 AliTRDtrackerV1::AliTRDtrackerV1(AliTRDrecoParam *p) 
82   :AliTRDtracker()
83   ,fSieveSeeding(0)
84   ,fTracklets(0x0)
85   ,fRecoParam(p)
86   ,fFitter(0x0)
87 {
88   //
89   // Default constructor. Nothing is initialized.
90   //
91
92 }
93
94 //____________________________________________________________________
95 AliTRDtrackerV1::AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p) 
96   :AliTRDtracker(in)
97   ,fSieveSeeding(0)
98   ,fTracklets(0x0)
99   ,fRecoParam(p)
100   ,fFitter(0x0)
101 {
102   //
103   // Standard constructor.
104   // Setting of the geometry file, debug output (if enabled)
105   // and initilize fitter helper.
106   //
107
108         fFitter = new AliTRDtrackerFitter();
109
110 #ifdef DEBUG
111         fFitter->SetDebugStream(fDebugStreamer);
112 #endif
113
114 }
115   
116 //____________________________________________________________________
117 AliTRDtrackerV1::~AliTRDtrackerV1()
118
119   //
120   // Destructor
121   //
122
123         if(fFitter) delete fFitter;
124         if(fRecoParam) delete fRecoParam;
125         if(fTracklets) {fTracklets->Delete(); delete fTracklets;}
126 }
127
128 //____________________________________________________________________
129 Int_t AliTRDtrackerV1::Clusters2Tracks(AliESDEvent *esd)
130 {
131   //
132   // Steering stand alone tracking for full TRD detector
133   //
134   // Parameters :
135   //   esd     : The ESD event. On output it contains 
136   //             the ESD tracks found in TRD.
137   //
138   // Output :
139   //   Number of tracks found in the TRD detector.
140   // 
141   // Detailed description
142   // 1. Launch individual SM trackers. 
143   //    See AliTRDtrackerV1::Clusters2TracksSM() for details.
144   //
145
146         if(!fRecoParam){
147                 AliError("Reconstruction configuration not initialized. Call first AliTRDtrackerV1::SetRecoParam().");
148                 return 0;
149         }
150         
151         //AliInfo("Start Track Finder ...");
152         Int_t ntracks = 0;
153         for(int ism=0; ism<AliTRDtracker::kTrackingSectors; ism++){
154                         //AliInfo(Form("Processing supermodule %i ...", ism));
155                         ntracks += Clusters2TracksSM(fTrSec[ism], esd);
156         }
157         AliInfo(Form("Found %d TRD tracks.", ntracks));
158         return ntracks;
159 }
160
161
162 //_____________________________________________________________________________
163 Bool_t AliTRDtrackerV1::GetTrackPoint(Int_t /*index*/, AliTrackPoint &/*p*/) const
164 {
165         //AliInfo(Form("Asking for tracklet %d", index));
166         
167         if(index<0) return kFALSE;
168         //AliTRDseedV1 *tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index);
169         // etc
170         return kTRUE;
171 }
172
173
174 //_____________________________________________________________________________
175 Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) 
176 {
177   //
178   // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
179   // backpropagated by the TPC tracker. Each seed is first propagated 
180   // to the TRD, and then its prolongation is searched in the TRD.
181   // If sufficiently long continuation of the track is found in the TRD
182   // the track is updated, otherwise it's stored as originaly defined 
183   // by the TPC tracker.   
184   //  
185
186         Int_t   found    = 0;     // number of tracks found
187         Float_t foundMin = 20.0;
188         
189         AliTRDseed::SetNTimeBins(fTimeBinsPerPlane);
190         Int_t    nSeed   = event->GetNumberOfTracks();
191         if(!nSeed){
192                 // run stand alone tracking
193                 if (AliTRDReconstructor::SeedingOn()) Clusters2Tracks(event);
194                 return 0;
195         }
196         
197         Float_t *quality = new Float_t[nSeed];
198         Int_t   *index   = new Int_t[nSeed];
199         for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
200                 AliESDtrack *seed = event->GetTrack(iSeed);
201                 Double_t covariance[15];
202                 seed->GetExternalCovariance(covariance);
203                 quality[iSeed] = covariance[0] + covariance[2];
204         }
205         // Sort tracks according to covariance of local Y and Z
206         TMath::Sort(nSeed,quality,index,kFALSE);
207         
208         // Backpropagate all seeds
209         for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
210         
211                 // Get the seeds in sorted sequence
212                 AliESDtrack *seed = event->GetTrack(index[iSeed]);
213         
214                 // Check the seed status
215                 ULong_t status = seed->GetStatus();
216                 if ((status & AliESDtrack::kTPCout) == 0) continue;
217                 if ((status & AliESDtrack::kTRDout) != 0) continue;
218         
219                 // Do the back prolongation
220                 Int_t   lbl         = seed->GetLabel();
221                 AliTRDtrackV1 *track  = new AliTRDtrackV1(*seed);
222                 //track->Print();
223                 track->SetSeedLabel(lbl);
224                 seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); // Make backup
225                 fNseeds++;
226                 Float_t p4          = track->GetC();
227                 Int_t   expectedClr = FollowBackProlongation(*track);
228                 //AliInfo(Form("\nTRACK %d Clusters %d [%d] in chi2 %f", index[iSeed], expectedClr, track->GetNumberOfClusters(), track->GetChi2()));
229                 //track->Print();
230
231                 //Double_t cov[15];
232                 //seed->GetExternalCovariance(cov);
233                 //AliInfo(Form("track %d cov[%f %f] 0", index[iSeed], cov[0], cov[2]));
234
235                 if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
236                                 (track->Pt() > 0.8)) {
237                         //
238                         // Make backup for back propagation
239                         //
240                         Int_t foundClr = track->GetNumberOfClusters();
241                         if (foundClr >= foundMin) {
242                                 //AliInfo(Form("Making backup track ncls [%d]...", foundClr));
243                                 track->CookdEdx();
244                                 track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
245                                 CookLabel(track,1 - fgkLabelFraction);
246                                 if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
247                                 
248                                 
249                 //seed->GetExternalCovariance(cov);
250                 //AliInfo(Form("track %d cov[%f %f] 0 test", index[iSeed], cov[0], cov[2]));
251
252                                 // Sign only gold tracks
253                                 if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
254                                         if ((seed->GetKinkIndex(0)      ==   0) &&
255                                                         (track->Pt()                <  1.5)) UseClusters(track);
256                                 }
257                                 Bool_t isGold = kFALSE;
258         
259                                 // Full gold track
260                                 if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
261                                         if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
262
263                                         isGold = kTRUE;
264                                 }
265                 //seed->GetExternalCovariance(cov);
266                 //AliInfo(Form("track %d cov[%f %f] 00", index[iSeed], cov[0], cov[2]));
267         
268                                 // Almost gold track
269                                 if ((!isGold)  && (track->GetNCross() == 0) &&
270                                                 (track->GetChi2() / track->GetNumberOfClusters()  < 7)) {
271                                         //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
272                                         if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
273                                         
274                                         isGold = kTRUE;
275                                 }
276                 //seed->GetExternalCovariance(cov);
277                 //AliInfo(Form("track %d cov[%f %f] 01", index[iSeed], cov[0], cov[2]));
278                                 
279                                 if ((!isGold) && (track->GetBackupTrack())) {
280                                         if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) && ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
281                                                 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
282                                                 isGold = kTRUE;
283                                         }
284                                 }
285                 //seed->GetExternalCovariance(cov);
286                 //AliInfo(Form("track %d cov[%f %f] 02", index[iSeed], cov[0], cov[2]));
287         
288                                 //if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected())  > 0.4)) {
289                                         //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
290                                 //}
291                         }
292                 }
293                 /**/
294         
295                 /**/
296                 // Debug part of tracking
297 /*              TTreeSRedirector &cstream = *fDebugStreamer;
298                 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.
299                 if (AliTRDReconstructor::StreamLevel() > 0) {
300                         if (track->GetBackupTrack()) {
301                                 cstream << "Tracks"
302                                 << "EventNrInFile="  << eventNrInFile
303                                 << "ESD.="     << seed
304                                 << "trd.="     << track
305                                 << "trdback.=" << track->GetBackupTrack()
306                                 << "\n";
307                         }
308                         else {
309                                 cstream << "Tracks"
310                                 << "EventNrInFile="  << eventNrInFile
311                                 << "ESD.="     << seed
312                                 << "trd.="     << track
313                                 << "trdback.=" << track
314                                 << "\n";
315                         }
316                 }*/
317                 /**/
318         
319                 //seed->GetExternalCovariance(cov);
320                 //AliInfo(Form("track %d cov[%f %f] 1", index[iSeed], cov[0], cov[2]));
321
322                 // Propagation to the TOF (I.Belikov)
323                 if (track->GetStop() == kFALSE) {
324                         //AliInfo("Track not stopped in TRD ...");
325                         Double_t xtof  = 371.0;
326                         Double_t xTOF0 = 370.0;
327                 
328                         Double_t c2    = track->GetSnp() + track->GetC() * (xtof - track->GetX());
329                         if (TMath::Abs(c2) >= 0.99) {
330                                 delete track;
331                                 continue;
332                         }
333                         
334                         PropagateToX(*track,xTOF0,fgkMaxStep);
335         
336                         // Energy losses taken to the account - check one more time
337                         c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
338                         if (TMath::Abs(c2) >= 0.99) {
339                                 delete track;
340                                 continue;
341                         }
342                         
343                         //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
344                         //      fHBackfit->Fill(7);
345                         //delete track;
346                         //      continue;
347                         //}
348         
349                         Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
350                         Double_t y;
351                         track->GetYAt(xtof,GetBz(),y);
352                         if (y >  ymax) {
353                                 if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
354                                         delete track;
355                                         continue;
356                                 }
357                         }else if (y < -ymax) {
358                                 if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
359                                         delete track;
360                                         continue;
361                                 }
362                         }
363                                         
364                         if (track->PropagateTo(xtof)) {
365                                 //AliInfo("set kTRDout");
366                                 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
367         
368                                 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
369                                         for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
370                                                 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
371                                         }
372                                         seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
373                                 }
374                                 //seed->SetTRDtrack(new AliTRDtrack(*track));
375                                 if (track->GetNumberOfClusters() > foundMin) found++;
376                         }
377                 } else {
378                         //AliInfo("Track stopped in TRD ...");
379                         
380                         if ((track->GetNumberOfClusters() >              15) &&
381                                         (track->GetNumberOfClusters() > 0.5*expectedClr)) {
382                                 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
383         
384                                 //seed->SetStatus(AliESDtrack::kTRDStop);
385                                 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
386                                         for (Int_t j = 0; j <AliESDtrack::kNSlice; j++) {
387                                                 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
388                                         }
389                                         seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
390                                 }
391                                 //seed->SetTRDtrack(new AliTRDtrack(*track));
392                                 found++;
393                         }
394                 }
395
396                 //if (((t->GetStatus()&AliESDtrack::kTRDout)!=0 )
397         
398                 seed->SetTRDQuality(track->StatusForTOF());
399                 seed->SetTRDBudget(track->GetBudget(0));
400                 
401 //              if ((seed->GetStatus()&AliESDtrack::kTRDin)!=0 ) printf("TRDin ");      
402 //              if ((seed->GetStatus()&AliESDtrack::kTRDbackup)!=0 ) printf("TRDbackup ");      
403 //              if ((seed->GetStatus()&AliESDtrack::kTRDStop)!=0 ) printf("TRDstop ");  
404 //              if ((seed->GetStatus()&AliESDtrack::kTRDout)!=0 ) printf("TRDout ");    
405 //              printf("\n");
406                 delete track;
407
408                 //seed->GetExternalCovariance(cov);
409                 //AliInfo(Form("track %d cov[%f %f] 2", index[iSeed], cov[0], cov[2]));
410         }
411         
412
413         AliInfo(Form("Number of seeds: %d",fNseeds));
414         AliInfo(Form("Number of back propagated TRD tracks: %d",found));
415                 
416         //fSeeds->Clear(); 
417         fNseeds = 0;
418         
419         delete [] index;
420         delete [] quality;
421         
422   return 0;
423 }
424
425
426 //____________________________________________________________________
427 Int_t AliTRDtrackerV1::RefitInward(AliESDEvent *event)
428 {
429   //
430   // Refits tracks within the TRD. The ESD event is expected to contain seeds 
431   // at the outer part of the TRD. 
432   // The tracks are propagated to the innermost time bin 
433   // of the TRD and the ESD event is updated
434   // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
435   //
436
437   Int_t   nseed    = 0; // contor for loaded seeds
438   Int_t   found    = 0; // contor for updated TRD tracks
439   
440   // Calibration monitor
441   AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
442   if (!calibra) AliInfo("Could not get Calibra instance\n");
443   
444   
445   AliTRDtrackV1 track;
446   for (Int_t itrack = 0; itrack < event->GetNumberOfTracks(); itrack++) {
447     AliESDtrack *seed = event->GetTrack(itrack);
448                 new(&track) AliTRDtrackV1(*seed);
449
450     if (track.GetX() < 270.0) {
451       seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup);
452       //AliInfo(Form("Remove for X = %7.3f [270.]\n", track.GetX()));
453                         continue;
454     }
455
456     ULong_t status = seed->GetStatus();
457     if((status & AliESDtrack::kTRDout) == 0) continue;
458     if((status & AliESDtrack::kTRDin)  != 0) continue;
459     nseed++; 
460
461     track.ResetCovariance(50.0);
462
463                 // do the propagation and processing
464     Bool_t kUPDATE = kFALSE;
465                 Double_t xTPC = 250.0;
466     if(FollowProlongation(track)){
467                         // computes PID for track
468                         track.CookPID();
469                         // update calibration references using this track
470                         if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(&track);
471         
472                         // Prolongate to TPC
473                         if (PropagateToX(track, xTPC, fgkMaxStep)) { //  -with update
474                                 seed->UpdateTrackParams(&track, AliESDtrack::kTRDrefit);
475                                 track.UpdateESDtrack(seed);
476                                 // Add TRD track to ESDfriendTrack
477                                 if (AliTRDReconstructor::StreamLevel() > 0 /*&& quality TODO*/){ 
478                                         AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
479                                         calibTrack->SetOwner();
480                                         seed->AddCalibObject(calibTrack);
481                                 }
482                                 found++;
483                                 kUPDATE = kTRUE;
484                         }
485                 }        
486                 
487                 // Prolongate to TPC without update
488                 if(!kUPDATE) {
489       AliTRDtrackV1 tt(*seed);
490       if (PropagateToX(tt, xTPC, fgkMaxStep)) seed->UpdateTrackParams(&tt, AliESDtrack::kTRDrefit);
491     }
492   }
493   AliInfo(Form("Number of loaded seeds: %d",nseed));
494   AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
495   
496         return 0;
497 }
498
499
500 //____________________________________________________________________
501 Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
502 {
503 // Extrapolates the TRD track in the TPC direction.
504 //
505 // Parameters
506 //   t : the TRD track which has to be extrapolated
507 // 
508 // Output
509 //   number of clusters attached to the track
510 //
511 // Detailed description
512 //
513 // Starting from current radial position of track <t> this function
514 // extrapolates the track through the 6 TRD layers. The following steps
515 // are being performed for each plane:
516 // 1. prepare track:
517 //   a. get plane limits in the local x direction
518 //   b. check crossing sectors 
519 //   c. check track inclination
520 // 2. search tracklet in the tracker list (see GetTracklet() for details)
521 // 3. evaluate material budget using the geo manager
522 // 4. propagate and update track using the tracklet information.
523 //
524 // Debug level 2
525 //
526   
527         //AliInfo("");
528         Int_t    nClustersExpected = 0;
529         Int_t lastplane = 5; //GetLastPlane(&t);
530         for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
531                 //AliInfo(Form("plane %d", iplane));
532     Int_t row1 = GetGlobalTimeBin(0, iplane, 0); // to be modified to the true time bin in the geometrical acceptance
533                 //AliInfo(Form("row1 %d", row1));
534
535     // Propagate track close to the plane if neccessary
536                 AliTRDpropagationLayer *layer = fTrSec[0]->GetLayer(row1);
537                 Double_t currentx  = layer->GetX();
538     if (currentx < (-fgkMaxStep + t.GetX())) 
539       if (!PropagateToX(t, currentx+fgkMaxStep, fgkMaxStep)) break;
540
541     if (!AdjustSector(&t)) break;
542      
543                 Int_t row0    = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
544                 //AliInfo(Form("row0 %d", row0));
545
546     // Start global position
547     Double_t xyz0[3];
548     t.GetXYZ(xyz0);
549
550                 // End global position
551     Double_t x = fTrSec[0]->GetLayer(row0)->GetX(), y, z;
552     if (!t.GetProlongation(x,y,z)) break;    
553     Double_t xyz1[3];
554     xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
555     xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
556     xyz1[2] =  z;
557                 
558     // Get material budget
559     Double_t param[7];
560     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
561     Double_t xrho= param[0]*param[4];
562     Double_t xx0 = param[1]; // Get mean propagation parameters
563
564     // Propagate and update
565     //Int_t sector = t.GetSector();
566     Int_t   index   = 0;
567                 //AliInfo(Form("sector %d", sector));
568     AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
569           //AliInfo(Form("tracklet %p @ %d", tracklet, index));
570                 if(!tracklet) continue;
571                 //AliInfo(Form("reco %p", tracklet->GetRecoParam()));
572                 t.SetTracklet(tracklet, iplane, index);
573                 
574                 t.PropagateTo(tracklet->GetX0(), xx0, xrho); // not correct
575           if (!AdjustSector(&t)) break;
576           
577     Double_t maxChi2 = t.GetPredictedChi2(tracklet);
578           if (maxChi2 < 1e+10 && t.Update(tracklet, maxChi2)){ 
579                 nClustersExpected += tracklet->GetN();
580         }
581   }
582
583 #ifdef DEBUG
584         if(AliTRDReconstructor::StreamLevel() > 1){
585                 Int_t index;
586                 for(int iplane=0; iplane<6; iplane++){
587                         AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
588                         if(!tracklet) continue;
589                         t.SetTracklet(tracklet, iplane, index);
590                 }
591
592                 TTreeSRedirector &cstreamer = *fDebugStreamer;
593                 cstreamer << "FollowProlongation"
594                         << "ncl="      << nClustersExpected
595                         << "track.="   << &t
596                         << "\n";
597         }
598 #endif
599
600   return nClustersExpected;
601
602 }
603
604 //_____________________________________________________________________________
605 Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
606 {
607 // Extrapolates the TRD track in the TOF direction.
608 //
609 // Parameters
610 //   t : the TRD track which has to be extrapolated
611 // 
612 // Output
613 //   number of clusters attached to the track
614 //
615 // Detailed description
616 //
617 // Starting from current radial position of track <t> this function
618 // extrapolates the track through the 6 TRD layers. The following steps
619 // are being performed for each plane:
620 // 1. prepare track:
621 //   a. get plane limits in the local x direction
622 //   b. check crossing sectors 
623 //   c. check track inclination
624 // 2. build tracklet (see AliTRDseed::AttachClusters() for details)
625 // 3. evaluate material budget using the geo manager
626 // 4. propagate and update track using the tracklet information.
627 //
628 // Debug level 2
629 //
630
631         Int_t nClustersExpected = 0;
632
633   // Loop through the TRD planes
634   for (Int_t iplane = 0; iplane < AliTRDgeometry::Nplan(); iplane++) {
635                 //AliInfo(Form("Processing plane %d ...", iplane));
636                 // Get the global time bin for the first local time bin of the given plane
637     Int_t row0 = GetGlobalTimeBin(0, iplane, fTimeBinsPerPlane-1);
638
639                 // Retrive first propagation layer in the chamber
640                 AliTRDpropagationLayer *layer = fTrSec[0]->GetLayer(row0);
641
642                 // Get the X coordinates of the propagation layer for the first time bin
643     Double_t currentx = layer->GetX(); // what if X is not defined ???
644     if (currentx < t.GetX()) continue;
645     
646                 // Get the global time bin for the last local time bin of the given plane
647                 Int_t row1 = GetGlobalTimeBin(0, iplane, 0);
648     
649                 // Propagate closer to the current chamber if neccessary 
650     if (currentx > (fgkMaxStep + t.GetX()) && !PropagateToX(t, currentx-fgkMaxStep, fgkMaxStep)) break;
651    
652                 // Rotate track to adjacent sector if neccessary
653     if (!AdjustSector(&t)) break;
654                 Int_t sector =  Int_t(TMath::Abs(t.GetAlpha()/AliTRDgeometry::GetAlpha()));
655                 if(t.GetAlpha() < 0) sector = AliTRDgeometry::Nsect() - sector-1;
656                 
657                 //AliInfo(Form("sector %d [%f]", sector, t.GetAlpha()));
658
659                 // Check whether azimuthal angle is getting too large
660     if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) break;
661    
662     //Calculate global entry and exit positions of the track in chamber (only track prolongation)
663     Double_t xyz0[3]; // entry point 
664                 t.GetXYZ(xyz0);   
665                 //printf("Entry global x[%7.3f] y[%7.3f] z[%7.3f]\n", xyz0[0], xyz0[1], xyz0[2]);
666                                 
667                 // Get local Y and Z at the X-position of the end of the chamber
668                 Double_t x0 = fTrSec[sector]->GetLayer(row1)->GetX(), y, z;
669                 if (!t.GetProlongation(x0, y, z)) break;
670                 //printf("Exit  local x[%7.3f] y[%7.3f] z[%7.3f]\n", x0, y, z);
671
672                 Double_t xyz1[3]; // exit point
673                 xyz1[0] =  x0 * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha()); 
674     xyz1[1] = +x0 * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
675     xyz1[2] =  z;
676
677                 //printf("Exit  global x[%7.3f] y[%7.3f] z[%7.3f]\n", xyz1[0], xyz1[1], xyz1[2]);
678                 // Find tracklet along the path inside the chamber
679     AliTRDseedV1 tracklet(*t.GetTracklet(iplane));
680                 // if the track is not already build (e.g. stand alone tracker) we build it now.
681                 if(!tracklet.GetN()){ // a better check has to be implemented TODO!!!!!!!
682                 
683                         //AliInfo(Form("Building tracklet for plane %d ...", iplane));
684                         // check if we are inside detection volume
685                         Int_t ichmb = fGeom->GetChamber(xyz0[2], iplane); 
686                         if(ichmb<0) ichmb = fGeom->GetChamber(xyz1[2], iplane);
687                         if(ichmb<0){
688                                 // here we should decide what to do with the track. The space between the pads in 2 chambers is 4cm+. Is it making sense to continue building the tracklet here TODO????
689                                 AliWarning(Form("Track prolongated in the interspace between TRD detectors in plane %d. Skip plane. To be fixed !", iplane));
690                                 continue;
691                         }
692         
693                         // temporary until the functionalities of AliTRDpropagationLayer and AliTRDstackLayer are merged TODO
694                         AliTRDpadPlane *pp = fGeom->GetPadPlane(iplane, ichmb);
695                         Int_t nrows = pp->GetNrows();
696                         Double_t stacklength = pp->GetRow0ROC() - pp->GetRowEndROC();/*(nrows - 2) * pp->GetLengthIPad()        + 2 * pp->GetLengthOPad() + (nrows - 1) * pp->GetRowSpacing();*/
697                         Double_t z0  = fGeom->GetRow0(iplane, ichmb, 0);
698
699                         Int_t nClustersChmb = 0;
700                         AliTRDstackLayer stackLayer[35];
701                         for(int itb=0; itb<fTimeBinsPerPlane; itb++){
702                                 const AliTRDpropagationLayer ksmLayer(*(fTrSec[sector]->GetLayer(row1 - itb)));
703                                 stackLayer[itb] = ksmLayer;
704 #ifdef DEBUG
705                                 stackLayer[itb].SetDebugStream(fDebugStreamer);
706 #endif                  
707                                 stackLayer[itb].SetRange(z0 - stacklength, stacklength);
708                                 stackLayer[itb].SetSector(sector);
709                                 stackLayer[itb].SetStackNr(ichmb);
710                                 stackLayer[itb].SetNRows(nrows);
711                                 stackLayer[itb].SetRecoParam(fRecoParam);
712                                 stackLayer[itb].BuildIndices();
713                                 nClustersChmb += stackLayer[itb].GetNClusters();
714                         }
715                         if(!nClustersChmb) continue;
716                         //AliInfo(Form("Detector p[%d] c[%d]. Building tracklet from %d clusters ... ", iplane, ichmb, nClustersChmb));
717
718                         tracklet.SetRecoParam(fRecoParam);
719                         tracklet.SetTilt(TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle()));
720                         tracklet.SetPadLength(pp->GetLengthIPad());
721                         tracklet.SetPlane(iplane);
722                         //Int_t tbRange   = fTimeBinsPerPlane; //Int_t(AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght() * AliTRDCommonParam::Instance()->GetSamplingFrequency()/AliTRDcalibDB::Instance()->GetVdriftDet()->GetValue(det));
723                         //printf("%d hl[%f] pl[%f] tb[%d]\n", il, hL[il], padlength[il], tbRange[il]);
724                         //tracklet.SetNTimeBinsRange(tbRange);
725                         tracklet.SetX0(x0);
726                         tracklet.Init(&t);
727                         if(!tracklet.AttachClustersIter(stackLayer, 1000.)) continue;
728                         tracklet.Init(&t);
729
730                         //if(!tracklet.AttachClusters(stackLayer, kTRUE)) continue;
731                         //if(!tracklet.Fit()) continue;
732                 }
733                 Int_t ncl = tracklet.GetN();
734                 //AliInfo(Form("N clusters %d", ncl));
735                 
736                 // Discard tracklet if bad quality.
737                 //Check if this condition is not already checked during building of the tracklet
738     if(ncl < fTimeBinsPerPlane * fRecoParam->GetFindableClusters()){
739                         //AliInfo(Form("Discard tracklet for %d nclusters", ncl));
740                         continue;
741                 }
742                 
743                 // load tracklet to the tracker and the track
744                 Int_t index = SetTracklet(&tracklet);
745                 t.SetTracklet(&tracklet, iplane, index);
746                 
747                 // Calculate the mean material budget along the path inside the chamber
748     Double_t param[7];
749                 AliTracker::MeanMaterialBudget(xyz0, xyz1, param);      
750     // The mean propagation parameters
751     Double_t xrho = param[0]*param[4]; // density*length
752     Double_t xx0  = param[1]; // radiation length
753                 
754                 // Propagate and update track
755                 t.PropagateTo(tracklet.GetX0(), xx0, xrho);
756           if (!AdjustSector(&t)) break;
757                 Double_t maxChi2 = t.GetPredictedChi2(&tracklet);
758                 if (maxChi2<1e+10 && t.Update(&tracklet, maxChi2)){ 
759                         nClustersExpected += ncl;
760                 }
761                 // Reset material budget if 2 consecutive gold
762                 if(iplane>0 && ncl + t.GetTracklet(iplane-1)->GetN() > 20) t.SetBudget(2, 0.);
763
764                 // Make backup of the track until is gold
765                 // TO DO update quality check of the track.
766                 // consider comparison with fTimeBinsRange
767                 Float_t ratio0 = ncl / Float_t(fTimeBinsPerPlane);
768                 //Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);    
769     //printf("tracklet.GetChi2() %f     [< 18.0]\n", tracklet.GetChi2()); 
770                 //printf("ratio0    %f              [>   0.8]\n", ratio0);
771                 //printf("ratio1     %f             [>   0.6]\n", ratio1); 
772                 //printf("ratio0+ratio1 %f          [>   1.5]\n", ratio0+ratio1); 
773                 //printf("t.GetNCross()  %d         [==    0]\n", t.GetNCross()); 
774                 //printf("TMath::Abs(t.GetSnp()) %f [<  0.85]\n", TMath::Abs(t.GetSnp()));
775                 //printf("t.GetNumberOfClusters() %d [>    20]\n", t.GetNumberOfClusters());
776     
777                 if (//(tracklet.GetChi2()      <  18.0) && TO DO check with FindClusters and move it to AliTRDseed::Update 
778         (ratio0                  >   0.8) && 
779         //(ratio1                  >   0.6) && 
780         //(ratio0+ratio1           >   1.5) && 
781         (t.GetNCross()           ==    0) && 
782         (TMath::Abs(t.GetSnp())  <  0.85) &&
783         (t.GetNumberOfClusters() >    20)) t.MakeBackupTrack();
784                 
785         } // end planes loop
786
787 #ifdef DEBUG
788         if(AliTRDReconstructor::StreamLevel() > 1){
789                 TTreeSRedirector &cstreamer = *fDebugStreamer;
790                 cstreamer << "FollowBackProlongation"
791                         << "ncl="      << nClustersExpected
792                         << "track.="   << &t
793                         << "\n";
794         }
795 #endif
796         
797         return nClustersExpected;
798 }
799
800 //____________________________________________________________________
801 void AliTRDtrackerV1::UnloadClusters() 
802
803   //
804   // Clears the arrays of clusters and tracks. Resets sectors and timebins 
805   //
806
807   Int_t i;
808   Int_t nentr;
809
810   nentr = fClusters->GetEntriesFast();
811         //AliInfo(Form("clearing %d clusters", nentr));
812   for (i = 0; i < nentr; i++) {
813     delete fClusters->RemoveAt(i);
814   }
815   fNclusters = 0;
816
817   nentr = fTracklets->GetEntriesFast();
818         //AliInfo(Form("clearing %d tracklets", nentr));
819   for (i = 0; i < nentr; i++) {
820     delete fTracklets->RemoveAt(i);
821   }
822
823   nentr = fSeeds->GetEntriesFast();
824         //AliInfo(Form("clearing %d seeds", nentr));
825   for (i = 0; i < nentr; i++) {
826     delete fSeeds->RemoveAt(i);
827   }
828
829   nentr = fTracks->GetEntriesFast();
830   //AliInfo(Form("clearing %d tracks", nentr));
831         for (i = 0; i < nentr; i++) {
832     delete fTracks->RemoveAt(i);
833   }
834
835   Int_t nsec = AliTRDgeometry::kNsect;
836   for (i = 0; i < nsec; i++) {    
837     for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
838       fTrSec[i]->GetLayer(pl)->Clear();
839     }
840   }
841
842 }
843
844 //____________________________________________________________________
845 AliTRDseedV1* AliTRDtrackerV1::GetTracklet(AliTRDtrackV1 *track, Int_t p, Int_t &idx)
846 {
847 // Find tracklet for TRD track <track>
848 // Parameters
849 // - track
850 // - sector
851 // - plane
852 // - index
853 // Output
854 // tracklet
855 // index
856 // Detailed description
857 //
858         idx = track->GetTrackletIndex(p);
859         //AliInfo(Form("looking for tracklet in plane %d idx %d [%d]", p, idx, track->GetTrackletIndex(p)));
860         AliTRDseedV1 *tracklet = idx<0 ? 0x0 : (AliTRDseedV1*)fTracklets->UncheckedAt(idx);
861         //AliInfo(Form("found 0x%x @ %d", tracklet, idx));
862
863 //   Int_t *index = track->GetTrackletIndexes();
864 //   for (UInt_t i = 0; i < 6; i++) AliInfo(Form("index[%d] = %d", i, index[i]));
865 // 
866 //      for (UInt_t i = 0; i < 6/*kMaxTimeBinIndex*/; i++) {
867 //     if (index[i] < 0) continue;
868 // 
869 //              tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index[i]);
870 //              if(!tracklet) break;
871 // 
872 //              if(tracklet->GetPlane() != p) continue;
873 // 
874 //              idx = index[i];
875 //      }
876
877         return tracklet;
878 }
879
880 //____________________________________________________________________
881 Int_t AliTRDtrackerV1::SetTracklet(AliTRDseedV1 *tracklet)
882 {
883 // Add this tracklet to the list of tracklets stored in the tracker
884 //
885 // Parameters
886 //   - tracklet : pointer to the tracklet to be added to the list
887 //
888 // Output
889 //   - the index of the new tracklet in the tracker tracklets list
890 //
891 // Detailed description
892 // Build the tracklets list if it is not yet created (late initialization)
893 // and adds the new tracklet to the list.
894 //
895         if(!fTracklets){
896                 fTracklets = new TClonesArray("AliTRDseedV1", AliTRDgeometry::Nsect()*kMaxTracksStack);
897                 fTracklets->SetOwner(kTRUE);
898         }
899         Int_t nentries = fTracklets->GetEntriesFast();
900         new ((*fTracklets)[nentries]) AliTRDseedV1(*tracklet);
901         return nentries;
902 }
903
904 //____________________________________________________________________
905 Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *sector
906                                        , AliESDEvent *esd)
907 {
908   //
909   // Steer tracking for one SM.
910   //
911   // Parameters :
912   //   sector  : Array of (SM) propagation layers containing clusters
913   //   esd     : The current ESD event. On output it contains the also
914   //             the ESD (TRD) tracks found in this SM. 
915   //
916   // Output :
917   //   Number of tracks found in this TRD supermodule.
918   // 
919   // Detailed description
920   //
921   // 1. Unpack AliTRDpropagationLayers objects for each stack.
922   // 2. Launch stack tracking. 
923   //    See AliTRDtrackerV1::Clusters2TracksStack() for details.
924   // 3. Pack results in the ESD event.
925   //
926         
927         AliTRDpadPlane *pp = 0x0;
928         
929         // allocate space for esd tracks in this SM
930         TClonesArray esdTrackList("AliESDtrack", 2*kMaxTracksStack);
931         esdTrackList.SetOwner();
932         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
933         Int_t nTimeBins = cal->GetNumberOfTimeBins();
934         const Int_t kFindable = Int_t(fRecoParam->GetFindableClusters()*6.*nTimeBins);
935         
936         Int_t ntracks = 0;
937         Int_t nClStack = 0;
938         for(int istack = 0; istack<AliTRDpropagationLayer::kZones; istack++){
939                 AliTRDstackLayer stackLayer[kNPlanes*kNTimeBins];
940                 
941                 nClStack = 0;
942                 //AliInfo(Form("Processing stack %i ...",istack));
943                 //AliInfo("Building stack propagation layers ...");
944                 for(int ilayer=0; ilayer<kNPlanes*nTimeBins; ilayer++){
945                         pp = fGeom->GetPadPlane((Int_t)(ilayer/nTimeBins), istack);
946                         Double_t stacklength = (pp->GetNrows() - 2) * pp->GetLengthIPad() 
947                                              + 2 * pp->GetLengthOPad() + 2 * pp->GetLengthRim();
948                         //Debug
949                         Double_t z0  = fGeom->GetRow0((Int_t)(ilayer/nTimeBins),istack,0);
950                         const AliTRDpropagationLayer ksmLayer(*(sector->GetLayer(ilayer)));
951                         stackLayer[ilayer] = ksmLayer;
952 #ifdef DEBUG
953                         stackLayer[ilayer].SetDebugStream(fDebugStreamer);
954 #endif                  
955                         stackLayer[ilayer].SetRange(z0 - stacklength, stacklength);
956                         stackLayer[ilayer].SetSector(sector->GetSector());
957                         stackLayer[ilayer].SetStackNr(istack);
958                         stackLayer[ilayer].SetNRows(pp->GetNrows());
959                         stackLayer[ilayer].SetRecoParam(fRecoParam);
960                         stackLayer[ilayer].BuildIndices();
961                         nClStack += stackLayer[ilayer].GetNClusters();
962                 }
963                 //AliInfo(Form("Finish building stack propagation layers. nClusters %d.", nClStack));
964                 if(nClStack < kFindable) continue;
965                 ntracks += Clusters2TracksStack(&stackLayer[0], &esdTrackList);
966         }
967         //AliInfo(Form("Found %d tracks in SM", ntracks));
968         
969         for(int itrack=0; itrack<ntracks; itrack++) 
970           esd->AddTrack((AliESDtrack*)esdTrackList[itrack]);
971
972         return ntracks;
973 }
974
975 //____________________________________________________________________
976 Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
977                                           , TClonesArray *esdTrackList)
978 {
979   //
980   // Make tracks in one TRD stack.
981   //
982   // Parameters :
983   //   layer  : Array of stack propagation layers containing clusters
984   //   esdTrackList  : Array of ESD tracks found by the stand alone tracker. 
985   //                   On exit the tracks found in this stack are appended.
986   //
987   // Output :
988   //   Number of tracks found in this stack.
989   // 
990   // Detailed description
991   //
992   // 1. Find the 3 most useful seeding chambers. See BuildSeedingConfigs() for details.
993   // 2. Steer AliTRDtrackerV1::MakeSeeds() for 3 seeding layer configurations. 
994   //    See AliTRDtrackerV1::MakeSeeds() for more details.
995   // 3. Arrange track candidates in decreasing order of their quality
996   // 4. Classify tracks in 5 categories according to:
997   //    a) number of layers crossed
998   //    b) track quality 
999   // 5. Sign clusters by tracks in decreasing order of track quality
1000   // 6. Build AliTRDtrack out of seeding tracklets
1001   // 7. Cook MC label
1002   // 8. Build ESD track and register it to the output list
1003   //
1004
1005         AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized
1006         Int_t pars[4]; // MakeSeeds parameters
1007
1008         //Double_t alpha = AliTRDgeometry::GetAlpha();
1009         //Double_t shift = .5 * alpha;
1010         Int_t configs[kNConfigs];
1011         
1012         // Build initial seeding configurations
1013         Double_t quality = BuildSeedingConfigs(layer, configs);
1014 #ifdef DEBUG
1015                 if(AliTRDReconstructor::StreamLevel() > 1) 
1016                   AliInfo(Form("Plane config %d %d %d Quality %f"
1017                               , configs[0], configs[1], configs[2], quality));
1018 #endif
1019
1020         // Initialize contors
1021         Int_t ntracks,      // number of TRD track candidates
1022               ntracks1,     // number of registered TRD tracks/iter
1023               ntracks2 = 0; // number of all registered TRD tracks in stack
1024         fSieveSeeding = 0;
1025         do{
1026                 // Loop over seeding configurations
1027                 ntracks = 0; ntracks1 = 0;
1028                 for (Int_t iconf = 0; iconf<3; iconf++) {
1029                         pars[0] = configs[iconf];
1030                         pars[1] = layer->GetStackNr();
1031                         pars[2] = ntracks;
1032                         ntracks = MakeSeeds(layer, &sseed[6*ntracks], pars);
1033                         if(ntracks == kMaxTracksStack) break;
1034                 }
1035 #ifdef DEBUG
1036                 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Candidate TRD tracks %d in stack %d iteration %d.", ntracks, pars[1], fSieveSeeding));
1037 #endif          
1038                 if(!ntracks) break;
1039                 
1040                 // Sort the seeds according to their quality
1041                 Int_t sort[kMaxTracksStack];
1042                 TMath::Sort(ntracks, fTrackQuality, sort, kTRUE);
1043         
1044                 // Initialize number of tracks so far and logic switches
1045                 Int_t ntracks0 = esdTrackList->GetEntriesFast();
1046                 Bool_t signedTrack[kMaxTracksStack];
1047                 Bool_t fakeTrack[kMaxTracksStack];
1048                 for (Int_t i=0; i<ntracks; i++){
1049                         signedTrack[i] = kFALSE;
1050                         fakeTrack[i] = kFALSE;
1051                 }
1052                 //AliInfo("Selecting track candidates ...");
1053                 
1054                 // Sieve clusters in decreasing order of track quality
1055                 Double_t trackParams[7];
1056 //              AliTRDseedV1 *lseed = 0x0;
1057                 Int_t jSieve = 0, candidates;
1058                 do{
1059                         //AliInfo(Form("\t\tITER = %i ", jSieve));
1060
1061                         // Check track candidates
1062                         candidates = 0;
1063                         for (Int_t itrack = 0; itrack < ntracks; itrack++) {
1064                                 Int_t trackIndex = sort[itrack];
1065                                 if (signedTrack[trackIndex] || fakeTrack[trackIndex]) continue;
1066         
1067                                 
1068                                 // Calculate track parameters from tracklets seeds
1069                                 Int_t labelsall[1000];
1070                                 Int_t nlabelsall = 0;
1071                                 Int_t naccepted  = 0;
1072                                 Int_t ncl        = 0;
1073                                 Int_t nused      = 0;
1074                                 Int_t nlayers    = 0;
1075                                 Int_t findable   = 0;
1076                                 for (Int_t jLayer = 0; jLayer < kNPlanes; jLayer++) {
1077                                         Int_t jseed = kNPlanes*trackIndex+jLayer;
1078                                         if (TMath::Abs(sseed[jseed].GetYref(0) / sseed[jseed].GetX0()) < 0.15) 
1079                                           findable++;
1080         
1081                                         if(!sseed[jseed].IsOK()) continue;
1082                                         sseed[jseed].UpdateUsed();
1083                                         ncl   += sseed[jseed].GetN2();
1084                                         nused += sseed[jseed].GetNUsed();
1085                                         nlayers++;
1086         
1087                                         // Cooking label
1088                                         for (Int_t itime = 0; itime < fTimeBinsPerPlane; itime++) {
1089                                                 if(!sseed[jseed].IsUsable(itime)) continue;
1090                                                 naccepted++;
1091                                                 Int_t tindex = 0, ilab = 0;
1092                                                 while(ilab<3 && (tindex = sseed[jseed].GetClusters(itime)->GetLabel(ilab)) >= 0){
1093                                                         labelsall[nlabelsall++] = tindex;
1094                                                         ilab++;
1095                                                 }
1096                                         }
1097                                 }
1098                                 // Filter duplicated tracks
1099                                 if (nused > 30){
1100                                         printf("Skip %d nused %d\n", trackIndex, nused);
1101                                         fakeTrack[trackIndex] = kTRUE;
1102                                         continue;
1103                                 }
1104                                 if (Float_t(nused)/ncl >= .25){
1105                                         printf("Skip %d nused/ncl >= .25\n", trackIndex);
1106                                         fakeTrack[trackIndex] = kTRUE;
1107                                         continue;
1108                                 }
1109                                 
1110                                 // Classify tracks
1111                                 Bool_t skip = kFALSE;
1112                                 switch(jSieve){
1113                                 case 0:
1114                                         if(nlayers < 6) {skip = kTRUE; break;}
1115                                         if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
1116                                         break;
1117         
1118                                 case 1:
1119                                         if(nlayers < findable){skip = kTRUE; break;}
1120                                         if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -4.){skip = kTRUE; break;}
1121                                         break;
1122         
1123                                 case 2:
1124                                         if ((nlayers == findable) || (nlayers == 6)) { skip = kTRUE; break;}
1125                                         if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -6.0){skip = kTRUE; break;}
1126                                         break;
1127         
1128                                 case 3:
1129                                         if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
1130                                         break;
1131         
1132                                 case 4:
1133                                         if (nlayers == 3){skip = kTRUE; break;}
1134                                         //if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) - nused/(nlayers-3.0) < -15.0){skip = kTRUE; break;}
1135                                         break;
1136                                 }
1137                                 if(skip){
1138                                         candidates++;
1139                                         printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused);
1140                                         continue;
1141                                 }
1142                                 signedTrack[trackIndex] = kTRUE;
1143                                                 
1144
1145                                 // Build track label - what happens if measured data ???
1146                                 Int_t labels[1000];
1147                                 Int_t outlab[1000];
1148                                 Int_t nlab = 0;
1149                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1150                                         Int_t jseed = kNPlanes*trackIndex+iLayer;
1151                                         if(!sseed[jseed].IsOK()) continue;
1152                                         for(int ilab=0; ilab<2; ilab++){
1153                                                 if(sseed[jseed].GetLabels(ilab) < 0) continue;
1154                                                 labels[nlab] = sseed[jseed].GetLabels(ilab);
1155                                                 nlab++;
1156                                         }
1157                                 }
1158                                 Freq(nlab,labels,outlab,kFALSE);
1159                                 Int_t   label     = outlab[0];
1160                                 Int_t   frequency = outlab[1];
1161                                 Freq(nlabelsall,labelsall,outlab,kFALSE);
1162                                 Int_t   label1    = outlab[0];
1163                                 Int_t   label2    = outlab[2];
1164                                 Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
1165         
1166                                 
1167                                 // Sign clusters
1168                                 AliTRDcluster *cl = 0x0; Int_t clusterIndex = -1;
1169                                 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1170                                         Int_t jseed = kNPlanes*trackIndex+jLayer;
1171                                         if(!sseed[jseed].IsOK()) continue;
1172                                         if(TMath::Abs(sseed[jseed].GetYfit(1) - sseed[jseed].GetYfit(1)) >= .2) continue; // check this condition with Marian
1173                                         sseed[jseed].UseClusters();
1174                                         if(!cl){
1175                                                 Int_t ic = 0;
1176                                                 while(!(cl = sseed[jseed].GetClusters(ic))) ic++;
1177                                                 clusterIndex =  sseed[jseed].GetIndexes(ic);
1178                                         }
1179                                 }
1180                                 if(!cl) continue;
1181
1182                                 
1183                                 // Build track parameters
1184                                 AliTRDseedV1 *lseed =&sseed[trackIndex*6];
1185                                 Int_t idx = 0;
1186                                 while(idx<3 && !lseed->IsOK()) {
1187                                         idx++;
1188                                         lseed++;
1189                                 }
1190                                 Double_t cR = lseed->GetC();
1191                                 trackParams[1] = lseed->GetYref(0);
1192                                 trackParams[2] = lseed->GetZref(0);
1193                                 trackParams[3] = lseed->GetX0() * cR - TMath::Sin(TMath::ATan(lseed->GetYref(1)));
1194                                 trackParams[4] = lseed->GetZref(1) / TMath::Sqrt(1. + lseed->GetYref(1) * lseed->GetYref(1));
1195                                 trackParams[5] = cR;
1196                                 trackParams[0] = lseed->GetX0();
1197                                 trackParams[6] = layer[0].GetSector();/* *alpha+shift;  // Supermodule*/
1198
1199 #ifdef DEBUG
1200                                 if(AliTRDReconstructor::StreamLevel() > 1) printf("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]);
1201                                 
1202                                 if(AliTRDReconstructor::StreamLevel() >= 1){
1203                                         Int_t sector = layer[0].GetSector();
1204                                         Int_t nclusters = 0;
1205                                         AliTRDseedV1 *dseed[6];
1206                                         for(int is=0; is<6; is++){
1207                                                 dseed[is] = new AliTRDseedV1(sseed[trackIndex*6+is]);
1208                                                 dseed[is]->SetOwner();
1209                                                 nclusters += sseed[is].GetN2();
1210                                                 //for(int ic=0; ic<30; ic++) if(sseed[trackIndex*6+is].GetClusters(ic)) printf("l[%d] tb[%d] cptr[%p]\n", is, ic, sseed[trackIndex*6+is].GetClusters(ic));
1211                                         }
1212                                         //Int_t eventNrInFile = esd->GetEventNumberInFile();
1213                                         //AliInfo(Form("Number of clusters %d.", nclusters));
1214                                         TTreeSRedirector &cstreamer = *fDebugStreamer;
1215                                         cstreamer << "Clusters2TracksStack"
1216                                                 << "Iter="      << fSieveSeeding
1217                                                 << "Like="      << fTrackQuality[trackIndex]
1218                                                 << "S0.="       << dseed[0]
1219                                                 << "S1.="       << dseed[1]
1220                                                 << "S2.="       << dseed[2]
1221                                                 << "S3.="       << dseed[3]
1222                                                 << "S4.="       << dseed[4]
1223                                                 << "S5.="       << dseed[5]
1224                                                 << "p0=" << trackParams[0]
1225                                                 << "p1=" << trackParams[1]
1226                                                 << "p2=" << trackParams[2]
1227                                                 << "p3=" << trackParams[3]
1228                                                 << "p4=" << trackParams[4]
1229                                                 << "p5=" << trackParams[5]
1230                                                 << "p6=" << trackParams[6]
1231                                                 << "Sector="    << sector
1232                                                 << "Stack="     << pars[1]
1233                                                 << "Label="     << label
1234                                                 << "Label1="    << label1
1235                                                 << "Label2="    << label2
1236                                                 << "FakeRatio=" << fakeratio
1237                                                 << "Freq="      << frequency
1238                                                 << "Ncl="       << ncl
1239                                                 << "NLayers="   << nlayers
1240                                                 << "Findable="  << findable
1241                                                 << "NUsed="     << nused
1242                                                 << "\n";
1243                                         //???for(int is=0; is<6; is++) delete dseed[is];
1244                                 }
1245 #endif
1246                         
1247                                 AliTRDtrackV1 *track = AliTRDtrackerV1::MakeTrack(&sseed[trackIndex*kNPlanes], trackParams);
1248                                 if(!track){
1249                                         AliWarning("Fail to build a TRD Track.");
1250                                         continue;
1251                                 }
1252                                 AliInfo("End of MakeTrack()");
1253                                 AliESDtrack esdTrack;
1254                                 esdTrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
1255                                 esdTrack.SetLabel(track->GetLabel());
1256                                 new ((*esdTrackList)[ntracks0++]) AliESDtrack(esdTrack);
1257                                 ntracks1++;
1258                         }
1259
1260                         jSieve++;
1261                 } while(jSieve<5 && candidates); // end track candidates sieve
1262                 if(!ntracks1) break;
1263
1264                 // increment counters
1265                 ntracks2 += ntracks1;
1266                 fSieveSeeding++;
1267
1268                 // Rebuild plane configurations and indices taking only unused clusters into account
1269                 quality = BuildSeedingConfigs(layer, configs);
1270                 //if(quality < fRecoParam->GetPlaneQualityThreshold()) break;
1271                 
1272                 for(Int_t il = 0; il < kNPlanes * fTimeBinsPerPlane; il++) layer[il].BuildIndices(fSieveSeeding);
1273
1274 #ifdef DEBUG
1275                                 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Sieve level %d Plane config %d %d %d Quality %f", fSieveSeeding, configs[0], configs[1], configs[2], quality));
1276 #endif
1277         } while(fSieveSeeding<10); // end stack clusters sieve
1278         
1279
1280
1281         //AliInfo(Form("Registered TRD tracks %d in stack %d.", ntracks2, pars[1]));
1282
1283         return ntracks2;
1284 }
1285
1286 //___________________________________________________________________
1287 Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers
1288                                             , Int_t *configs)
1289 {
1290   //
1291   // Assign probabilities to chambers according to their
1292   // capability of producing seeds.
1293   // 
1294   // Parameters :
1295   //
1296   //   layers : Array of stack propagation layers for all 6 chambers in one stack
1297   //   configs : On exit array of configuration indexes (see GetSeedingConfig()
1298   // for details) in the decreasing order of their seeding probabilities. 
1299   //
1300   // Output :
1301   //
1302   //  Return top configuration quality 
1303   //
1304   // Detailed description:
1305   //
1306   // To each chamber seeding configuration (see GetSeedingConfig() for
1307   // the list of all configurations) one defines 2 quality factors:
1308   //  - an apriori topological quality (see GetSeedingConfig() for details) and
1309   //  - a data quality based on the uniformity of the distribution of
1310   //    clusters over the x range (time bins population). See CookChamberQA() for details.
1311   // The overall chamber quality is given by the product of this 2 contributions.
1312   // 
1313
1314         Double_t chamberQA[kNPlanes];
1315         for(int iplane=0; iplane<kNPlanes; iplane++){
1316                 chamberQA[iplane] = MakeSeedingPlanes(&layers[iplane*fTimeBinsPerPlane]);
1317                 //printf("chamberQA[%d] = %f\n", iplane, chamberQA[iplane]);
1318         }
1319
1320         Double_t tconfig[kNConfigs];
1321         Int_t planes[4];
1322         for(int iconf=0; iconf<kNConfigs; iconf++){
1323                 GetSeedingConfig(iconf, planes);
1324                 tconfig[iconf] = fgTopologicQA[iconf];
1325                 for(int iplane=0; iplane<4; iplane++) tconfig[iconf] *= chamberQA[planes[iplane]]; 
1326         }
1327         
1328         TMath::Sort(kNConfigs, tconfig, configs, kTRUE);
1329         return tconfig[configs[0]];
1330 }
1331
1332 //____________________________________________________________________
1333 Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
1334                                , AliTRDseedV1 *sseed
1335                                , Int_t *ipar)
1336 {
1337   //
1338   // Make tracklet seeds in the TRD stack.
1339   //
1340   // Parameters :
1341   //   layers : Array of stack propagation layers containing clusters
1342   //   sseed  : Array of empty tracklet seeds. On exit they are filled.
1343   //   ipar   : Control parameters:
1344   //       ipar[0] -> seeding chambers configuration
1345   //       ipar[1] -> stack index
1346   //       ipar[2] -> number of track candidates found so far
1347   //
1348   // Output :
1349   //   Number of tracks candidates found.
1350   // 
1351   // Detailed description
1352   //
1353   // The following steps are performed:
1354   // 1. Select seeding layers from seeding chambers
1355   // 2. Select seeding clusters from the seeding AliTRDpropagationLayerStack.
1356   //   The clusters are taken from layer 3, layer 0, layer 1 and layer 2, in
1357   //   this order. The parameters controling the range of accepted clusters in
1358   //   layer 0, 1, and 2 are defined in AliTRDstackLayer::BuildCond().
1359   // 3. Helix fit of the cluster set. (see AliTRDtrackerFitter::FitRieman(AliTRDcluster**))
1360   // 4. Initialize seeding tracklets in the seeding chambers.
1361   // 5. Filter 0.
1362   //   Chi2 in the Y direction less than threshold ... (1./(3. - sLayer))
1363   //   Chi2 in the Z direction less than threshold ... (1./(3. - sLayer))
1364   // 6. Attach clusters to seeding tracklets and find linear approximation of
1365   //   the tracklet (see AliTRDseedV1::AttachClustersIter()). The number of used
1366   //   clusters used by current seeds should not exceed ... (25).
1367   // 7. Filter 1.
1368   //   All 4 seeding tracklets should be correctly constructed (see
1369   //   AliTRDseedV1::AttachClustersIter())
1370   // 8. Helix fit of the seeding tracklets
1371   // 9. Filter 2.
1372   //   Likelihood calculation of the fit. (See AliTRDtrackerV1::CookLikelihood() for details)
1373   // 10. Extrapolation of the helix fit to the other 2 chambers:
1374   //    a) Initialization of extrapolation tracklet with fit parameters
1375   //    b) Helix fit of tracklets
1376   //    c) Attach clusters and linear interpolation to extrapolated tracklets
1377   //    d) Helix fit of tracklets
1378   // 11. Improve seeding tracklets quality by reassigning clusters.
1379   //      See AliTRDtrackerV1::ImproveSeedQuality() for details.
1380   // 12. Helix fit of all 6 seeding tracklets and chi2 calculation
1381   // 13. Hyperplane fit and track quality calculation. See AliTRDtrackerFitter::FitHyperplane() for details.
1382   // 14. Cooking labels for tracklets. Should be done only for MC
1383   // 15. Register seeds.
1384   //
1385
1386         AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0}; // initilize seeding clusters
1387         AliTRDseedV1 *cseed = &sseed[0]; // initialize tracklets for first track
1388         Int_t ncl, mcl; // working variable for looping over clusters
1389         Int_t index[AliTRDstackLayer::kMaxClustersLayer], jndex[AliTRDstackLayer::kMaxClustersLayer];
1390         // chi2 storage
1391         // chi2[0] = tracklet chi2 on the Z direction
1392         // chi2[1] = tracklet chi2 on the R direction
1393         Double_t chi2[4];
1394
1395
1396         // this should be data member of AliTRDtrack
1397         Double_t seedQuality[kMaxTracksStack];
1398         
1399         // unpack control parameters
1400         Int_t config  = ipar[0];
1401         Int_t istack  = ipar[1];
1402         Int_t ntracks = ipar[2];
1403         Int_t planes[kNSeedPlanes]; GetSeedingConfig(config, planes);   
1404 #ifdef DEBUG
1405                 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks));
1406 #endif
1407         
1408         // Init chambers geometry
1409         Int_t det/*, tbRange[6]*/; // time bins inside the detector geometry
1410         Double_t hL[kNPlanes];       // Tilting angle
1411         Float_t padlength[kNPlanes]; // pad lenghts
1412         AliTRDpadPlane *pp;
1413         for(int il=0; il<kNPlanes; il++){
1414                 pp = fGeom->GetPadPlane(il, istack); // istack has to be imported
1415                 hL[il]        = TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle());
1416                 padlength[il] = pp->GetLengthIPad();
1417                 det           = il; // to be fixed !!!!!
1418                 //tbRange[il]   = fTimeBinsPerPlane; //Int_t(AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght() * AliTRDCommonParam::Instance()->GetSamplingFrequency()/AliTRDcalibDB::Instance()->GetVdriftDet()->GetValue(det));
1419                 //printf("%d hl[%f] pl[%f] tb[%d]\n", il, hL[il], padlength[il], tbRange[il]);
1420         }
1421
1422         Double_t cond0[4], cond1[4], cond2[4];
1423         // make seeding layers (to be moved in Clusters2TracksStack)
1424         AliTRDstackLayer *layer[] = {0x0, 0x0, 0x0, 0x0};
1425         for(int isl=0; isl<kNSeedPlanes; isl++) layer[isl] = MakeSeedingLayer(&layers[planes[isl] * fTimeBinsPerPlane], planes[isl]);
1426
1427
1428         // Start finding seeds
1429         Int_t icl = 0;
1430         while((c[3] = (*layer[3])[icl++])){
1431                 if(!c[3]) continue;
1432                 layer[0]->BuildCond(c[3], cond0, 0);
1433                 layer[0]->GetClusters(cond0, index, ncl);
1434                 Int_t jcl = 0;
1435                 while(jcl<ncl) {
1436                         c[0] = (*layer[0])[index[jcl++]];
1437                         if(!c[0]) continue;
1438                         Double_t dx    = c[3]->GetX() - c[0]->GetX();
1439                         Double_t theta = (c[3]->GetZ() - c[0]->GetZ())/dx;
1440                         Double_t phi   = (c[3]->GetY() - c[0]->GetY())/dx;
1441                         layer[1]->BuildCond(c[0], cond1, 1, theta, phi);
1442                         layer[1]->GetClusters(cond1, jndex, mcl);
1443
1444                         Int_t kcl = 0;
1445                         while(kcl<mcl) {
1446                                 c[1] = (*layer[1])[jndex[kcl++]];
1447                                 if(!c[1]) continue;
1448                                 layer[2]->BuildCond(c[1], cond2, 2, theta, phi);
1449                                 c[2] = layer[2]->GetNearestCluster(cond2);
1450                                 if(!c[2]) continue;
1451                                 
1452                                 //AliInfo("Seeding clusters found. Building seeds ...");
1453                                 //for(Int_t i = 0; i < kNSeedPlanes; i++) printf("%i. coordinates: x = %3.3f, y = %3.3f, z = %3.3f\n", i, c[i]->GetX(), c[i]->GetY(), c[i]->GetZ());
1454                                 for (Int_t il = 0; il < 6; il++) cseed[il].Reset();
1455
1456                                 fFitter->Reset();
1457
1458                                 fFitter->FitRieman(c, kNSeedPlanes);
1459
1460                                 chi2[0] = 0.; chi2[1] = 0.;
1461                                 AliTRDseedV1 *tseed = 0x0;
1462                                 for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
1463                                         Int_t jLayer = planes[iLayer];
1464                                         tseed = &cseed[jLayer];
1465                                         tseed->SetRecoParam(fRecoParam);
1466                                         tseed->SetPlane(jLayer);
1467                                         tseed->SetTilt(hL[jLayer]);
1468                                         tseed->SetPadLength(padlength[jLayer]);
1469                                         //tseed->SetNTimeBinsRange(tbRange[jLayer]);
1470                                         tseed->SetX0(layer[iLayer]->GetX());//layers[jLayer*fTimeBinsPerPlane].GetX());
1471
1472                                         tseed->Init(fFitter->GetRiemanFitter());
1473                                         // temporary until new AttachClusters()
1474                                         tseed->SetX0(layers[(jLayer+1)*fTimeBinsPerPlane-1].GetX());
1475                                         chi2[0] += tseed->GetChi2Z(c[iLayer]->GetZ());
1476                                         chi2[1] += tseed->GetChi2Y(c[iLayer]->GetY());
1477                                 }
1478
1479                                 Bool_t isFake = kFALSE;
1480                                 if (c[0]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
1481                                 if (c[1]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
1482                                 if (c[2]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
1483 #ifdef DEBUG
1484                                 if(AliTRDReconstructor::StreamLevel() >= 2){
1485                                         Float_t yref[4], ycluster[4];
1486                                         for(int il=0; il<4; il++){
1487                                                 tseed = &cseed[planes[il]];
1488                                                 yref[il] = tseed->GetYref(0);
1489                                                 ycluster[il] = c[il]->GetY();
1490                                         }
1491                                         Float_t threshold = .5;//1./(3. - sLayer);
1492                                         Int_t ll = c[3]->GetLabel(0);
1493                                         TTreeSRedirector &cs0 = *fDebugStreamer;
1494                                                         cs0 << "MakeSeeds0"
1495                                                         <<"isFake=" << isFake
1496                                                         <<"label=" << ll
1497                                                         <<"threshold=" << threshold
1498                                                         <<"chi2=" << chi2[1]
1499                                                         <<"yref0="<<yref[0]
1500                                                         <<"yref1="<<yref[1]
1501                                                         <<"yref2="<<yref[2]
1502                                                         <<"yref3="<<yref[3]
1503                                                         <<"ycluster0="<<ycluster[0]
1504                                                         <<"ycluster1="<<ycluster[1]
1505                                                         <<"ycluster2="<<ycluster[2]
1506                                                         <<"ycluster3="<<ycluster[3]
1507                                                         <<"\n";
1508                                 }
1509 #endif
1510
1511                                 if(chi2[0] > fRecoParam->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){
1512                                         //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0]));
1513                                         continue;
1514                                 }
1515                                 if(chi2[1] > fRecoParam->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){
1516                                         //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1]));
1517                                         continue;
1518                                 }
1519                                 //AliInfo("Passed chi2 filter.");
1520
1521 #ifdef DEBUG
1522                                 if(AliTRDReconstructor::StreamLevel() >= 2){
1523                                         Float_t minmax[2] = { -100.0,  100.0 };
1524                                         for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1525                                                 Float_t max = c[iLayer]->GetZ() + cseed[planes[iLayer]].GetPadLength() * 0.5 + 1.0 - cseed[planes[iLayer]].GetZref(0);
1526                                                 if (max < minmax[1]) minmax[1] = max;
1527                                                 Float_t min = c[iLayer]->GetZ()-cseed[planes[iLayer]].GetPadLength() * 0.5 - 1.0 - cseed[planes[iLayer]].GetZref(0);
1528                                                 if (min > minmax[0]) minmax[0] = min;
1529                                         }
1530                                         Double_t xpos[4];
1531                                         for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
1532                                         TTreeSRedirector &cstreamer = *fDebugStreamer;
1533                                                         cstreamer << "MakeSeeds1"
1534                                                 << "isFake=" << isFake
1535                                                 << "config="   << config
1536                                                 << "Cl0.="   << c[0]
1537                                                 << "Cl1.="   << c[1]
1538                                                 << "Cl2.="   << c[2]
1539                                                 << "Cl3.="   << c[3]
1540                                                 << "X0="     << xpos[0] //layer[sLayer]->GetX()
1541                                                 << "X1="     << xpos[1] //layer[sLayer + 1]->GetX()
1542                                                 << "X2="     << xpos[2] //layer[sLayer + 2]->GetX()
1543                                                 << "X3="     << xpos[3] //layer[sLayer + 3]->GetX()
1544                                                 << "Y2exp="  << cond2[0]
1545                                                 << "Z2exp="  << cond2[1]
1546                                                 << "Chi2R="  << chi2[0]
1547                                                 << "Chi2Z="  << chi2[1]
1548                                                 << "Seed0.=" << &cseed[planes[0]]
1549                                                 << "Seed1.=" << &cseed[planes[1]]
1550                                                 << "Seed2.=" << &cseed[planes[2]]
1551                                                 << "Seed3.=" << &cseed[planes[3]]
1552                                                 << "Zmin="   << minmax[0]
1553                                                 << "Zmax="   << minmax[1]
1554                                                 << "\n" ;
1555                                 }               
1556 #endif
1557                                 // try attaching clusters to tracklets
1558                                 Int_t nUsedCl = 0;
1559                                 Int_t nlayers = 0;
1560                                 for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
1561                                         Int_t jLayer = planes[iLayer];
1562                                         if(!cseed[jLayer].AttachClustersIter(&layers[jLayer*fTimeBinsPerPlane], 5., kFALSE, c[iLayer])) continue;
1563                                         nUsedCl += cseed[jLayer].GetNUsed();
1564                                         if(nUsedCl > 25) break;
1565                                         nlayers++;
1566                                 }
1567                                 if(nlayers < kNSeedPlanes){ 
1568                                         //AliInfo("Failed updating all seeds.");
1569                                         continue;
1570                                 }
1571                                 // fit tracklets and cook likelihood
1572                                 chi2[0] = 0.; chi2[1] = 0.;
1573                                 fFitter->FitRieman(&cseed[0], &planes[0]);
1574                                 AliRieman *rim = fFitter->GetRiemanFitter();
1575                                 for(int iLayer=0; iLayer<4; iLayer++){
1576                                         cseed[planes[iLayer]].Init(rim);
1577                                         chi2[0] += (Float_t)cseed[planes[iLayer]].GetChi2Z();
1578                                         chi2[1] += cseed[planes[iLayer]].GetChi2Y();
1579                                 }
1580                                 Double_t chi2r = chi2[1], chi2z = chi2[0];
1581                                 Double_t like = CookLikelihood(&cseed[0], planes, chi2); // to be checked
1582                                 if (TMath::Log(1.E-9 + like) < fRecoParam->GetTrackLikelihood()){
1583                                         //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
1584                                         continue;
1585                                 }
1586                                 //AliInfo(Form("Passed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
1587
1588
1589                                 // book preliminary results
1590                                 seedQuality[ntracks] = like;
1591                                 fSeedLayer[ntracks]  = config;/*sLayer;*/
1592
1593                                 // attach clusters to the extrapolation seeds
1594                                 Int_t lextrap[2];
1595                                 GetExtrapolationConfig(config, lextrap);
1596                                 Int_t nusedf   = 0; // debug value
1597                                 for(int iLayer=0; iLayer<2; iLayer++){
1598                                         Int_t jLayer = lextrap[iLayer];
1599                                         
1600                                         // prepare extrapolated seed
1601                                         cseed[jLayer].Reset();
1602                                         cseed[jLayer].SetRecoParam(fRecoParam);
1603                                         cseed[jLayer].SetPlane(jLayer);
1604                                         cseed[jLayer].SetTilt(hL[jLayer]);
1605                                         cseed[jLayer].SetX0(layers[(jLayer +1) * fTimeBinsPerPlane-1].GetX());
1606                                         cseed[jLayer].SetPadLength(padlength[jLayer]);
1607                                         //cseed[jLayer].SetNTimeBinsRange(tbRange[jLayer]);
1608                                         cseed[jLayer].Init(rim);
1609 //                                      AliTRDcluster *cd = FindSeedingCluster(&layers[jLayer*fTimeBinsPerPlane], &cseed[jLayer]);
1610 //                                      if(cd == 0x0) continue;
1611
1612                                         // fit extrapolated seed
1613                                         AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
1614                                         if ((jLayer == 0) && !(cseed[1].IsOK())) continue;
1615                                         if ((jLayer == 5) && !(cseed[4].IsOK())) continue;
1616                                         AliTRDseedV1 tseed = cseed[jLayer];
1617                                         if(!tseed.AttachClustersIter(&layers[jLayer*fTimeBinsPerPlane], 1000.)) continue;
1618                                         cseed[jLayer] = tseed;
1619                                         nusedf += cseed[jLayer].GetNUsed(); // debug value
1620                                         AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
1621                                 }
1622                                 //AliInfo("Extrapolation done.");
1623
1624                                 ImproveSeedQuality(layers, cseed);
1625                                 //AliInfo("Improve seed quality done.");
1626
1627                                 nlayers   = 0;
1628                                 Int_t nclusters = 0;
1629                                 Int_t findable  = 0;
1630                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1631                                         if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) findable++;
1632                                         if (!cseed[iLayer].IsOK()) continue;
1633                                         nclusters += cseed[iLayer].GetN2();
1634                                         nlayers++;
1635                                 }
1636                                 if (nlayers < 3){ 
1637                                         //AliInfo("Failed quality check on seeds.");
1638                                         continue;
1639                                 }
1640
1641                                 // fit full track and cook likelihoods
1642                                 fFitter->FitRieman(&cseed[0]);
1643                                 Double_t chi2ZF = 0., chi2RF = 0.;
1644                                 for(int ilayer=0; ilayer<6; ilayer++){
1645                                         cseed[ilayer].Init(fFitter->GetRiemanFitter());
1646                                         if (!cseed[ilayer].IsOK()) continue;
1647                                         //tchi2 = cseed[ilayer].GetChi2Z();
1648                                         //printf("layer %d chi2 %e\n", ilayer, tchi2);
1649                                         chi2ZF += cseed[ilayer].GetChi2Z();
1650                                         chi2RF += cseed[ilayer].GetChi2Y();
1651                                 }
1652                                 chi2ZF /= TMath::Max((nlayers - 3.), 1.);
1653                                 chi2RF /= TMath::Max((nlayers - 3.), 1.);
1654
1655                                 // do the final track fitting
1656                                 fFitter->SetLayers(nlayers);
1657 #ifdef DEBUG
1658                                 fFitter->SetDebugStream(fDebugStreamer);
1659 #endif
1660                                 fTrackQuality[ntracks] = fFitter->FitHyperplane(&cseed[0], chi2ZF, GetZ());
1661                                 Double_t param[3];
1662                                 Double_t chi2[2];
1663                                 fFitter->GetHyperplaneFitResults(param);
1664                                 fFitter->GetHyperplaneFitChi2(chi2);
1665                                 //AliInfo("Hyperplane fit done\n");
1666
1667                                 // finalize tracklets
1668                                 Int_t labels[12];
1669                                 Int_t outlab[24];
1670                                 Int_t nlab = 0;
1671                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1672                                         if (!cseed[iLayer].IsOK()) continue;
1673
1674                                         if (cseed[iLayer].GetLabels(0) >= 0) {
1675                                                 labels[nlab] = cseed[iLayer].GetLabels(0);
1676                                                 nlab++;
1677                                         }
1678
1679                                         if (cseed[iLayer].GetLabels(1) >= 0) {
1680                                                 labels[nlab] = cseed[iLayer].GetLabels(1);
1681                                                 nlab++;
1682                                         }
1683                                 }
1684                                 Freq(nlab,labels,outlab,kFALSE);
1685                                 Int_t label     = outlab[0];
1686                                 Int_t frequency = outlab[1];
1687                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1688                                         cseed[iLayer].SetFreq(frequency);
1689                                         cseed[iLayer].SetC(param[1]/*cR*/);
1690                                         cseed[iLayer].SetCC(param[0]/*cC*/);
1691                                         cseed[iLayer].SetChi2(chi2[0]);
1692                                         cseed[iLayer].SetChi2Z(chi2ZF);
1693                                 }
1694             
1695 #ifdef DEBUG
1696                                 if(AliTRDReconstructor::StreamLevel() >= 2){
1697                                         Double_t curv = (fFitter->GetRiemanFitter())->GetC();
1698                                         TTreeSRedirector &cstreamer = *fDebugStreamer;
1699                                         cstreamer << "MakeSeeds2"
1700                                                 << "C="       << curv
1701                                                 << "Chi2R="   << chi2r
1702                                                 << "Chi2Z="   << chi2z
1703                                                 << "Chi2TR="  << chi2[0]
1704                                                 << "Chi2TC="  << chi2[1]
1705                                                 << "Chi2RF="  << chi2RF
1706                                                 << "Chi2ZF="  << chi2ZF
1707                                                 << "Ncl="     << nclusters
1708                                                 << "Nlayers=" << nlayers
1709                                                 << "NUsedS="  << nUsedCl
1710                                                 << "NUsed="   << nusedf
1711                                                 << "Findable" << findable
1712                                                 << "Like="    << like
1713                                                 << "S0.="     << &cseed[0]
1714                                                 << "S1.="     << &cseed[1]
1715                                                 << "S2.="     << &cseed[2]
1716                                                 << "S3.="     << &cseed[3]
1717                                                 << "S4.="     << &cseed[4]
1718                                                 << "S5.="     << &cseed[5]
1719                                                 << "Label="   << label
1720                                                 << "Freq="    << frequency
1721                                                 << "\n";
1722                                 }
1723 #endif
1724                                 
1725                                 ntracks++;
1726                                 if(ntracks == kMaxTracksStack){
1727                                         AliWarning(Form("Number of seeds reached maximum allowed (%d) in stack.", kMaxTracksStack));
1728                                         return ntracks;
1729                                 }
1730                                 cseed += 6;
1731                         }
1732                 }
1733         }
1734         for(int isl=0; isl<4; isl++) delete layer[isl];
1735         
1736         return ntracks;
1737 }
1738
1739 //_____________________________________________________________________________
1740 AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params)
1741 {
1742   //
1743   // Build a TRD track out of tracklet candidates
1744   //
1745   // Parameters :
1746   //   seeds  : array of tracklets
1747   //   params : track parameters (see MakeSeeds() function body for a detailed description)
1748   //
1749   // Output :
1750   //   The TRD track.
1751   //
1752   // Detailed description
1753   //
1754   // To be discussed with Marian !!
1755   //
1756
1757   Double_t alpha = AliTRDgeometry::GetAlpha();
1758   Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
1759   Double_t c[15];
1760
1761   c[ 0] = 0.2;
1762   c[ 1] = 0.0; c[ 2] = 2.0;
1763   c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
1764   c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0;  c[ 9] = 0.1;
1765   c[10] = 0.0; c[11] = 0.0; c[12] = 0.0;  c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
1766
1767   AliTRDtrackV1 *track = new AliTRDtrackV1(seeds, &params[1], c, params[0], params[6]*alpha+shift);
1768         track->PropagateTo(params[0]-5.0);
1769   track->ResetCovariance(1);
1770   Int_t nc = FollowBackProlongation(*track);
1771         AliInfo(Form("N clusters for track %d", nc));
1772         if (nc < 30) {
1773     delete track;
1774     track = 0x0;
1775   } else {
1776 //     track->CookdEdx();
1777 //     track->CookdEdxTimBin(-1);
1778 //     CookLabel(track, 0.9);
1779   }
1780
1781   return track;
1782 }
1783
1784 //____________________________________________________________________
1785 void AliTRDtrackerV1::CookLabel(AliKalmanTrack */*pt*/, Float_t /*wrong*/) const
1786 {
1787         // to be implemented, preferably at the level of TRD tracklet. !!!!!!!
1788 }
1789
1790 //____________________________________________________________________
1791 void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers
1792                                        , AliTRDseedV1 *cseed)
1793 {
1794   //
1795   // Sort tracklets according to "quality" and try to "improve" the first 4 worst
1796   //
1797   // Parameters :
1798   //  layers : Array of propagation layers for a stack/supermodule
1799   //  cseed  : Array of 6 seeding tracklets which has to be improved
1800   // 
1801   // Output :
1802   //   cssed : Improved seeds
1803   // 
1804   // Detailed description
1805   //
1806   // Iterative procedure in which new clusters are searched for each
1807   // tracklet seed such that the seed quality (see AliTRDseed::GetQuality())
1808   // can be maximized. If some optimization is found the old seeds are replaced.
1809   //
1810         
1811         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1812         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1813         
1814         // make a local working copy
1815         AliTRDseedV1 bseed[6];
1816         for (Int_t jLayer = 0; jLayer < 6; jLayer++) bseed[jLayer] = cseed[jLayer];
1817
1818
1819         Float_t lastquality = 10000.0;
1820         Float_t lastchi2    = 10000.0;
1821         Float_t chi2        =  1000.0;
1822
1823         for (Int_t iter = 0; iter < 4; iter++) {
1824                 Float_t sumquality = 0.0;
1825                 Float_t squality[6];
1826                 Int_t   sortindexes[6];
1827
1828                 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1829                         squality[jLayer]  = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : -1.;
1830                         sumquality +=squality[jLayer];
1831                 }
1832                 if ((sumquality >= lastquality) || (chi2       >     lastchi2)) break;
1833
1834                 
1835                 lastquality = sumquality;
1836                 lastchi2    = chi2;
1837                 if (iter > 0) for (Int_t jLayer = 0; jLayer < 6; jLayer++) cseed[jLayer] = bseed[jLayer];
1838
1839                 
1840                 TMath::Sort(6, squality, sortindexes, kFALSE);
1841                 for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
1842                         Int_t bLayer = sortindexes[jLayer];
1843                         bseed[bLayer].AttachClustersIter(&layers[bLayer*nTimeBins], squality[bLayer], kTRUE);
1844                 }
1845
1846                 chi2 = AliTRDseedV1::FitRiemanTilt(bseed,kTRUE);
1847         } // Loop: iter
1848 }
1849
1850 //____________________________________________________________________
1851 Double_t  AliTRDtrackerV1::MakeSeedingPlanes(AliTRDstackLayer *layers)
1852 {
1853   //
1854   // Calculate plane quality for seeding.
1855   // 
1856   //
1857   // Parameters :
1858   //   layers : Array of propagation layers for this plane.
1859   //
1860   // Output :
1861   //   plane quality factor for seeding
1862   // 
1863   // Detailed description
1864   //
1865   // The quality of the plane for seeding is higher if:
1866   //  1. the average timebin population is closer to an integer number
1867   //  2. the distribution of clusters/timebin is closer to a uniform distribution.
1868   //    - the slope of the first derivative of a parabolic fit is small or
1869   //    - the slope of a linear fit is small
1870   //
1871
1872         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1873         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1874
1875 //      Double_t x;
1876 //      TLinearFitter fitter(1, "pol1");
1877 //      fitter.ClearPoints();
1878         Int_t ncl = 0;
1879         Int_t nused = 0;
1880         Int_t nClLayer;
1881         for(int itb=0; itb<nTimeBins; itb++){
1882                 //x = layer[itb].GetX();
1883                 //printf("x[%d] = %f nCls %d\n", itb, x, layer[itb].GetNClusters());
1884                 //if(!layer[itb].GetNClusters()) continue;
1885                 //fitter.AddPoint(&x, layer[itb].GetNClusters(), 1.);
1886                 nClLayer = layers[itb].GetNClusters();
1887                 ncl += nClLayer;
1888                 for(Int_t incl = 0; incl < nClLayer; incl++)
1889                         if((layers[itb].GetCluster(incl))->IsUsed()) nused++;
1890         }
1891         
1892         // calculate the deviation of the mean number of clusters from the
1893         // closest integer values
1894         Float_t nclMed = float(ncl-nused)/nTimeBins;
1895         Int_t ncli = Int_t(nclMed);
1896         Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
1897         nclDev -= (nclDev>.5) && ncli ? 0. : 1.; 
1898         /*Double_t quality = */ return TMath::Exp(2.*nclDev);
1899
1900 //      // get slope of the derivative
1901 //      if(!fitter.Eval()) return quality;
1902 //      fitter.PrintResults(3);
1903 //      Double_t a = fitter.GetParameter(1);
1904 // 
1905 //      printf("ncl_dev(%f)  a(%f)\n", ncl_dev, a);
1906 //      return quality*TMath::Exp(-a);
1907 }
1908
1909 //____________________________________________________________________
1910 Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed
1911                                        , Int_t planes[4]
1912                                        , Double_t *chi2)
1913 {
1914   //
1915   // Calculate the probability of this track candidate.
1916   //
1917   // Parameters :
1918   //   cseeds : array of candidate tracklets
1919   //   planes : array of seeding planes (see seeding configuration)
1920   //   chi2   : chi2 values (on the Z and Y direction) from the rieman fit of the track.
1921   //
1922   // Output :
1923   //   likelihood value
1924   // 
1925   // Detailed description
1926   //
1927   // The track quality is estimated based on the following 4 criteria:
1928   //  1. precision of the rieman fit on the Y direction (likea)
1929   //  2. chi2 on the Y direction (likechi2y)
1930   //  3. chi2 on the Z direction (likechi2z)
1931   //  4. number of attached clusters compared to a reference value 
1932   //     (see AliTRDrecoParam::fkFindable) (likeN)
1933   //
1934   // The distributions for each type of probabilities are given below as of
1935   // (date). They have to be checked to assure consistency of estimation.
1936   //
1937  
1938         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1939         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1940         // ratio of the total number of clusters/track which are expected to be found by the tracker.
1941         Float_t fgFindable = fRecoParam->GetFindableClusters();
1942
1943         
1944         Int_t nclusters = 0;
1945         Double_t sumda = 0.;
1946         for(UChar_t ilayer = 0; ilayer < 4; ilayer++){
1947                 Int_t jlayer = planes[ilayer];
1948                 nclusters += cseed[jlayer].GetN2();
1949                 sumda += TMath::Abs(cseed[jlayer].GetYfitR(1) - cseed[jlayer].GetYref(1));
1950         }
1951         Double_t likea     = TMath::Exp(-sumda*10.6);
1952         Double_t likechi2y  = 0.0000000001;
1953         if (chi2[1] < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2[1]) * 7.73);
1954         Double_t likechi2z = TMath::Exp(-chi2[0] * 0.088) / TMath::Exp(-chi2[0] * 0.019);
1955         Int_t enc = Int_t(fgFindable*4.*nTimeBins);     // Expected Number Of Clusters, normally 72
1956         Double_t likeN     = TMath::Exp(-(enc - nclusters) * 0.19);
1957         
1958         Double_t like      = likea * likechi2y * likechi2z * likeN;
1959
1960 #ifdef DEBUG
1961         //AliInfo(Form("sumda(%f) chi2[0](%f) chi2[1](%f) likea(%f) likechi2y(%f) likechi2z(%f) nclusters(%d) likeN(%f)", sumda, chi2[0], chi2[1], likea, likechi2y, likechi2z, nclusters, likeN));
1962         if(AliTRDReconstructor::StreamLevel() >= 2){
1963                 TTreeSRedirector &cstreamer = *fDebugStreamer;
1964                 cstreamer << "CookLikelihood"
1965                         << "sumda="     << sumda
1966                         << "chi0="      << chi2[0]
1967                         << "chi1="      << chi2[1]
1968                         << "likea="     << likea
1969                         << "likechi2y=" << likechi2y
1970                         << "likechi2z=" << likechi2z
1971                         << "nclusters=" << nclusters
1972                         << "likeN="     << likeN
1973                         << "like="      << like
1974                         << "\n";
1975         }
1976 #endif
1977
1978         return like;
1979 }
1980
1981 //___________________________________________________________________
1982 void AliTRDtrackerV1::GetMeanCLStack(AliTRDstackLayer *layers
1983                                    , Int_t *planes
1984                                    , Double_t *params)
1985 {
1986   //
1987   // Determines the Mean number of clusters per layer.
1988   // Needed to determine good Seeding Layers
1989   //
1990   // Parameters:
1991   //    - Array of AliTRDstackLayers
1992   //    - Container for the params
1993   //
1994   // Detailed description
1995   //
1996   // Two Iterations:
1997   // In the first Iteration the mean is calculted using all layers.
1998   // After this, all layers outside the 1-sigma-region are rejected.
1999   // Then the mean value and the standard-deviation are calculted a second
2000   // time in order to select all layers in the 1-sigma-region as good-candidates.
2001   //
2002
2003         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2004         Int_t nTimeBins = cal->GetNumberOfTimeBins();
2005         
2006         Float_t mean = 0, stdev = 0;
2007         Double_t ncl[kNTimeBins*kNSeedPlanes], mcl[kNTimeBins*kNSeedPlanes];
2008         Int_t position = 0;
2009         memset(ncl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
2010         memset(mcl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
2011         Int_t nused = 0;
2012         for(Int_t ipl = 0; ipl < kNSeedPlanes; ipl++){
2013                 for(Int_t ils = 0; ils < nTimeBins; ils++){
2014                         position = planes[ipl]*nTimeBins + ils;
2015                         ncl[ipl * nTimeBins + ils] = layers[position].GetNClusters();
2016                         nused = 0;
2017                         for(Int_t icl = 0; icl < ncl[ipl * nTimeBins + ils]; icl++)
2018                                 if((layers[position].GetCluster(icl))->IsUsed()) nused++;
2019                         ncl[ipl * nTimeBins + ils] -= nused;
2020                 }
2021         }
2022         // Declaration of quartils:
2023         //Double_t qvals[3] = {0.0, 0.0, 0.0};
2024         //Double_t qprop[3] = {0.16667, 0.5, 0.83333};
2025         // Iterations
2026         Int_t counter;
2027         Double_t *array;
2028         Int_t *limit;
2029         Int_t nLayers = nTimeBins * kNSeedPlanes;
2030         for(Int_t iter = 0; iter < 2; iter++){
2031                 array = (iter == 0) ? &ncl[0] : &mcl[0];
2032                 limit = (iter == 0) ? &nLayers : &counter;
2033                 counter = 0;
2034                 if(iter == 1){
2035                         for(Int_t i = 0; i < nTimeBins *kNSeedPlanes; i++){
2036                                 if((ncl[i] >  mean + stdev) || (ncl[i] <  mean - stdev)) continue; // Outside 1-sigma region
2037 //                              if((ncl[i] >  qvals[2]) || (ncl[i] <  qvals[0])) continue; // Outside 1-sigma region
2038                                 if(ncl[i] == 0) continue;                                                // 0-Layers also rejected
2039                                 mcl[counter] = ncl[i];
2040                                 counter++;
2041                         }
2042                 }
2043                 if(*limit == 0) break;
2044                 printf("Limit = %d\n", *limit);
2045                 //using quartils instead of mean and RMS 
2046 //              TMath::Quantiles(*limit,3,array,qvals,qprop,kFALSE);
2047                 mean = TMath::Median(*limit, array, 0x0);
2048                 stdev  = TMath::RMS(*limit, array);
2049         }
2050 //      printf("Quantiles: 0.16667 = %3.3f, 0.5 = %3.3f, 0.83333 = %3.3f\n", qvals[0],qvals[1],qvals[2]);
2051 //      memcpy(params,qvals,sizeof(Double_t)*3);
2052         params[1] = (Double_t)TMath::Nint(mean);
2053         params[0] = (Double_t)TMath::Nint(mean - stdev);
2054         params[2] = (Double_t)TMath::Nint(mean + stdev);
2055
2056 }
2057
2058 //___________________________________________________________________
2059 Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
2060                                       , Double_t *params)
2061 {
2062   //
2063   // Algorithm to find optimal seeding layer
2064   // Layers inside one sigma region (given by Quantiles) are sorted
2065   // according to their difference.
2066   // All layers outside are sorted according t
2067   //
2068   // Parameters:
2069   //     - Array of AliTRDstackLayers (in the current plane !!!)
2070   //     - Container for the Indices of the seeding Layer candidates
2071   //
2072   // Output:
2073   //     - Number of Layers inside the 1-sigma-region
2074   //
2075   // The optimal seeding layer should contain the mean number of
2076   // custers in the layers in one chamber.
2077   //
2078
2079         //printf("Params: %3.3f, %3.3f, %3.3f\n", params[0], params[1], params[2]);
2080         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2081         const Int_t kMaxClustersLayer = AliTRDstackLayer::kMaxClustersLayer;
2082         Int_t nTimeBins = cal->GetNumberOfTimeBins();
2083         Int_t ncl[kNTimeBins], indices[kNTimeBins], bins[kMaxClustersLayer];
2084         memset(ncl, 0, sizeof(Int_t)*kNTimeBins);
2085         memset(indices, 0, sizeof(Int_t)*kNTimeBins);
2086         memset(bins, 0, sizeof(Int_t)*kMaxClustersLayer);
2087         Int_t nused = 0;
2088         for(Int_t ils = 0; ils < nTimeBins; ils++){
2089                 ncl[ils] = layers[ils].GetNClusters();
2090                 nused = 0;
2091                 for(Int_t icl = 0; icl < ncl[ils]; icl++)
2092                         if((layers[ils].GetCluster(icl))->IsUsed()) nused++;
2093                 ncl[ils] -= nused;
2094         }
2095         
2096         Float_t mean = params[1];
2097         for(Int_t ils = 0; ils < nTimeBins; ils++){
2098                 memmove(indices + bins[ncl[ils]+1] + 1, indices + bins[ncl[ils]+1], sizeof(Int_t)*(nTimeBins - ils));
2099                 indices[bins[ncl[ils]+1]] = ils;
2100                 for(Int_t i = ncl[ils]+1; i < kMaxClustersLayer; i++)
2101                         bins[i]++;
2102         }
2103         
2104         //for(Int_t i = 0; i < nTimeBins; i++) printf("Bin %d = %d\n", i, bins[i]);
2105         Int_t sbin = -1;
2106         Int_t nElements;
2107         Int_t position = 0;
2108         TRandom *r = new TRandom();
2109         Int_t iter = 0;
2110         while(1){
2111                 while(sbin < (Int_t)params[0] || sbin > (Int_t)params[2]){
2112                         // Randomly selecting one bin
2113                         sbin = (Int_t)r->Poisson(mean);
2114                 }
2115                 printf("Bin = %d\n",sbin);
2116                 //Randomly selecting one Layer in the bin
2117                 nElements = bins[sbin + 1] - bins[sbin];
2118                 printf("nElements = %d\n", nElements);
2119                 if(iter == 5){
2120                         position = (Int_t)(gRandom->Rndm()*(nTimeBins-1));
2121                         break;
2122                 }
2123                 else if(nElements==0){
2124                         iter++;
2125                         continue;
2126                 }
2127                 position = (Int_t)(gRandom->Rndm()*(nElements-1)) + bins[sbin];
2128                 break;
2129         }
2130         delete r;
2131         return indices[position];
2132 }
2133
2134 //____________________________________________________________________
2135 AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDstackLayer *layers
2136                                                  , AliTRDseedV1/*AliRieman*/ *reference)
2137 {
2138   //
2139   // Finds a seeding Cluster for the extrapolation chamber.
2140   //
2141   // The seeding cluster should be as close as possible to the assumed
2142   // track which is represented by a Rieman fit.
2143   // Therefore the selecting criterion is the minimum distance between
2144   // the best fitting cluster and the Reference which is derived from
2145   // the AliTRDseed. Because all layers are assumed to be equally good
2146   // a linear search is performed.
2147   //
2148   // Imput parameters: - layers: array of AliTRDstackLayers (in one chamber!!!)
2149   //                   - sfit: the reference
2150   //
2151   // Output:           - the best seeding cluster
2152   //
2153
2154         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2155         Int_t nTimeBins = cal->GetNumberOfTimeBins();
2156         
2157         // distances as squared distances
2158         Int_t index = 0;
2159         Float_t ypos = 0.0, zpos = 0.0, distance = 0.0, nearestDistance =100000.0; 
2160         ypos = reference->GetYref(0);
2161         zpos = reference->GetZref(0);
2162         AliTRDcluster *currentBest = 0x0, *temp = 0x0;
2163         for(Int_t ils = 0; ils < nTimeBins; ils++){
2164                 // Reference positions
2165 //              ypos = reference->GetYat(layers[ils].GetX());
2166 //              zpos = reference->GetZat(layers[ils].GetX());
2167                 index = layers[ils].SearchNearestCluster(ypos, zpos, fRecoParam->GetRoad2y(), fRecoParam->GetRoad2z());
2168                 if(index == -1) continue;
2169                 temp = layers[ils].GetCluster(index);
2170                 if(!temp) continue;
2171                 distance = (temp->GetY() - ypos) * (temp->GetY() - ypos) + (temp->GetZ() - zpos) * (temp->GetZ() - zpos);
2172                 if(distance < nearestDistance){
2173                         nearestDistance = distance;
2174                         currentBest = temp;
2175                 }
2176         }
2177         return currentBest;
2178 }
2179
2180 //____________________________________________________________________
2181 AliTRDstackLayer *AliTRDtrackerV1::MakeSeedingLayer(AliTRDstackLayer *layers
2182                                                   , Int_t plane)
2183 {
2184   //
2185   // Creates a seeding layer
2186   //
2187         
2188         // constants
2189         const Int_t kMaxRows = 16;
2190         const Int_t kMaxCols = 144;
2191         const Int_t kMaxPads = 2304;
2192         
2193         // Get the calculation
2194         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2195         Int_t nTimeBins = cal->GetNumberOfTimeBins();
2196         
2197         // Get the geometrical data of the chamber
2198         AliTRDpadPlane *pp = fGeom->GetPadPlane(plane, layers[0].GetStackNr());
2199         Int_t nCols = pp->GetNcols();
2200         Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
2201         Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
2202         Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
2203         Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
2204         Int_t nRows = pp->GetNrows();
2205         Float_t binlength = (ymax - ymin)/nCols; 
2206         //AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
2207         
2208         // Fill the histogram
2209         Int_t arrpos;
2210         Float_t ypos;
2211         Int_t irow, nClusters;  
2212         Int_t *histogram[kMaxRows];                                                                                     // 2D-Histogram
2213         Int_t hvals[kMaxPads];  memset(hvals, 0, sizeof(Int_t)*kMaxPads);       
2214         Float_t *sigmas[kMaxRows];
2215         Float_t svals[kMaxPads];        memset(svals, 0, sizeof(Float_t)*kMaxPads);     
2216         AliTRDcluster *c = 0x0;
2217         for(Int_t irs = 0; irs < kMaxRows; irs++){
2218                 histogram[irs] = &hvals[irs*kMaxCols];
2219                 sigmas[irs] = &svals[irs*kMaxCols];
2220         }
2221         for(Int_t iTime = 0; iTime < nTimeBins; iTime++){
2222                 nClusters = layers[iTime].GetNClusters();
2223                 for(Int_t incl = 0; incl < nClusters; incl++){
2224                         c = layers[iTime].GetCluster(incl);     
2225                         ypos = c->GetY();
2226                         if(ypos > ymax && ypos < ymin) continue;
2227                         irow = pp->GetPadRowNumber(c->GetZ());                          // Zbin
2228                         if(irow < 0)continue;
2229                         arrpos = static_cast<Int_t>((ypos - ymin)/binlength);
2230                         if(ypos == ymax) arrpos = nCols - 1;
2231                         histogram[irow][arrpos]++;
2232                         sigmas[irow][arrpos] += c->GetSigmaZ2();
2233                 }
2234         }
2235         
2236 // Now I have everything in the histogram, do the selection
2237 //      printf("Starting the analysis\n");
2238         //Int_t nPads = nCols * nRows;
2239         // This is what we are interested in: The center of gravity of the best candidates
2240         Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
2241         Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
2242         Float_t *cogy[kMaxRows];
2243         Float_t *cogz[kMaxRows];
2244         // Lookup-Table storing coordinates according ti the bins
2245         Float_t yLengths[kMaxCols];
2246         Float_t zLengths[kMaxRows];
2247         for(Int_t icnt = 0; icnt < nCols; icnt++){
2248                 yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
2249         }
2250         for(Int_t icnt = 0; icnt < nRows; icnt++){
2251                 zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
2252         }
2253
2254         // A bitfield is used to mask the pads as usable
2255         Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
2256         for(UChar_t icount = 0; icount < nRows; icount++){
2257                 cogy[icount] = &cogyvals[icount*kMaxCols];
2258                 cogz[icount] = &cogzvals[icount*kMaxCols];
2259         }
2260         // In this array the array position of the best candidates will be stored
2261         Int_t cand[kMaxTracksStack];
2262         Float_t sigcands[kMaxTracksStack];
2263         
2264         // helper variables
2265         Int_t indices[kMaxPads]; memset(indices, 0, sizeof(Int_t)*kMaxPads);
2266         Int_t nCandidates = 0;
2267         Float_t norm, cogv;
2268         // histogram filled -> Select best bins
2269         TMath::Sort(kMaxPads, hvals, indices);                  // bins storing a 0 should not matter
2270         // Set Threshold
2271         Int_t maximum = hvals[indices[0]];      // best
2272         Int_t threshold = static_cast<UChar_t>(maximum * fRecoParam->GetFindableClusters());
2273         Int_t col, row, lower, lower1, upper, upper1;
2274         for(Int_t ib = 0; ib < kMaxPads; ib++){
2275                 if(nCandidates >= kMaxTracksStack){
2276                         AliWarning(Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, kMaxTracksStack));
2277                         break;
2278                 }
2279                 // Positions
2280                 row = indices[ib]/nCols;
2281                 col = indices[ib]%nCols;
2282                 // here will be the threshold condition:
2283                 if((mask[col] & (1 << row)) != 0) continue;             // Pad is masked: continue
2284                 if(histogram[row][col] < TMath::Max(threshold, 1)){     // of course at least one cluster is needed
2285                         break;                  // number of clusters below threshold: break;
2286                 } 
2287                 // passing: Mark the neighbors
2288                 lower  = TMath::Max(col - 1, 0); upper  = TMath::Min(col + 2, nCols);
2289                 lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
2290                 for(Int_t ic = lower; ic < upper; ++ic)
2291                         for(Int_t ir = lower1; ir < upper1; ++ir){
2292                                 if(ic == col && ir == row) continue;
2293                                 mask[ic] |= (1 << ir);
2294                         }
2295                 // Storing the position in an array
2296                 // testing for neigboring
2297                 cogv = 0;
2298                 norm = 0;
2299                 lower = TMath::Max(col - 1,0);
2300                 upper = TMath::Min(col + 2, nCols);
2301                 for(Int_t inb = lower; inb < upper; ++inb){
2302                         cogv += yLengths[inb] * histogram[row][inb];
2303                         norm += histogram[row][inb];
2304                 }
2305                 cogy[row][col] = cogv / norm;
2306                 cogv = 0; norm = 0;
2307                 lower = TMath::Max(row - 1, 0);
2308                 upper = TMath::Min(row + 2, nRows);
2309                 for(Int_t inb = lower; inb < upper; ++inb){
2310                         cogv += zLengths[inb] * histogram[inb][col];
2311                         norm += histogram[inb][col];
2312                 }
2313                 cogz[row][col] = cogv /  norm;
2314                 // passed the filter
2315                 cand[nCandidates] = row*kMaxCols + col; // store the position of a passig candidate into an Array
2316                 sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
2317                 // Analysis output
2318                 nCandidates++;
2319         }
2320         AliTRDstackLayer *fakeLayer = new AliTRDstackLayer(layers[0].GetZ0(), layers[0].GetDZ0(), layers[0].GetStackNr());
2321         fakeLayer->SetX((TMath::Abs(layers[nTimeBins-1].GetX() + layers[0].GetX()))/2);
2322         fakeLayer->SetSector(layers[0].GetSector());
2323         AliTRDcluster **fakeClusters = 0x0;
2324         UInt_t *fakeIndices = 0x0;
2325         if(nCandidates){
2326                 fakeClusters = new AliTRDcluster*[nCandidates];
2327                 fakeIndices = new UInt_t[nCandidates];
2328                 UInt_t fakeIndex = 0;
2329                 for(Int_t ican = 0; ican < nCandidates; ican++){
2330                         fakeClusters[ican] = new AliTRDcluster();
2331                         fakeClusters[ican]->SetX(fakeLayer->GetX());
2332                         fakeClusters[ican]->SetY(cogyvals[cand[ican]]);
2333                         fakeClusters[ican]->SetZ(cogzvals[cand[ican]]);
2334                         fakeClusters[ican]->SetSigmaZ2(sigcands[ican]);
2335                         fakeIndices[ican] = fakeIndex++;// fantasy number
2336                 }
2337         }
2338         fakeLayer->SetRecoParam(fRecoParam);
2339         fakeLayer->SetClustersArray(fakeClusters, nCandidates);
2340         fakeLayer->SetIndexArray(fakeIndices);
2341         fakeLayer->SetNRows(nRows);
2342         fakeLayer->BuildIndices();
2343         //fakeLayer->PrintClusters();
2344         
2345 #ifdef DEBUG
2346         if(AliTRDReconstructor::StreamLevel() >= 3){
2347                 TMatrixD hist(nRows, nCols);
2348                 for(Int_t i = 0; i < nRows; i++)
2349                         for(Int_t j = 0; j < nCols; j++)
2350                                 hist(i,j) = histogram[i][j];
2351                 TTreeSRedirector &cstreamer = *fDebugStreamer;
2352                 cstreamer << "MakeSeedingLayer"
2353                         << "Iteration="  << fSieveSeeding
2354                         << "plane="      << plane
2355                         << "ymin="       << ymin
2356                         << "ymax="       << ymax
2357                         << "zmin="       << zmin
2358                         << "zmax="       << zmax
2359                         << "L.="         << fakeLayer
2360                         << "Histogram.=" << &hist
2361                         << "\n";
2362         }
2363 #endif
2364         return fakeLayer;
2365 }
2366
2367 //____________________________________________________________________
2368 void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig, Int_t planes[4])
2369 {
2370   //
2371   // Map seeding configurations to detector planes.
2372   //
2373   // Parameters :
2374   //   iconfig : configuration index
2375   //   planes  : member planes of this configuration. On input empty.
2376   //
2377   // Output :
2378   //   planes : contains the planes which are defining the configuration
2379   // 
2380   // Detailed description
2381   //
2382   // Here is the list of seeding planes configurations together with
2383   // their topological classification:
2384   //
2385   //  0 - 5432 TQ 0
2386   //  1 - 4321 TQ 0
2387   //  2 - 3210 TQ 0
2388   //  3 - 5321 TQ 1
2389   //  4 - 4210 TQ 1
2390   //  5 - 5431 TQ 1
2391   //  6 - 4320 TQ 1
2392   //  7 - 5430 TQ 2
2393   //  8 - 5210 TQ 2
2394   //  9 - 5421 TQ 3
2395   // 10 - 4310 TQ 3
2396   // 11 - 5410 TQ 4
2397   // 12 - 5420 TQ 5
2398   // 13 - 5320 TQ 5
2399   // 14 - 5310 TQ 5
2400   //
2401   // The topologic quality is modeled as follows:
2402   // 1. The general model is define by the equation:
2403   //  p(conf) = exp(-conf/2)
2404   // 2. According to the topologic classification, configurations from the same
2405   //    class are assigned the agerage value over the model values.
2406   // 3. Quality values are normalized.
2407   // 
2408   // The topologic quality distribution as function of configuration is given below:
2409   //Begin_Html
2410   // <img src="gif/topologicQA.gif">
2411   //End_Html
2412   //
2413
2414         switch(iconfig){
2415         case 0: // 5432 TQ 0
2416                 planes[0] = 2;
2417                 planes[1] = 3;
2418                 planes[2] = 4;
2419                 planes[3] = 5;
2420                 break;
2421         case 1: // 4321 TQ 0
2422                 planes[0] = 1;
2423                 planes[1] = 2;
2424                 planes[2] = 3;
2425                 planes[3] = 4;
2426                 break;
2427         case 2: // 3210 TQ 0
2428                 planes[0] = 0;
2429                 planes[1] = 1;
2430                 planes[2] = 2;
2431                 planes[3] = 3;
2432                 break;
2433         case 3: // 5321 TQ 1
2434                 planes[0] = 1;
2435                 planes[1] = 2;
2436                 planes[2] = 3;
2437                 planes[3] = 5;
2438                 break;
2439         case 4: // 4210 TQ 1
2440                 planes[0] = 0;
2441                 planes[1] = 1;
2442                 planes[2] = 2;
2443                 planes[3] = 4;
2444                 break;
2445         case 5: // 5431 TQ 1
2446                 planes[0] = 1;
2447                 planes[1] = 3;
2448                 planes[2] = 4;
2449                 planes[3] = 5;
2450                 break;
2451         case 6: // 4320 TQ 1
2452                 planes[0] = 0;
2453                 planes[1] = 2;
2454                 planes[2] = 3;
2455                 planes[3] = 4;
2456                 break;
2457         case 7: // 5430 TQ 2
2458                 planes[0] = 0;
2459                 planes[1] = 3;
2460                 planes[2] = 4;
2461                 planes[3] = 5;
2462                 break;
2463         case 8: // 5210 TQ 2
2464                 planes[0] = 0;
2465                 planes[1] = 1;
2466                 planes[2] = 2;
2467                 planes[3] = 5;
2468                 break;
2469         case 9: // 5421 TQ 3
2470                 planes[0] = 1;
2471                 planes[1] = 2;
2472                 planes[2] = 4;
2473                 planes[3] = 5;
2474                 break;
2475         case 10: // 4310 TQ 3
2476                 planes[0] = 0;
2477                 planes[1] = 1;
2478                 planes[2] = 3;
2479                 planes[3] = 4;
2480                 break;
2481         case 11: // 5410 TQ 4
2482                 planes[0] = 0;
2483                 planes[1] = 1;
2484                 planes[2] = 4;
2485                 planes[3] = 5;
2486                 break;
2487         case 12: // 5420 TQ 5
2488                 planes[0] = 0;
2489                 planes[1] = 2;
2490                 planes[2] = 4;
2491                 planes[3] = 5;
2492                 break;
2493         case 13: // 5320 TQ 5
2494                 planes[0] = 0;
2495                 planes[1] = 2;
2496                 planes[2] = 3;
2497                 planes[3] = 5;
2498                 break;
2499         case 14: // 5310 TQ 5
2500                 planes[0] = 0;
2501                 planes[1] = 1;
2502                 planes[2] = 3;
2503                 planes[3] = 5;
2504                 break;
2505         }
2506 }
2507
2508 //____________________________________________________________________
2509 void AliTRDtrackerV1::GetExtrapolationConfig(Int_t iconfig, Int_t planes[2])
2510 {
2511   //
2512   // Returns the extrapolation planes for a seeding configuration.
2513   //
2514   // Parameters :
2515   //   iconfig : configuration index
2516   //   planes  : planes which are not in this configuration. On input empty.
2517   //
2518   // Output :
2519   //   planes : contains the planes which are not in the configuration
2520   // 
2521   // Detailed description
2522   //
2523
2524         switch(iconfig){
2525         case 0: // 5432 TQ 0
2526                 planes[0] = 1;
2527                 planes[1] = 0;
2528                 break;
2529         case 1: // 4321 TQ 0
2530                 planes[0] = 5;
2531                 planes[1] = 0;
2532                 break;
2533         case 2: // 3210 TQ 0
2534                 planes[0] = 4;
2535                 planes[1] = 5;
2536                 break;
2537         case 3: // 5321 TQ 1
2538                 planes[0] = 4;
2539                 planes[1] = 0;
2540                 break;
2541         case 4: // 4210 TQ 1
2542                 planes[0] = 5;
2543                 planes[1] = 3;
2544                 break;
2545         case 5: // 5431 TQ 1
2546                 planes[0] = 2;
2547                 planes[1] = 0;
2548                 break;
2549         case 6: // 4320 TQ 1
2550                 planes[0] = 5;
2551                 planes[1] = 1;
2552                 break;
2553         case 7: // 5430 TQ 2
2554                 planes[0] = 2;
2555                 planes[1] = 1;
2556                 break;
2557         case 8: // 5210 TQ 2
2558                 planes[0] = 4;
2559                 planes[1] = 3;
2560                 break;
2561         case 9: // 5421 TQ 3
2562                 planes[0] = 3;
2563                 planes[1] = 0;
2564                 break;
2565         case 10: // 4310 TQ 3
2566                 planes[0] = 5;
2567                 planes[1] = 2;
2568                 break;
2569         case 11: // 5410 TQ 4
2570                 planes[0] = 3;
2571                 planes[1] = 2;
2572                 break;
2573         case 12: // 5420 TQ 5
2574                 planes[0] = 3;
2575                 planes[1] = 1;
2576                 break;
2577         case 13: // 5320 TQ 5
2578                 planes[0] = 4;
2579                 planes[1] = 1;
2580                 break;
2581         case 14: // 5310 TQ 5
2582                 planes[0] = 4;
2583                 planes[1] = 2;
2584                 break;
2585         }
2586 }