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