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