e896d954ce8fa725ac67817ad909658b77e35f9a
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUTrackerGlo.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17   TO CHECK
18   GetBz usage
19
20  */
21
22 //-------------------------------------------------------------------------
23 //  Implementation of the ITS Upgrade tracker mother class.              
24 //-------------------------------------------------------------------------
25 #include <TTree.h>
26 #include <Riostream.h> 
27 #include <TMath.h>
28 #include <TFile.h>
29
30 #include "AliITSUTrackerGlo.h"
31 #include "AliESDEvent.h"
32 #include "AliESDtrack.h"
33 #include "AliITSURecoDet.h"
34 #include "AliITSURecoSens.h"
35 #include "AliITSUReconstructor.h"
36 #include "AliITSReconstructor.h"
37 #include "AliITSUSeed.h"
38 #include "AliITSUAux.h"
39 #include "AliITSUClusterPix.h"
40 #include "AliITSUGeomTGeo.h"
41 using namespace AliITSUAux;
42 using namespace TMath;
43
44 //----------------- tmp stuff -----------------
45
46 ClassImp(AliITSUTrackerGlo)
47
48 const Double_t AliITSUTrackerGlo::fgkToler =  1e-6;// tolerance for layer on-surface check
49
50
51 //_________________________________________________________________________
52 AliITSUTrackerGlo::AliITSUTrackerGlo(AliITSUReconstructor* rec)
53 :  fReconstructor(rec)
54   ,fITS(0)
55   ,fCurrESDtrack(0)
56   ,fCurrESDtrMClb(kDummyLabel)
57   ,fCurrMass(kPionMass)
58   ,fCountProlongationTrials(0)
59   ,fCountITSin(0)
60   ,fCountITSout(0)
61   ,fCountITSrefit(0)
62   ,fHypStore(100)
63   ,fLayerMaxCandidates(1000)
64   ,fLayerCandidates(0)
65   ,fNBranchesAdded(0)
66   ,fNCandidatesAdded(0)
67   ,fCurrHyp(0)
68   ,fWorkHyp(0)
69   ,fSeedsPool("AliITSUSeed",0)
70   ,fFreeSeedsID(0)
71   ,fNFreeSeeds(0)
72   ,fLastSeedID(0)
73   ,fNLrActive(0)
74   ,fDefTrackConds(0)
75   ,fCurrTrackCond(0)
76   ,fCurrActLrID(-1)
77   ,fCurrLayer(0)
78   ,fTrackPhase(-1)
79   ,fClInfo(0)
80 #ifdef  _FILL_CONTROL_HISTOS_
81   ,fCHistoArr(0)
82 #endif
83 {
84   // Default constructor
85   if (rec) Init(rec);
86 }
87
88 //_________________________________________________________________________
89 AliITSUTrackerGlo::~AliITSUTrackerGlo()
90 {
91  // Default destructor
92  //  
93   delete[] fClInfo;
94   //
95 #ifdef  _FILL_CONTROL_HISTOS_
96   if (fCHistoArr) {
97     TFile* ctrOut = TFile::Open("itsuControlHistos.root","recreate");
98     ctrOut->cd();
99     AliInfo("Storing control histos");
100     fCHistoArr->Print();
101     //    ctrOut->WriteObject(fCHistoArr,"controlH","kSingleKey");
102     fCHistoArr->Write();
103     ctrOut->Close();
104     delete ctrOut;
105     fCHistoArr = 0;
106   }
107 #endif 
108   //
109 }
110
111 //_________________________________________________________________________
112 void AliITSUTrackerGlo::Init(AliITSUReconstructor* rec)
113 {
114   // init with external reconstructor
115   //
116   fITS = rec->GetITSInterface();
117   fNLrActive = fITS->GetNLayersActive();
118   fClInfo = new Int_t[fNLrActive*2];
119   fWorkHyp = new AliITSUTrackHyp(fNLrActive);
120   //
121   if (fLayerMaxCandidates<1) fLayerMaxCandidates = 1000;
122   fLayerCandidates = new AliITSUSeed*[fLayerMaxCandidates];
123   fSeedsPool.ExpandCreateFast(1000); // RS TOCHECK
124   fFreeSeedsID.Set(1000);
125   //
126 }
127
128 //_________________________________________________________________________
129 void AliITSUTrackerGlo::CreateDefaultTrackCond()
130 {
131   // creates default tracking conditions to be used when recoparam does not provide them
132
133   AliITSUTrackCond* cond = new AliITSUTrackCond();
134   //
135   cond->SetNLayers(fNLrActive);
136   cond->AddNewCondition(fNLrActive);
137   cond->AddGroupPattern( 0xffff ); // require all layers hit
138   cond->Init();
139   //
140   fDefTrackConds.AddLast(cond);
141   //
142   AliInfo("Created conditions: ");
143   for (int i=0;i<fDefTrackConds.GetEntriesFast();i++) fDefTrackConds[i]->Print();
144   //
145 }
146
147
148 //_________________________________________________________________________
149 Int_t AliITSUTrackerGlo::Clusters2Tracks(AliESDEvent *esdEv)
150 {
151   //
152   SetTrackingPhase(kClus2Tracks);
153   //
154 #ifdef  _FILL_CONTROL_HISTOS_
155   if (!fCHistoArr) BookControlHistos();
156 #endif
157   static int evID = 0;
158   TObjArray *trackConds = 0;
159   //
160   fCountProlongationTrials = 0;
161   fCountITSin = 0;
162   fCountITSout = 0;
163   fCountITSrefit = 0;
164   //
165   ResetSeedsPool();
166   int nTrESD = esdEv->GetNumberOfTracks();
167   AliInfo(Form("Will try to find prolongations for %d tracks",nTrESD));
168   int nTrackCond = AliITSUReconstructor::GetRecoParam()->GetNTrackingConditions();
169   if (nTrackCond<1) {
170     if (!fDefTrackConds.GetEntriesFast()) {
171       AliInfo("No tracking conditions found in recoparams, creating default one requesting all layers hit");
172       CreateDefaultTrackCond();
173     }
174     trackConds = &fDefTrackConds;
175     nTrackCond = trackConds->GetEntriesFast();
176   } 
177   else trackConds = AliITSUReconstructor::GetRecoParam()->GetTrackingConditions();
178   //
179   static Bool_t first = kTRUE;
180   if (first) {
181     AliITSUReconstructor::GetRecoParam()->Print();
182     first = kFALSE;
183   }
184   fHypStore.Delete();
185   if (fHypStore.GetSize()<nTrESD) fHypStore.Expand(nTrESD+100);
186   //
187   fITS->ProcessClusters();
188   //
189   FlagSplitClusters(); // tmp RS
190   //
191   for (int icnd=0;icnd<nTrackCond;icnd++) {
192     fCurrTrackCond = (AliITSUTrackCond*)trackConds->UncheckedAt(icnd);
193     if (!fCurrTrackCond->IsInitDone()) fCurrTrackCond->Init();
194     // select ESD tracks to propagate
195     for (int itr=0;itr<nTrESD;itr++) {
196       fCurrESDtrack = esdEv->GetTrack(itr);
197       fCurrESDtrMClb = fCurrESDtrack->GetLabel();
198       //
199       if (!NeedToProlong(fCurrESDtrack)) continue;  // are we interested in this track?
200       /*
201         // if speceific tracks should be checked, create a global array int wtc[] = {evITS*10000+trID, ...};
202       Bool_t dbg = kFALSE;
203       int nwtc = sizeof(wtc)/sizeof(int);
204
205       for (int k=0;k<nwtc;k++) {
206         if (wtc[k]==evID*10000+itr) {
207           dbg = kTRUE;
208           printf("\n\n\n\n\n\n\n");
209           printf("At esdTr: %d MC: %d\n",wtc[k],fCurrESDtrMClb);
210           printf("\n\n\n\n\n\n\n");
211           break;
212         }
213       }
214       AliLog::SetClassDebugLevel("AliITSUTrackerGlo",dbg ? 10:0);
215       */
216
217       AliDebug(-1,Form("Processing track %d | M=%.3f Pt=%.3f | MCLabel: %d",itr,fCurrESDtrack->GetMass(kTRUE),fCurrESDtrack->Pt(),fCurrESDtrMClb));//RS
218       FindTrack(fCurrESDtrack, itr);
219     }   
220     //
221     if (AliDebugLevelClass()>+2) {
222       AliInfo(Form("SeedsPool: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
223       fHypStore.Print();
224     }
225     FinalizeHypotheses();
226     //
227     //AliLog::SetClassDebugLevel("AliITSUTrackerGlo",0);     // in case wtc array was defined, uncomment this
228   }
229   //
230   AliInfo(Form("%d ITSin for %d tried TPC seeds out of %d ESD tracks\n",fCountITSin,fCountProlongationTrials,nTrESD));
231   evID++;
232   return 0;
233 }
234
235 //_________________________________________________________________________
236 Int_t AliITSUTrackerGlo::PropagateBack(AliESDEvent *esdEv)
237 {
238   //
239   // Do outward fits in ITS
240   //
241   SetTrackingPhase(kPropBack);
242   int nTrESD = esdEv->GetNumberOfTracks();
243   AliDebug(1,Form("Will propagate back %d tracks",nTrESD));
244   //
245   double bz0 = GetBz();
246   Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
247   AliITSUTrackHyp dummyTr,*currTr=0;
248   const double kWatchStep=10.; // for larger steps watch arc vs segment difference
249   Double_t times[AliPID::kSPECIES];
250   //
251   //
252   for (int itr=0;itr<nTrESD;itr++) {
253     fCurrESDtrack = esdEv->GetTrack(itr);
254     fCurrESDtrMClb = fCurrESDtrack->GetLabel();
255     // Start time integral and add distance from current position to vertex 
256     if (fCurrESDtrack->IsOn(AliESDtrack::kITSout)) continue; // already done
257     //
258     fCurrESDtrack->GetXYZ(xyzTrk); 
259     Double_t dst = 0.;     // set initial track lenght, tof
260     {
261       double dxs = xyzTrk[0] - xyzVtx[0];
262       double dys = xyzTrk[1] - xyzVtx[1];
263       double dzs = xyzTrk[2] - xyzVtx[2];
264       // RS: for large segment steps d use approximation of cicrular arc s by
265       // s = 2R*asin(d/2R) = d/p asin(p) \approx d/p (p + 1/6 p^3) = d (1+1/6 p^2)
266       // where R is the track radius, p = d/2R = 2C*d (C is the abs curvature)
267       // Hence s^2/d^2 = (1+1/6 p^2)^2
268       dst = dxs*dxs + dys*dys;
269       if (dst > kWatchStep*kWatchStep) { // correct circular part for arc/segment factor
270         double crv = Abs(fCurrESDtrack->GetC(bz0));
271         double fcarc = 1.+crv*crv*dst/6.;
272         dst *= fcarc*fcarc;
273       }
274       dst += dzs*dzs;
275       dst = Sqrt(dst); 
276     }
277     //    
278     fCurrESDtrack->SetStatus(AliESDtrack::kTIME);
279     //
280     if (!fCurrESDtrack->IsOn(AliESDtrack::kITSin)) { // case of tracks w/o ITS prolongation: just set the time integration
281       dummyTr.AliExternalTrackParam::operator=(*fCurrESDtrack);
282       dummyTr.StartTimeIntegral();
283       dummyTr.AddTimeStep(dst);
284       dummyTr.GetIntegratedTimes(times); 
285       fCurrESDtrack->SetIntegratedTimes(times);
286       fCurrESDtrack->SetIntegratedLength(dummyTr.GetIntegratedLength());
287       continue;
288     }
289     //
290     currTr = GetTrackHyp(itr);
291     currTr->StartTimeIntegral();
292     currTr->AddTimeStep(dst);
293     //    printf("Before resetCov: "); currTr->AliExternalTrackParam::Print();
294     currTr->ResetCovariance(10000);
295     if (RefitTrack(currTr,fITS->GetRMax())) { // propagate to exit from the ITS/TPC screen
296       UpdateESDTrack(currTr,AliESDtrack::kITSout);
297       fCountITSout++;
298     }
299     else {
300       AliDebug(2,Form("Refit Failed for track %d | ESDtrack#%d (MClb:%d)",itr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
301       //currTr->AliExternalTrackParam::Print();
302       //currTr->GetWinner()->Print();
303     }
304     //
305   }
306   //
307   AliInfo(Form("%d ITSout in %d ITSin tracks for %d tried TPC seeds out of %d ESD tracks\n",
308                fCountITSout,fCountITSin,fCountProlongationTrials,nTrESD));
309   //
310   return 0;
311 }
312
313 //_________________________________________________________________________
314 Int_t AliITSUTrackerGlo::RefitInward(AliESDEvent *esdEv)
315 {
316   //
317   // refit the tracks inward, using current cov.matrix
318   //
319   SetTrackingPhase(kRefitInw);
320   Int_t nTrESD = esdEv->GetNumberOfTracks();
321   //  AliLog::SetClassDebugLevel("AliITSUTrackerGlo",10);
322
323   AliDebug(1,Form("Will refit inward %d tracks",nTrESD));
324   AliITSUTrackHyp *currTr=0;
325   //
326   for (int itr=0;itr<nTrESD;itr++) {
327     fCurrESDtrack = esdEv->GetTrack(itr);
328     fCurrESDtrMClb = fCurrESDtrack->GetLabel();
329     // Start time integral and add distance from current position to vertex 
330     UInt_t trStat = fCurrESDtrack->GetStatus();
331     if ( !(trStat & AliESDtrack::kITSout) ) continue;
332     if (   trStat & AliESDtrack::kITSrefit ) continue; // already done
333     if (  (trStat & AliESDtrack::kTPCout) && !(trStat & AliESDtrack::kTPCrefit) ) continue;
334     //
335     currTr = GetTrackHyp(itr);
336     currTr->AliExternalTrackParam::operator=(*fCurrESDtrack);  // fetch current ESDtrack kinematics
337     if (RefitTrack(currTr,fITS->GetRMin())) { // propagate up to inside radius of the beam pipe
338       UpdateESDTrack(currTr,AliESDtrack::kITSrefit);
339       fCountITSrefit++;
340     }
341     else {
342       AliDebug(2,Form("Refit Failed for track %d |ESDtrack#%d (MClb:%d)",itr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
343     }
344   }    
345   //
346   AliInfo(Form("%d ITSrefit in %d ITSout in %d ITSin tracks for %d tried TPC seeds out of %d ESD tracks\n",
347                fCountITSrefit,fCountITSout,fCountITSin,fCountProlongationTrials,nTrESD));
348   //
349   //  AliLog::SetClassDebugLevel("AliITSUTrackerGlo",0);
350   return 0;
351 }
352
353 //_________________________________________________________________________
354 Int_t AliITSUTrackerGlo::LoadClusters(TTree * treeRP)
355 {
356   // read from tree (if pointer provided) or directly from the ITS reco interface
357   //
358   return fReconstructor->LoadClusters(treeRP);
359
360
361 //_________________________________________________________________________
362 void AliITSUTrackerGlo::UnloadClusters()
363 {
364   //
365   // To be implemented 
366   //
367   
368   Info("UnloadClusters","To be implemented");
369
370 //_________________________________________________________________________
371 AliCluster * AliITSUTrackerGlo::GetCluster(Int_t /*index*/) const
372 {
373   //
374   // To be implemented 
375   //
376   
377   Info("GetCluster","To be implemented");
378   return 0x0;
379
380
381 //_________________________________________________________________________
382 Bool_t AliITSUTrackerGlo::NeedToProlong(AliESDtrack* esdTr)
383 {
384   // do we need to match this track to ITS?
385   //
386   static double bz = GetBz();
387   if (!esdTr->IsOn(AliESDtrack::kTPCin) ||
388       esdTr->IsOn(AliESDtrack::kTPCout) ||
389       esdTr->IsOn(AliESDtrack::kITSin)  ||
390       esdTr->GetKinkIndex(0)>0) return kFALSE;
391   //
392   if (esdTr->Pt()<AliITSUReconstructor::GetRecoParam()->GetMinPtForProlongation()) return kFALSE;
393   //
394   float dtz[2];
395   esdTr->GetDZ(GetX(),GetY(),GetZ(),bz,dtz); 
396   // if track is not V0 candidata but has large offset wrt IP, reject it. RS TOCHECK
397   if ( !(esdTr->GetV0Index(0)>0 && dtz[0]>AliITSUReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()) 
398        && (Abs(dtz[0])>AliITSUReconstructor::GetRecoParam()->GetMaxDForProlongation() ||
399            Abs(dtz[1])>AliITSUReconstructor::GetRecoParam()->GetMaxDZForProlongation())) return kFALSE;
400   //
401   //  if (esdTr->Pt()<3) return kFALSE;//RS
402   return kTRUE;
403 }
404
405 //_________________________________________________________________________
406 void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr, Int_t esdID)
407 {
408   // find prolongaion candidates finding for single seed
409   //
410   AliITSUSeed seedUC;  // copy of the seed from the upper layer
411   AliITSUSeed seedT;   // transient seed between the seedUC and new prolongation hypothesis
412   //
413   AliITSUTrackHyp* hypTr = InitHypothesis(esdTr,esdID);  // initialize prolongations hypotheses tree
414   fCurrHyp = fWorkHyp;
415   fCurrHyp->InitFrom(hypTr);
416   //
417   AliITSURecoSens *hitSens[AliITSURecoSens::kNNeighbors+1];
418   //
419   TObjArray clArr; // container for transfer of clusters matching to seed
420   //
421   for (int ila=fNLrActive;ila--;) {
422     fCurrActLrID = ila;
423     fCurrLayer = fITS->GetLayerActive(ila);
424     Bool_t noClSharing = fCurrTrackCond->GetClSharing(ila)==0;
425     int ilaUp = ila+1;                         // prolong seeds from layer above
426     //
427     // for the outermost layer the seed is created from the ESD track
428     int nSeedsUp = (ilaUp==fNLrActive) ? 1 : fCurrHyp->GetNSeeds(ilaUp);
429     int maxNBranches   = fCurrTrackCond->GetMaxBranches(ila);
430     int maxNCandidates = fCurrTrackCond->GetMaxCandidates(ila);
431     //
432     for (int isd=0;isd<nSeedsUp;isd++) {
433       AliITSUSeed* seedU;
434       if (ilaUp==fNLrActive) {
435         seedU = 0;
436         seedUC.InitFromESDTrack(esdTr);
437       }
438       else {
439         seedU = fCurrHyp->GetSeed(ilaUp,isd);  // seed on prev.active layer to prolong  
440         if (seedU->IsKilled()) continue;
441         seedUC = *seedU;                       // its copy will be prolonged
442         seedUC.SetParent(seedU);        
443       }
444       seedUC.ResetFMatrix();                    // reset the matrix for propagation to next layer
445       // go till next active layer
446       AliDebug(2,Form("working on Lr:%d Seed:%d of %d for esdID=%d (MClb:%d) | pT=%.3f",ila,isd,nSeedsUp,esdID,fCurrESDtrMClb,seedUC.Pt()));
447       if (!TransportToLayer(&seedUC, fITS->GetLrIDActive(ilaUp), fITS->GetLrIDActive(ila)) ) {
448         //
449         AliDebug(2,Form("Transport failed | esdID=%d (MClb:%d)",esdID,fCurrESDtrMClb));
450         // Check if the seed satisfies to track definition
451         if (NeedToKill(&seedUC,kTransportFailed) && seedU) KillSeed(seedU,kTRUE);
452         continue; // RS TODO: decide what to do with tracks stopped on higher layers w/o killing
453       }
454       if (!GetRoadWidth(&seedUC, ila)) { // failed to find road width on the layer
455         if (NeedToKill(&seedUC,kRWCheckFailed) && seedU) KillSeed(seedU,kTRUE);
456         continue;
457       }
458       /*
459       //RS toremove
460       int mcquest = -1;
461       if (!seedUC.ContainsFake() && AliDebugLevelClass()>2) {
462         mcquest = fCurrESDtrMClb;
463         seedUC.Print("etp");
464         printf("FSParams: "); for (int ip=0;ip<kNTrImpData;ip++) printf("%+e ",fTrImpData[ip]); printf("\n");
465       }
466       //
467       int nsens = fCurrLayer->FindSensors(&fTrImpData[kTrPhi0], hitSens, mcquest);  // find detectors which may be hit by the track
468       */
469       int nsens = fCurrLayer->FindSensors(&fTrImpData[kTrPhi0], hitSens);  // find detectors which may be hit by the track
470       AliDebug(2,Form("Will check %d sensors on lr:%d | esdID=%d (MClb:%d)",nsens,ila,esdID,fCurrESDtrMClb));
471       //
472       seedUC.SetLr(ila);
473       //
474       double bz = GetBz();
475       for (int isn=nsens;isn--;) {
476         AliITSURecoSens* sens = hitSens[isn];
477         int ncl = sens->GetNClusters();
478         if (!ncl) continue;
479         seedT = seedUC;
480         //
481         // We need to propagate the seed to sensor on fCurrLayer staying in the frame of the sensor from prev.layer,
482         // since the transport matrix should be defined in this frame.
483         double xs; // X in the TF of current seed, corresponding to intersection with sensor plane
484         if (!seedT.GetTrackingXAtXAlpha(sens->GetXTF(),sens->GetPhiTF(),bz, xs)) {
485           if (AliDebugLevelClass()>2) {
486             printf("Failed on GetTrackingXAtXAlpha: X=%.4f alp=%.3f\n",sens->GetXTF(),sens->GetPhiTF());
487             seedT.Print("etp");
488           }
489           continue;
490         }
491         if (xs<seedT.GetX()) {
492           if (!PropagateSeed(&seedT,xs,fCurrMass)) continue;
493         }
494         else { // some low precision tracks may hit the sensor plane outside of the layer radius
495           if (AliDebugLevelClass()>2) {
496             if (!seedT.ContainsFake()) {
497               printf("WRONG PROP on L%d, sens%d of %d: %.4f > %.4f\n",ila,isn,nsens,xs,seedT.GetX());
498               sens->Print();
499               seedT.Print("etp");           
500             }
501           }
502           if (!seedT.PropagateParamOnlyTo(xs,bz)) continue;
503         }
504         //      if (!seedT.PropagateToX(xs,bz)) continue;
505         //      if (!seedT.Rotate(sens->GetPhiTF())) continue;
506         if (!seedT.RotateToAlpha(sens->GetPhiTF())) continue;
507         //
508         int clID0 = sens->GetFirstClusterId();
509         for (int icl=ncl;icl--;) { // don't change the order, clusters are sorted
510           //AliLog::SetClassDebugLevel("AliITSUTrackerGlo",10);
511           int clID = clID0 + icl;
512           if (noClSharing && fCurrLayer->GetCluster(clID)->IsClusterUsed()) continue;
513           int res = CheckCluster(&seedT,ila,clID0+icl);
514           //AliLog::SetClassDebugLevel("AliITSUTrackerGlo", 0);
515           //
516           if (res==kStopSearchOnSensor) break;     // stop looking on this sensor
517           if (res==kClusterNotMatching) continue;  // cluster does not match
518           // cluster is matching and it was added to the hypotheses tree
519         }
520       }
521       // cluster search is done. Do we need to have a version of this seed skipping current layer
522       if (!NeedToKill(&seedUC,kMissingCluster)) {
523         AliITSUSeed* seedSkp = NewSeedFromPool(&seedUC);
524         double penalty = -fCurrTrackCond->GetMissPenalty(ila);
525         // to do: make penalty to account for probability to miss the cluster for good reason
526         seedSkp->SetChi2Cl(penalty);
527         AddSeedBranch(seedSkp);
528         //      AddProlongationHypothesis(seedSkp,ila);      
529       }
530       // transfer the new branches of the seed to the hypothesis container
531       if (fNBranchesAdded) ValidateAllowedBranches(maxNBranches);
532       //
533     }
534     if (fNCandidatesAdded) ValidateAllowedCandidates(ila,maxNCandidates);
535     //    ((TObjArray*)fCurrHyp->GetLayerSeeds(ila))->Sort();
536     if (AliDebugLevelClass()>2) { //RS
537       printf(">>> All hypotheses on lr %d: \n",ila);
538       for (int ih=0;ih<fCurrHyp->GetNSeeds(ila);ih++) {printf(" #%3d ",ih); fCurrHyp->GetSeed(ila,ih)->Print();}
539     }
540     //
541     /*
542     for (int ih=0;ih<fCurrHyp->GetNSeeds(ila);ih++) {
543       if (ila!=0) continue;
544       double vecL[5] = {0};
545       double matL[15] = {0};
546       AliITSUSeed* sp = fCurrHyp->GetSeed(ila,ih);
547       while(sp->GetParent()) {
548         sp->Smooth(vecL,matL);
549         if (sp->GetLayerID()>=fNLrActive-1) break;
550         sp = (AliITSUSeed*)sp->GetParent();
551       }
552       */
553   }
554   //
555   SaveReducedHypothesesTree(hypTr);
556   if (AliDebugLevelClass()>1) {
557     printf("\nSaved hypotheses for esdTrack %d (MCLab:%d)\n",esdID,fCurrESDtrMClb);
558     hypTr->Print("l");
559   }
560   fCurrHyp = 0;
561   //
562 }
563
564 //_________________________________________________________________________
565 AliITSUTrackHyp* AliITSUTrackerGlo::InitHypothesis(AliESDtrack *esdTr, Int_t esdID)
566 {
567   // init prolongaion candidates finding for single seed
568   AliITSUTrackHyp* hyp = GetTrackHyp(esdID);
569   if (hyp) return hyp;
570   //
571   fCountProlongationTrials++;
572   //
573   fCurrMass = esdTr->GetMass();
574   if (fCurrMass<kPionMass*0.9) fCurrMass = kPionMass; // don't trust to mu, e identification from TPCin
575   //
576   hyp = new AliITSUTrackHyp(fNLrActive);
577   hyp->SetESDTrack(esdTr);
578   hyp->SetUniqueID(esdID);
579   hyp->SetMass(fCurrMass);
580   SetTrackHyp(hyp,esdID);
581   //
582   return hyp;
583   // TO DO
584 }
585
586 //_________________________________________________________________________
587 Bool_t AliITSUTrackerGlo::TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t lTo)
588 {
589   // transport seed from layerFrom to the entrance (along the direction of the track) of layerTo
590   //  
591   //
592   if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
593   //
594   int dir = lTo > lFrom ? 1:-1;
595   AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
596   Bool_t checkFirst = kTRUE;
597   while(lFrom!=lTo) {
598     if (lrFr) {
599       if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) {
600         //printf("FailHere0 Dir=%d\n",dir);
601         //seed->Print("etp");
602         return kFALSE; // go till the end of current layer
603       }
604       checkFirst = kFALSE;
605     }
606     AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
607     if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
608     //
609     // go the entrance of the layer, assuming no materials in between
610     double xToGo = lrTo->GetR(-dir);
611     //    double xts = xToGo;
612     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
613       //printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
614       //seed->Print("etp");
615       //
616       
617       return kFALSE;
618     }
619     AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
620     if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
621       //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
622       //seed->Print("etp");
623       return kFALSE;
624     }
625     lrFr = lrTo;
626   }
627   return kTRUE;
628   //
629 }
630
631 //_________________________________________________________________________
632 Bool_t AliITSUTrackerGlo::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo)
633 {
634   // transport track from layerFrom to the entrance of layerTo
635   //  
636   if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
637   //
638   int dir = lTo > lFrom ? 1:-1;
639   AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
640   Bool_t checkFirst = kTRUE;
641   while(lFrom!=lTo) {
642     if (lrFr) {
643       if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
644       checkFirst = kFALSE;
645     }
646     AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
647     if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
648     //
649     // go the entrance of the layer, assuming no materials in between
650     double xToGo = lrTo->GetR(-dir);
651     //    double xts = xToGo;
652     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
653       //      printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
654       //      seed->Print("etp");
655       return kFALSE;
656     }
657     AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
658     if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
659       //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
660       //seed->Print("etp");
661       return kFALSE;
662     }
663     lrFr = lrTo;
664   }
665   return kTRUE;
666   //
667 }
668
669 //_________________________________________________________________________
670 Bool_t AliITSUTrackerGlo::TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t xStop)
671 {
672   // transport track from layerFrom to the entrance of layerTo but do not pass control parameter X
673   //  
674   if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
675   //
676   int dir = lTo > lFrom ? 1:-1;
677   AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
678   Bool_t checkFirst = kTRUE;
679   while(lFrom!=lTo) {
680     if (lrFr) {
681       if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
682       checkFirst = kFALSE;
683     }
684     AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
685     if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
686     //
687     // go the entrance of the layer, assuming no materials in between
688     double xToGo = lrTo->GetR(-dir); // R of the entrance to layer
689     //
690     //    double xts = xToGo;
691     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
692       //      printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
693       //      seed->Print("etp");
694       return kFALSE;
695     }
696     if ( (dir>0&&xToGo>xStop) || (dir<0&&xToGo<xStop) ) xToGo = xStop;
697     //
698     AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
699     if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
700       //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
701       //seed->Print("etp");
702       return kFALSE;
703     }
704     lrFr = lrTo;
705   }
706   return kTRUE;
707   //
708 }
709
710 //_________________________________________________________________________
711 Bool_t AliITSUTrackerGlo::GoToExitFromLayer(AliITSUSeed* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
712 {
713   // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
714   // If check is requested, do this only provided the track has not exited the layer already
715   double xToGo = lr->GetR(dir);
716   if (check) { // do we need to track till the surface of the current layer ?
717     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
718     AliDebug(3,Form(" dir:%d Cur: %e Tgt: %e",dir,Sqrt(curR2),xToGo));
719     if      (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
720     else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
721   }
722   if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
723   // go via layer to its boundary, applying material correction.
724   if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
725   //
726   return kTRUE;
727   //
728 }
729
730 //_________________________________________________________________________
731 Bool_t AliITSUTrackerGlo::GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
732 {
733   // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
734   // If check is requested, do this only provided the track has not exited the layer already
735   double xToGo = lr->GetR(dir);
736   if (check) { // do we need to track till the surface of the current layer ?
737     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
738     if      (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
739     else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
740   }
741   AliDebug(2,Form(" dir=%d : from R=%.4f -> R=%.4f",dir,Sqrt(seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY()), xToGo));
742   if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
743   // go via layer to its boundary, applying material correction.
744   if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
745   //
746   return kTRUE;
747   //
748 }
749
750
751 //_________________________________________________________________________
752 Bool_t AliITSUTrackerGlo::GoToEntranceToLayer(AliITSUSeed* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
753 {
754   // go to the entrance of lr in direction dir, w/o applying material corrections.
755   // If check is requested, do this only provided the track did not reach the layer already
756   double xToGo = lr->GetR(-dir);
757   if (check) { // do we need to track till the surface of the current layer ?
758     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
759     if      (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
760     else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
761   }
762   if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
763   // go via layer to its boundary, applying material correction.
764   if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
765   return kTRUE;
766   //
767 }
768
769 //_________________________________________________________________________
770 Bool_t AliITSUTrackerGlo::GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
771 {
772   // go to the entrance of lr in direction dir, w/o applying material corrections.
773   // If check is requested, do this only provided the track did not reach the layer already
774   double xToGo = lr->GetR(-dir);
775   if (check) { // do we need to track till the surface of the current layer ?
776     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
777     if      (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
778     else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
779   }
780   if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
781   // go via layer to its boundary, applying material correction.
782   if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
783   return kTRUE;
784   //
785 }
786
787 //_________________________________________________________________________
788 Bool_t AliITSUTrackerGlo::GetRoadWidth(AliITSUSeed* seed, int ilrA)
789 {
790   // calculate road width in terms of phi and z for the track which MUST be on the external radius of the layer
791   // as well as some aux info
792   double bz = GetBz();
793   AliITSURecoLayer* lrA = fITS->GetLayerActive(ilrA);
794   seed->GetXYZ(&fTrImpData[kTrXIn]);    // lab position at the entrance from above
795   static AliExternalTrackParam sc;   // seed copy for manipulations
796   sc = *seed;
797   //
798   fTrImpData[kTrPhiIn] = ATan2(fTrImpData[kTrYIn],fTrImpData[kTrXIn]);
799   if (!sc.Rotate(fTrImpData[kTrPhiIn])) return kFALSE; // go to the frame of the entry point into the layer
800   double dr  = lrA->GetDR();                              // approximate X dist at the inner radius
801   if (!sc.GetXYZAt(sc.GetX()-dr, bz, fTrImpData + kTrXOut)) {
802     // special case: track does not reach inner radius, might be tangential
803     double r = sc.GetD(0,0,bz);
804     double x;
805     if (!sc.GetXatLabR(r,x,bz,-1)) {
806       sc.Print();
807       AliFatal(Form("This should not happen: r=%f",r));
808     }
809     dr = Abs(sc.GetX() - x);
810     if (!sc.GetXYZAt(x, bz, fTrImpData + kTrXOut)) {
811       sc.Print();
812       AliFatal(Form("This should not happen: x=%f",x));
813     }
814   }
815   //
816   fTrImpData[kTrPhiOut] = ATan2(fTrImpData[kTrYOut],fTrImpData[kTrXOut]);
817   double sgy = sc.GetSigmaY2() + dr*dr*sc.GetSigmaSnp2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(ilrA);
818   double sgz = sc.GetSigmaZ2() + dr*dr*sc.GetSigmaTgl2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(ilrA);
819   sgy = Sqrt(sgy)*fCurrTrackCond->GetNSigmaRoadY(ilrA);
820   sgz = Sqrt(sgz)*fCurrTrackCond->GetNSigmaRoadZ(ilrA);
821   double dphi0 = 0.5*Abs(fTrImpData[kTrPhiOut]-fTrImpData[kTrPhiIn]);
822   double phi0  = 0.5*(fTrImpData[kTrPhiOut]+fTrImpData[kTrPhiIn]);
823   if ( dphi0>(0.5*Pi()) ) {
824     // special case of angles around pi 
825     dphi0 = Abs(phi0);    
826     phi0  += Pi();
827   }
828
829   fTrImpData[kTrPhi0] = phi0;
830   fTrImpData[kTrZ0]   = 0.5*(fTrImpData[kTrZOut]+fTrImpData[kTrZIn]);
831   fTrImpData[kTrDPhi] = dphi0 + sgy/lrA->GetR();
832   fTrImpData[kTrDZ]   = 0.5*Abs(fTrImpData[kTrZOut]-fTrImpData[kTrZIn])   + sgz;
833   //  
834   return kTRUE;
835 }
836
837 //________________________________________
838 void AliITSUTrackerGlo::ResetSeedsPool()
839 {
840   // mark all seeds in the pool as unused
841   AliDebug(-1,Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
842   fNFreeSeeds = 0;
843   fSeedsPool.Clear(); // seeds don't allocate memory
844 }
845
846
847 //________________________________________
848 void AliITSUTrackerGlo::MarkSeedFree(AliITSUSeed *sd) 
849 {
850   // account that this seed is "deleted" 
851   int id = sd->GetPoolID();
852   if (id<0) {
853     AliError(Form("Freeing of seed %p NOT from the pool is requested",sd)); 
854     return;
855   }
856   //  AliInfo(Form("%d %p",id, seed));
857   fSeedsPool.RemoveAt(id);
858   if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 );
859   fFreeSeedsID.GetArray()[fNFreeSeeds++] = id;
860 }
861
862 //_________________________________________________________________________
863 Int_t AliITSUTrackerGlo::CheckCluster(AliITSUSeed* track, Int_t lr, Int_t clID) 
864 {
865   // Check if the cluster (in tracking frame!) is matching to track. 
866   // The track must be already propagated to sensor tracking frame.
867   // Returns:  kStopSearchOnSensor if the search on given sensor should be stopped, 
868   //           kClusterMatching    if the cluster is matching
869   //           kClusterMatching    otherwise
870   //
871   const double kTolerX = 5e-4;
872   //  AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
873   AliCluster *cl = fCurrLayer->GetCluster(clID);
874   //
875   Bool_t goodCl = kFALSE;
876   int currLabel = Abs(fCurrESDtrMClb);
877   //
878   if (cl->GetLabel(0)>=0) {
879     for (int i=0;i<3;i++) if (cl->GetLabel(i)>=0 && cl->GetLabel(i)==currLabel) {goodCl = kTRUE; break;}
880   }
881   Bool_t goodSeed = !track->ContainsFake();
882   //
883   if (TMath::Abs(cl->GetX())>kTolerX) { // if due to the misalingment X is large, propagate track only
884     if (!track->PropagateParamOnlyTo(track->GetX()+cl->GetX(),GetBz())) {
885       if (goodCl&&goodSeed && AliDebugLevelClass()>2 ) {
886         AliDebug(2,Form("Lost good cl on L:%d failed propagation. |ESDtrack#%d (MClb:%d)",lr,fCurrESDtrack->GetID(),fCurrESDtrMClb)); 
887         track->Print("etp");
888         cl->Print();
889       }
890       return kStopSearchOnSensor; // propagation failed, seedT is intact
891     }
892   }
893   double dy = cl->GetY()-track->GetY();
894   double dz = cl->GetZ()-track->GetZ();
895   //
896 #ifdef  _FILL_CONTROL_HISTOS_
897   int hcOffs = fTrackPhase*kHistosPhase + lr;
898   double htrPt=-1;
899   if (goodCl && (((AliITSUClusterPix*)cl)->GetNPix()>1 || !((AliITSUClusterPix*)cl)->IsSplit()) && goodSeed && fCHistoArr /* && track->GetChi2Penalty()<1e-5*/) {
900     htrPt = track->Pt();
901     ((TH2*)fCHistoArr->At(kHResY+hcOffs))->Fill(htrPt,dy);
902     ((TH2*)fCHistoArr->At(kHResZ+hcOffs))->Fill(htrPt,dz);
903     double errY = track->GetSigmaY2();
904     double errZ = track->GetSigmaZ2();
905     if (errY>0) ((TH2*)fCHistoArr->At(kHResYP+hcOffs))->Fill(htrPt,dy/Sqrt(errY));
906     if (errZ>0) ((TH2*)fCHistoArr->At(kHResZP+hcOffs))->Fill(htrPt,dz/Sqrt(errZ));
907   }
908 #endif  
909   //
910   double dy2 = dy*dy;
911   double tol2 = (track->GetSigmaY2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(lr))*
912     fCurrTrackCond->GetNSigmaRoadY(lr)*fCurrTrackCond->GetNSigmaRoadY(lr); // RS TOOPTIMIZE
913   if (dy2>tol2) {                          // the clusters are sorted in Z(col) then in Y(row). 
914     if (goodCl&&goodSeed && AliDebugLevelClass()>2) {    
915       AliDebug(2,Form("Lost good cl: dy2=%e > tol2=%e |ESDtrack#%d (MClb:%d)",dy2,tol2,fCurrESDtrack->GetID(),fCurrESDtrMClb)); 
916       track->Print("etp");
917       cl->Print();
918       PrintSeedClusters(track);
919     }
920     // clusters are sorted in increasing Y and the loop goes from last (highers Y) to 1st.
921     // if the track has too large y for this cluster (dy<0) it will be so for all other clusters of the sensor
922     if (dy<0) return kStopSearchOnSensor;  // No chance that other cluster of this sensor will match (all Y's will be even larger)
923     else      return kClusterNotMatching;   // Other clusters may match
924   }
925   double dz2 = dz*dz;
926   tol2 = (track->GetSigmaZ2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(lr))*
927     fCurrTrackCond->GetNSigmaRoadZ(lr)*fCurrTrackCond->GetNSigmaRoadZ(lr); // RS TOOPTIMIZE
928   if (dz2>tol2) {
929     if (goodCl&&goodSeed && AliDebugLevelClass()>2) {
930       AliDebug(2,Form("Lost good cl on L:%d : dz2=%e > tol2=%e |ESDtrack#%d (MClb:%d)",lr,dz2,tol2,fCurrESDtrack->GetID(),fCurrESDtrMClb)); 
931       track->Print("etp");
932       cl->Print();
933       PrintSeedClusters(track);
934     }
935     return kClusterNotMatching; // Other clusters may match
936   }
937   //
938   // check chi2
939   Double_t p[2]={cl->GetY(), cl->GetZ()};
940   Double_t cov[3]={cl->GetSigmaY2(), cl->GetSigmaYZ(), cl->GetSigmaZ2()};
941   double chi2 = track->GetPredictedChi2(p,cov);
942   //
943 #ifdef  _FILL_CONTROL_HISTOS_
944   if (htrPt>0) {
945     ((TH2*)fCHistoArr->At(kHChi2Cl+hcOffs))->Fill(htrPt,chi2);
946   }
947 #endif
948   //
949   if (chi2>fCurrTrackCond->GetMaxTr2ClChi2(lr)) {
950     if (goodCl&&goodSeed && AliDebugLevelClass()>2) {
951       AliDebug(2,Form("Lost good cl on L:%d : Chi2=%e > Chi2Max=%e |dy: %+.3e dz: %+.3e |ESDtrack#%d (MClb:%d)\n",
952                       lr,chi2,fCurrTrackCond->GetMaxTr2ClChi2(lr),dy,dz,fCurrESDtrack->GetID(),fCurrESDtrMClb)); 
953       track->Print("etp");
954       cl->Print("");
955       PrintSeedClusters(track);
956     }
957     return kClusterNotMatching;
958   }
959   //
960   track = NewSeedFromPool(track);  // input track will be reused, use its clone for updates
961   if (!track->Update()) {
962     if (goodCl&&goodSeed && AliDebugLevelClass()>2) {
963       AliDebug(2,Form("Lost good cluster on L:%d : failed to update",lr));
964       track->Print("etp");
965       cl->Print("");
966       PrintSeedClusters(track);
967     }
968     MarkSeedFree(track);
969     return kClusterNotMatching;
970   }
971   track->SetChi2Cl(chi2);
972   track->SetLrClusterID(lr,clID);
973   //  cl->IncreaseClusterUsage(); // do this only for winners
974   //
975   track->SetFake(!goodCl);
976   //
977   AliDebug(2,Form("AddCl(%d) Cl%d lr:%d: dY:%+8.4f dZ:%+8.4f (MC: %5d %5d %5d) |Chi2=%f(%c)| ",
978                goodCl,clID,lr,dy,dz2,cl->GetLabel(0),cl->GetLabel(1),cl->GetLabel(2), chi2, track->IsFake() ? '-':'+'));
979   //
980   AddSeedBranch(track);
981 #ifdef  _FILL_CONTROL_HISTOS_
982   if (htrPt>0) {
983     ((TH2*)fCHistoArr->At(kHChi2Nrm+hcOffs))->Fill(htrPt,track->GetChi2GloNrm());
984   }
985 #endif
986   //
987   return kClusterMatching;
988 }
989
990 //_________________________________________________________________________
991 Bool_t AliITSUTrackerGlo::NeedToKill(AliITSUSeed *seed, Int_t flag)
992 {
993   // check if the seed should not be discarded
994   const UShort_t kMask = 0xffff;
995   if (flag==kMissingCluster) {
996     int lastChecked = seed->GetLayerID();
997     UShort_t patt = seed->GetHitsPattern();
998     if (lastChecked) patt |= ~(kMask<<lastChecked); // not all layers were checked, complete unchecked ones by potential hits
999     Bool_t seedOK = fCurrTrackCond->CheckPattern(patt);
1000     return !seedOK;
1001   }
1002   return kTRUE;
1003 }
1004
1005 //______________________________________________________________________________
1006 Bool_t AliITSUTrackerGlo::PropagateSeed(AliITSUSeed *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr) 
1007 {
1008   // propagate seed to given x applying material correction if requested
1009   const Double_t kEpsilon = 1e-5;
1010   Double_t xpos     = seed->GetX();
1011   Int_t dir         = (xpos<xToGo) ? 1:-1;
1012   Double_t xyz0[3],xyz1[3],param[7];
1013   //
1014   Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
1015   if (matCorr || updTime) seed->GetXYZ(xyz1);   //starting global position
1016   while ( (xToGo-xpos)*dir > kEpsilon){
1017     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
1018     Double_t x    = xpos+step;
1019     Double_t bz=GetBz();   // getting the local Bz
1020     if (!seed->PropagateToX(x,bz)) return kFALSE;
1021     double ds = 0;
1022     if (matCorr || updTime) {
1023       xyz0[0]=xyz1[0]; // global pos at the beginning of step
1024       xyz0[1]=xyz1[1];
1025       xyz0[2]=xyz1[2];
1026       seed->GetXYZ(xyz1);    //  // global pos at the end of step
1027       if (matCorr) {
1028         MeanMaterialBudget(xyz0,xyz1,param);    
1029         Double_t xrho=param[0]*param[4], xx0=param[1];
1030         if (dir>0) xrho = -xrho; // outward should be negative
1031         if (!seed->ApplyMaterialCorrection(xx0,xrho,mass,kFALSE)) {AliInfo("Fail2"); seed->Print("etp"); return kFALSE;}
1032         ds = param[4];
1033       }
1034        else { // matCorr is not requested but time integral is
1035         double d0 = xyz1[0]-xyz0[0];
1036         double d1 = xyz1[1]-xyz0[1];
1037         double d2 = xyz1[2]-xyz0[2];    
1038         ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
1039       }     
1040     }
1041     if (updTime) seed->AddTimeStep(ds);
1042     xpos = seed->GetX();
1043   }
1044   return kTRUE;
1045 }
1046
1047 //______________________________________________________________________________
1048 Bool_t AliITSUTrackerGlo::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr) 
1049 {
1050   // propagate seed to given x applying material correction if requested
1051   const Double_t kEpsilon = 1e-5;
1052   Double_t xpos     = seed->GetX();
1053   Int_t dir         = (xpos<xToGo) ? 1:-1;
1054   Double_t xyz0[3],xyz1[3],param[7];
1055   //
1056   if (AliDebugLevelClass()>1) {
1057     AliDebug(2,Form("Before Propagate to X=%f with M=%.3f MaxStep=%.4f MatCorr=%d",xToGo,mass,maxStep,matCorr));
1058     seed->AliExternalTrackParam::Print();
1059   }
1060   Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
1061   if (matCorr || updTime) seed->GetXYZ(xyz1);   //starting global position
1062   while ( (xToGo-xpos)*dir > kEpsilon){
1063     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
1064     Double_t x    = xpos+step;
1065     Double_t bz=GetBz();   // getting the local Bz
1066     if (!seed->PropagateTo(x,bz))  return kFALSE;
1067     double ds = 0;
1068     if (matCorr || updTime) {
1069       xyz0[0]=xyz1[0]; // global pos at the beginning of step
1070       xyz0[1]=xyz1[1];
1071       xyz0[2]=xyz1[2];
1072       seed->GetXYZ(xyz1);    //  // global pos at the end of step
1073       //
1074       if (matCorr) {
1075         MeanMaterialBudget(xyz0,xyz1,param);    
1076         Double_t xrho=param[0]*param[4], xx0=param[1];
1077         if (dir>0) xrho = -xrho; // outward should be negative
1078         if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
1079         ds = param[4];
1080       }
1081       else { // matCorr is not requested but time integral is
1082         double d0 = xyz1[0]-xyz0[0];
1083         double d1 = xyz1[1]-xyz0[1];
1084         double d2 = xyz1[2]-xyz0[2];    
1085         ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
1086       }
1087     }
1088     if (updTime) seed->AddTimeStep(ds);
1089     //
1090     xpos = seed->GetX();
1091   }
1092   if (AliDebugLevelClass()>1) {
1093     AliDebug(2,Form("After Propagate to X=%f",xToGo));
1094     seed->AliExternalTrackParam::Print();
1095   }
1096   return kTRUE;
1097 }
1098
1099 //______________________________________________________________________________
1100 void AliITSUTrackerGlo::FinalizeHypotheses()
1101 {
1102   // select winner for each hypothesis, remove cl. sharing conflicts
1103   AliDebug(2,"TODO");
1104   //
1105   int nh = fHypStore.GetEntriesFast();
1106   AliITSUSeed* winner = 0;
1107   for (int ih=0;ih<nh;ih++) {
1108     AliITSUTrackHyp* hyp = (AliITSUTrackHyp*) fHypStore.UncheckedAt(ih); 
1109     if (!hyp) continue; // removed
1110     if (hyp->GetESDTrack()->IsOn(AliESDtrack::kITSin)) continue; // winner already defined
1111     //
1112     if (!(winner=hyp->DefineWinner())) continue;   // failed to find winner
1113     FlagSeedClusters(winner,kTRUE);
1114     fCountITSin++;
1115     CookMCLabel(hyp);
1116     UpdateESDTrack(hyp,AliESDtrack::kITSin);
1117   }
1118   //
1119 }
1120
1121 //______________________________________________________________________________
1122 void AliITSUTrackerGlo::UpdateESDTrack(AliITSUTrackHyp* hyp,Int_t flag)
1123 {
1124   // update ESD track with current best hypothesis
1125   AliESDtrack* esdTr = hyp->GetESDTrack();
1126   if (!esdTr) return;
1127   AliITSUSeed* win = hyp->GetWinner();
1128   if (!win) return;
1129   switch (flag) {
1130   case AliESDtrack::kITSin: 
1131     esdTr->UpdateTrackParams(hyp,flag); // update kinematics
1132     // TODO: set cluster info
1133     break;
1134     //
1135   case AliESDtrack::kITSout: 
1136     esdTr->UpdateTrackParams(hyp,flag); // update kinematics
1137     // TODO: avoid setting friend
1138     break;
1139     //
1140   case AliESDtrack::kITSrefit: 
1141     esdTr->UpdateTrackParams(hyp,flag); // update kinematics
1142     // TODO: avoid setting cluster info
1143     break;
1144   default:
1145     AliFatal(Form("Unknown flag %d",flag));
1146   }
1147   //
1148   // transfer module indices
1149   // TODO
1150 }
1151
1152 //______________________________________________________________________________
1153 Bool_t AliITSUTrackerGlo::RefitTrack(AliITSUTrackHyp* trc, Double_t rDest, Bool_t stopAtLastCl)
1154 {
1155   // refit track till radius rDest
1156   AliITSUTrackHyp tmpTr;
1157   //
1158   double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY());
1159   int dir,lrStart,lrStop;
1160   //
1161   dir = rCurr<rDest ? 1 : -1;
1162   lrStart = fITS->FindFirstLayerID(rCurr,dir);
1163   lrStop  = fITS->FindLastLayerID(rDest,dir);
1164   if (AliDebugLevelClass()>2) {
1165     printf("Refit %d: Lr:%d (%f) -> Lr:%d (%f)\n",dir,lrStart,rCurr, lrStop,rDest);
1166     printf("Before refit: "); trc->AliExternalTrackParam::Print();
1167   }
1168   if (lrStop<0 || lrStart<0) AliFatal(Form("Failed to find start(%d) or last(%d) layers. Track from %.3f to %.3f",lrStart,lrStop,rCurr,rDest));
1169   //
1170   int ncl = trc->FetchClusterInfo(fClInfo);
1171   fCurrMass = trc->GetMass();
1172   tmpTr.AliKalmanTrack::operator=(*trc);
1173   tmpTr.SetChi2(0);
1174   int iclLr[2],nclLr,clCount=0;
1175   //
1176   for (int ilr=lrStart;ilr!=lrStop;ilr+=dir) {
1177     AliITSURecoLayer* lr = fITS->GetLayer(ilr);
1178     if ( dir*(rCurr-lr->GetR(dir))>0) continue; // this layer is already passed
1179     int ilrA2,ilrA = lr->GetActiveID();
1180     // passive layer or active w/o hits will be traversed on the way to next cluster
1181     if (!lr->IsActive() || fClInfo[ilrA2=(ilrA<<1)]<0) continue; 
1182     //
1183     // select the order in which possible 2 clusters (in case of the overlap) will be traversed and fitted
1184     nclLr=0;
1185     if (dir>0) { // clusters are stored in increasing radius order
1186       iclLr[nclLr++]=fClInfo[ilrA2++];
1187       if (fClInfo[ilrA2]>=0) iclLr[nclLr++]=fClInfo[ilrA2];
1188     }
1189     else {
1190       if ( fClInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=fClInfo[ilrA2+1];
1191       iclLr[nclLr++]=fClInfo[ilrA2];
1192     }
1193     //
1194     Bool_t transportedToLayer = kFALSE;
1195     for (int icl=0;icl<nclLr;icl++) {
1196       AliITSUClusterPix* clus =  (AliITSUClusterPix*)lr->GetCluster(iclLr[icl]);
1197       AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
1198       if (!tmpTr.Rotate(sens->GetPhiTF())) {
1199         AliDebug(2,Form("Failed on rotate to %f | ESDtrack#%d (MClb:%d)",sens->GetPhiTF(),fCurrESDtrack->GetID(),fCurrESDtrMClb));
1200         return kFALSE;
1201       }
1202       //
1203       double xClus = sens->GetXTF()+clus->GetX();
1204       if (!transportedToLayer) {
1205         if (ilr!=lrStart && !TransportToLayerX(&tmpTr,lrStart,ilr,xClus)) {
1206           AliDebug(2,Form("Failed to transport %d -> %d | ESDtrack#%d (MClb:%d)\n",lrStart,ilr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1207           //tmpTr.AliExternalTrackParam::Print(); trc->GetWinner()->Print("etp");
1208           return kFALSE; // go to the entrance to the layer
1209         }
1210         lrStart = ilr;
1211         transportedToLayer = kTRUE;
1212       }
1213       //
1214       if (AliDebugLevelClass()>1) {
1215         AliDebug(2,Form("Propagate to cl:%d on lr %d Need to go %.4f -> %.4f",icl,ilrA,tmpTr.GetX(),xClus));
1216         //      printf("Before: "); tmpTr.AliExternalTrackParam::Print();
1217       }
1218       //
1219       if (!PropagateSeed(&tmpTr,xClus,fCurrMass)) {
1220         AliDebug(2,Form("Failed on propagate to %f (dir=%d) | ESDtrack#%d (MClb:%d)",xClus,dir,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1221         //tmpTr.AliExternalTrackParam::Print(""); trc->GetWinner()->Print("etp");
1222         return kFALSE;
1223       }
1224       //
1225 #ifdef  _FILL_CONTROL_HISTOS_
1226       int hcOffs = fTrackPhase*kHistosPhase + ilrA;
1227       double htrPt=-1;
1228       if (fCHistoArr && trc->GetLabel()>=0/* && trc->Charge()>0*/) {
1229         htrPt = tmpTr.Pt();
1230         double dy = clus->GetY()-tmpTr.GetY();
1231         double dz = clus->GetZ()-tmpTr.GetZ();
1232         ((TH2*)fCHistoArr->At(kHResY+hcOffs))->Fill(htrPt,dy);
1233         ((TH2*)fCHistoArr->At(kHResZ+hcOffs))->Fill(htrPt,dz);
1234         double errY = tmpTr.GetSigmaY2();
1235         double errZ = tmpTr.GetSigmaZ2();
1236         if (errY>0) ((TH2*)fCHistoArr->At(kHResYP+hcOffs))->Fill(htrPt,dy/Sqrt(errY));
1237         if (errZ>0) ((TH2*)fCHistoArr->At(kHResZP+hcOffs))->Fill(htrPt,dz/Sqrt(errZ));
1238       }
1239 #endif  
1240       //
1241       double chi2;
1242       if ( (chi2=tmpTr.Update(clus))<0 ) {
1243         AliDebug(2,Form("Failed on Update | ESDtrack#%d (MClb:%d)",fCurrESDtrack->GetID(),fCurrESDtrMClb));             
1244         return kFALSE;
1245       }
1246 #ifdef  _FILL_CONTROL_HISTOS_
1247       if (htrPt>0) {
1248         ((TH2*)fCHistoArr->At(kHChi2Cl+hcOffs))->Fill(htrPt,chi2);
1249       }
1250 #endif 
1251       if (AliDebugLevelClass()>1) {
1252         printf("AfterRefit: "); tmpTr.AliExternalTrackParam::Print();
1253       }
1254       if (stopAtLastCl && ++clCount==ncl) return kTRUE; // it was requested to not propagate after last update
1255     }
1256     //
1257   }
1258   // All clusters were succesfully fitted. Even if the track does not reach rDest, this is enough to validate it.
1259   // Still, try to go as close as possible to rDest.
1260   //
1261   if (lrStart!=lrStop) {
1262     //printf("Going to last layer %d -> %d\n",lrStart,lrStop);
1263     if (!TransportToLayer(&tmpTr,lrStart,lrStop)) {
1264       AliDebug(1,Form("Failed on TransportToLayer %d->%d | ESDtrack#%d (MClb:%d)",lrStart,lrStop,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1265       //tmpTr.AliExternalTrackParam::Print();
1266       //trc->GetWinner()->Print("etp");
1267       return kTRUE;
1268     }    
1269     if (!GoToExitFromLayer(&tmpTr,fITS->GetLayer(lrStop),dir)) {
1270       AliDebug(3,Form("Failed on GoToExitFromLayer %d | ESDtrack#%d (MClb:%d)",lrStop,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1271       return kTRUE; // go till the exit from layer
1272     }
1273     //
1274     //printf("On exit from last layer\n");
1275     //tmpTr.AliExternalTrackParam::Print();
1276     // go to the destination radius
1277     if (!tmpTr.GetXatLabR(rDest,rDest,GetBz(),dir)) return kTRUE;
1278     if (!PropagateSeed(&tmpTr,rDest,fCurrMass, 100, kFALSE)) return kTRUE;
1279   }
1280   trc->AliKalmanTrack::operator=(tmpTr);
1281   if (AliDebugLevelClass()>2) {
1282     printf("After refit (now at lr %d): ",lrStop); trc->AliExternalTrackParam::Print();
1283   }
1284   return kTRUE;
1285 }
1286
1287 //__________________________________________________________________
1288 void AliITSUTrackerGlo::CookMCLabel(AliITSUTrackHyp* hyp) 
1289 {
1290   // build MC label
1291   //
1292   const int kMaxLbPerCl = 3;
1293   int lbID[kMaxLayers*6],lbStat[kMaxLayers*6];
1294   Int_t lr,nLab=0,nCl=0;
1295   AliITSUSeed *seed = hyp->GetWinner();
1296   while(seed) {
1297     int clID = seed->GetLrCluster(lr);
1298     if (clID>=0) {
1299       AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
1300       nCl++;
1301       for (int imc=0;imc<kMaxLbPerCl;imc++) { // labels within single cluster
1302         int trLb = cl->GetLabel(imc);
1303         if (trLb<0) break;
1304         // search this mc track in already accounted ones
1305         int iLab;
1306         for (iLab=0;iLab<nLab;iLab++) if (lbID[iLab]==trLb) break;
1307         if (iLab<nLab) lbStat[iLab]++;
1308         else {
1309           lbID[nLab] = trLb;
1310           lbStat[nLab++] = 1;
1311         }
1312       } // loop over given cluster's labels
1313     }
1314     seed = (AliITSUSeed*)seed->GetParent();
1315   } // loop over clusters
1316   // 
1317   AliESDtrack* esdTr = hyp->GetESDTrack();
1318   int tpcLab = esdTr ? Abs(esdTr->GetTPCLabel()) : -kDummyLabel;
1319   if (nCl && nLab) {
1320     int maxLab=0,nTPCok=0;
1321     for (int ilb=nLab;ilb--;) {
1322       int st = lbStat[ilb];
1323       if (lbStat[maxLab]<st) maxLab=ilb;
1324       if (tpcLab==lbID[ilb]) nTPCok += st;
1325     }
1326     hyp->SetFakeRatio(1.-float(nTPCok)/nCl);
1327     hyp->SetLabel( nTPCok==nCl ? tpcLab : -tpcLab);
1328     hyp->SetITSLabel( lbStat[maxLab]==nCl ? lbStat[maxLab] : -lbStat[maxLab]); // winner label
1329     return;
1330   }
1331   //
1332   hyp->SetFakeRatio(-1.);
1333   hyp->SetLabel( -tpcLab );
1334   hyp->SetITSLabel( kDummyLabel );
1335 }
1336
1337 //__________________________________________________________________
1338 Bool_t AliITSUTrackerGlo::AddSeedBranch(AliITSUSeed* seed)
1339 {
1340   // add new prolongation branch for the seed to temporary array. The seeds are sorted according to Compare f-n
1341   if (fNCandidatesAdded+fNBranchesAdded>=fLayerMaxCandidates ) {
1342     AliInfo(Form("Number of candidates at layer %d for current seed reached %d, increasing buffer",fCurrActLrID,fLayerMaxCandidates));
1343     fLayerMaxCandidates = 2*(fLayerMaxCandidates+1);
1344     AliITSUSeed** tmpArr = fLayerCandidates;
1345     fLayerCandidates = new AliITSUSeed*[fLayerMaxCandidates];
1346     memcpy(fLayerCandidates,tmpArr,(fNCandidatesAdded+fNBranchesAdded)*sizeof(AliITSUSeed*));
1347     delete tmpArr; // delete only array, not objects
1348   }
1349   AliITSUSeed** branches = &fLayerCandidates[fNCandidatesAdded]; // note: fNCandidatesAdded is incremented after adding all branches of current seed
1350   int slot=fNBranchesAdded++;
1351   for (int slotF=slot;slotF--;) { // slotF is always slot-1
1352     AliITSUSeed* si = branches[slotF];
1353     if (si->Compare(seed)<0) break; // found the last seed with better quality
1354     // otherwise, shift the worse seed to the next slot
1355     branches[slot] = si;
1356     slot = slotF; // slot should be slotF+1
1357   }
1358   // if needed, move worse seeds
1359   branches[slot] = seed;
1360   return kTRUE;
1361   //
1362 }
1363
1364 //__________________________________________________________________
1365 void AliITSUTrackerGlo::ValidateAllowedBranches(Int_t acceptMax)
1366 {
1367   // keep allowed number of branches for current seed and disable extras
1368   int nb = Min(fNBranchesAdded,acceptMax);
1369   //  if (nb<fNBranchesAdded) printf("ValidateAllowedBranches: %d of %d (%d) on lr %d\n",nb,fNBranchesAdded,fNCandidatesAdded,ilr);
1370   // disable unused branches
1371   AliITSUSeed** branches = &fLayerCandidates[fNCandidatesAdded];
1372 #ifdef  _FILL_CONTROL_HISTOS_
1373   int bestID = -1;
1374   for (int ib=0;ib<fNBranchesAdded;ib++) {
1375     AliITSUSeed* sd = branches[ib];
1376     if (!sd->ContainsFake() && (bestID<0 || sd->Compare(branches[bestID])<0) ) bestID = ib;
1377   }
1378   if (bestID>=0) {
1379     TH2* hb = (TH2*)fCHistoArr->At(kHBestInBranch + fTrackPhase*kHistosPhase + fCurrActLrID);
1380     if (hb) hb->Fill(branches[bestID]->Pt(), bestID);
1381   }
1382 #endif
1383   for (int ib=nb;ib<fNBranchesAdded;ib++) {
1384     /*
1385     if (AliDebugLevelClass()>+2 && !branches[ib]->ContainsFake()) {
1386       printf("Suppress good branch as %d of %d |",ib,fNBranchesAdded); branches[ib]->Print();
1387       printf("Survivors : \n");
1388       for (int j=0;j<nb;j++) branches[j]->Print();
1389     }
1390     */
1391     MarkSeedFree(branches[ib]);
1392   }
1393   fNCandidatesAdded += nb; // update total candidates counter
1394   fNBranchesAdded = 0; // reset branches counter
1395   //
1396 }
1397
1398 //__________________________________________________________________
1399 void AliITSUTrackerGlo::ValidateAllowedCandidates(Int_t ilr, Int_t acceptMax)
1400 {
1401   // transfer allowed number of branches to hypothesis container
1402   //
1403   // sort candidates in increasing order of chi2
1404   static int lastSize = 0;
1405   static int *index = 0;
1406   static float *chi2 = 0;
1407   if (fLayerMaxCandidates>lastSize) {
1408     lastSize = fLayerMaxCandidates;
1409     delete[] index;
1410     delete[] chi2;
1411     index = new int[lastSize];
1412     chi2  = new float[lastSize];
1413   }
1414   for (int i=fNCandidatesAdded;i--;) chi2[i] = fLayerCandidates[i]->GetChi2GloNrm();
1415   Sort(fNCandidatesAdded,chi2,index,kFALSE);
1416   //
1417 #ifdef  _FILL_CONTROL_HISTOS_
1418   int bestID = -1;
1419   for (int ib=0;ib<fNCandidatesAdded;ib++) {
1420     AliITSUSeed* sd = fLayerCandidates[index[ib]];
1421     if (!sd->ContainsFake() && (bestID<0 || sd->Compare(fLayerCandidates[index[bestID]])<0) ) bestID = ib;
1422   }
1423   if (bestID>=0) {
1424     TH2* hb = (TH2*)fCHistoArr->At(kHBestInCand + fTrackPhase*kHistosPhase + fCurrActLrID);
1425     if (hb) hb->Fill(fLayerCandidates[index[bestID]]->Pt(), bestID);
1426   }
1427 #endif
1428   //
1429   int nb = Min(fNCandidatesAdded,acceptMax);
1430   //
1431   for (int ib=0;ib<nb;ib++) AddProlongationHypothesis(fLayerCandidates[index[ib]],ilr);
1432   // disable unused candidates
1433   for (int ib=nb;ib<fNCandidatesAdded;ib++) {
1434     /*
1435     if (AliDebugLevelClass()>+2 && !fLayerCandidates[index[ib]]->ContainsFake()) {
1436       printf("Suppress good candidate as %d of %d |",index[ib],fNCandidatesAdded); fLayerCandidates[index[ib]]->Print();
1437     }
1438     */
1439     MarkSeedFree(fLayerCandidates[index[ib]]);    
1440   }
1441   fNCandidatesAdded = 0; // reset candidates counter
1442   //
1443 }
1444
1445 //__________________________________________________________________
1446 void AliITSUTrackerGlo::SaveReducedHypothesesTree(AliITSUTrackHyp* dest)
1447 {
1448   // remove those hypothesis seeds which dont lead to candidates at final layer
1449   //
1450   // 1st, flag the seeds to save
1451   int lr0 = 0;
1452   for (int isd=0;isd<fCurrHyp->GetNSeeds(lr0);isd++) {
1453     AliITSUSeed* seed = fCurrHyp->RemoveSeed(lr0,isd);
1454     if (!seed) continue;
1455     seed->FlagTree(AliITSUSeed::kSave);
1456     dest->AddSeed(seed,lr0);
1457   }
1458   for (int ilr=1;ilr<fNLrActive;ilr++) {
1459     int nsd = fCurrHyp->GetNSeeds(ilr);
1460     for (int isd=0;isd<nsd;isd++) {
1461       AliITSUSeed* seed = fCurrHyp->RemoveSeed(ilr,isd);
1462       if (!seed) continue; // already discarded or saved
1463       if (seed->IsSaved()) dest->AddSeed(seed,ilr);
1464       else MarkSeedFree(seed);
1465     }
1466   }
1467   //
1468   //  AliInfo(Form("SeedsPool: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
1469 }
1470
1471 //__________________________________________________________________
1472 void AliITSUTrackerGlo::FlagSplitClusters()
1473 {
1474   // set special bit on split clusters using MC info
1475   for (int ilr=fNLrActive;ilr--;) {
1476     int nsplit=0;
1477     AliITSURecoLayer* lr = fITS->GetLayerActive(ilr);
1478     for (int isn=lr->GetNSensors();isn--;) {
1479       AliITSURecoSens* sens = lr->GetSensor(isn);
1480       int nCl = sens->GetNClusters();
1481       if (!nCl) continue;
1482       int cl0 = sens->GetFirstClusterId();
1483       for (int icl=nCl;icl--;) {
1484         AliITSUClusterPix *cl = (AliITSUClusterPix*)lr->GetCluster(cl0+icl);
1485         for (int icl1=icl;icl1--;) {
1486           AliITSUClusterPix *cl1 = (AliITSUClusterPix*)lr->GetCluster(cl0+icl1);
1487           if (cl->HasCommonTrack(cl1)) {
1488             if (!cl->IsSplit())  nsplit++;
1489             if (!cl1->IsSplit()) nsplit++;
1490             cl->SetSplit();
1491             cl1->SetSplit();
1492           }
1493         }
1494       }
1495     }
1496     AliInfo(Form("%4d out of %4d clusters are split",nsplit,lr->GetNClusters()));
1497   }
1498   //
1499 }
1500
1501 //__________________________________________________________________
1502 Bool_t AliITSUTrackerGlo::ContainsSplitCluster(const AliITSUSeed* seed, Int_t maxSize)
1503 {
1504   // check if the seed contains split cluster with size < maxSize
1505   int lrID,clID;
1506   if ( (clID=seed->GetLrCluster(lrID))>=0 ) {
1507     AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1508     if (cl->IsSplit() && cl->GetNPix()<maxSize ) return kTRUE;
1509   }
1510   return seed->GetParent() ? ContainsSplitCluster((AliITSUSeed*)seed->GetParent(),maxSize) : kFALSE;
1511   //
1512 }
1513
1514 //__________________________________________________________________
1515 void AliITSUTrackerGlo::PrintSeedClusters(const AliITSUSeed* seed, Option_t* option)
1516 {
1517   // print seeds clusters
1518   int lrID,clID;
1519   if ( (clID=seed->GetLrCluster(lrID))>=0 ) {
1520     AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1521     cl->Print(option);
1522   }
1523   if (seed->GetParent()) PrintSeedClusters((AliITSUSeed*)seed->GetParent(), option);
1524   //
1525 }
1526
1527 //__________________________________________________________________
1528 void AliITSUTrackerGlo::FlagSeedClusters(const AliITSUSeed* seed, Bool_t flg)
1529 {
1530   // mark used clusters
1531   int lrID,clID;
1532   while (seed) {
1533     if ( (clID=seed->GetLrCluster(lrID))>=0 ) fITS->GetLayerActive(lrID)->GetCluster(clID)->SetBit(AliCluster::kUsed,flg);
1534     seed = (AliITSUSeed*)seed->GetParent();
1535   }
1536   //
1537 }
1538
1539
1540
1541 #ifdef  _FILL_CONTROL_HISTOS_
1542 //__________________________________________________________________
1543 void AliITSUTrackerGlo::BookControlHistos()
1544 {
1545   // book special control histos
1546   if (!fCHistoArr) { // create control histos
1547     const int kNResDef=7;
1548     const double kResDef[kNResDef]={0.05,0.05,0.3, 0.05,1,0.5,1.5};
1549     fCHistoArr = new TObjArray();
1550     fCHistoArr->SetOwner(kTRUE);
1551     const double ptMax=10;
1552     const double plMax=10;
1553     const double chiMax=100;
1554     const int nptbins=50;
1555     const int nresbins=400;
1556     const int nplbins=50;
1557     const int nchbins=200;
1558     const int maxBr  = 15;
1559     const int maxCand = 200;
1560     TString ttl;
1561     for (int stp=0;stp<kNTrackingPhases;stp++) {
1562       for (int ilr=0;ilr<fNLrActive;ilr++) {
1563         int hoffs = stp*kHistosPhase + ilr;
1564         double mxdf = ilr>=kNResDef ? kResDef[kNResDef-1] : kResDef[ilr];
1565         ttl = Form("S%d_residY%d",stp,ilr);
1566         TH2F* hdy = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nresbins,-mxdf,mxdf);
1567         fCHistoArr->AddAtAndExpand(hdy,hoffs + kHResY);
1568         hdy->SetDirectory(0);
1569         //
1570         ttl = Form("S%d_residYPull%d",stp,ilr); 
1571         TH2F* hdyp = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nplbins,-plMax,plMax);
1572         fCHistoArr->AddAtAndExpand(hdyp,hoffs + kHResYP);
1573         hdyp->SetDirectory(0);
1574         //
1575         ttl = Form("S%d_residZ%d",stp,ilr);     
1576         TH2F* hdz = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nresbins,-mxdf,mxdf);
1577         fCHistoArr->AddAtAndExpand(hdz,hoffs + kHResZ);
1578         hdz->SetDirectory(0);
1579         //
1580         ttl = Form("S%d_residZPull%d",stp,ilr);         
1581         TH2F* hdzp = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nplbins,-plMax,plMax);
1582         hdzp->SetDirectory(0);
1583         fCHistoArr->AddAtAndExpand(hdzp,hoffs + kHResZP);
1584         //
1585         ttl = Form("S%d_chi2Cl%d",stp,ilr);             
1586         TH2F* hchi = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins,0.,chiMax);
1587         hchi->SetDirectory(0);
1588         fCHistoArr->AddAtAndExpand(hchi,hoffs + kHChi2Cl);
1589         //
1590         ttl = Form("S%d_chi2Nrm%d",stp,ilr);            
1591         TH2F* hchiN = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins,0.,chiMax);
1592         hchiN->SetDirectory(0);
1593         fCHistoArr->AddAtAndExpand(hchiN,hoffs + kHChi2Nrm);
1594         //
1595         if (stp==0) { // these histos make sense only for clusters2tracks stage
1596           ttl = Form("S%d_bestInBranch%d",stp,ilr);             
1597           TH2* hnbr = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, maxBr,-0.5,maxBr-0.5);
1598           hnbr->SetDirectory(0);
1599           fCHistoArr->AddAtAndExpand(hnbr,hoffs + kHBestInBranch);
1600           //
1601           ttl = Form("S%d_bestInCands%d",stp,ilr);              
1602           TH2* hncn = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, maxCand,-0.5,maxCand-0.5);
1603           hncn->SetDirectory(0);
1604           fCHistoArr->AddAtAndExpand(hncn,hoffs + kHBestInCand);
1605           //
1606         }
1607       }
1608     }
1609   }
1610   //  
1611 }
1612 #endif