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