All modifications to track reconstruction removed with the exception of the flag...
[u/mrichter/AliRoot.git] / ITS / AliITStrackerMI.cxx
1 /**************************************************************************
2  * Copyright(c) 2007-2009, 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 /* $Id$ */
16
17 //-------------------------------------------------------------------------
18 //               Implementation of the ITS tracker class
19 //    It reads AliITSRecPoint clusters and creates AliITStrackMI tracks
20 //                   and fills with them the ESD
21 //          Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch 
22 //          Current support and development: 
23 //                     Andrea Dainese, andrea.dainese@lnl.infn.it
24 //     dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
25 //     Params moved to AliITSRecoParam by: Andrea Dainese, INFN
26 //     Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
27 //-------------------------------------------------------------------------
28
29 #include <TMatrixD.h>
30 #include <TTree.h>
31 #include <TDatabasePDG.h>
32 #include <TString.h>
33 #include <TRandom.h>
34 #include <TTreeStream.h>
35 #include <TVector3.h>
36
37 #include "AliLog.h"
38 #include "AliGeomManager.h"
39 #include "AliITSPlaneEff.h"
40 #include "AliTrackPointArray.h"
41 #include "AliESDEvent.h"
42 #include "AliESDtrack.h"
43 #include "AliV0.h"
44 #include "AliITSChannelStatus.h"
45 #include "AliITSDetTypeRec.h"
46 #include "AliITSRecPoint.h"
47 #include "AliITSRecPointContainer.h"
48 #include "AliITSgeomTGeo.h"
49 #include "AliITSReconstructor.h"
50 #include "AliITSClusterParam.h"
51 #include "AliITSsegmentation.h"
52 #include "AliITSCalibration.h"
53 #include "AliITSPlaneEffSPD.h"
54 #include "AliITSPlaneEffSDD.h"
55 #include "AliITSPlaneEffSSD.h"
56 #include "AliITSV0Finder.h"
57 #include "AliITStrackerMI.h"
58 #include "AliMathBase.h"
59
60
61 ClassImp(AliITStrackerMI)
62
63 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
64
65 AliITStrackerMI::AliITStrackerMI():AliTracker(),
66 fI(0),
67 fBestTrack(),
68 fTrackToFollow(),
69 fTrackHypothesys(),
70 fBestHypothesys(),
71 fOriginal(),
72 fCurrentEsdTrack(),
73 fPass(0),
74 fAfterV0(kFALSE),
75 fLastLayerToTrackTo(0),
76 fCoefficients(0),
77 fEsd(0),
78 fTrackingPhase("Default"),
79 fUseTGeo(3),
80 fNtracks(0),
81 fFlagFakes(kFALSE),
82 fSelectBestMIP03(kFALSE),
83 fUseImproveKalman(kFALSE),
84 fxOverX0Pipe(-1.),
85 fxTimesRhoPipe(-1.),
86 fxOverX0PipeTrks(0),
87 fxTimesRhoPipeTrks(0),
88 fxOverX0ShieldTrks(0),
89 fxTimesRhoShieldTrks(0),
90 fxOverX0LayerTrks(0),
91 fxTimesRhoLayerTrks(0),
92 fDebugStreamer(0),
93 fITSChannelStatus(0),
94 fkDetTypeRec(0),
95 fPlaneEff(0),
96 fITSPid(0)
97
98  {
99   //Default constructor
100   Int_t i;
101   for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
102   for(i=0;i<2;i++) {
103     fxOverX0Shield[i]=-1.;
104     fxTimesRhoShield[i]=-1.;
105     fConstraint[i]=0;
106   }
107   for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
108   fOriginal.SetOwner();
109   for(i=0;i<AliITSgeomTGeo::kNLayers;i++)fForceSkippingOfLayer[i]=0;
110   for(i=0;i<100000;i++)fBestTrackIndex[i]=0;
111   fITSPid=new AliITSPIDResponse();
112
113 }
114 //------------------------------------------------------------------------
115 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
116 fI(AliITSgeomTGeo::GetNLayers()),
117 fBestTrack(),
118 fTrackToFollow(),
119 fTrackHypothesys(),
120 fBestHypothesys(),
121 fOriginal(),
122 fCurrentEsdTrack(),
123 fPass(0),
124 fAfterV0(kFALSE),
125 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
126 fCoefficients(0),
127 fEsd(0),
128 fTrackingPhase("Default"),
129 fUseTGeo(3),
130 fNtracks(0),
131 fFlagFakes(kFALSE),
132 fSelectBestMIP03(kFALSE),
133 fUseImproveKalman(kFALSE),
134 fxOverX0Pipe(-1.),
135 fxTimesRhoPipe(-1.),
136 fxOverX0PipeTrks(0),
137 fxTimesRhoPipeTrks(0),
138 fxOverX0ShieldTrks(0),
139 fxTimesRhoShieldTrks(0),
140 fxOverX0LayerTrks(0),
141 fxTimesRhoLayerTrks(0),
142 fDebugStreamer(0),
143 fITSChannelStatus(0),
144 fkDetTypeRec(0),
145 fPlaneEff(0),
146 fITSPid(0) {
147   //--------------------------------------------------------------------
148   //This is the AliITStrackerMI constructor
149   //--------------------------------------------------------------------
150   if (geom) {
151     AliWarning("\"geom\" is actually a dummy argument !");
152   }
153
154   fOriginal.SetOwner();
155   fCoefficients = 0;
156   fAfterV0     = kFALSE;
157
158   for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
159     Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
160     Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
161
162     Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
163     AliITSgeomTGeo::GetTranslation(i,1,1,xyz); 
164     Double_t poff=TMath::ATan2(y,x);
165     Double_t zoff=z;
166
167     AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
168     Double_t r=TMath::Sqrt(x*x + y*y);
169
170     AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
171     r += TMath::Sqrt(x*x + y*y);
172     AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
173     r += TMath::Sqrt(x*x + y*y);
174     AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
175     r += TMath::Sqrt(x*x + y*y);
176     r*=0.25;
177
178     new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
179
180     for (Int_t j=1; j<nlad+1; j++) {
181       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
182         TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
183         const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
184         m.Multiply(tm);
185         Double_t txyz[3]={0.};
186         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
187         m.LocalToMaster(txyz,xyz);
188         r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
189         Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
190
191         if (phi<0) phi+=TMath::TwoPi();
192         else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
193
194         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1); 
195         new(&det) AliITSdetector(r,phi); 
196         // compute the real radius (with misalignment)
197         TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
198         mmisal.Multiply(tm);
199         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
200         mmisal.LocalToMaster(txyz,xyz);
201         Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
202         det.SetRmisal(rmisal);
203         
204       } // end loop on detectors
205     } // end loop on ladders
206     fForceSkippingOfLayer[i-1] = 0;
207   } // end loop on layers
208
209
210   fI=AliITSgeomTGeo::GetNLayers();
211
212   fPass=0;
213   fConstraint[0]=1; fConstraint[1]=0;
214
215   Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
216                      AliITSReconstructor::GetRecoParam()->GetYVdef(),
217                      AliITSReconstructor::GetRecoParam()->GetZVdef()}; 
218   Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
219                      AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
220                      AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()}; 
221   SetVertex(xyzVtx,ersVtx);
222
223   fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
224   for (Int_t i=0;i<100000;i++){
225     fBestTrackIndex[i]=0;
226   }
227
228   // store positions of centre of SPD modules (in z)
229   //  The convetion is that fSPDdetzcentre is ordered from -z to +z
230   Double_t tr[3];
231   AliITSgeomTGeo::GetTranslation(1,1,1,tr);
232   if (tr[2]<0) { // old geom (up to v5asymmPPR)
233     AliITSgeomTGeo::GetTranslation(1,1,1,tr);
234     fSPDdetzcentre[0] = tr[2];
235     AliITSgeomTGeo::GetTranslation(1,1,2,tr);
236     fSPDdetzcentre[1] = tr[2];
237     AliITSgeomTGeo::GetTranslation(1,1,3,tr);
238     fSPDdetzcentre[2] = tr[2];
239     AliITSgeomTGeo::GetTranslation(1,1,4,tr);
240     fSPDdetzcentre[3] = tr[2];
241   } else { // new geom (from v11Hybrid)
242     AliITSgeomTGeo::GetTranslation(1,1,4,tr);
243     fSPDdetzcentre[0] = tr[2];
244     AliITSgeomTGeo::GetTranslation(1,1,3,tr);
245     fSPDdetzcentre[1] = tr[2];
246     AliITSgeomTGeo::GetTranslation(1,1,2,tr);
247     fSPDdetzcentre[2] = tr[2];
248     AliITSgeomTGeo::GetTranslation(1,1,1,tr);
249     fSPDdetzcentre[3] = tr[2];
250   }
251
252   fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
253   if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
254     AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
255     fUseTGeo = 3;
256   }
257
258   for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
259   for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
260   
261   if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0)
262     fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
263
264   // only for plane efficiency evaluation
265   if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
266       AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
267     Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
268     if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
269       AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
270     if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
271     else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
272     else fPlaneEff = new AliITSPlaneEffSSD();
273     if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
274        if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
275     if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
276   }
277   //
278   // RS
279   fSelectBestMIP03  = kFALSE;//AliITSReconstructor::GetRecoParam()->GetSelectBestMIP03();
280   fFlagFakes        = AliITSReconstructor::GetRecoParam()->GetFlagFakes();
281   fUseImproveKalman = AliITSReconstructor::GetRecoParam()->GetUseImproveKalman();
282   //
283   fITSPid=new AliITSPIDResponse();
284 }
285 /*
286 //------------------------------------------------------------------------
287 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
288 fI(tracker.fI),
289 fBestTrack(tracker.fBestTrack),
290 fTrackToFollow(tracker.fTrackToFollow),
291 fTrackHypothesys(tracker.fTrackHypothesys),
292 fBestHypothesys(tracker.fBestHypothesys),
293 fOriginal(tracker.fOriginal),
294 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
295 fPass(tracker.fPass),
296 fAfterV0(tracker.fAfterV0),
297 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
298 fCoefficients(tracker.fCoefficients),
299 fEsd(tracker.fEsd),
300 fTrackingPhase(tracker.fTrackingPhase),
301 fUseTGeo(tracker.fUseTGeo),
302 fNtracks(tracker.fNtracks),
303 fFlagFakes(tracker.fFlagFakes),
304 fSelectBestMIP03(tracker.fSelectBestMIP03),
305 fxOverX0Pipe(tracker.fxOverX0Pipe),
306 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
307 fxOverX0PipeTrks(0),
308 fxTimesRhoPipeTrks(0),
309 fxOverX0ShieldTrks(0),
310 fxTimesRhoShieldTrks(0),
311 fxOverX0LayerTrks(0),
312 fxTimesRhoLayerTrks(0),
313 fDebugStreamer(tracker.fDebugStreamer),
314 fITSChannelStatus(tracker.fITSChannelStatus),
315 fkDetTypeRec(tracker.fkDetTypeRec),
316 fPlaneEff(tracker.fPlaneEff) {
317   //Copy constructor
318   fOriginal.SetOwner();
319   Int_t i;
320   for(i=0;i<4;i++) {
321     fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
322   }
323   for(i=0;i<6;i++) {
324     fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
325     fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
326   }
327   for(i=0;i<2;i++) {
328     fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
329     fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
330   }
331 }
332
333 //------------------------------------------------------------------------
334 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
335   //Assignment operator
336   this->~AliITStrackerMI();
337   new(this) AliITStrackerMI(tracker);
338   return *this;
339 }
340 */
341 //------------------------------------------------------------------------
342 AliITStrackerMI::~AliITStrackerMI()
343 {
344   //
345   //destructor
346   //
347   if (fCoefficients) delete [] fCoefficients;
348   DeleteTrksMaterialLUT();
349   if (fDebugStreamer) {
350     //fDebugStreamer->Close();
351     delete fDebugStreamer;
352   }
353   if(fITSChannelStatus) delete fITSChannelStatus;
354   if(fPlaneEff) delete fPlaneEff;
355   if(fITSPid) delete fITSPid;
356
357 }
358 //------------------------------------------------------------------------
359 void AliITStrackerMI::ReadBadFromDetTypeRec() {
360   //--------------------------------------------------------------------
361   //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
362   //i.e. from OCDB
363   //--------------------------------------------------------------------
364
365   if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
366
367   Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
368
369   if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
370
371   // ITS channels map
372   if(fITSChannelStatus) delete fITSChannelStatus;
373   fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
374
375   // ITS detectors and chips
376   Int_t i=0,j=0,k=0,ndet=0;
377   for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
378     Int_t nBadDetsPerLayer=0;
379     ndet=AliITSgeomTGeo::GetNDetectors(i);    
380     for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
381       for (k=1; k<ndet+1; k++) {
382         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);  
383         det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
384         if(det.IsBad()) {nBadDetsPerLayer++;}
385       } // end loop on detectors
386     } // end loop on ladders
387     AliInfo(Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
388   } // end loop on layers
389   
390   return;
391 }
392 //------------------------------------------------------------------------
393 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
394   //--------------------------------------------------------------------
395   //This function loads ITS clusters
396   //--------------------------------------------------------------------
397  
398   TClonesArray *clusters = NULL;
399   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
400   clusters=rpcont->FetchClusters(0,cTree);
401   if(!clusters) return 1;
402
403   if(!(rpcont->IsSPDActive() || rpcont->IsSDDActive() || rpcont->IsSSDActive())){
404       AliError("ITS is not in a known running configuration: SPD, SDD and SSD are not active");
405       return 1;
406   }
407   Int_t i=0,j=0,ndet=0;
408   Int_t detector=0;
409   for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
410     ndet=fgLayers[i].GetNdetectors();
411     Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
412     for (; j<jmax; j++) {           
413       //      if (!cTree->GetEvent(j)) continue;
414       clusters = rpcont->UncheckedGetClusters(j);
415       if(!clusters)continue;
416       Int_t ncl=clusters->GetEntriesFast();
417       SignDeltas(clusters,GetZ());
418  
419       while (ncl--) {
420         AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
421         detector=c->GetDetectorIndex();
422
423         if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
424         
425         Int_t retval = fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
426         if(retval) {
427           AliWarning(Form("Too many clusters on layer %d!",i));
428           break;  
429         } 
430       }
431
432       // add dead zone "virtual" cluster in SPD, if there is a cluster within 
433       // zwindow cm from the dead zone      
434       //  This method assumes that fSPDdetzcentre is ordered from -z to +z
435       if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
436         for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
437           Int_t lab[4]   = {0,0,0,detector};
438           Int_t info[3]  = {0,0,i};
439           Float_t q      = 0.; // this identifies virtual clusters
440           Float_t hit[6] = {xdead,
441                             0.,
442                             AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
443                             AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
444                             q,
445                             0.};
446           Bool_t local   = kTRUE;
447           Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
448           hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
449           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow) 
450             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
451           hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
452           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow) 
453             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
454           hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
455           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow) 
456             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
457           hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
458           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow) 
459             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
460           hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
461           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow) 
462             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
463           hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
464           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow) 
465             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
466         }
467       } // "virtual" clusters in SPD
468       
469     }
470     //
471     fgLayers[i].ResetRoad(); //road defined by the cluster density
472     fgLayers[i].SortClusters();
473   }
474
475   // check whether we have to skip some layers
476   SetForceSkippingOfLayer();
477
478   return 0;
479 }
480 //------------------------------------------------------------------------
481 void AliITStrackerMI::UnloadClusters() {
482   //--------------------------------------------------------------------
483   //This function unloads ITS clusters
484   //--------------------------------------------------------------------
485   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
486 }
487 //------------------------------------------------------------------------
488 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
489   //--------------------------------------------------------------------
490   // Publishes all pointers to clusters known to the tracker into the
491   // passed object array.
492   // The ownership is not transfered - the caller is not expected to delete
493   // the clusters.
494   //--------------------------------------------------------------------
495
496   for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
497     for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
498       AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
499       array->AddLast(cl);
500     }
501   }
502
503   return;
504 }
505 //------------------------------------------------------------------------
506 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
507   //--------------------------------------------------------------------
508   // Correction for the material between the TPC and the ITS
509   //--------------------------------------------------------------------
510   if (t->GetX() > AliITSRecoParam::Getriw()) {   // inward direction 
511       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
512       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
513       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))  return 0;// ITS screen
514   } else if (t->GetX() < AliITSRecoParam::Getrs()) {  // outward direction
515       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))        return 0;// ITS screen
516       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1))       return 0;// TPC central drum
517       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
518   } else {
519     printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
520     return 0;
521   }
522   
523   return 1;
524 }
525 //------------------------------------------------------------------------
526 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
527   //--------------------------------------------------------------------
528   // This functions reconstructs ITS tracks
529   // The clusters must be already loaded !
530   //--------------------------------------------------------------------
531
532   AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
533
534   fTrackingPhase="Clusters2Tracks";
535   //
536   // RS
537   fSelectBestMIP03  = kFALSE;//AliITSReconstructor::GetRecoParam()->GetSelectBestMIP03();
538   fFlagFakes        = AliITSReconstructor::GetRecoParam()->GetFlagFakes();
539   fUseImproveKalman = AliITSReconstructor::GetRecoParam()->GetUseImproveKalman();
540   //
541   TObjArray itsTracks(15000);
542   fOriginal.Clear();
543   fEsd = event;         // store pointer to the esd 
544
545   // temporary (for cosmics)
546   if(event->GetVertex()) {
547     TString title = event->GetVertex()->GetTitle();
548     if(title.Contains("cosmics")) {
549       Double_t xyz[3]={GetX(),GetY(),GetZ()};
550       Double_t exyz[3]={0.1,0.1,0.1};
551       SetVertex(xyz,exyz);
552     }
553   }
554   // temporary
555   Int_t noesd = 0;
556   {/* Read ESD tracks */
557     Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
558     Int_t nentr=event->GetNumberOfTracks();
559     noesd=nentr;
560     //    Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
561     while (nentr--) {
562       AliESDtrack *esd=event->GetTrack(nentr);
563       //  ---- for debugging:
564       //if(TMath::Abs(esd->GetX()-83.65)<0.1) { FILE *f=fopen("tpc.dat","a"); fprintf(f,"%f %f %f %f %f %f\n",(Float_t)event->GetEventNumberInFile(),(Float_t)TMath::Abs(esd->GetLabel()),(Float_t)esd->GetX(),(Float_t)esd->GetY(),(Float_t)esd->GetZ(),(Float_t)esd->Pt()); fclose(f); }
565
566       if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
567       if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
568       if (esd->GetStatus()&AliESDtrack::kITSin) continue;
569       if (esd->GetKinkIndex(0)>0) continue;   //kink daughter
570       AliITStrackMI *t = new AliITStrackMI(*esd);
571       t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP());              //I.B.
572       Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
573
574
575       // look at the ESD mass hypothesys !
576       if (t->GetMass()<0.9*pimass) t->SetMass(pimass); 
577       // write expected q
578       t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
579
580       if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
581         //track - can be  V0 according to TPC
582       } else {  
583         if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
584           delete t;
585           continue;
586         }       
587         if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
588           delete t;
589           continue;
590         }
591         if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
592           delete t;
593           continue;
594         }
595         if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
596           delete t;
597           continue;
598         }
599       }
600       t->SetReconstructed(kFALSE);
601       itsTracks.AddLast(t);
602       fOriginal.AddLast(t);
603     }
604   } /* End Read ESD tracks */
605
606   itsTracks.Sort();
607   fOriginal.Sort();
608   Int_t nentr=itsTracks.GetEntriesFast();
609   fTrackHypothesys.Expand(nentr);
610   fBestHypothesys.Expand(nentr);
611   MakeCoefficients(nentr);
612   if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
613   Int_t ntrk=0;
614   // THE TWO TRACKING PASSES
615   for (fPass=0; fPass<2; fPass++) {
616      Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
617      for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
618        AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
619        if (t==0) continue;              //this track has been already tracked
620        //cout<<"========== "<<fPass<<"    "<<fCurrentEsdTrack<<" =========\n";
621        if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue;  //this track was  already  "succesfully" reconstructed
622        Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz);              //I.B.
623        if (fConstraint[fPass]) { 
624          if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
625              TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
626        }
627
628        Int_t tpcLabel=t->GetLabel(); //save the TPC track label       
629        AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
630        fI = 6;
631        ResetTrackToFollow(*t);
632        ResetBestTrack();
633
634        FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
635  
636
637        SortTrackHypothesys(fCurrentEsdTrack,20,0);  //MI change
638        //
639        AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
640        if (!besttrack) continue;
641        besttrack->SetLabel(tpcLabel);
642        //       besttrack->CookdEdx();
643        CookdEdx(besttrack);
644        besttrack->SetFakeRatio(1.);
645        CookLabel(besttrack,0.); //For comparison only
646        UpdateESDtrack(besttrack,AliESDtrack::kITSin);
647        t->SetWinner(besttrack);
648
649        if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue;  //to be tracked also without vertex constrain 
650
651        t->SetReconstructed(kTRUE);
652        ntrk++;  
653        AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
654      }
655      GetBestHypothesysMIP(itsTracks); 
656   } // end loop on the two tracking passes
657   //
658   if (fFlagFakes) FlagFakes(itsTracks);
659   //
660   if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
661   if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
662   fAfterV0 = kTRUE;
663   //
664   itsTracks.Clear();
665   //
666   Int_t entries = fTrackHypothesys.GetEntriesFast();
667   for (Int_t ientry=0; ientry<entries; ientry++) {
668     TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
669     if (array) array->Delete();
670     delete fTrackHypothesys.RemoveAt(ientry); 
671   }
672
673   fTrackHypothesys.Delete();
674   entries = fBestHypothesys.GetEntriesFast();
675   for (Int_t ientry=0; ientry<entries; ientry++) {
676     TObjArray * array =(TObjArray*)fBestHypothesys.UncheckedAt(ientry);
677     if (array) array->Delete();
678     delete fBestHypothesys.RemoveAt(ientry);
679   }
680   fBestHypothesys.Delete();
681   fOriginal.Clear();
682   delete [] fCoefficients;
683   fCoefficients=0;
684   DeleteTrksMaterialLUT();
685
686   AliInfo(Form("Number of prolonged tracks: %d out of %d ESD tracks",ntrk,noesd));
687
688   fTrackingPhase="Default";
689   
690   return 0;
691 }
692 //------------------------------------------------------------------------
693 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
694   //--------------------------------------------------------------------
695   // This functions propagates reconstructed ITS tracks back
696   // The clusters must be loaded !
697   //--------------------------------------------------------------------
698   fTrackingPhase="PropagateBack";
699   Int_t nentr=event->GetNumberOfTracks();
700   //  Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
701
702   Int_t ntrk=0;
703   for (Int_t i=0; i<nentr; i++) {
704      AliESDtrack *esd=event->GetTrack(i);
705
706      if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
707      if (esd->GetStatus()&AliESDtrack::kITSout) continue;
708
709      AliITStrackMI *t = new AliITStrackMI(*esd);
710
711      t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
712
713      ResetTrackToFollow(*t);
714
715      /*
716      // propagate to vertex [SR, GSI 17.02.2003]
717      // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
718      if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
719        if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
720          fTrackToFollow.StartTimeIntegral();
721        // from vertex to outside pipe
722        CorrectForPipeMaterial(&fTrackToFollow,"outward");
723        }*/
724      // Start time integral and add distance from current position to vertex 
725      Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
726      fTrackToFollow.GetXYZ(xyzTrk);
727      Double_t dst2 = 0.;
728      for (Int_t icoord=0; icoord<3; icoord++) {  
729        Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
730        dst2 += di*di; 
731      }
732      fTrackToFollow.StartTimeIntegral();
733      fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
734
735      fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
736      if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
737         if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
738           delete t;
739           continue;
740         }
741         fTrackToFollow.SetLabel(t->GetLabel());
742         //fTrackToFollow.CookdEdx();
743         CookLabel(&fTrackToFollow,0.); //For comparison only
744         fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
745         //UseClusters(&fTrackToFollow);
746         ntrk++;
747      }
748      delete t;
749   }
750
751   AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
752
753   fTrackingPhase="Default";
754
755   return 0;
756 }
757 //------------------------------------------------------------------------
758 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
759   //--------------------------------------------------------------------
760   // This functions refits ITS tracks using the 
761   // "inward propagated" TPC tracks
762   // The clusters must be loaded !
763   //--------------------------------------------------------------------
764   fTrackingPhase="RefitInward";
765
766   if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
767
768   Bool_t doExtra=AliITSReconstructor::GetRecoParam()->GetSearchForExtraClusters();
769   if(!doExtra) AliDebug(2,"Do not search for extra clusters");
770
771   Int_t nentr=event->GetNumberOfTracks();
772   //  Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
773
774   Int_t ntrk=0;
775   for (Int_t i=0; i<nentr; i++) {
776     AliESDtrack *esd=event->GetTrack(i);
777
778     if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
779     if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
780     if (esd->GetStatus()&AliESDtrack::kTPCout)
781       if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
782
783     AliITStrackMI *t = new AliITStrackMI(*esd);
784
785     t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
786     if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
787        delete t;
788        continue;
789     }
790
791     ResetTrackToFollow(*t);
792     fTrackToFollow.ResetClusters();
793
794     // ITS standalone tracks
795     if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) {
796       fTrackToFollow.ResetCovariance(10.);
797       // protection for loopers that can have parameters screwed up
798       if(TMath::Abs(fTrackToFollow.GetY())>1000. ||
799          TMath::Abs(fTrackToFollow.GetZ())>1000.) {
800         delete t;
801         continue;
802       }
803     }
804
805     //Refitting...
806     Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
807                AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
808
809     AliDebug(2,Form("Refit LABEL %d  %d",t->GetLabel(),t->GetNumberOfClusters()));
810     if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,doExtra,pe)) {
811        AliDebug(2,"  refit OK");
812        fTrackToFollow.SetLabel(t->GetLabel());
813        //       fTrackToFollow.CookdEdx();
814        CookdEdx(&fTrackToFollow);
815
816        CookLabel(&fTrackToFollow,0.0); //For comparison only
817
818        //The beam pipe
819        if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
820          fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
821          AliESDtrack  *esdTrack =fTrackToFollow.GetESDtrack();
822          //printf("                                       %d\n",esdTrack->GetITSModuleIndex(0));
823          //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
824          Double_t r[3]={0.,0.,0.};
825          Double_t maxD=3.;
826          esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
827          ntrk++;
828        }
829     }
830     delete t;
831   }
832
833   AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
834
835   fTrackingPhase="Default";
836
837   return 0;
838 }
839 //------------------------------------------------------------------------
840 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
841   //--------------------------------------------------------------------
842   //       Return pointer to a given cluster
843   //--------------------------------------------------------------------
844   Int_t l=(index & 0xf0000000) >> 28;
845   Int_t c=(index & 0x0fffffff) >> 00;
846   return fgLayers[l].GetCluster(c);
847 }
848 //------------------------------------------------------------------------
849 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
850   //--------------------------------------------------------------------
851   // Get track space point with index i
852   //--------------------------------------------------------------------
853
854   Int_t l=(index & 0xf0000000) >> 28;
855   Int_t c=(index & 0x0fffffff) >> 00;
856   AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
857   Int_t idet = cl->GetDetectorIndex();
858
859   Float_t xyz[3];
860   Float_t cov[6];
861   cl->GetGlobalXYZ(xyz);
862   cl->GetGlobalCov(cov);
863   p.SetXYZ(xyz, cov);
864   p.SetCharge(cl->GetQ());
865   p.SetDriftTime(cl->GetDriftTime());
866   p.SetChargeRatio(cl->GetChargeRatio());
867   p.SetClusterType(cl->GetClusterType());
868   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer; 
869   switch (l) {
870   case 0:
871     iLayer = AliGeomManager::kSPD1;
872     break;
873   case 1:
874     iLayer = AliGeomManager::kSPD2;
875     break;
876   case 2:
877     iLayer = AliGeomManager::kSDD1;
878     break;
879   case 3:
880     iLayer = AliGeomManager::kSDD2;
881     break;
882   case 4:
883     iLayer = AliGeomManager::kSSD1;
884     break;
885   case 5:
886     iLayer = AliGeomManager::kSSD2;
887     break;
888   default:
889     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
890     break;
891   };
892   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
893   p.SetVolumeID((UShort_t)volid);
894   return kTRUE;
895 }
896 //------------------------------------------------------------------------
897 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index, 
898                         AliTrackPoint& p, const AliESDtrack *t) {
899   //--------------------------------------------------------------------
900   // Get track space point with index i
901   // (assign error estimated during the tracking)
902   //--------------------------------------------------------------------
903
904   Int_t l=(index & 0xf0000000) >> 28;
905   Int_t c=(index & 0x0fffffff) >> 00;
906   const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
907   Int_t idet = cl->GetDetectorIndex();
908
909   const AliITSdetector &det=fgLayers[l].GetDetector(idet);
910
911   // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
912   Float_t detxy[2];
913   detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
914   detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
915   Double_t alpha = t->GetAlpha();
916   Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
917   Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe+cl->GetX(),GetBz()));
918   phi += alpha-det.GetPhi();
919   Float_t tgphi = TMath::Tan(phi);
920
921   Float_t tgl = t->GetTgl(); // tgl about const along track
922   Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
923
924   Float_t errtrky,errtrkz,covyz;
925   Bool_t addMisalErr=kFALSE;
926   AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
927
928   Float_t xyz[3];
929   Float_t cov[6];
930   cl->GetGlobalXYZ(xyz);
931   //  cl->GetGlobalCov(cov);
932   Float_t pos[3] = {0.,0.,0.};
933   AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
934   tmpcl.GetGlobalCov(cov);
935
936   p.SetXYZ(xyz, cov);
937   p.SetCharge(cl->GetQ());
938   p.SetDriftTime(cl->GetDriftTime());
939   p.SetChargeRatio(cl->GetChargeRatio());
940   p.SetClusterType(cl->GetClusterType());
941
942   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer; 
943   switch (l) {
944   case 0:
945     iLayer = AliGeomManager::kSPD1;
946     break;
947   case 1:
948     iLayer = AliGeomManager::kSPD2;
949     break;
950   case 2:
951     iLayer = AliGeomManager::kSDD1;
952     break;
953   case 3:
954     iLayer = AliGeomManager::kSDD2;
955     break;
956   case 4:
957     iLayer = AliGeomManager::kSSD1;
958     break;
959   case 5:
960     iLayer = AliGeomManager::kSSD2;
961     break;
962   default:
963     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
964     break;
965   };
966   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
967
968   p.SetVolumeID((UShort_t)volid);
969   return kTRUE;
970 }
971 //------------------------------------------------------------------------
972 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain) 
973 {
974   //--------------------------------------------------------------------
975   // Follow prolongation tree
976   //--------------------------------------------------------------------
977   //
978   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
979   Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
980
981
982   AliESDtrack * esd = otrack->GetESDtrack();
983   if (esd->GetV0Index(0)>0) {
984     // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
985     //                      mapping of ESD track is different as ITS track in Containers
986     //                      Need something more stable
987     //                      Indexes are set back again to the ESD track indexes in UpdateTPCV0
988     for (Int_t i=0;i<3;i++){
989       Int_t  index = esd->GetV0Index(i);
990       if (index==0) break;
991       AliESDv0 * vertex = fEsd->GetV0(index);
992       if (vertex->GetStatus()<0) continue;     // rejected V0
993       //
994       if (esd->GetSign()>0) {
995         vertex->SetIndex(0,esdindex);
996       } else {
997         vertex->SetIndex(1,esdindex);
998       }
999     }
1000   }
1001   TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
1002   if (!bestarray){
1003     bestarray = new TObjArray(5);
1004     bestarray->SetOwner();
1005     fBestHypothesys.AddAt(bestarray,esdindex);
1006   }
1007
1008   //
1009   //setup tree of the prolongations
1010   //
1011   const int kMaxTr = 100; //RS
1012   static AliITStrackMI tracks[7][kMaxTr];
1013   AliITStrackMI *currenttrack;
1014   static AliITStrackMI currenttrack1;
1015   static AliITStrackMI currenttrack2;  
1016   static AliITStrackMI backuptrack;
1017   Int_t ntracks[7];
1018   Int_t nindexes[7][kMaxTr];
1019   Float_t normalizedchi2[kMaxTr];
1020   for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
1021   otrack->SetNSkipped(0);
1022   new (&(tracks[6][0])) AliITStrackMI(*otrack);
1023   ntracks[6]=1;
1024   for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
1025   Int_t modstatus = 1; // found 
1026   Float_t xloc,zloc;
1027   // 
1028   //
1029   // follow prolongations
1030   for (Int_t ilayer=5; ilayer>=0; ilayer--) {
1031     AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
1032     fI = ilayer;
1033     //
1034     AliITSlayer &layer=fgLayers[ilayer];
1035     Double_t r = layer.GetR(); 
1036     ntracks[ilayer]=0;
1037     //
1038     //
1039     Int_t nskipped=0;
1040     Float_t nused =0;
1041     for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
1042       //      printf("LR %d Tr:%d NSeeds: %d\n",ilayer,itrack,ntracks[ilayer+1]);
1043       //set current track
1044       if (ntracks[ilayer]>=kMaxTr) break;  
1045       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
1046       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
1047       if (ntracks[ilayer]>15+ilayer){
1048         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
1049         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
1050       }
1051
1052       new(&currenttrack1)  AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
1053   
1054       // material between SSD and SDD, SDD and SPD
1055       if (ilayer==3) 
1056         if(!CorrectForShieldMaterial(&currenttrack1,"SDD","inward")) continue;
1057       if (ilayer==1) 
1058         if(!CorrectForShieldMaterial(&currenttrack1,"SPD","inward")) continue;
1059
1060       // detector number
1061       Double_t phi,z;
1062       if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
1063       Int_t idet=layer.FindDetectorIndex(phi,z);
1064
1065       Double_t trackGlobXYZ1[3];
1066       if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1067
1068       // Get the budget to the primary vertex for the current track being prolonged
1069       Double_t budgetToPrimVertex = 0;
1070       double xMSLrs[9],x2X0MSLrs[9]; // needed for ImproveKalman
1071       int nMSLrs = 0;
1072       //
1073       if (fUseImproveKalman) nMSLrs = GetEffectiveThicknessLbyL(xMSLrs,x2X0MSLrs);
1074       else budgetToPrimVertex = GetEffectiveThickness();
1075       //
1076       // check if we allow a prolongation without point
1077       Int_t skip = CheckSkipLayer(&currenttrack1,ilayer,idet);
1078       if (skip) {
1079         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1080         // propagate to the layer radius
1081         Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1082         if(!vtrack->Propagate(xToGo)) continue;
1083         // apply correction for material of the current layer
1084         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1085         vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1086         vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1087         vtrack->SetClIndex(ilayer,-1);
1088         modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1089         if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1090           vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1091         }
1092         if(constrain && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1093           fUseImproveKalman ? 
1094             vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) : 
1095             vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1096         }
1097         ntracks[ilayer]++;
1098         continue;
1099       }
1100
1101       // track outside layer acceptance in z
1102       if (idet<0) continue;
1103       
1104       //propagate to the intersection with the detector plane
1105       const AliITSdetector &det=layer.GetDetector(idet);
1106       new(&currenttrack2)  AliITStrackMI(currenttrack1);
1107       if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1108       if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1109       currenttrack1.SetDetectorIndex(idet);
1110       currenttrack2.SetDetectorIndex(idet);
1111       if(!LocalModuleCoord(ilayer,idet,&currenttrack1,xloc,zloc)) continue; // local module coords
1112
1113       //***************
1114       // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1115       //
1116       // road in global (rphi,z) [i.e. in tracking ref. system]
1117       Double_t zmin,zmax,ymin,ymax;
1118       if (!ComputeRoad(&currenttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1119
1120       // select clusters in road
1121       layer.SelectClusters(zmin,zmax,ymin,ymax); 
1122       //********************
1123
1124       // Define criteria for track-cluster association
1125       Double_t msz = currenttrack1.GetSigmaZ2() + 
1126         AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1127         AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1128         AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1129       Double_t msy = currenttrack1.GetSigmaY2() + 
1130         AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1131         AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1132         AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1133       if (constrain) {
1134         msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1135         msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC(); 
1136       }  else {
1137         msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1138         msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC(); 
1139       }
1140       msz = 1./msz; // 1/RoadZ^2
1141       msy = 1./msy; // 1/RoadY^2
1142
1143       //
1144       //
1145       // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1146       //
1147       const AliITSRecPoint *cl=0; 
1148       Int_t clidx=-1;
1149       Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1150       Bool_t deadzoneSPD=kFALSE;
1151       currenttrack = &currenttrack1;
1152
1153       // check if the road contains a dead zone 
1154       Bool_t noClusters = kFALSE;
1155       if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1156       if (noClusters) AliDebug(2,"no clusters in road");
1157       Double_t dz=0.5*(zmax-zmin);
1158       Double_t dy=0.5*(ymax-ymin);
1159       Int_t dead = CheckDeadZone(&currenttrack1,ilayer,idet,dz,dy,noClusters); 
1160       if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1161       // create a prolongation without clusters (check also if there are no clusters in the road)
1162       if (dead || 
1163           (noClusters && 
1164            AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1165         AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1166         updatetrack->SetClIndex(ilayer,-1);
1167         if (dead==0) {
1168           modstatus = 5; // no cls in road
1169         } else if (dead==1) {
1170           modstatus = 7; // holes in z in SPD
1171         } else if (dead==2 || dead==3 || dead==4) {
1172           modstatus = 2; // dead from OCDB
1173         }
1174         updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1175         // apply correction for material of the current layer
1176         CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1177         if (constrain) { // apply vertex constrain
1178           updatetrack->SetConstrain(constrain);
1179           Bool_t isPrim = kTRUE;
1180           if (ilayer<4) { // check that it's close to the vertex
1181             updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1182             if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1183                 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() || 
1184                 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1185                 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1186           }
1187           if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1188             fUseImproveKalman ? 
1189               updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) : 
1190               updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1191           }
1192         }
1193         updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1194         if (dead) {
1195           if (dead==1) { // dead zone at z=0,+-7cm in SPD
1196             updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1197             deadzoneSPD=kTRUE;
1198           } else if (dead==2 || dead==3) { // dead module or chip from OCDB  
1199             updatetrack->SetDeadZoneProbability(ilayer,1.); 
1200           } else if (dead==4) { // at least a single dead channel from OCDB  
1201             updatetrack->SetDeadZoneProbability(ilayer,0.); 
1202           } 
1203         }
1204         ntracks[ilayer]++;
1205       }
1206
1207       clidx=-1;
1208       // loop over clusters in the road
1209       while ((cl=layer.GetNextCluster(clidx))!=0) { 
1210         if (ntracks[ilayer]>int(0.95*kMaxTr)) break; //space for skipped clusters  
1211         Bool_t changedet =kFALSE;  
1212         if (TMath::Abs(cl->GetQ())<1.e-13 && deadzoneSPD==kTRUE) continue;
1213         Int_t idetc=cl->GetDetectorIndex();
1214
1215         if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1216           // take into account misalignment (bring track to real detector plane)
1217           Double_t xTrOrig = currenttrack->GetX();
1218           if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1219           // a first cut on track-cluster distance
1220           if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz + 
1221                (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. ) 
1222             {  // cluster not associated to track
1223               AliDebug(2,"not associated");
1224               // MvL: added here as well
1225               // bring track back to ideal detector plane
1226               currenttrack->Propagate(xTrOrig);
1227               continue;
1228             }
1229           // bring track back to ideal detector plane
1230           if (!currenttrack->Propagate(xTrOrig)) continue;
1231         } else {                                      // have to move track to cluster's detector
1232           const AliITSdetector &detc=layer.GetDetector(idetc);
1233           // a first cut on track-cluster distance
1234           Double_t y;
1235           if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1236           if ( (z-cl->GetZ())*(z-cl->GetZ())*msz + 
1237                (y-cl->GetY())*(y-cl->GetY())*msy > 1. ) 
1238             continue; // cluster not associated to track
1239           //
1240           new (&backuptrack) AliITStrackMI(currenttrack2);
1241           changedet = kTRUE;
1242           currenttrack =&currenttrack2;
1243           if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1244             new (currenttrack) AliITStrackMI(backuptrack);
1245             changedet = kFALSE;
1246             continue;
1247           }
1248           currenttrack->SetDetectorIndex(idetc);
1249           // Get again the budget to the primary vertex 
1250           // for the current track being prolonged, if had to change detector 
1251           //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1252         }
1253
1254         // calculate track-clusters chi2
1255         chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer); 
1256         // chi2 cut
1257         AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1258         if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1259           if (TMath::Abs(cl->GetQ())<1.e-13) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster       
1260           if (ntracks[ilayer]>=kMaxTr) continue;
1261           AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1262           updatetrack->SetClIndex(ilayer,-1);
1263           if (changedet) new (&currenttrack2) AliITStrackMI(backuptrack);
1264
1265           if (TMath::Abs(cl->GetQ())>1.e-13) { // real cluster
1266             if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1267               AliDebug(2,"update failed");
1268               continue;
1269             } 
1270             updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2); 
1271             modstatus = 1; // found
1272           } else {             // virtual cluster in dead zone
1273             updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1274             updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1275             modstatus = 7; // holes in z in SPD
1276           }
1277
1278           if (changedet) {
1279             Float_t xlocnewdet,zlocnewdet;
1280             if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1281               updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1282             }
1283           } else {
1284             updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1285           }
1286           if (cl->IsUsed()) updatetrack->IncrementNUsed();
1287
1288           // apply correction for material of the current layer
1289           CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1290
1291           if (constrain) { // apply vertex constrain
1292             updatetrack->SetConstrain(constrain);
1293             Bool_t isPrim = kTRUE;
1294             if (ilayer<4) { // check that it's close to the vertex
1295               updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1296               if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1297                   AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() || 
1298                   TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1299                   AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1300             }
1301             if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1302               fUseImproveKalman ? 
1303                 updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) : 
1304                 updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1305             }
1306           } //apply vertex constrain              
1307           ntracks[ilayer]++;
1308         }  // create new hypothesis
1309         else {
1310           AliDebug(2,"chi2 too large");
1311         }
1312
1313       } // loop over possible prolongations 
1314      
1315       // allow one prolongation without clusters
1316       if (constrain && itrack<=1 && TMath::Abs(currenttrack1.GetNSkipped())<1.e-13 && deadzoneSPD==kFALSE && ntracks[ilayer]<kMaxTr) {
1317         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1318         // apply correction for material of the current layer
1319         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1320         vtrack->SetClIndex(ilayer,-1);
1321         modstatus = 3; // skipped 
1322         vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1323         if(AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1324           fUseImproveKalman ? 
1325             vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) : 
1326             vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1327         }
1328         vtrack->IncrementNSkipped();
1329         ntracks[ilayer]++;
1330       }
1331       
1332
1333     } // loop over tracks in layer ilayer+1
1334
1335     //loop over track candidates for the current layer
1336     //
1337     //
1338     Int_t accepted=0;
1339     
1340     Int_t golden=0;
1341     for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1342       normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer); 
1343       if (normalizedchi2[itrack] < 
1344           AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1345       if (ilayer>4) {
1346         accepted++;
1347       } else {
1348         if (constrain) { // constrain
1349           if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1) 
1350             accepted++;
1351         } else {        // !constrain
1352           if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1) 
1353             accepted++;
1354         }
1355       }
1356     }
1357     // sort tracks by increasing normalized chi2
1358     TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE); 
1359     ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1360     if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1361     //    if (ntracks[ilayer]>90) ntracks[ilayer]=90; 
1362     if (ntracks[ilayer]>int(kMaxTr*0.9)) ntracks[ilayer]=int(kMaxTr*0.9); 
1363   } // end loop over layers
1364
1365
1366   //
1367   // Now select tracks to be kept
1368   //
1369   Int_t max = constrain ? 20 : 5;
1370
1371   // tracks that reach layer 0 (SPD inner)
1372   for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1373     AliITStrackMI & track= tracks[0][nindexes[0][i]];
1374     if (track.GetNumberOfClusters()<2) continue;
1375     if (!constrain && track.GetNormChi2(0) >
1376         AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1377       continue;
1378     }
1379     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1380   }
1381
1382   // tracks that reach layer 1 (SPD outer)
1383   for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1384     AliITStrackMI & track= tracks[1][nindexes[1][i]];
1385     if (track.GetNumberOfClusters()<4) continue;
1386     if (!constrain && track.GetNormChi2(1) >
1387         AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1388     if (constrain) track.IncrementNSkipped();
1389     if (!constrain) {
1390       track.SetD(0,track.GetD(GetX(),GetY()));   
1391       track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1392       if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1393         track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1394       }
1395     }
1396     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1397   }
1398
1399   // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1400   if (!constrain){  
1401     for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1402       AliITStrackMI & track= tracks[2][nindexes[2][i]];
1403       if (track.GetNumberOfClusters()<3) continue;
1404       if (track.GetNormChi2(2) >
1405           AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1406       track.SetD(0,track.GetD(GetX(),GetY()));
1407       track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1408       if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1409         track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1410       }
1411       AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1412     }
1413   }
1414   
1415   if (!constrain) {
1416     //
1417     // register best track of each layer - important for V0 finder
1418     //
1419     for (Int_t ilayer=0;ilayer<5;ilayer++){
1420       if (ntracks[ilayer]==0) continue;
1421       AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1422       if (track.GetNumberOfClusters()<1) continue;
1423       CookLabel(&track,0);
1424       bestarray->AddAt(new AliITStrackMI(track),ilayer);
1425     }
1426   }
1427   //
1428   // update TPC V0 information
1429   //
1430   if (otrack->GetESDtrack()->GetV0Index(0)>0){    
1431     Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1432     for (Int_t i=0;i<3;i++){
1433       Int_t  index = otrack->GetESDtrack()->GetV0Index(i); 
1434       if (index==0) break;
1435       AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1436       if (vertex->GetStatus()<0) continue;     // rejected V0
1437       //
1438       if (otrack->GetSign()>0) {
1439         vertex->SetIndex(0,esdindex);
1440       }
1441       else{
1442         vertex->SetIndex(1,esdindex);
1443       }
1444       //find nearest layer with track info
1445       Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
1446       Int_t nearestold  = GetNearestLayer(xrp);               //I.B.
1447       Int_t nearest     = nearestold; 
1448       for (Int_t ilayer =nearest;ilayer<7;ilayer++){
1449         if (ntracks[nearest]==0){
1450           nearest = ilayer;
1451         }
1452       }
1453       //
1454       AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1455       if (nearestold<5&&nearest<5){
1456         Bool_t accept = track.GetNormChi2(nearest)<10; 
1457         if (accept){
1458           if (track.GetSign()>0) {
1459             vertex->SetParamP(track);
1460             vertex->Update(fprimvertex);
1461             //vertex->SetIndex(0,track.fESDtrack->GetID()); 
1462             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1463           }else{
1464             vertex->SetParamN(track);
1465             vertex->Update(fprimvertex);
1466             //vertex->SetIndex(1,track.fESDtrack->GetID());
1467             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1468           }
1469           vertex->SetStatus(vertex->GetStatus()+1);
1470         }else{
1471           //vertex->SetStatus(-2);  // reject V0  - not enough clusters
1472         }
1473       }
1474     }
1475   } 
1476   
1477 }
1478 //------------------------------------------------------------------------
1479 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1480 {
1481   //--------------------------------------------------------------------
1482   //
1483   //
1484   return fgLayers[layer];
1485 }
1486 //------------------------------------------------------------------------
1487 AliITStrackerMI::AliITSlayer::AliITSlayer():
1488 fR(0),
1489 fPhiOffset(0),
1490 fNladders(0),
1491 fZOffset(0),
1492 fNdetectors(0),
1493 fDetectors(0),
1494 fN(0),
1495 fDy5(0),
1496 fDy10(0),
1497 fDy20(0),
1498 fClustersCs(0),
1499 fClusterIndexCs(0),
1500 fYcs(0),
1501 fZcs(0),
1502 fNcs(0),
1503 fCurrentSlice(-1),
1504 fZmin(0),
1505 fZmax(0),
1506 fYmin(0),
1507 fYmax(0),
1508 fI(0),
1509 fImax(0),
1510 fSkip(0),
1511 fAccepted(0),
1512 fRoad(0),
1513 fMaxSigmaClY(0),
1514 fMaxSigmaClZ(0),
1515 fNMaxSigmaCl(3)
1516 {
1517   //--------------------------------------------------------------------
1518   //default AliITSlayer constructor
1519   //--------------------------------------------------------------------
1520   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1521     fClusterWeight[i]=0;
1522     fClusterTracks[0][i]=-1;
1523     fClusterTracks[1][i]=-1;
1524     fClusterTracks[2][i]=-1;    
1525     fClusterTracks[3][i]=-1;
1526     fY[i]=0;    
1527     fZ[i]=0;    
1528   }
1529   fYB[0]=0;
1530   fYB[1]=0;
1531
1532   for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer5; j++) {
1533     for (Int_t j1=0; j1<6; j1++) {
1534       fClusters5[j1][j]=0;
1535       fClusterIndex5[j1][j]=-1;
1536       fY5[j1][j]=0;
1537       fZ5[j1][j]=0;
1538       fN5[j1]=0;
1539       fBy5[j1][0]=0;
1540       fBy5[j1][1]=0;
1541     }
1542   }
1543
1544   for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer10; j++) {
1545     for (Int_t j1=0; j1<11; j1++) {
1546       fClusters10[j1][j]=0;
1547       fClusterIndex10[j1][j]=-1;
1548       fY10[j1][j]=0;
1549       fZ10[j1][j]=0;
1550       fN10[j1]=0;
1551       fBy10[j1][0]=0;
1552       fBy10[j1][1]=0;
1553     }
1554   }
1555
1556   for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer20; j++) {
1557     for (Int_t j1=0; j1<21; j1++) {
1558       fClusters20[j1][j]=0;
1559       fClusterIndex20[j1][j]=-1;
1560       fY20[j1][j]=0;
1561       fZ20[j1][j]=0;
1562       fN20[j1]=0;
1563       fBy20[j1][0]=0;
1564       fBy20[j1][1]=0;
1565     }
1566   }
1567   for(Int_t i=0;i<AliITSRecoParam::fgkMaxClusterPerLayer;i++){
1568     fClusters[i]=NULL;
1569     fClusterIndex[i]=0;
1570   }
1571 }
1572 //------------------------------------------------------------------------
1573 AliITStrackerMI::AliITSlayer::
1574 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1575 fR(r),
1576 fPhiOffset(p),
1577 fNladders(nl),
1578 fZOffset(z),
1579 fNdetectors(nd),
1580 fDetectors(0),
1581 fN(0),
1582 fDy5(0),
1583 fDy10(0),
1584 fDy20(0),
1585 fClustersCs(0),
1586 fClusterIndexCs(0),
1587 fYcs(0),
1588 fZcs(0),
1589 fNcs(0),
1590 fCurrentSlice(-1),
1591 fZmin(0),
1592 fZmax(0),
1593 fYmin(0),
1594 fYmax(0),
1595 fI(0),
1596 fImax(0),
1597 fSkip(0),
1598 fAccepted(0),
1599 fRoad(0),
1600 fMaxSigmaClY(0),
1601 fMaxSigmaClZ(0),
1602 fNMaxSigmaCl(3) {
1603   //--------------------------------------------------------------------
1604   //main AliITSlayer constructor
1605   //--------------------------------------------------------------------
1606   fDetectors=new AliITSdetector[fNladders*fNdetectors];
1607   fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1608
1609   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1610     fClusterWeight[i]=0;
1611     fClusterTracks[0][i]=-1;
1612     fClusterTracks[1][i]=-1;
1613     fClusterTracks[2][i]=-1;    
1614     fClusterTracks[3][i]=-1;    
1615     fY[i]=0;    
1616     fZ[i]=0;    
1617   }
1618
1619   fYB[0]=0;
1620   fYB[1]=0;
1621
1622   for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer5; j++) {
1623     for (Int_t j1=0; j1<6; j1++) {
1624       fClusters5[j1][j]=0;
1625       fClusterIndex5[j1][j]=-1;
1626       fY5[j1][j]=0;
1627       fZ5[j1][j]=0;
1628       fN5[j1]=0;
1629       fBy5[j1][0]=0;
1630       fBy5[j1][1]=0;
1631     }
1632   }
1633
1634   for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer10; j++) {
1635     for (Int_t j1=0; j1<11; j1++) {
1636       fClusters10[j1][j]=0;
1637       fClusterIndex10[j1][j]=-1;
1638       fY10[j1][j]=0;
1639       fZ10[j1][j]=0;
1640       fN10[j1]=0;
1641       fBy10[j1][0]=0;
1642       fBy10[j1][1]=0;
1643     }
1644   }
1645
1646   for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer20; j++) {
1647     for (Int_t j1=0; j1<21; j1++) {
1648       fClusters20[j1][j]=0;
1649       fClusterIndex20[j1][j]=-1;
1650       fY20[j1][j]=0;
1651       fZ20[j1][j]=0;
1652       fN20[j1]=0;
1653       fBy20[j1][0]=0;
1654       fBy20[j1][1]=0;
1655     }
1656   }
1657   for(Int_t i=0;i<AliITSRecoParam::fgkMaxClusterPerLayer;i++){
1658     fClusters[i]=NULL;
1659     fClusterIndex[i]=0;
1660   }
1661 }
1662 /*
1663 //------------------------------------------------------------------------
1664 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1665 fR(layer.fR),
1666 fPhiOffset(layer.fPhiOffset),
1667 fNladders(layer.fNladders),
1668 fZOffset(layer.fZOffset),
1669 fNdetectors(layer.fNdetectors),
1670 fDetectors(layer.fDetectors),
1671 fN(layer.fN),
1672 fDy5(layer.fDy5),
1673 fDy10(layer.fDy10),
1674 fDy20(layer.fDy20),
1675 fClustersCs(layer.fClustersCs),
1676 fClusterIndexCs(layer.fClusterIndexCs),
1677 fYcs(layer.fYcs),
1678 fZcs(layer.fZcs),
1679 fNcs(layer.fNcs),
1680 fCurrentSlice(layer.fCurrentSlice),
1681 fZmin(layer.fZmin),
1682 fZmax(layer.fZmax),
1683 fYmin(layer.fYmin),
1684 fYmax(layer.fYmax),
1685 fI(layer.fI),
1686 fImax(layer.fImax),
1687 fSkip(layer.fSkip),
1688 fAccepted(layer.fAccepted),
1689 fRoad(layer.fRoad),
1690 fMaxSigmaClY(layer.fMaxSigmaClY),
1691 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1692 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1693 {
1694   //Copy constructor
1695 }
1696 */
1697 //------------------------------------------------------------------------
1698 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1699   //--------------------------------------------------------------------
1700   // AliITSlayer destructor
1701   //--------------------------------------------------------------------
1702   delete [] fDetectors;
1703   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1704   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1705     fClusterWeight[i]=0;
1706     fClusterTracks[0][i]=-1;
1707     fClusterTracks[1][i]=-1;
1708     fClusterTracks[2][i]=-1;    
1709     fClusterTracks[3][i]=-1;    
1710   }
1711 }
1712 //------------------------------------------------------------------------
1713 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1714   //--------------------------------------------------------------------
1715   // This function removes loaded clusters
1716   //--------------------------------------------------------------------
1717   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1718   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1719     fClusterWeight[i]=0;
1720     fClusterTracks[0][i]=-1;
1721     fClusterTracks[1][i]=-1;
1722     fClusterTracks[2][i]=-1;    
1723     fClusterTracks[3][i]=-1;  
1724   }
1725   
1726   fN=0;
1727   fI=0;
1728 }
1729 //------------------------------------------------------------------------
1730 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1731   //--------------------------------------------------------------------
1732   // This function reset weights of the clusters
1733   //--------------------------------------------------------------------
1734   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1735     fClusterWeight[i]=0;
1736     fClusterTracks[0][i]=-1;
1737     fClusterTracks[1][i]=-1;
1738     fClusterTracks[2][i]=-1;    
1739     fClusterTracks[3][i]=-1;  
1740   }
1741   for (Int_t i=0; i<fN;i++) {
1742     AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1743     if (cl&&cl->IsUsed()) cl->Use();
1744   }
1745
1746 }
1747 //------------------------------------------------------------------------
1748 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1749   //--------------------------------------------------------------------
1750   // This function calculates the road defined by the cluster density
1751   //--------------------------------------------------------------------
1752   Int_t n=0;
1753   for (Int_t i=0; i<fN; i++) {
1754      if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1755   }
1756   if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1757 }
1758 //------------------------------------------------------------------------
1759 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1760   //--------------------------------------------------------------------
1761   //This function adds a cluster to this layer
1762   //--------------------------------------------------------------------
1763   if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1764     return 1;
1765   }
1766   fCurrentSlice=-1;
1767   fClusters[fN]=cl;
1768   fN++;
1769   AliITSdetector &det=GetDetector(cl->GetDetectorIndex());    
1770   //AD
1771   Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1772   Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2()); 
1773   if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1774   if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1775   if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1776   if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1777   //AD               
1778   /*
1779   if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1780   if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1781   if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1782   if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1783   */                 
1784   return 0;
1785 }
1786 //------------------------------------------------------------------------
1787 void  AliITStrackerMI::AliITSlayer::SortClusters()
1788 {
1789   //
1790   //sort clusters
1791   //
1792   AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1793   Float_t *z                = new Float_t[fN];
1794   Int_t   * index           = new Int_t[fN];
1795   //
1796   fMaxSigmaClY=0.; //AD
1797   fMaxSigmaClZ=0.; //AD
1798
1799   for (Int_t i=0;i<fN;i++){
1800     z[i] = fClusters[i]->GetZ();
1801     // save largest errors in y and z for this layer
1802     fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1803     fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1804   }
1805   TMath::Sort(fN,z,index,kFALSE);
1806   for (Int_t i=0;i<fN;i++){
1807     clusters[i] = fClusters[index[i]];
1808   }
1809   //
1810   for (Int_t i=0;i<fN;i++){
1811     fClusters[i] = clusters[i];
1812     fZ[i]        = fClusters[i]->GetZ();
1813     AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());    
1814     Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1815     if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1816     fY[i] = y;
1817   }
1818   delete[] index;
1819   delete[] z;
1820   delete[] clusters;
1821   //
1822
1823   fYB[0]=10000000;
1824   fYB[1]=-10000000;
1825   for (Int_t i=0;i<fN;i++){
1826     if (fY[i]<fYB[0]) fYB[0]=fY[i];
1827     if (fY[i]>fYB[1]) fYB[1]=fY[i];
1828     fClusterIndex[i] = i;
1829   }
1830   //
1831   // fill slices
1832   fDy5 = (fYB[1]-fYB[0])/5.;
1833   fDy10 = (fYB[1]-fYB[0])/10.;
1834   fDy20 = (fYB[1]-fYB[0])/20.;
1835   for (Int_t i=0;i<6;i++)  fN5[i] =0;  
1836   for (Int_t i=0;i<11;i++) fN10[i]=0;  
1837   for (Int_t i=0;i<21;i++) fN20[i]=0;
1838   //  
1839   for (Int_t i=0;i<6;i++) {fBy5[i][0] =  fYB[0]+(i-0.75)*fDy5; fBy5[i][1] =  fYB[0]+(i+0.75)*fDy5;}
1840   for (Int_t i=0;i<11;i++) {fBy10[i][0] =  fYB[0]+(i-0.75)*fDy10; fBy10[i][1] =  fYB[0]+(i+0.75)*fDy10;} 
1841   for (Int_t i=0;i<21;i++) {fBy20[i][0] =  fYB[0]+(i-0.75)*fDy20; fBy20[i][1] =  fYB[0]+(i+0.75)*fDy20;}
1842   //
1843   //
1844   for (Int_t i=0;i<fN;i++)
1845     for (Int_t irot=-1;irot<=1;irot++){
1846       Float_t curY = fY[i]+irot*TMath::TwoPi()*fR; 
1847       // slice 5
1848       for (Int_t slice=0; slice<6;slice++){
1849         if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1850           fClusters5[slice][fN5[slice]] = fClusters[i];
1851           fY5[slice][fN5[slice]] = curY;
1852           fZ5[slice][fN5[slice]] = fZ[i];
1853           fClusterIndex5[slice][fN5[slice]]=i;
1854           fN5[slice]++;
1855         }
1856       }
1857       // slice 10
1858       for (Int_t slice=0; slice<11;slice++){
1859         if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1860           fClusters10[slice][fN10[slice]] = fClusters[i];
1861           fY10[slice][fN10[slice]] = curY;
1862           fZ10[slice][fN10[slice]] = fZ[i];
1863           fClusterIndex10[slice][fN10[slice]]=i;
1864           fN10[slice]++;
1865         }
1866       }
1867       // slice 20
1868       for (Int_t slice=0; slice<21;slice++){
1869         if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1870           fClusters20[slice][fN20[slice]] = fClusters[i];
1871           fY20[slice][fN20[slice]] = curY;
1872           fZ20[slice][fN20[slice]] = fZ[i];
1873           fClusterIndex20[slice][fN20[slice]]=i;
1874           fN20[slice]++;
1875         }
1876       }      
1877     }
1878
1879   //
1880   // consistency check
1881   //
1882   for (Int_t i=0;i<fN-1;i++){
1883     if (fZ[i]>fZ[i+1]){
1884       printf("Bug\n");
1885     }
1886   }
1887   //
1888   for (Int_t slice=0;slice<21;slice++)
1889   for (Int_t i=0;i<fN20[slice]-1;i++){
1890     if (fZ20[slice][i]>fZ20[slice][i+1]){
1891       printf("Bug\n");
1892     }
1893   }
1894
1895
1896 }
1897 //------------------------------------------------------------------------
1898 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1899   //--------------------------------------------------------------------
1900   // This function returns the index of the nearest cluster 
1901   //--------------------------------------------------------------------
1902   Int_t ncl=0;
1903   const Float_t *zcl;  
1904   if (fCurrentSlice<0) {
1905     ncl = fN;
1906     zcl   = fZ;
1907   }
1908   else{
1909     ncl   = fNcs;
1910     zcl   = fZcs;;
1911   }
1912   
1913   if (ncl==0) return 0;
1914   Int_t b=0, e=ncl-1, m=(b+e)/2;
1915   for (; b<e; m=(b+e)/2) {
1916     //    if (z > fClusters[m]->GetZ()) b=m+1;
1917     if (z > zcl[m]) b=m+1;
1918     else e=m; 
1919   }
1920   return m;
1921 }
1922 //------------------------------------------------------------------------
1923 Bool_t AliITStrackerMI::ComputeRoad(AliITStrackMI* track,Int_t ilayer,Int_t idet,Double_t &zmin,Double_t &zmax,Double_t &ymin,Double_t &ymax) const {
1924   //--------------------------------------------------------------------
1925   // This function computes the rectangular road for this track
1926   //--------------------------------------------------------------------
1927
1928
1929   AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1930   // take into account the misalignment: propagate track to misaligned detector plane
1931   if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1932
1933   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1934                     TMath::Sqrt(track->GetSigmaZ2() + 
1935                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1936                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1937                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1938   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1939                     TMath::Sqrt(track->GetSigmaY2() + 
1940                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1941                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1942                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1943       
1944   // track at boundary between detectors, enlarge road
1945   Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1946   if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) || 
1947        (track->GetY()+dy > det.GetYmax()-boundaryWidth) || 
1948        (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1949        (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1950     Float_t tgl = TMath::Abs(track->GetTgl());
1951     if (tgl > 1.) tgl=1.;
1952     Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1953     dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1954     Float_t snp = TMath::Abs(track->GetSnp());
1955     if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1956     dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1957   } // boundary
1958   
1959   // add to the road a term (up to 2-3 mm) to deal with misalignments
1960   dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1961   dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1962
1963   Double_t r = fgLayers[ilayer].GetR();
1964   zmin = track->GetZ() - dz; 
1965   zmax = track->GetZ() + dz;
1966   ymin = track->GetY() + r*det.GetPhi() - dy;
1967   ymax = track->GetY() + r*det.GetPhi() + dy;
1968
1969   // bring track back to idead detector plane
1970   if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1971
1972   return kTRUE;
1973 }
1974 //------------------------------------------------------------------------
1975 void AliITStrackerMI::AliITSlayer::
1976 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1977   //--------------------------------------------------------------------
1978   // This function sets the "window"
1979   //--------------------------------------------------------------------
1980  
1981   Double_t circle=2*TMath::Pi()*fR;
1982   fYmin = ymin; 
1983   fYmax = ymax;
1984   fZmin = zmin;
1985   fZmax = zmax;
1986   // AD
1987   // enlarge road in y by maximum cluster error on this layer (3 sigma)
1988   fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
1989   fYmax += fNMaxSigmaCl*fMaxSigmaClY;
1990   fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
1991   fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
1992
1993   Float_t ymiddle = (fYmin+fYmax)*0.5;
1994   if (ymiddle<fYB[0]) {
1995     fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1996   } else if (ymiddle>fYB[1]) {
1997     fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1998   }
1999   
2000   //
2001   fCurrentSlice =-1;
2002   // defualt take all
2003   fClustersCs = fClusters;
2004   fClusterIndexCs = fClusterIndex;
2005   fYcs  = fY;
2006   fZcs  = fZ;
2007   fNcs  = fN;
2008   //
2009   //is in 20 slice?
2010   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
2011     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
2012     if (slice<0) slice=0;
2013     if (slice>20) slice=20;
2014     Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
2015     if (isOK) {
2016       fCurrentSlice=slice;
2017       fClustersCs = fClusters20[fCurrentSlice];
2018       fClusterIndexCs = fClusterIndex20[fCurrentSlice];
2019       fYcs  = fY20[fCurrentSlice];
2020       fZcs  = fZ20[fCurrentSlice];
2021       fNcs  = fN20[fCurrentSlice];
2022     }
2023   }  
2024   //
2025   //is in 10 slice?
2026   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
2027     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
2028     if (slice<0) slice=0;
2029     if (slice>10) slice=10;
2030     Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
2031     if (isOK) {
2032       fCurrentSlice=slice;
2033       fClustersCs = fClusters10[fCurrentSlice];
2034       fClusterIndexCs = fClusterIndex10[fCurrentSlice];
2035       fYcs  = fY10[fCurrentSlice];
2036       fZcs  = fZ10[fCurrentSlice];
2037       fNcs  = fN10[fCurrentSlice];
2038     }
2039   }  
2040   //
2041   //is in 5 slice?
2042   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
2043     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
2044     if (slice<0) slice=0;
2045     if (slice>5) slice=5;
2046     Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
2047     if (isOK) {
2048       fCurrentSlice=slice;
2049       fClustersCs = fClusters5[fCurrentSlice];
2050       fClusterIndexCs = fClusterIndex5[fCurrentSlice];
2051       fYcs  = fY5[fCurrentSlice];
2052       fZcs  = fZ5[fCurrentSlice];
2053       fNcs  = fN5[fCurrentSlice];
2054     }
2055   }  
2056   //  
2057   fI        = FindClusterIndex(fZmin); 
2058   fImax     = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
2059   fSkip     = 0;
2060   fAccepted = 0;
2061
2062   return;
2063 }
2064 //------------------------------------------------------------------------
2065 Int_t AliITStrackerMI::AliITSlayer::
2066 FindDetectorIndex(Double_t phi, Double_t z) const {
2067   //--------------------------------------------------------------------
2068   //This function finds the detector crossed by the track
2069   //--------------------------------------------------------------------
2070   Double_t dphi;
2071   if (fZOffset<0)            // old geometry
2072     dphi = -(phi-fPhiOffset);
2073   else                       // new geometry
2074     dphi = phi-fPhiOffset;
2075
2076
2077   if      (dphi <  0) dphi += 2*TMath::Pi();
2078   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
2079   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
2080   if (np>=fNladders) np-=fNladders;
2081   if (np<0)          np+=fNladders;
2082
2083
2084   Double_t dz=fZOffset-z;
2085   Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
2086   Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
2087   if (nz>=fNdetectors || nz<0) {
2088     //printf("ndet %d phi %f z %f  np %d nz %d\n",fNdetectors,phi,z,np,nz);
2089     return -1;
2090   }
2091
2092   // ad hoc correction for 3rd ladder of SDD inner layer,
2093   // which is reversed (rotated by pi around local y)
2094   // this correction is OK only from AliITSv11Hybrid onwards
2095   if (GetR()>12. && GetR()<20.) { // SDD inner
2096     if(np==2) { // 3rd ladder
2097       Double_t posMod252[3];
2098       AliITSgeomTGeo::GetTranslation(252,posMod252);
2099       // check the Z coordinate of Mod 252: if negative 
2100       // (old SDD geometry in AliITSv11Hybrid)
2101       // the swap of numeration whould be applied
2102       if(posMod252[2]<0.){
2103         nz = (fNdetectors-1) - nz;
2104       } 
2105     }
2106   }
2107   //printf("ndet %d phi %f z %f  np %d nz %d\n",fNdetectors,phi,z,np,nz);
2108
2109
2110   return np*fNdetectors + nz;
2111 }
2112 //------------------------------------------------------------------------
2113 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
2114 {
2115   //--------------------------------------------------------------------
2116   // This function returns clusters within the "window" 
2117   //--------------------------------------------------------------------
2118
2119   if (fCurrentSlice<0) {
2120     Double_t rpi2 = 2.*fR*TMath::Pi();
2121     for (Int_t i=fI; i<fImax; i++) {
2122       Double_t y = fY[i];
2123       Double_t z = fZ[i];
2124       if (fYmax<y) y -= rpi2;
2125       if (fYmin>y) y += rpi2;
2126       if (y<fYmin) continue;
2127       if (y>fYmax) continue;
2128       // AD
2129       // skip clusters that are in "extended" road but they 
2130       // 3sigma error does not touch the original road
2131       if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
2132       if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
2133       //
2134       if (TMath::Abs(fClusters[i]->GetQ())<1.e-13 && fSkip==2) continue;
2135       ci=i;
2136       if (!test) fI=i+1;
2137       return fClusters[i];
2138     }
2139   } else {
2140     for (Int_t i=fI; i<fImax; i++) {
2141       if (fYcs[i]<fYmin) continue;
2142       if (fYcs[i]>fYmax) continue;
2143       if (TMath::Abs(fClustersCs[i]->GetQ())<1.e-13 && fSkip==2) continue;
2144       ci=fClusterIndexCs[i];
2145       if (!test) fI=i+1;
2146       return fClustersCs[i];
2147     }
2148   }
2149   return 0;
2150 }
2151 //------------------------------------------------------------------------
2152 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
2153 const {
2154   //--------------------------------------------------------------------
2155   // This function returns the layer thickness at this point (units X0)
2156   //--------------------------------------------------------------------
2157   Double_t d=0.0085;
2158   x0=AliITSRecoParam::GetX0Air();
2159   if (43<fR&&fR<45) { //SSD2
2160      Double_t dd=0.0034;
2161      d=dd;
2162      if (TMath::Abs(y-0.00)>3.40) d+=dd;
2163      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2164      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2165      for (Int_t i=0; i<12; i++) {
2166        if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
2167           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2168           d+=0.0034; 
2169           break;
2170        }
2171        if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
2172           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2173           d+=0.0034; 
2174           break;
2175        }         
2176        if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2177        if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2178      }
2179   } else 
2180   if (37<fR&&fR<41) { //SSD1
2181      Double_t dd=0.0034;
2182      d=dd;
2183      if (TMath::Abs(y-0.00)>3.40) d+=dd;
2184      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2185      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2186      for (Int_t i=0; i<11; i++) {
2187        if (TMath::Abs(z-3.9*i)<0.15) {
2188           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2189           d+=dd; 
2190           break;
2191        }
2192        if (TMath::Abs(z+3.9*i)<0.15) {
2193           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2194           d+=dd; 
2195           break;
2196        }         
2197        if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2198        if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}         
2199      }
2200   } else
2201   if (13<fR&&fR<26) { //SDD
2202      Double_t dd=0.0033;
2203      d=dd;
2204      if (TMath::Abs(y-0.00)>3.30) d+=dd;
2205
2206      if (TMath::Abs(y-1.80)<0.55) {
2207         d+=0.016;
2208         for (Int_t j=0; j<20; j++) {
2209           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2210           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2211         } 
2212      }
2213      if (TMath::Abs(y+1.80)<0.55) {
2214         d+=0.016;
2215         for (Int_t j=0; j<20; j++) {
2216           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2217           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2218         } 
2219      }
2220
2221      for (Int_t i=0; i<4; i++) {
2222        if (TMath::Abs(z-7.3*i)<0.60) {
2223           d+=dd;
2224           if (TMath::Abs(y-0.00)>3.30) d+=dd; 
2225           break;
2226        }
2227        if (TMath::Abs(z+7.3*i)<0.60) {
2228           d+=dd; 
2229           if (TMath::Abs(y-0.00)>3.30) d+=dd; 
2230           break;
2231        }
2232      }
2233   } else
2234   if (6<fR&&fR<8) {   //SPD2
2235      Double_t dd=0.0063; x0=21.5;
2236      d=dd;
2237      if (TMath::Abs(y-3.08)>0.5) d+=dd;
2238      if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2239   } else
2240   if (3<fR&&fR<5) {   //SPD1
2241      Double_t dd=0.0063; x0=21.5;
2242      d=dd;
2243      if (TMath::Abs(y+0.21)>0.6) d+=dd;
2244      if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2245   }
2246
2247   return d;
2248 }
2249 //------------------------------------------------------------------------
2250 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2251 fR(det.fR),
2252 fRmisal(det.fRmisal),
2253 fPhi(det.fPhi),
2254 fSinPhi(det.fSinPhi),
2255 fCosPhi(det.fCosPhi),
2256 fYmin(det.fYmin),
2257 fYmax(det.fYmax),
2258 fZmin(det.fZmin),
2259 fZmax(det.fZmax),
2260 fIsBad(det.fIsBad),
2261 fNChips(det.fNChips),
2262 fChipIsBad(det.fChipIsBad)
2263 {
2264   //Copy constructor
2265 }
2266 //------------------------------------------------------------------------
2267 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2268                                                const AliITSDetTypeRec *detTypeRec)
2269 {
2270   //--------------------------------------------------------------------
2271   // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2272   //--------------------------------------------------------------------
2273
2274   // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2275   // while in the tracker they start from 0 for each layer
2276   for(Int_t il=0; il<ilayer; il++) 
2277     idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2278
2279   Int_t detType;
2280   if (ilayer==0 || ilayer==1) {        // ----------  SPD
2281     detType = 0;
2282   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
2283     detType = 1;
2284   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
2285     detType = 2;
2286   } else {
2287     printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2288     return;
2289   }
2290
2291   // Get calibration from AliITSDetTypeRec
2292   AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2293   calib->SetModuleIndex(idet);
2294   AliITSCalibration *calibSPDdead = 0;
2295   if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2296   if (calib->IsBad() ||
2297       (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2298     {
2299       SetBad();
2300       //      printf("lay %d bad %d\n",ilayer,idet);
2301     }
2302
2303   // Get segmentation from AliITSDetTypeRec
2304   AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2305
2306   // Read info about bad chips
2307   fNChips = segm->GetMaximumChipIndex()+1;
2308   //printf("ilayer %d  detType %d idet %d fNChips %d %d  GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2309   if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2310   fChipIsBad = new Bool_t[fNChips];
2311   for (Int_t iCh=0;iCh<fNChips;iCh++) {
2312     fChipIsBad[iCh] = calib->IsChipBad(iCh);
2313     if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2314     //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2315   }
2316
2317   return;
2318 }
2319 //------------------------------------------------------------------------
2320 Double_t AliITStrackerMI::GetEffectiveThickness()
2321 {
2322   //--------------------------------------------------------------------
2323   // Returns the thickness between the current layer and the vertex (units X0)
2324   //--------------------------------------------------------------------
2325
2326   if(fUseTGeo!=0) {
2327     if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2328     if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2329     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2330   }
2331
2332   // beam pipe
2333   Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2334   Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2335
2336   // layers
2337   Double_t x0=0;
2338   Double_t xn=fgLayers[fI].GetR();
2339   for (Int_t i=0; i<fI; i++) {
2340     Double_t xi=fgLayers[i].GetR();
2341     Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2342     d+=dLayer*xi*xi;
2343   }
2344
2345   // shields
2346   if (fI>1) {
2347     Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2348     d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2349   }
2350   if (fI>3) {
2351     Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2352     d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2353   }
2354   return d/(xn*xn);
2355 }
2356
2357 //------------------------------------------------------------------------
2358 Int_t AliITStrackerMI::GetEffectiveThicknessLbyL(Double_t* xMS, Double_t* x2x0MS)
2359 {
2360   //--------------------------------------------------------------------
2361   // Returns the array of layers between the current layer and the vertex
2362   //--------------------------------------------------------------------
2363   //
2364   if(fUseTGeo!=0) {
2365     if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2366     if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2367     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2368   }
2369
2370   int nl = 0;
2371   double x0 = 0;
2372   for (int il=fI;il--;) {
2373     //
2374     if (il==3) {
2375       x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2376       xMS[nl++]  = AliITSRecoParam::GetrInsideShield(1);
2377     }
2378     else if (il==1) {
2379       x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2380       xMS[nl++]  = AliITSRecoParam::GetrInsideShield(0);
2381     }
2382     //
2383     x2x0MS[nl] = (fUseTGeo==0 ? fgLayers[il].GetThickness(0,0,x0) : fxOverX0Layer[il]);
2384     xMS[nl++]  = fgLayers[il].GetR();
2385     //
2386   }
2387   //
2388   // beam pipe
2389   x2x0MS[nl]  = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2390   xMS[nl++]  = AliITSRecoParam::GetrPipe();
2391   //
2392   return nl;
2393 }
2394
2395
2396 //------------------------------------------------------------------------
2397 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2398   //-------------------------------------------------------------------
2399   // This function returns number of clusters within the "window" 
2400   //--------------------------------------------------------------------
2401   Int_t ncl=0;
2402   for (Int_t i=fI; i<fN; i++) {
2403     const AliITSRecPoint *c=fClusters[i];
2404     if (c->GetZ() > fZmax) break;
2405     if (c->IsUsed()) continue;
2406     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());    
2407     Double_t y=fR*det.GetPhi() + c->GetY();
2408
2409     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2410     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2411
2412     if (y<fYmin) continue;
2413     if (y>fYmax) continue;
2414     ncl++;
2415   }
2416   return ncl;
2417 }
2418 //------------------------------------------------------------------------
2419 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2420                                 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff) 
2421 {
2422   //--------------------------------------------------------------------
2423   // This function refits the track "track" at the position "x" using
2424   // the clusters from "clusters"
2425   // If "extra"==kTRUE, 
2426   //    the clusters from overlapped modules get attached to "track" 
2427   // If "planeff"==kTRUE,
2428   //    special approach for plane efficiency evaluation is applyed
2429   //--------------------------------------------------------------------
2430
2431   Int_t index[AliITSgeomTGeo::kNLayers];
2432   Int_t k;
2433   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2434   Int_t nc=clusters->GetNumberOfClusters();
2435   for (k=0; k<nc; k++) { 
2436     Int_t idx=clusters->GetClusterIndex(k);
2437     Int_t ilayer=(idx&0xf0000000)>>28;
2438     index[ilayer]=idx; 
2439   }
2440
2441   return RefitAt(xx,track,index,extra,planeeff); // call the method below
2442 }
2443 //------------------------------------------------------------------------
2444 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2445                                 const Int_t *clusters,Bool_t extra, Bool_t planeeff) 
2446 {
2447   //--------------------------------------------------------------------
2448   // This function refits the track "track" at the position "x" using
2449   // the clusters from array
2450   // If "extra"==kTRUE, 
2451   //    the clusters from overlapped modules get attached to "track" 
2452   // If "planeff"==kTRUE,
2453   //    special approach for plane efficiency evaluation is applyed
2454   //--------------------------------------------------------------------
2455   Int_t index[AliITSgeomTGeo::kNLayers];
2456   Int_t k;
2457   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2458   //
2459   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) { 
2460     index[k]=clusters[k]; 
2461   }
2462
2463   // special for cosmics and TPC prolonged tracks: 
2464   // propagate to the innermost of:
2465   // - innermost layer crossed by the track
2466   // - innermost layer where a cluster was associated to the track
2467   static AliITSRecoParam *repa = NULL;
2468   if(!repa){
2469     repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
2470     if(!repa){
2471       repa = AliITSRecoParam::GetHighFluxParam();
2472       AliWarning("Using default AliITSRecoParam class");
2473     }
2474   }
2475   Int_t evsp=repa->GetEventSpecie();
2476   ULong_t trStatus=0;
2477   if(track->GetESDtrack()) trStatus=track->GetStatus();
2478   Int_t innermostlayer=0;
2479   if((evsp&AliRecoParam::kCosmic) || (trStatus&AliESDtrack::kTPCin))  {
2480     innermostlayer=5;
2481     Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2482     for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2483       if( (drphi < (fgLayers[innermostlayer].GetR()+1.)) ||
2484           index[innermostlayer] >= 0 ) break;
2485     }
2486
2487     AliDebug(2,Form(" drphi  %f  innermost %d",drphi,innermostlayer));
2488   }
2489
2490   Int_t modstatus=1; // found
2491   Float_t xloc,zloc;
2492   Int_t from, to, step;
2493   if (xx > track->GetX()) {
2494       from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2495       step=+1;
2496   } else {
2497       from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2498       step=-1;
2499   }
2500   TString dir = (step>0 ? "outward" : "inward");
2501
2502   for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2503      AliITSlayer &layer=fgLayers[ilayer];
2504      Double_t r=layer.GetR();
2505
2506      if (step<0 && xx>r) break;
2507
2508      // material between SSD and SDD, SDD and SPD
2509      Double_t hI=ilayer-0.5*step; 
2510      if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2511        if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2512      if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2513        if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2514
2515
2516      Double_t oldGlobXYZ[3];
2517      if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2518
2519      // continue if we are already beyond this layer
2520      Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2521      if(step>0 && oldGlobR > r) continue; // going outward
2522      if(step<0 && oldGlobR < r) continue; // going inward
2523
2524      Double_t phi,z;
2525      if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2526
2527      Int_t idet=layer.FindDetectorIndex(phi,z);
2528
2529      // check if we allow a prolongation without point for large-eta tracks
2530      Int_t skip = CheckSkipLayer(track,ilayer,idet);
2531      if (skip==2) {
2532        modstatus = 4; // out in z
2533        if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2534          track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2535        }
2536        // cross layer
2537        // apply correction for material of the current layer
2538        // add time if going outward
2539        CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2540        continue;
2541      }
2542
2543      if (idet<0) return kFALSE;
2544
2545
2546      const AliITSdetector &det=layer.GetDetector(idet);
2547      // only for ITS-SA tracks refit
2548      if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2549      // 
2550      if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2551
2552      track->SetDetectorIndex(idet);
2553      if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2554
2555      Double_t dz,zmin,zmax,dy,ymin,ymax;
2556
2557      const AliITSRecPoint *clAcc=0;
2558      Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2559
2560      Int_t idx=index[ilayer];
2561      if (idx>=0) { // cluster in this layer
2562        modstatus = 6; // no refit
2563        const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx); 
2564        if (cl) {
2565          if (idet != cl->GetDetectorIndex()) {
2566            idet=cl->GetDetectorIndex();
2567            const AliITSdetector &detc=layer.GetDetector(idet);
2568            if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2569            track->SetDetectorIndex(idet);
2570            if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2571          }
2572          Int_t cllayer = (idx & 0xf0000000) >> 28;;
2573          Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2574          if (chi2<maxchi2) { 
2575            clAcc=cl; 
2576            maxchi2=chi2; 
2577            modstatus = 1; // found
2578          } else {
2579             return kFALSE; //
2580          }
2581        }
2582      } else { // no cluster in this layer
2583        if (skip==1) {
2584          modstatus = 3; // skipped
2585          // Plane Eff determination:
2586          if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2587            if (IsOKForPlaneEff(track,clusters,ilayer))  // only adequate track for plane eff. evaluation
2588               UseTrackForPlaneEff(track,ilayer);
2589          }
2590        } else {
2591          modstatus = 5; // no cls in road
2592          // check dead
2593          if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2594          dz = 0.5*(zmax-zmin);
2595          dy = 0.5*(ymax-ymin);
2596          Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2597          if (dead==1) modstatus = 7; // holes in z in SPD
2598          if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2599        }
2600      }
2601      
2602      if (clAcc) {
2603        if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2604        track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2605      }
2606      track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2607
2608
2609      if (extra && clAcc) { // search for extra clusters in overlapped modules
2610        AliITStrackV2 tmp(*track);
2611        if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2612        layer.SelectClusters(zmin,zmax,ymin,ymax);
2613        
2614        const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2615        Int_t idetExtra=-1;  
2616        maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2617        Double_t tolerance=0.1;
2618        while ((clExtra=layer.GetNextCluster(ci))!=0) {
2619          // only clusters in another module! (overlaps)
2620          idetExtra = clExtra->GetDetectorIndex();
2621          if (idet == idetExtra) continue;
2622          
2623          const AliITSdetector &detx=layer.GetDetector(idetExtra);
2624          
2625          if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2626          if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2627          if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2628          if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2629          
2630          Double_t chi2=tmp.GetPredictedChi2(clExtra);
2631          if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2632        }
2633        if (cci>=0) {
2634          track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2635          track->SetExtraModule(ilayer,idetExtra);
2636        }
2637      } // end search for extra clusters in overlapped modules
2638      
2639      // Correct for material of the current layer
2640      // cross material
2641      // add time if going outward
2642      if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2643      track->SetCheckInvariant(kTRUE);
2644   } // end loop on layers
2645
2646   if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2647
2648   return kTRUE;
2649 }
2650 //------------------------------------------------------------------------
2651 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2652 {
2653   //
2654   // calculate normalized chi2
2655   //  return NormalizedChi2(track,0);
2656   Float_t chi2 = 0;
2657   Float_t sum=0;
2658   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2659   //  track->fdEdxMismatch=0;
2660   Float_t dedxmismatch =0;
2661   Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack); 
2662   if (mode<100){
2663     for (Int_t i = 0;i<6;i++){
2664       if (track->GetClIndex(i)>=0){
2665         Float_t cerry, cerrz;
2666         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2667         else 
2668           { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2669         cerry*=cerry;
2670         cerrz*=cerrz;   
2671         Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2672         if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2673           Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2674           if (ratio<0.5) {
2675             cchi2+=(0.5-ratio)*10.;
2676             //track->fdEdxMismatch+=(0.5-ratio)*10.;
2677             dedxmismatch+=(0.5-ratio)*10.;          
2678           }
2679         }
2680         if (i<2 ||i>3){
2681           AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));  
2682           Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2683           if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.); 
2684           if (i<2) chi2+=2*cl->GetDeltaProbability();
2685         }
2686         chi2+=cchi2;
2687         sum++;
2688       }
2689     }
2690     if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2691       track->SetdEdxMismatch(dedxmismatch);
2692     }
2693   }
2694   else{
2695     for (Int_t i = 0;i<4;i++){
2696       if (track->GetClIndex(i)>=0){
2697         Float_t cerry, cerrz;
2698         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2699         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2700         cerry*=cerry;
2701         cerrz*=cerrz;
2702         chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2703         chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);      
2704         sum++;
2705       }
2706     }
2707     for (Int_t i = 4;i<6;i++){
2708       if (track->GetClIndex(i)>=0){     
2709         Float_t cerry, cerrz;
2710         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2711         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2712         cerry*=cerry;
2713         cerrz*=cerrz;   
2714         Float_t cerryb, cerrzb;
2715         if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2716         else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2717         cerryb*=cerryb;
2718         cerrzb*=cerrzb;
2719         chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2720         chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);      
2721         sum++;
2722       }
2723     }
2724   }
2725   if (track->GetESDtrack()->GetTPCsignal()>85){
2726     Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2727     if (ratio<0.5) {
2728       chi2+=(0.5-ratio)*5.;      
2729     }
2730     if (ratio>2){
2731       chi2+=(ratio-2.0)*3; 
2732     }
2733   }
2734   //
2735   Double_t match = TMath::Sqrt(track->GetChi22());
2736   if (track->GetConstrain())  match/=track->GetNumberOfClusters();
2737   if (!track->GetConstrain()) { 
2738     if (track->GetNumberOfClusters()>2) {
2739       match/=track->GetNumberOfClusters()-2.;
2740     } else {
2741       match=0;
2742     }
2743   }
2744   if (match<0) match=0;
2745
2746   // penalty factor for missing points (NDeadZone>0), but no penalty
2747   // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2748   // or there is a dead from OCDB)
2749   Float_t deadzonefactor = 0.; 
2750   if (track->GetNDeadZone()>0.) {    
2751     Int_t sumDeadZoneProbability=0; 
2752     for(Int_t ilay=0;ilay<6;ilay++) {
2753       if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2754     }
2755     Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2756     if(nDeadZoneWithProbNot1>0) {
2757       Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
2758       AliDebug(2,Form("nDeadZone %f sumDZProbability %d nDZWithProbNot1 %d deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
2759       deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2760       Float_t one = 1.;
2761       deadZoneProbability = TMath::Min(deadZoneProbability,one);
2762       deadzonefactor = 3.*(1.1-deadZoneProbability);  
2763     }
2764   }  
2765
2766   Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2767     (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2768                                 1./(1.+track->GetNSkipped()));     
2769   AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped %f\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
2770   AliDebug(2,Form("NormChi2 %f  cls %d\n",normchi2,track->GetNumberOfClusters()));
2771   return normchi2;
2772 }
2773 //------------------------------------------------------------------------
2774 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2775 {
2776   //
2777   // return matching chi2 between two tracks
2778   Double_t largeChi2=1000.;
2779
2780   AliITStrackMI track3(*track2);
2781   if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2782   TMatrixD vec(5,1);
2783   vec(0,0)=track1->GetY()   - track3.GetY();
2784   vec(1,0)=track1->GetZ()   - track3.GetZ();
2785   vec(2,0)=track1->GetSnp() - track3.GetSnp();
2786   vec(3,0)=track1->GetTgl() - track3.GetTgl();
2787   vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2788   //
2789   TMatrixD cov(5,5);
2790   cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2791   cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2792   cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2793   cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2794   cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2795   
2796   cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2797   cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2798   cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2799   cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2800   //
2801   cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2802   cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2803   cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2804   //
2805   cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2806   cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2807   //
2808   cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2809   
2810   cov.Invert();
2811   TMatrixD vec2(cov,TMatrixD::kMult,vec);
2812   TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2813   return chi2(0,0);
2814 }
2815 //------------------------------------------------------------------------
2816 Double_t  AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2817 {
2818   //
2819   //  return probability that given point (characterized by z position and error) 
2820   //  is in SPD dead zone
2821   //     This method assumes that fSPDdetzcentre is ordered from -z to +z
2822   //
2823   Double_t probability = 0.;
2824   Double_t nearestz = 0.,distz=0.;
2825   Int_t    nearestzone = -1;
2826   Double_t mindistz = 1000.;
2827
2828   // find closest dead zone
2829   for (Int_t i=0; i<3; i++) {
2830     distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2831     if (distz<mindistz) {
2832       nearestzone=i;
2833       nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2834       mindistz=distz;
2835     }
2836   }
2837
2838   // too far from dead zone
2839   if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2840
2841
2842   Double_t zmin, zmax;   
2843   if (nearestzone==0) { // dead zone at z = -7
2844     zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2845     zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2846   } else if (nearestzone==1) { // dead zone at z = 0
2847     zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2848     zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2849   } else if (nearestzone==2) { // dead zone at z = +7
2850     zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2851     zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2852   } else {
2853     zmin = 0.;
2854     zmax = 0.;
2855   }
2856   // probability that the true z is in the range [zmin,zmax] (i.e. inside 
2857   // dead zone)
2858   probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) - 
2859                       AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2860   AliDebug(2,Form("zpos %f +- %f nearestzone %d  zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
2861   return probability;
2862 }
2863 //------------------------------------------------------------------------
2864 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2865 {
2866   //
2867   // calculate normalized chi2
2868   Float_t chi2[6];
2869   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2870   Float_t ncl = 0;
2871   for (Int_t i = 0;i<6;i++){
2872     if (TMath::Abs(track->GetDy(i))>0){      
2873       chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2874       chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2875       ncl++;
2876     }
2877     else{chi2[i]=10000;}
2878   }
2879   Int_t index[6];
2880   TMath::Sort(6,chi2,index,kFALSE);
2881   Float_t max = float(ncl)*fac-1.;
2882   Float_t sumchi=0, sumweight=0; 
2883   for (Int_t i=0;i<max+1;i++){
2884     Float_t weight = (i<max)?1.:(max+1.-i);
2885     sumchi+=weight*chi2[index[i]];
2886     sumweight+=weight;
2887   }
2888   Double_t normchi2 = sumchi/sumweight;
2889   return normchi2;
2890 }
2891 //------------------------------------------------------------------------
2892 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2893 {
2894   //
2895   // calculate normalized chi2
2896   //  if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2897   Int_t npoints = 0;
2898   Double_t res =0;
2899   for (Int_t i=0;i<6;i++){
2900     if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2901     Double_t sy1 = forwardtrack->GetSigmaY(i);
2902     Double_t sz1 = forwardtrack->GetSigmaZ(i);
2903     Double_t sy2 = backtrack->GetSigmaY(i);
2904     Double_t sz2 = backtrack->GetSigmaZ(i);
2905     if (i<2){ sy2=1000.;sz2=1000;}
2906     //    
2907     Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2908     Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2909     // 
2910     Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2911     Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2912     //
2913     res+= nz0*nz0+ny0*ny0;
2914     npoints++;
2915   }
2916   if (npoints>1) return 
2917                    TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2918                    //2*forwardtrack->fNUsed+
2919                    res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2920                                   1./(1.+forwardtrack->GetNSkipped()));
2921   return 1000;
2922 }
2923 //------------------------------------------------------------------------
2924 Float_t  *AliITStrackerMI::GetWeight(Int_t index) {
2925   //--------------------------------------------------------------------
2926   //       Return pointer to a given cluster
2927   //--------------------------------------------------------------------
2928   Int_t l=(index & 0xf0000000) >> 28;
2929   Int_t c=(index & 0x0fffffff) >> 00;
2930   return fgLayers[l].GetWeight(c);
2931 }
2932 //------------------------------------------------------------------------
2933 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2934 {
2935   //---------------------------------------------
2936   // register track to the list
2937   //
2938   if (track->GetESDtrack()->GetKinkIndex(0)!=0) return;  //don't register kink tracks
2939   //
2940   //
2941   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2942     Int_t index = track->GetClusterIndex(icluster);
2943     Int_t l=(index & 0xf0000000) >> 28;
2944     Int_t c=(index & 0x0fffffff) >> 00;
2945     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2946     for (Int_t itrack=0;itrack<4;itrack++){
2947       if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2948         fgLayers[l].SetClusterTracks(itrack,c,id);
2949         break;
2950       }
2951     }
2952   }
2953 }
2954 //------------------------------------------------------------------------
2955 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2956 {
2957   //---------------------------------------------
2958   // unregister track from the list
2959   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2960     Int_t index = track->GetClusterIndex(icluster);
2961     Int_t l=(index & 0xf0000000) >> 28;
2962     Int_t c=(index & 0x0fffffff) >> 00;
2963     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2964     for (Int_t itrack=0;itrack<4;itrack++){
2965       if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2966         fgLayers[l].SetClusterTracks(itrack,c,-1);
2967       }
2968     }
2969   }
2970 }
2971 //------------------------------------------------------------------------
2972 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2973 {
2974   //-------------------------------------------------------------
2975   //get number of shared clusters
2976   //-------------------------------------------------------------
2977   Float_t shared=0;
2978   for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2979   // mean  number of clusters
2980   Float_t *ny = GetNy(id), *nz = GetNz(id); 
2981
2982  
2983   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2984     Int_t index = track->GetClusterIndex(icluster);
2985     Int_t l=(index & 0xf0000000) >> 28;
2986     Int_t c=(index & 0x0fffffff) >> 00;
2987     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2988     // if (ny[l]<1.e-13){
2989     //   printf("problem\n");
2990     // }
2991     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2992     Float_t weight=1;
2993     //
2994     Float_t deltan = 0;
2995     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2996     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2997       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2998     if (l<2 || l>3){      
2999       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3000     }
3001     else{
3002       deltan = (cl->GetNz()-nz[l]);
3003     }
3004     if (deltan>2.0) continue;  // extended - highly probable shared cluster
3005     weight = 2./TMath::Max(3.+deltan,2.);
3006     //
3007     for (Int_t itrack=0;itrack<4;itrack++){
3008       if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
3009         list[l]=index;
3010         clist[l] = (AliITSRecPoint*)GetCluster(index);
3011         track->SetSharedWeight(l,weight);
3012         shared+=weight; 
3013         break;
3014       }
3015     }
3016   }
3017   track->SetNUsed(shared);
3018   return shared;
3019 }
3020 //------------------------------------------------------------------------
3021 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
3022 {
3023   //
3024   // find first shared track 
3025   //
3026   // mean  number of clusters
3027   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID); 
3028   //
3029   for (Int_t i=0;i<6;i++) overlist[i]=-1;
3030   Int_t sharedtrack=100000;
3031   Int_t tracks[24],trackindex=0;
3032   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
3033   //
3034   for (Int_t icluster=0;icluster<6;icluster++){
3035     if (clusterlist[icluster]<0) continue;
3036     Int_t index = clusterlist[icluster];
3037     Int_t l=(index & 0xf0000000) >> 28;
3038     Int_t c=(index & 0x0fffffff) >> 00;
3039     // if (ny[l]<1.e-13){
3040     //   printf("problem\n");
3041     // }
3042     if (c>fgLayers[l].GetNumberOfClusters()) continue;
3043     //if (l>3) continue;
3044     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3045     //
3046     Float_t deltan = 0;
3047     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
3048     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
3049       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3050     if (l<2 || l>3){      
3051       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3052     }
3053     else{
3054       deltan = (cl->GetNz()-nz[l]);
3055     }
3056     if (deltan>2.0) continue;  // extended - highly probable shared cluster
3057     //
3058     for (Int_t itrack=3;itrack>=0;itrack--){
3059       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3060       if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
3061        tracks[trackindex]  = fgLayers[l].GetClusterTracks(itrack,c);
3062        trackindex++;
3063       }
3064     }
3065   }
3066   if (trackindex==0) return -1;
3067   if (trackindex==1){    
3068     sharedtrack = tracks[0];
3069   }else{
3070     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
3071     else{
3072       //
3073       Int_t tracks2[24], cluster[24];
3074       for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
3075       Int_t index =0;
3076       //
3077       for (Int_t i=0;i<trackindex;i++){
3078         if (tracks[i]<0) continue;
3079         tracks2[index] = tracks[i];
3080         cluster[index]++;       
3081         for (Int_t j=i+1;j<trackindex;j++){
3082           if (tracks[j]<0) continue;
3083           if (tracks[j]==tracks[i]){
3084             cluster[index]++;
3085             tracks[j]=-1;
3086           }
3087         }
3088         index++;
3089       }
3090       Int_t max=0;
3091       for (Int_t i=0;i<index;i++){
3092         if (cluster[index]>max) {
3093           sharedtrack=tracks2[index];
3094           max=cluster[index];
3095         }
3096       }
3097     }
3098   }
3099   
3100   if (sharedtrack>=100000) return -1;
3101   //
3102   // make list of overlaps
3103   shared =0;
3104   for (Int_t icluster=0;icluster<6;icluster++){
3105     if (clusterlist[icluster]<0) continue;
3106     Int_t index = clusterlist[icluster];
3107     Int_t l=(index & 0xf0000000) >> 28;
3108     Int_t c=(index & 0x0fffffff) >> 00;
3109     if (c>fgLayers[l].GetNumberOfClusters()) continue;
3110     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3111     if (l==0 || l==1){
3112       if (cl->GetNy()>2) continue;
3113       if (cl->GetNz()>2) continue;
3114     }
3115     if (l==4 || l==5){
3116       if (cl->GetNy()>3) continue;
3117       if (cl->GetNz()>3) continue;
3118     }
3119     //
3120     for (Int_t itrack=3;itrack>=0;itrack--){
3121       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3122       if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
3123         overlist[l]=index;
3124         shared++;      
3125       }
3126     }
3127   }
3128   return sharedtrack;
3129 }
3130 //------------------------------------------------------------------------
3131 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1,AliITStrackMI* original){
3132   //
3133   // try to find track hypothesys without conflicts
3134   // with minimal chi2;
3135   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
3136   Int_t entries1 = arr1->GetEntriesFast();
3137   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
3138   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
3139   Int_t entries2 = arr2->GetEntriesFast();
3140   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
3141   //
3142   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
3143   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
3144   if (track10->Pt()>0.5+track20->Pt()) return track10;
3145   //
3146   for (Int_t itrack=0;itrack<entries1;itrack++){
3147     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3148     UnRegisterClusterTracks(track,trackID1);
3149   }
3150   //
3151   for (Int_t itrack=0;itrack<entries2;itrack++){
3152     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3153     UnRegisterClusterTracks(track,trackID2);
3154   }
3155   Int_t index1=0;
3156   Int_t index2=0;
3157   Float_t maxconflicts=6;
3158   Double_t maxchi2 =1000.;
3159   //
3160   // get weight of hypothesys - using best hypothesy
3161   Double_t w1,w2;
3162  
3163   Int_t list1[6],list2[6];
3164   AliITSRecPoint *clist1[6], *clist2[6] ;
3165   RegisterClusterTracks(track10,trackID1);
3166   RegisterClusterTracks(track20,trackID2);
3167   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
3168   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
3169   UnRegisterClusterTracks(track10,trackID1);
3170   UnRegisterClusterTracks(track20,trackID2);
3171   //
3172   // normalized chi2
3173   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
3174   Float_t nerry[6],nerrz[6];
3175   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
3176   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
3177   for (Int_t i=0;i<6;i++){
3178      if ( (erry1[i]>0) && (erry2[i]>0)) {
3179        nerry[i] = TMath::Min(erry1[i],erry2[i]);
3180        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
3181      }else{
3182        nerry[i] = TMath::Max(erry1[i],erry2[i]);
3183        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
3184      }
3185      if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
3186        chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
3187        chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
3188        ncl1++;
3189      }
3190      if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
3191        chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
3192        chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
3193        ncl2++;
3194      }
3195   }
3196   chi21/=ncl1;
3197   chi22/=ncl2;
3198   //
3199   // 
3200   Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
3201   Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
3202   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
3203   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
3204   //
3205   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
3206         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
3207         +1.*track10->Pt()/(track10->Pt()+track20->Pt())
3208         );
3209   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
3210         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
3211         +1.*track20->Pt()/(track10->Pt()+track20->Pt())
3212         );
3213
3214   Double_t sumw = w1+w2;
3215   w1/=sumw;
3216   w2/=sumw;
3217   if (w1<w2*0.5) {
3218     w1 = (d2+0.5)/(d1+d2+1.);
3219     w2 = (d1+0.5)/(d1+d2+1.);
3220   }
3221   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
3222   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
3223   //
3224   // get pair of "best" hypothesys
3225   //  
3226   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1); 
3227   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2); 
3228
3229   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3230     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3231     //if (track1->fFakeRatio>0) continue;
3232     RegisterClusterTracks(track1,trackID1);
3233     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3234       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3235
3236       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3237       //if (track2->fFakeRatio>0) continue;
3238       Float_t nskipped=0;            
3239       RegisterClusterTracks(track2,trackID2);
3240       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3241       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3242       UnRegisterClusterTracks(track2,trackID2);
3243       //
3244       if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3245       if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3246       if (nskipped>0.5) continue;
3247       //
3248       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3249       if (conflict1+1<cconflict1) continue;
3250       if (conflict2+1<cconflict2) continue;
3251       Float_t conflict=0;
3252       Float_t sumchi2=0;
3253       Float_t sum=0;
3254       for (Int_t i=0;i<6;i++){
3255         //
3256         Float_t c1 =0.; // conflict coeficients
3257         Float_t c2 =0.; 
3258         if (clist1[i]&&clist2[i]){
3259           Float_t deltan = 0;
3260           if (i<2 || i>3){      
3261             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3262           }
3263           else{
3264             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3265           }
3266           c1  = 2./TMath::Max(3.+deltan,2.);      
3267           c2  = 2./TMath::Max(3.+deltan,2.);      
3268         }
3269         else{
3270           if (clist1[i]){
3271             Float_t deltan = 0;
3272             if (i<2 || i>3){      
3273               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3274             }
3275             else{
3276               deltan = (clist1[i]->GetNz()-nz1[i]);
3277             }
3278             c1  = 2./TMath::Max(3.+deltan,2.);    
3279             c2  = 0;
3280           }
3281
3282           if (clist2[i]){
3283             Float_t deltan = 0;
3284             if (i<2 || i>3){      
3285               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3286             }
3287             else{
3288               deltan = (clist2[i]->GetNz()-nz2[i]);
3289             }
3290             c2  = 2./TMath::Max(3.+deltan,2.);    
3291             c1  = 0;
3292           }       
3293         }
3294         //
3295         chi21=0;chi22=0;
3296         if (TMath::Abs(track1->GetDy(i))>0.) {
3297           chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3298             (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3299           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3300           //  (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3301         }else{
3302           if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3303         }
3304         //
3305         if (TMath::Abs(track2->GetDy(i))>0.) {
3306           chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3307             (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3308           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3309           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3310         }
3311         else{
3312           if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3313         }
3314         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3315         if (chi21>0) sum+=w1;
3316         if (chi22>0) sum+=w2;
3317         conflict+=(c1+c2);
3318       }
3319       Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3320       if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3321       Double_t normchi2 = 2*conflict+sumchi2/sum;
3322       if ( normchi2 <maxchi2 ){   
3323         index1 = itrack1;
3324         index2 = itrack2;
3325         maxconflicts = conflict;
3326         maxchi2 = normchi2;
3327       }      
3328     }
3329     UnRegisterClusterTracks(track1,trackID1);
3330   }
3331   //
3332   //  if (maxconflicts<4 && maxchi2<th0){   
3333   if (maxchi2<th0*2.){   
3334     Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3335     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3336     track1->SetChi2MIP(5,maxconflicts);
3337     track1->SetChi2MIP(6,maxchi2);
3338     track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3339     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
3340     track1->SetChi2MIP(8,index1);
3341     fBestTrackIndex[trackID1] =index1;
3342     UpdateESDtrack(track1, AliESDtrack::kITSin);
3343     original->SetWinner(track1);
3344   }  
3345   else if (track10->GetChi2MIP(0)<th1){
3346     track10->SetChi2MIP(5,maxconflicts);
3347     track10->SetChi2MIP(6,maxchi2);    
3348     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
3349     UpdateESDtrack(track10,AliESDtrack::kITSin);
3350     original->SetWinner(track10);
3351   }   
3352   
3353   for (Int_t itrack=0;itrack<entries1;itrack++){
3354     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3355     UnRegisterClusterTracks(track,trackID1);
3356   }
3357   //
3358   for (Int_t itrack=0;itrack<entries2;itrack++){
3359     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3360     UnRegisterClusterTracks(track,trackID2);
3361   }
3362
3363   if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3364       &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){ 
3365     //  if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3366   //    &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){ 
3367     RegisterClusterTracks(track10,trackID1);
3368   }
3369   if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3370       &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){ 
3371     //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3372     //  &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){ 
3373     RegisterClusterTracks(track20,trackID2);  
3374   }
3375   return track10; 
3376  
3377 }
3378 //------------------------------------------------------------------------
3379 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3380   //--------------------------------------------------------------------
3381   // This function marks clusters assigned to the track
3382   //--------------------------------------------------------------------
3383   AliTracker::UseClusters(t,from);
3384
3385   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3386   //if (c->GetQ()>2) c->Use();
3387   if (c->GetSigmaZ2()>0.1) c->Use();
3388   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3389   //if (c->GetQ()>2) c->Use();
3390   if (c->GetSigmaZ2()>0.1) c->Use();
3391
3392 }
3393 //------------------------------------------------------------------------
3394 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3395 {
3396   //------------------------------------------------------------------
3397   // add track to the list of hypothesys
3398   //------------------------------------------------------------------
3399
3400   if (esdindex>=fTrackHypothesys.GetEntriesFast()) 
3401     fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3402   //
3403   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3404   if (!array) {
3405     array = new TObjArray(10);
3406     fTrackHypothesys.AddAt(array,esdindex);
3407   }
3408   array->AddLast(track);
3409 }
3410 //------------------------------------------------------------------------
3411 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3412 {
3413   //-------------------------------------------------------------------
3414   // compress array of track hypothesys
3415   // keep only maxsize best hypothesys
3416   //-------------------------------------------------------------------
3417   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3418   if (! (fTrackHypothesys.At(esdindex)) ) return;
3419   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3420   Int_t entries = array->GetEntriesFast();
3421   //
3422   //- find preliminary besttrack as a reference
3423   Float_t minchi2=10000;
3424   Int_t maxn=0;
3425   AliITStrackMI * besttrack=0;
3426   //
3427   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3428     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3429     if (!track) continue;
3430     Float_t chi2 = NormalizedChi2(track,0);
3431     //
3432     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3433     track->SetLabel(tpcLabel);
3434     CookdEdx(track);
3435     track->SetFakeRatio(1.);
3436     CookLabel(track,0.); //For comparison only
3437     //
3438     //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3439     if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3440       if (track->GetNumberOfClusters()<maxn) continue;
3441       maxn = track->GetNumberOfClusters();
3442       //      if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2 *= track->GetChi2MIP(3); // RS    
3443       if (chi2<minchi2){
3444         minchi2=chi2;
3445         besttrack=track;
3446       }
3447     }
3448     else{
3449       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3450         delete array->RemoveAt(itrack);
3451       }  
3452     }
3453   }
3454   if (!besttrack) return;
3455   //
3456   //
3457   //take errors of best track as a reference
3458   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3459   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3460   for (Int_t j=0;j<6;j++) {
3461     if (besttrack->GetClIndex(j)>=0){
3462       erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3463       errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3464       ny[j]   = besttrack->GetNy(j);
3465       nz[j]   = besttrack->GetNz(j);
3466     }
3467   }
3468   //
3469   // calculate normalized chi2
3470   //
3471   Float_t * chi2        = new Float_t[entries];
3472   Int_t * index         = new Int_t[entries];  
3473   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3474   for (Int_t itrack=0;itrack<entries;itrack++){
3475     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3476     if (track){
3477       AliDebug(2,Form("track %d  ncls %d\n",itrack,track->GetNumberOfClusters()));      
3478       double chi2t = GetNormalizedChi2(track, mode);
3479       track->SetChi2MIP(0,chi2t);
3480       if (chi2t<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3481         if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3482         chi2[itrack] = chi2t;
3483       }
3484       else{
3485         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3486           delete array->RemoveAt(itrack);            
3487         }
3488       }
3489     }
3490   }
3491   //
3492   TMath::Sort(entries,chi2,index,kFALSE);
3493   besttrack = (AliITStrackMI*)array->At(index[0]);
3494   if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3495   if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3496     for (Int_t j=0;j<6;j++){
3497       if (besttrack->GetClIndex(j)>=0){
3498         erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3499         errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3500         ny[j]   = besttrack->GetNy(j);
3501         nz[j]   = besttrack->GetNz(j);
3502       }
3503     }
3504   }
3505   //
3506   // calculate one more time with updated normalized errors
3507   for (Int_t i=0;i<entries;i++) chi2[i] =10000;  
3508   for (Int_t itrack=0;itrack<entries;itrack++){
3509     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3510     if (track){      
3511       double chi2t = GetNormalizedChi2(track, mode);
3512       track->SetChi2MIP(0,chi2t);
3513       AliDebug(2,Form("track %d  ncls %d\n",itrack,track->GetNumberOfClusters()));            
3514       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3515         if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3516         chi2[itrack] = chi2t;  //-0*(track->GetNumberOfClusters()+track->GetNDeadZone()); 
3517       }
3518       else {
3519         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3520           delete array->RemoveAt(itrack);       
3521         }
3522       }
3523     }   
3524   }
3525   entries = array->GetEntriesFast();  
3526   //
3527   //
3528   if (entries>0){
3529     TObjArray * newarray = new TObjArray();  
3530     TMath::Sort(entries,chi2,index,kFALSE);
3531     besttrack = (AliITStrackMI*)array->At(index[0]);
3532     if (besttrack){
3533       AliDebug(2,Form("ncls best track %d     %f   %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3534       //
3535       for (Int_t j=0;j<6;j++){
3536         if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3537           erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3538           errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3539           ny[j]   = besttrack->GetNy(j);
3540           nz[j]   = besttrack->GetNz(j);
3541         }
3542       }
3543       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3544       minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3545       Float_t minn = besttrack->GetNumberOfClusters()-3;
3546       Int_t accepted=0;
3547       for (Int_t i=0;i<entries;i++){
3548         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);    
3549         if (!track) continue;
3550         if (accepted>maxcut) break;
3551         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3552         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3553           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3554             delete array->RemoveAt(index[i]);
3555             continue;
3556           }
3557         }
3558         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3559         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3560           if (!shortbest) accepted++;
3561           //
3562           newarray->AddLast(array->RemoveAt(index[i]));      
3563           for (Int_t j=0;j<6;j++){
3564             if (nz[j]==0){
3565               erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3566               errz[j] = track->GetSigmaZ(j); errz[j]   = track->GetSigmaZ(j+6);
3567               ny[j]   = track->GetNy(j);
3568               nz[j]   = track->GetNz(j);
3569             }
3570           }
3571         }
3572         else{
3573           delete array->RemoveAt(index[i]);
3574         }
3575       }
3576       array->Delete();
3577       delete fTrackHypothesys.RemoveAt(esdindex);
3578       fTrackHypothesys.AddAt(newarray,esdindex);
3579     }
3580     else{
3581       array->Delete();
3582       delete fTrackHypothesys.RemoveAt(esdindex);
3583     }
3584   }
3585   delete [] chi2;
3586   delete [] index;
3587 }
3588 //------------------------------------------------------------------------
3589 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3590 {
3591   //-------------------------------------------------------------
3592   // try to find best hypothesy
3593   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3594   // RS: optionally changing this to product of Chi2MIP(0)*Chi2MIP(3) == (chi2*chi2_interpolated)
3595   //-------------------------------------------------------------
3596   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3597   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3598   if (!array) return 0;
3599   Int_t entries = array->GetEntriesFast();
3600   if (!entries) return 0;  
3601   Float_t minchi2 = 100000;
3602   AliITStrackMI * besttrack=0;
3603   //
3604   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
3605   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3606   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};     
3607   Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3608   //
3609   for (Int_t i=0;i<entries;i++){    
3610     AliITStrackMI * track = (AliITStrackMI*)array->At(i);    
3611     if (!track) continue;
3612     Float_t sigmarfi,sigmaz;
3613     GetDCASigma(track,sigmarfi,sigmaz);
3614     track->SetDnorm(0,sigmarfi);
3615     track->SetDnorm(1,sigmaz);
3616     //
3617     track->SetChi2MIP(1,1000000);
3618     track->SetChi2MIP(2,1000000);
3619     track->SetChi2MIP(3,1000000);
3620     //
3621     // backtrack
3622     backtrack = new(backtrack) AliITStrackMI(*track); 
3623     if (track->GetConstrain()) {
3624       if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3625       if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3626         if (fUseImproveKalman) {if (!backtrack->ImproveKalman(xyzVtx,ersVtx,0,0,0)) continue;}
3627         else                   {if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;}
3628       }
3629       backtrack->ResetCovariance(10.);      
3630     }else{
3631       backtrack->ResetCovariance(10.);
3632     }
3633     backtrack->ResetClusters();
3634
3635     Double_t x = original->GetX();
3636     if (!RefitAt(x,backtrack,track)) continue;
3637     //
3638     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3639     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3640     if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.)  continue;
3641     track->SetChi22(GetMatchingChi2(backtrack,original));
3642
3643     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
3644     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
3645     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
3646
3647
3648     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1))  continue;
3649     //
3650     //forward track - without constraint
3651     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3652     forwardtrack->ResetClusters();
3653     x = track->GetX();
3654     if (!RefitAt(x,forwardtrack,track) && fSelectBestMIP03) continue;  // w/o fwd track MIP03 is meaningless
3655     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));    
3656     if  (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0)  continue;
3657     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2))  continue;
3658     
3659     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3660     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3661     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
3662     forwardtrack->SetD(0,track->GetD(0));
3663     forwardtrack->SetD(1,track->GetD(1));    
3664     {
3665       Int_t list[6];
3666       AliITSRecPoint* clist[6];
3667       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));      
3668       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3669     }
3670     
3671     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3672     if  ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;    
3673     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3674       track->SetChi2MIP(3,1000);
3675       continue; 
3676     }
3677     Double_t chi2 = track->GetChi2MIP(0); // +track->GetNUsed();    //RS
3678     if (fSelectBestMIP03) chi2 *= track->GetChi2MIP(3);
3679     else chi2 += track->GetNUsed();
3680     //
3681     for (Int_t ichi=0;ichi<5;ichi++){
3682       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3683     }
3684     if (chi2 < minchi2){
3685       //besttrack = new AliITStrackMI(*forwardtrack);
3686       besttrack = track;
3687       besttrack->SetLabel(track->GetLabel());
3688       besttrack->SetFakeRatio(track->GetFakeRatio());
3689       minchi2   = chi2;
3690       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3691       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3692       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
3693     }    
3694   }
3695   delete backtrack;
3696   delete forwardtrack;
3697
3698   if (!besttrack)  return 0;
3699
3700   Int_t accepted=0;
3701   for (Int_t i=0;i<entries;i++){    
3702     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3703    
3704     if (!track) continue;
3705     
3706     if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. || 
3707         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)
3708         // RS: don't apply this cut when fSelectBestMIP03 is on
3709         || (!fSelectBestMIP03 && (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.))
3710         ){
3711       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3712         delete array->RemoveAt(i);    
3713         continue;
3714       }
3715     }
3716     else{
3717       accepted++;
3718     }
3719   }
3720   //
3721   array->Compress();
3722   SortTrackHypothesys(esdindex,checkmax,1);
3723
3724   array = (TObjArray*) fTrackHypothesys.At(esdindex);
3725   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3726   besttrack = (AliITStrackMI*)array->At(0);  
3727   if (!besttrack)  return 0;
3728   besttrack->SetChi2MIP(8,0);
3729   fBestTrackIndex[esdindex]=0;
3730   entries = array->GetEntriesFast();
3731   AliITStrackMI *longtrack =0;
3732   minchi2 =1000;
3733   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3734   for (Int_t itrack=entries-1;itrack>0;itrack--) {
3735     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3736     if (!track->GetConstrain()) continue;
3737     if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3738     if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3739     if (track->GetChi2MIP(0)>4.) continue;
3740     minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3741     longtrack =track;
3742   }
3743   //if (longtrack) besttrack=longtrack;
3744   //
3745   // RS do shared cluster analysis here only if the new sharing analysis is not requested
3746   //RRR if (fFlagFakes) return besttrack;
3747
3748   Int_t list[6];
3749   AliITSRecPoint * clist[6];
3750   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3751   if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3752       &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){ 
3753     RegisterClusterTracks(besttrack,esdindex);
3754   }
3755   //
3756   //
3757   if (shared>0.0){
3758     Int_t nshared;
3759     Int_t overlist[6];
3760     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3761     if (sharedtrack>=0){
3762       //
3763       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5,original);     
3764       if (besttrack){
3765         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3766       }
3767       else return 0;
3768     }
3769   }  
3770   
3771   if (shared>2.5) return 0;
3772   if (shared>1.0) return besttrack;
3773   //
3774   // Don't sign clusters if not gold track
3775   //
3776   if (!besttrack->IsGoldPrimary()) return besttrack;
3777   if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
3778   //
3779   if (fConstraint[fPass]){
3780     //
3781     // sign clusters
3782     //
3783     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3784     for (Int_t i=0;i<6;i++){
3785       Int_t index = besttrack->GetClIndex(i);
3786       if (index<0) continue; 
3787       Int_t ilayer =  (index & 0xf0000000) >> 28;
3788       if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3789       AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);     
3790       if (!c) continue;
3791       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3792       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3793       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
3794       if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer)) 
3795         if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3796       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
3797
3798       Bool_t cansign = kTRUE;
3799       for (Int_t itrack=0;itrack<entries; itrack++){
3800         AliITStrackMI * track = (AliITStrackMI*)array->At(i);   
3801         if (!track) continue;
3802         if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3803         if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3804           cansign = kFALSE;
3805           break;
3806         }
3807       }
3808       if (cansign){
3809         if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3810         if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;    
3811         if (!c->IsUsed()) c->Use();
3812       }
3813     }
3814   }
3815   return besttrack;
3816
3817 //------------------------------------------------------------------------
3818 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3819 {
3820   //
3821   // get "best" hypothesys
3822   //
3823
3824   Int_t nentries = itsTracks.GetEntriesFast();
3825   for (Int_t i=0;i<nentries;i++){
3826     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3827     if (!track) continue;
3828     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3829     if (!array) continue;
3830     if (array->GetEntriesFast()<=0) continue;
3831     //
3832     AliITStrackMI* longtrack=0;
3833     Float_t minn=0;
3834     Float_t maxchi2=1000;
3835     for (Int_t j=0;j<array->GetEntriesFast();j++){
3836       AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3837       if (!trackHyp) continue;
3838       if (trackHyp->GetGoldV0()) {
3839         longtrack = trackHyp;   //gold V0 track taken
3840         break;
3841       }
3842       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3843       Float_t chi2 = trackHyp->GetChi2MIP(0);
3844       if (fSelectBestMIP03) chi2 *= trackHyp->GetChi2MIP(3);
3845       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = chi2; //trackHyp->GetChi2MIP(0);
3846       //
3847       if (fAfterV0){ // ??? RS
3848         if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3849       }
3850       if (chi2 > maxchi2) continue;
3851       minn = trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3852       if (fSelectBestMIP03) minn++; // allow next to longest to win
3853       maxchi2 = chi2;
3854       longtrack=trackHyp;
3855     }    
3856     //
3857     //
3858     //
3859     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3860     if (!longtrack) {longtrack = besttrack;}
3861     else besttrack= longtrack;
3862     //
3863     if (besttrack) {
3864       Int_t list[6];
3865       AliITSRecPoint * clist[6];
3866       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3867       //
3868       track->SetNUsed(shared);      
3869       track->SetNSkipped(besttrack->GetNSkipped());
3870       track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3871       if (shared>0) {
3872         if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3873         Int_t nshared;
3874         Int_t overlist[6]; 
3875         //
3876         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3877         //if (sharedtrack==-1) sharedtrack=0;
3878         if (sharedtrack>=0) {       
3879           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5,track);                         
3880         }
3881       }   
3882       if (besttrack&&fAfterV0) {
3883         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3884         track->SetWinner(besttrack);
3885       }
3886       if (besttrack) {
3887         if (fConstraint[fPass]) {
3888           UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3889           track->SetWinner(besttrack);
3890         }
3891         if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3892           if ( TMath::Abs(besttrack->GetD(0))>0.1 || 
3893                TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);    
3894         }       
3895       }
3896     }
3897   }
3898
3899
3900 //------------------------------------------------------------------------
3901 void AliITStrackerMI::FlagFakes(TObjArray &itsTracks)
3902 {
3903   //
3904   // RS: flag those tracks which are suxpected to have fake clusters
3905   //
3906   const double kThreshPt = 0.5;
3907   AliRefArray *refArr[6];
3908   //
3909   for (int i=0;i<6;i++) {
3910     int ncl = fgLayers[i].GetNumberOfClusters();
3911     refArr[i] = new AliRefArray(ncl,TMath::Min(ncl,1000));
3912   }
3913   Int_t nentries = itsTracks.GetEntriesFast();
3914   //
3915   // fill cluster->track associations
3916   for (Int_t itr=0;itr<nentries;itr++){
3917     AliITStrackMI* track = (AliITStrackMI*)itsTracks.UncheckedAt(itr);   
3918     if (!track) continue;
3919     AliITStrackMI* trackITS = track->GetWinner();
3920     if (!trackITS) continue;
3921     for (int il=trackITS->GetNumberOfClusters();il--;) {
3922       int idx = trackITS->GetClusterIndex(il);
3923       Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3924       //      if (c>fgLayers[l].GetNumberOfClusters()) continue;
3925       refArr[l]->AddReference(c, itr);
3926     }
3927   }
3928   //
3929   const UInt_t kMaxRef = 100;
3930   UInt_t crefs[kMaxRef];
3931   Int_t ncrefs=0;
3932   // process tracks with shared clusters
3933   for (int itr=0;itr<nentries;itr++){
3934     AliITStrackMI* track0 = (AliITStrackMI*)itsTracks.UncheckedAt(itr);  
3935     AliITStrackMI* trackH0 = track0->GetWinner(); 
3936     if (!trackH0) continue;
3937     AliESDtrack* esd0 = track0->GetESDtrack();
3938     //
3939     for (int il=0;il<trackH0->GetNumberOfClusters();il++) {
3940       int idx = trackH0->GetClusterIndex(il);
3941       Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3942       ncrefs = refArr[l]->GetReferences(c,crefs,kMaxRef);                
3943       if (ncrefs<2) continue; // there will be always self-reference, for sharing needs at least 2
3944       esd0->SetITSSharedFlag(l); 
3945       for (int ir=ncrefs;ir--;) {
3946         if (int(crefs[ir]) <= itr) continue; // ==:selfreference, <: the same pair will be checked with >
3947         AliITStrackMI* track1 = (AliITStrackMI*)itsTracks.UncheckedAt(crefs[ir]);
3948         AliITStrackMI* trackH1 = track1->GetWinner(); 
3949         AliESDtrack* esd1 = track1->GetESDtrack();
3950         esd1->SetITSSharedFlag(l);
3951         //
3952         double pt0 = trackH0->Pt(), pt1 = trackH1->Pt(), res = 0.;      
3953         if      (pt0>kThreshPt && pt0-pt1>0.2+0.2*(pt0-kThreshPt) ) res = -100;
3954         else if (pt1>kThreshPt && pt1-pt0>0.2+0.2*(pt1-kThreshPt) ) res = 100;
3955
3956         // select the one with smallest chi2's product
3957         res += trackH0->GetChi2MIP(0)*trackH0->GetChi2MIP(3); 
3958         res -= trackH1->GetChi2MIP(0)*trackH1->GetChi2MIP(3);
3959         //
3960         if (res<0) esd1->SetITSFakeFlag();  // esd0 is winner
3961         else       esd0->SetITSFakeFlag();  // esd1 is winner
3962       }
3963       //
3964     }
3965     //
3966   }
3967   //
3968   for (int i=6;i--;) delete refArr[i];
3969 }
3970
3971 //------------------------------------------------------------------------
3972 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3973   //--------------------------------------------------------------------
3974   //This function "cooks" a track label. If label<0, this track is fake.
3975   //--------------------------------------------------------------------
3976   Int_t tpcLabel=-1; 
3977      
3978   if (track->GetESDtrack()){
3979     tpcLabel = track->GetESDtrack()->GetTPCLabel();
3980     ULong_t trStatus=track->GetESDtrack()->GetStatus();
3981     if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
3982   }
3983    track->SetChi2MIP(9,0);
3984    Int_t nwrong=0;
3985    for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3986      Int_t cindex = track->GetClusterIndex(i);
3987      Int_t l=(cindex & 0xf0000000) >> 28;
3988      AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3989      Int_t isWrong=1;
3990      for (Int_t ind=0;ind<3;ind++){
3991        if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
3992        //AliDebug(2,Form("icl %d  ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3993      }
3994      track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3995      nwrong+=isWrong;
3996    }
3997    Int_t nclusters = track->GetNumberOfClusters();
3998    if (nclusters > 0) //PH Some tracks don't have any cluster
3999      track->SetFakeRatio(double(nwrong)/double(nclusters));
4000    if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
4001      track->SetLabel(-tpcLabel);
4002    } else {
4003      track->SetLabel(tpcLabel);
4004    }
4005    AliDebug(2,Form(" nls %d wrong %d  label %d  tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
4006    
4007 }
4008 //------------------------------------------------------------------------
4009 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
4010   //
4011   // Fill the dE/dx in this track
4012   //
4013   track->SetChi2MIP(9,0);
4014   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
4015     Int_t cindex = track->GetClusterIndex(i);
4016     Int_t l=(cindex & 0xf0000000) >> 28;
4017     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4018     Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
4019     Int_t isWrong=1;
4020     for (Int_t ind=0;ind<3;ind++){
4021       if (cl->GetLabel(ind)==lab) isWrong=0;
4022     }
4023     track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4024   }
4025   Double_t low=0.;
4026   Double_t up=0.51;    
4027   track->CookdEdx(low,up);
4028 }
4029 //------------------------------------------------------------------------
4030 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
4031   //
4032   // Create some arrays
4033   //
4034   if (fCoefficients) delete []fCoefficients;
4035   fCoefficients = new Float_t[ntracks*48];
4036   for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
4037 }
4038 //------------------------------------------------------------------------
4039 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer) 
4040 {
4041   //
4042   // Compute predicted chi2
4043   //
4044   // Take into account the mis-alignment (bring track to cluster plane)
4045   Double_t xTrOrig=track->GetX();
4046   if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
4047   Float_t erry,errz,covyz;
4048   Float_t theta = track->GetTgl();
4049   Float_t phi   = track->GetSnp();
4050   phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
4051   AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
4052   AliDebug(2,Form(" chi2: tr-cl   %f  %f   tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
4053   AliDebug(2,Form(" chi2: tr-cl   %f  %f   tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
4054   Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
4055   // Bring the track back to detector plane in ideal geometry
4056   // [mis-alignment will be accounted for in UpdateMI()]
4057   if (!track->Propagate(xTrOrig)) return 1000.;
4058   Float_t ny,nz;
4059   AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);  
4060   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
4061   if (delta>1){
4062     chi2+=0.5*TMath::Min(delta/2,2.);
4063     chi2+=2.*cluster->GetDeltaProbability();
4064   }
4065   //
4066   track->SetNy(layer,ny);
4067   track->SetNz(layer,nz);
4068   track->SetSigmaY(layer,erry);
4069   track->SetSigmaZ(layer, errz);
4070   track->SetSigmaYZ(layer,covyz);
4071   //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
4072   track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
4073   return chi2;
4074
4075 }
4076 //------------------------------------------------------------------------
4077 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const 
4078 {
4079   //
4080   // Update ITS track
4081   //
4082   Int_t layer = (index & 0xf0000000) >> 28;
4083   track->SetClIndex(layer, index);
4084   if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
4085     if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
4086       chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
4087       track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
4088     }
4089   }
4090
4091   if (TMath::Abs(cl->GetQ())<1.e-13) return 0;  // ingore the "virtual" clusters
4092
4093
4094   // Take into account the mis-alignment (bring track to cluster plane)
4095   Double_t xTrOrig=track->GetX();
4096   Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
4097   AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
4098   AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
4099   AliDebug(2,Form(" xtr %f  xcl %f",track->GetX(),cl->GetX()));
4100
4101   if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
4102   
4103   AliCluster c(*cl);
4104   c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
4105   c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
4106   c.SetSigmaYZ(track->GetSigmaYZ(layer));
4107
4108
4109   Int_t updated = track->UpdateMI(&c,chi2,index);
4110   // Bring the track back to detector plane in ideal geometry
4111   if (!track->Propagate(xTrOrig)) return 0;
4112  
4113   if(!updated) AliDebug(2,"update failed");
4114   return updated;
4115 }
4116
4117 //------------------------------------------------------------------------
4118 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
4119 {