Fix for bug 71598 (Marco Van Leeuwen)
[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               // MvL: added here as well
1173               // bring track back to ideal detector plane
1174               currenttrack->Propagate(xTrOrig);
1175               continue;
1176             }
1177           // bring track back to ideal detector plane
1178           if (!currenttrack->Propagate(xTrOrig)) continue;
1179         } else {                                      // have to move track to cluster's detector
1180           const AliITSdetector &detc=layer.GetDetector(idetc);
1181           // a first cut on track-cluster distance
1182           Double_t y;
1183           if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1184           if ( (z-cl->GetZ())*(z-cl->GetZ())*msz + 
1185                (y-cl->GetY())*(y-cl->GetY())*msy > 1. ) 
1186             continue; // cluster not associated to track
1187           //
1188           new (&backuptrack) AliITStrackMI(currenttrack2);
1189           changedet = kTRUE;
1190           currenttrack =&currenttrack2;
1191           if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1192             new (currenttrack) AliITStrackMI(backuptrack);
1193             changedet = kFALSE;
1194             continue;
1195           }
1196           currenttrack->SetDetectorIndex(idetc);
1197           // Get again the budget to the primary vertex 
1198           // for the current track being prolonged, if had to change detector 
1199           //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1200         }
1201
1202         // calculate track-clusters chi2
1203         chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer); 
1204         // chi2 cut
1205         AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1206         if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1207           if (TMath::Abs(cl->GetQ())<1.e-13) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster       
1208           if (ntracks[ilayer]>=100) continue;
1209           AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1210           updatetrack->SetClIndex(ilayer,-1);
1211           if (changedet) new (&currenttrack2) AliITStrackMI(backuptrack);
1212
1213           if (TMath::Abs(cl->GetQ())>1.e-13) { // real cluster
1214             if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1215               AliDebug(2,"update failed");
1216               continue;
1217             } 
1218             updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2); 
1219             modstatus = 1; // found
1220           } else {             // virtual cluster in dead zone
1221             updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1222             updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1223             modstatus = 7; // holes in z in SPD
1224           }
1225
1226           if (changedet) {
1227             Float_t xlocnewdet,zlocnewdet;
1228             if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1229               updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1230             }
1231           } else {
1232             updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1233           }
1234           if (cl->IsUsed()) updatetrack->IncrementNUsed();
1235
1236           // apply correction for material of the current layer
1237           CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1238
1239           if (constrain) { // apply vertex constrain
1240             updatetrack->SetConstrain(constrain);
1241             Bool_t isPrim = kTRUE;
1242             if (ilayer<4) { // check that it's close to the vertex
1243               updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1244               if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1245                   AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() || 
1246                   TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1247                   AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1248             }
1249             if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1250           } //apply vertex constrain              
1251           ntracks[ilayer]++;
1252         }  // create new hypothesis
1253         else {
1254           AliDebug(2,"chi2 too large");
1255         }
1256
1257       } // loop over possible prolongations 
1258      
1259       // allow one prolongation without clusters
1260       if (constrain && itrack<=1 && TMath::Abs(currenttrack1.GetNSkipped())<1.e-13 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1261         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1262         // apply correction for material of the current layer
1263         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1264         vtrack->SetClIndex(ilayer,-1);
1265         modstatus = 3; // skipped 
1266         vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1267         vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1268         vtrack->IncrementNSkipped();
1269         ntracks[ilayer]++;
1270       }
1271       
1272
1273     } // loop over tracks in layer ilayer+1
1274
1275     //loop over track candidates for the current layer
1276     //
1277     //
1278     Int_t accepted=0;
1279     
1280     Int_t golden=0;
1281     for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1282       normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer); 
1283       if (normalizedchi2[itrack] < 
1284           AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1285       if (ilayer>4) {
1286         accepted++;
1287       } else {
1288         if (constrain) { // constrain
1289           if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1) 
1290             accepted++;
1291         } else {        // !constrain
1292           if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1) 
1293             accepted++;
1294         }
1295       }
1296     }
1297     // sort tracks by increasing normalized chi2
1298     TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE); 
1299     ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1300     if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1301     if (ntracks[ilayer]>90) ntracks[ilayer]=90; 
1302   } // end loop over layers
1303
1304
1305   //
1306   // Now select tracks to be kept
1307   //
1308   Int_t max = constrain ? 20 : 5;
1309
1310   // tracks that reach layer 0 (SPD inner)
1311   for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1312     AliITStrackMI & track= tracks[0][nindexes[0][i]];
1313     if (track.GetNumberOfClusters()<2) continue;
1314     if (!constrain && track.GetNormChi2(0) >
1315         AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1316       continue;
1317     }
1318     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1319   }
1320
1321   // tracks that reach layer 1 (SPD outer)
1322   for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1323     AliITStrackMI & track= tracks[1][nindexes[1][i]];
1324     if (track.GetNumberOfClusters()<4) continue;
1325     if (!constrain && track.GetNormChi2(1) >
1326         AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1327     if (constrain) track.IncrementNSkipped();
1328     if (!constrain) {
1329       track.SetD(0,track.GetD(GetX(),GetY()));   
1330       track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1331       if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1332         track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1333       }
1334     }
1335     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1336   }
1337
1338   // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1339   if (!constrain){  
1340     for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1341       AliITStrackMI & track= tracks[2][nindexes[2][i]];
1342       if (track.GetNumberOfClusters()<3) continue;
1343       if (!constrain && track.GetNormChi2(2) >
1344           AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1345       if (constrain) track.SetNSkipped(track.GetNSkipped()+2);      
1346       if (!constrain){
1347         track.SetD(0,track.GetD(GetX(),GetY()));
1348         track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1349         if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1350           track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1351         }
1352       }
1353       AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1354     }
1355   }
1356   
1357   if (!constrain) {
1358     //
1359     // register best track of each layer - important for V0 finder
1360     //
1361     for (Int_t ilayer=0;ilayer<5;ilayer++){
1362       if (ntracks[ilayer]==0) continue;
1363       AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1364       if (track.GetNumberOfClusters()<1) continue;
1365       CookLabel(&track,0);
1366       bestarray->AddAt(new AliITStrackMI(track),ilayer);
1367     }
1368   }
1369   //
1370   // update TPC V0 information
1371   //
1372   if (otrack->GetESDtrack()->GetV0Index(0)>0){    
1373     Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1374     for (Int_t i=0;i<3;i++){
1375       Int_t  index = otrack->GetESDtrack()->GetV0Index(i); 
1376       if (index==0) break;
1377       AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1378       if (vertex->GetStatus()<0) continue;     // rejected V0
1379       //
1380       if (otrack->GetSign()>0) {
1381         vertex->SetIndex(0,esdindex);
1382       }
1383       else{
1384         vertex->SetIndex(1,esdindex);
1385       }
1386       //find nearest layer with track info
1387       Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
1388       Int_t nearestold  = GetNearestLayer(xrp);               //I.B.
1389       Int_t nearest     = nearestold; 
1390       for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1391         if (ntracks[nearest]==0){
1392           nearest = ilayer;
1393         }
1394       }
1395       //
1396       AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1397       if (nearestold<5&&nearest<5){
1398         Bool_t accept = track.GetNormChi2(nearest)<10; 
1399         if (accept){
1400           if (track.GetSign()>0) {
1401             vertex->SetParamP(track);
1402             vertex->Update(fprimvertex);
1403             //vertex->SetIndex(0,track.fESDtrack->GetID()); 
1404             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1405           }else{
1406             vertex->SetParamN(track);
1407             vertex->Update(fprimvertex);
1408             //vertex->SetIndex(1,track.fESDtrack->GetID());
1409             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1410           }
1411           vertex->SetStatus(vertex->GetStatus()+1);
1412         }else{
1413           //vertex->SetStatus(-2);  // reject V0  - not enough clusters
1414         }
1415       }
1416     }
1417   } 
1418   
1419 }
1420 //------------------------------------------------------------------------
1421 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1422 {
1423   //--------------------------------------------------------------------
1424   //
1425   //
1426   return fgLayers[layer];
1427 }
1428 //------------------------------------------------------------------------
1429 AliITStrackerMI::AliITSlayer::AliITSlayer():
1430 fR(0),
1431 fPhiOffset(0),
1432 fNladders(0),
1433 fZOffset(0),
1434 fNdetectors(0),
1435 fDetectors(0),
1436 fN(0),
1437 fDy5(0),
1438 fDy10(0),
1439 fDy20(0),
1440 fClustersCs(0),
1441 fClusterIndexCs(0),
1442 fYcs(0),
1443 fZcs(0),
1444 fNcs(0),
1445 fCurrentSlice(-1),
1446 fZmin(0),
1447 fZmax(0),
1448 fYmin(0),
1449 fYmax(0),
1450 fI(0),
1451 fImax(0),
1452 fSkip(0),
1453 fAccepted(0),
1454 fRoad(0),
1455 fMaxSigmaClY(0),
1456 fMaxSigmaClZ(0),
1457 fNMaxSigmaCl(3)
1458 {
1459   //--------------------------------------------------------------------
1460   //default AliITSlayer constructor
1461   //--------------------------------------------------------------------
1462   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1463     fClusterWeight[i]=0;
1464     fClusterTracks[0][i]=-1;
1465     fClusterTracks[1][i]=-1;
1466     fClusterTracks[2][i]=-1;    
1467     fClusterTracks[3][i]=-1;    
1468   }
1469 }
1470 //------------------------------------------------------------------------
1471 AliITStrackerMI::AliITSlayer::
1472 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1473 fR(r),
1474 fPhiOffset(p),
1475 fNladders(nl),
1476 fZOffset(z),
1477 fNdetectors(nd),
1478 fDetectors(0),
1479 fN(0),
1480 fDy5(0),
1481 fDy10(0),
1482 fDy20(0),
1483 fClustersCs(0),
1484 fClusterIndexCs(0),
1485 fYcs(0),
1486 fZcs(0),
1487 fNcs(0),
1488 fCurrentSlice(-1),
1489 fZmin(0),
1490 fZmax(0),
1491 fYmin(0),
1492 fYmax(0),
1493 fI(0),
1494 fImax(0),
1495 fSkip(0),
1496 fAccepted(0),
1497 fRoad(0),
1498 fMaxSigmaClY(0),
1499 fMaxSigmaClZ(0),
1500 fNMaxSigmaCl(3) {
1501   //--------------------------------------------------------------------
1502   //main AliITSlayer constructor
1503   //--------------------------------------------------------------------
1504   fDetectors=new AliITSdetector[fNladders*fNdetectors];
1505   fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1506 }
1507 //------------------------------------------------------------------------
1508 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1509 fR(layer.fR),
1510 fPhiOffset(layer.fPhiOffset),
1511 fNladders(layer.fNladders),
1512 fZOffset(layer.fZOffset),
1513 fNdetectors(layer.fNdetectors),
1514 fDetectors(layer.fDetectors),
1515 fN(layer.fN),
1516 fDy5(layer.fDy5),
1517 fDy10(layer.fDy10),
1518 fDy20(layer.fDy20),
1519 fClustersCs(layer.fClustersCs),
1520 fClusterIndexCs(layer.fClusterIndexCs),
1521 fYcs(layer.fYcs),
1522 fZcs(layer.fZcs),
1523 fNcs(layer.fNcs),
1524 fCurrentSlice(layer.fCurrentSlice),
1525 fZmin(layer.fZmin),
1526 fZmax(layer.fZmax),
1527 fYmin(layer.fYmin),
1528 fYmax(layer.fYmax),
1529 fI(layer.fI),
1530 fImax(layer.fImax),
1531 fSkip(layer.fSkip),
1532 fAccepted(layer.fAccepted),
1533 fRoad(layer.fRoad),
1534 fMaxSigmaClY(layer.fMaxSigmaClY),
1535 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1536 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1537 {
1538   //Copy constructor
1539 }
1540 //------------------------------------------------------------------------
1541 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1542   //--------------------------------------------------------------------
1543   // AliITSlayer destructor
1544   //--------------------------------------------------------------------
1545   delete [] fDetectors;
1546   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1547   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1548     fClusterWeight[i]=0;
1549     fClusterTracks[0][i]=-1;
1550     fClusterTracks[1][i]=-1;
1551     fClusterTracks[2][i]=-1;    
1552     fClusterTracks[3][i]=-1;    
1553   }
1554 }
1555 //------------------------------------------------------------------------
1556 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1557   //--------------------------------------------------------------------
1558   // This function removes loaded clusters
1559   //--------------------------------------------------------------------
1560   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1561   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1562     fClusterWeight[i]=0;
1563     fClusterTracks[0][i]=-1;
1564     fClusterTracks[1][i]=-1;
1565     fClusterTracks[2][i]=-1;    
1566     fClusterTracks[3][i]=-1;  
1567   }
1568   
1569   fN=0;
1570   fI=0;
1571 }
1572 //------------------------------------------------------------------------
1573 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1574   //--------------------------------------------------------------------
1575   // This function reset weights of the clusters
1576   //--------------------------------------------------------------------
1577   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1578     fClusterWeight[i]=0;
1579     fClusterTracks[0][i]=-1;
1580     fClusterTracks[1][i]=-1;
1581     fClusterTracks[2][i]=-1;    
1582     fClusterTracks[3][i]=-1;  
1583   }
1584   for (Int_t i=0; i<fN;i++) {
1585     AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1586     if (cl&&cl->IsUsed()) cl->Use();
1587   }
1588
1589 }
1590 //------------------------------------------------------------------------
1591 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1592   //--------------------------------------------------------------------
1593   // This function calculates the road defined by the cluster density
1594   //--------------------------------------------------------------------
1595   Int_t n=0;
1596   for (Int_t i=0; i<fN; i++) {
1597      if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1598   }
1599   if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1600 }
1601 //------------------------------------------------------------------------
1602 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1603   //--------------------------------------------------------------------
1604   //This function adds a cluster to this layer
1605   //--------------------------------------------------------------------
1606   if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1607     return 1;
1608   }
1609   fCurrentSlice=-1;
1610   fClusters[fN]=cl;
1611   fN++;
1612   AliITSdetector &det=GetDetector(cl->GetDetectorIndex());    
1613   //AD
1614   Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1615   Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2()); 
1616   if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1617   if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1618   if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1619   if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1620   //AD               
1621   /*
1622   if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1623   if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1624   if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1625   if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1626   */                 
1627   return 0;
1628 }
1629 //------------------------------------------------------------------------
1630 void  AliITStrackerMI::AliITSlayer::SortClusters()
1631 {
1632   //
1633   //sort clusters
1634   //
1635   AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1636   Float_t *z                = new Float_t[fN];
1637   Int_t   * index           = new Int_t[fN];
1638   //
1639   fMaxSigmaClY=0.; //AD
1640   fMaxSigmaClZ=0.; //AD
1641
1642   for (Int_t i=0;i<fN;i++){
1643     z[i] = fClusters[i]->GetZ();
1644     // save largest errors in y and z for this layer
1645     fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1646     fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1647   }
1648   TMath::Sort(fN,z,index,kFALSE);
1649   for (Int_t i=0;i<fN;i++){
1650     clusters[i] = fClusters[index[i]];
1651   }
1652   //
1653   for (Int_t i=0;i<fN;i++){
1654     fClusters[i] = clusters[i];
1655     fZ[i]        = fClusters[i]->GetZ();
1656     AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());    
1657     Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1658     if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1659     fY[i] = y;
1660   }
1661   delete[] index;
1662   delete[] z;
1663   delete[] clusters;
1664   //
1665
1666   fYB[0]=10000000;
1667   fYB[1]=-10000000;
1668   for (Int_t i=0;i<fN;i++){
1669     if (fY[i]<fYB[0]) fYB[0]=fY[i];
1670     if (fY[i]>fYB[1]) fYB[1]=fY[i];
1671     fClusterIndex[i] = i;
1672   }
1673   //
1674   // fill slices
1675   fDy5 = (fYB[1]-fYB[0])/5.;
1676   fDy10 = (fYB[1]-fYB[0])/10.;
1677   fDy20 = (fYB[1]-fYB[0])/20.;
1678   for (Int_t i=0;i<6;i++)  fN5[i] =0;  
1679   for (Int_t i=0;i<11;i++) fN10[i]=0;  
1680   for (Int_t i=0;i<21;i++) fN20[i]=0;
1681   //  
1682   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;}
1683   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;} 
1684   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;}
1685   //
1686   //
1687   for (Int_t i=0;i<fN;i++)
1688     for (Int_t irot=-1;irot<=1;irot++){
1689       Float_t curY = fY[i]+irot*TMath::TwoPi()*fR; 
1690       // slice 5
1691       for (Int_t slice=0; slice<6;slice++){
1692         if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1693           fClusters5[slice][fN5[slice]] = fClusters[i];
1694           fY5[slice][fN5[slice]] = curY;
1695           fZ5[slice][fN5[slice]] = fZ[i];
1696           fClusterIndex5[slice][fN5[slice]]=i;
1697           fN5[slice]++;
1698         }
1699       }
1700       // slice 10
1701       for (Int_t slice=0; slice<11;slice++){
1702         if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1703           fClusters10[slice][fN10[slice]] = fClusters[i];
1704           fY10[slice][fN10[slice]] = curY;
1705           fZ10[slice][fN10[slice]] = fZ[i];
1706           fClusterIndex10[slice][fN10[slice]]=i;
1707           fN10[slice]++;
1708         }
1709       }
1710       // slice 20
1711       for (Int_t slice=0; slice<21;slice++){
1712         if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1713           fClusters20[slice][fN20[slice]] = fClusters[i];
1714           fY20[slice][fN20[slice]] = curY;
1715           fZ20[slice][fN20[slice]] = fZ[i];
1716           fClusterIndex20[slice][fN20[slice]]=i;
1717           fN20[slice]++;
1718         }
1719       }      
1720     }
1721
1722   //
1723   // consistency check
1724   //
1725   for (Int_t i=0;i<fN-1;i++){
1726     if (fZ[i]>fZ[i+1]){
1727       printf("Bug\n");
1728     }
1729   }
1730   //
1731   for (Int_t slice=0;slice<21;slice++)
1732   for (Int_t i=0;i<fN20[slice]-1;i++){
1733     if (fZ20[slice][i]>fZ20[slice][i+1]){
1734       printf("Bug\n");
1735     }
1736   }
1737
1738
1739 }
1740 //------------------------------------------------------------------------
1741 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1742   //--------------------------------------------------------------------
1743   // This function returns the index of the nearest cluster 
1744   //--------------------------------------------------------------------
1745   Int_t ncl=0;
1746   const Float_t *zcl;  
1747   if (fCurrentSlice<0) {
1748     ncl = fN;
1749     zcl   = fZ;
1750   }
1751   else{
1752     ncl   = fNcs;
1753     zcl   = fZcs;;
1754   }
1755   
1756   if (ncl==0) return 0;
1757   Int_t b=0, e=ncl-1, m=(b+e)/2;
1758   for (; b<e; m=(b+e)/2) {
1759     //    if (z > fClusters[m]->GetZ()) b=m+1;
1760     if (z > zcl[m]) b=m+1;
1761     else e=m; 
1762   }
1763   return m;
1764 }
1765 //------------------------------------------------------------------------
1766 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 {
1767   //--------------------------------------------------------------------
1768   // This function computes the rectangular road for this track
1769   //--------------------------------------------------------------------
1770
1771
1772   AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1773   // take into account the misalignment: propagate track to misaligned detector plane
1774   if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1775
1776   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1777                     TMath::Sqrt(track->GetSigmaZ2() + 
1778                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1779                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1780                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1781   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1782                     TMath::Sqrt(track->GetSigmaY2() + 
1783                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1784                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1785                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1786       
1787   // track at boundary between detectors, enlarge road
1788   Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1789   if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) || 
1790        (track->GetY()+dy > det.GetYmax()-boundaryWidth) || 
1791        (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1792        (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1793     Float_t tgl = TMath::Abs(track->GetTgl());
1794     if (tgl > 1.) tgl=1.;
1795     Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1796     dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1797     Float_t snp = TMath::Abs(track->GetSnp());
1798     if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1799     dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1800   } // boundary
1801   
1802   // add to the road a term (up to 2-3 mm) to deal with misalignments
1803   dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1804   dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1805
1806   Double_t r = fgLayers[ilayer].GetR();
1807   zmin = track->GetZ() - dz; 
1808   zmax = track->GetZ() + dz;
1809   ymin = track->GetY() + r*det.GetPhi() - dy;
1810   ymax = track->GetY() + r*det.GetPhi() + dy;
1811
1812   // bring track back to idead detector plane
1813   if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1814
1815   return kTRUE;
1816 }
1817 //------------------------------------------------------------------------
1818 void AliITStrackerMI::AliITSlayer::
1819 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1820   //--------------------------------------------------------------------
1821   // This function sets the "window"
1822   //--------------------------------------------------------------------
1823  
1824   Double_t circle=2*TMath::Pi()*fR;
1825   fYmin = ymin; 
1826   fYmax = ymax;
1827   fZmin = zmin;
1828   fZmax = zmax;
1829   // AD
1830   // enlarge road in y by maximum cluster error on this layer (3 sigma)
1831   fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
1832   fYmax += fNMaxSigmaCl*fMaxSigmaClY;
1833   fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
1834   fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
1835
1836   Float_t ymiddle = (fYmin+fYmax)*0.5;
1837   if (ymiddle<fYB[0]) {
1838     fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1839   } else if (ymiddle>fYB[1]) {
1840     fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1841   }
1842   
1843   //
1844   fCurrentSlice =-1;
1845   // defualt take all
1846   fClustersCs = fClusters;
1847   fClusterIndexCs = fClusterIndex;
1848   fYcs  = fY;
1849   fZcs  = fZ;
1850   fNcs  = fN;
1851   //
1852   //is in 20 slice?
1853   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1854     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1855     if (slice<0) slice=0;
1856     if (slice>20) slice=20;
1857     Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1858     if (isOK) {
1859       fCurrentSlice=slice;
1860       fClustersCs = fClusters20[fCurrentSlice];
1861       fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1862       fYcs  = fY20[fCurrentSlice];
1863       fZcs  = fZ20[fCurrentSlice];
1864       fNcs  = fN20[fCurrentSlice];
1865     }
1866   }  
1867   //
1868   //is in 10 slice?
1869   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1870     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1871     if (slice<0) slice=0;
1872     if (slice>10) slice=10;
1873     Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1874     if (isOK) {
1875       fCurrentSlice=slice;
1876       fClustersCs = fClusters10[fCurrentSlice];
1877       fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1878       fYcs  = fY10[fCurrentSlice];
1879       fZcs  = fZ10[fCurrentSlice];
1880       fNcs  = fN10[fCurrentSlice];
1881     }
1882   }  
1883   //
1884   //is in 5 slice?
1885   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1886     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1887     if (slice<0) slice=0;
1888     if (slice>5) slice=5;
1889     Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1890     if (isOK) {
1891       fCurrentSlice=slice;
1892       fClustersCs = fClusters5[fCurrentSlice];
1893       fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1894       fYcs  = fY5[fCurrentSlice];
1895       fZcs  = fZ5[fCurrentSlice];
1896       fNcs  = fN5[fCurrentSlice];
1897     }
1898   }  
1899   //  
1900   fI        = FindClusterIndex(fZmin); 
1901   fImax     = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
1902   fSkip     = 0;
1903   fAccepted = 0;
1904
1905   return;
1906 }
1907 //------------------------------------------------------------------------
1908 Int_t AliITStrackerMI::AliITSlayer::
1909 FindDetectorIndex(Double_t phi, Double_t z) const {
1910   //--------------------------------------------------------------------
1911   //This function finds the detector crossed by the track
1912   //--------------------------------------------------------------------
1913   Double_t dphi;
1914   if (fZOffset<0)            // old geometry
1915     dphi = -(phi-fPhiOffset);
1916   else                       // new geometry
1917     dphi = phi-fPhiOffset;
1918
1919
1920   if      (dphi <  0) dphi += 2*TMath::Pi();
1921   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1922   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1923   if (np>=fNladders) np-=fNladders;
1924   if (np<0)          np+=fNladders;
1925
1926
1927   Double_t dz=fZOffset-z;
1928   Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
1929   Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
1930   if (nz>=fNdetectors || nz<0) {
1931     //printf("ndet %d phi %f z %f  np %d nz %d\n",fNdetectors,phi,z,np,nz);
1932     return -1;
1933   }
1934
1935   // ad hoc correction for 3rd ladder of SDD inner layer,
1936   // which is reversed (rotated by pi around local y)
1937   // this correction is OK only from AliITSv11Hybrid onwards
1938   if (GetR()>12. && GetR()<20.) { // SDD inner
1939     if(np==2) { // 3rd ladder
1940       Double_t posMod252[3];
1941       AliITSgeomTGeo::GetTranslation(252,posMod252);
1942       // check the Z coordinate of Mod 252: if negative 
1943       // (old SDD geometry in AliITSv11Hybrid)
1944       // the swap of numeration whould be applied
1945       if(posMod252[2]<0.){
1946         nz = (fNdetectors-1) - nz;
1947       } 
1948     }
1949   }
1950   //printf("ndet %d phi %f z %f  np %d nz %d\n",fNdetectors,phi,z,np,nz);
1951
1952
1953   return np*fNdetectors + nz;
1954 }
1955 //------------------------------------------------------------------------
1956 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
1957 {
1958   //--------------------------------------------------------------------
1959   // This function returns clusters within the "window" 
1960   //--------------------------------------------------------------------
1961
1962   if (fCurrentSlice<0) {
1963     Double_t rpi2 = 2.*fR*TMath::Pi();
1964     for (Int_t i=fI; i<fImax; i++) {
1965       Double_t y = fY[i];
1966       Double_t z = fZ[i];
1967       if (fYmax<y) y -= rpi2;
1968       if (fYmin>y) y += rpi2;
1969       if (y<fYmin) continue;
1970       if (y>fYmax) continue;
1971       // AD
1972       // skip clusters that are in "extended" road but they 
1973       // 3sigma error does not touch the original road
1974       if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
1975       if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
1976       //
1977       if (TMath::Abs(fClusters[i]->GetQ())<1.e-13 && fSkip==2) continue;
1978       ci=i;
1979       if (!test) fI=i+1;
1980       return fClusters[i];
1981     }
1982   } else {
1983     for (Int_t i=fI; i<fImax; i++) {
1984       if (fYcs[i]<fYmin) continue;
1985       if (fYcs[i]>fYmax) continue;
1986       if (TMath::Abs(fClustersCs[i]->GetQ())<1.e-13 && fSkip==2) continue;
1987       ci=fClusterIndexCs[i];
1988       if (!test) fI=i+1;
1989       return fClustersCs[i];
1990     }
1991   }
1992   return 0;
1993 }
1994 //------------------------------------------------------------------------
1995 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1996 const {
1997   //--------------------------------------------------------------------
1998   // This function returns the layer thickness at this point (units X0)
1999   //--------------------------------------------------------------------
2000   Double_t d=0.0085;
2001   x0=AliITSRecoParam::GetX0Air();
2002   if (43<fR&&fR<45) { //SSD2
2003      Double_t dd=0.0034;
2004      d=dd;
2005      if (TMath::Abs(y-0.00)>3.40) d+=dd;
2006      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2007      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2008      for (Int_t i=0; i<12; i++) {
2009        if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
2010           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2011           d+=0.0034; 
2012           break;
2013        }
2014        if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
2015           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2016           d+=0.0034; 
2017           break;
2018        }         
2019        if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2020        if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2021      }
2022   } else 
2023   if (37<fR&&fR<41) { //SSD1
2024      Double_t dd=0.0034;
2025      d=dd;
2026      if (TMath::Abs(y-0.00)>3.40) d+=dd;
2027      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2028      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2029      for (Int_t i=0; i<11; i++) {
2030        if (TMath::Abs(z-3.9*i)<0.15) {
2031           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2032           d+=dd; 
2033           break;
2034        }
2035        if (TMath::Abs(z+3.9*i)<0.15) {
2036           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2037           d+=dd; 
2038           break;
2039        }         
2040        if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2041        if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}         
2042      }
2043   } else
2044   if (13<fR&&fR<26) { //SDD
2045      Double_t dd=0.0033;
2046      d=dd;
2047      if (TMath::Abs(y-0.00)>3.30) d+=dd;
2048
2049      if (TMath::Abs(y-1.80)<0.55) {
2050         d+=0.016;
2051         for (Int_t j=0; j<20; j++) {
2052           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2053           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2054         } 
2055      }
2056      if (TMath::Abs(y+1.80)<0.55) {
2057         d+=0.016;
2058         for (Int_t j=0; j<20; j++) {
2059           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2060           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2061         } 
2062      }
2063
2064      for (Int_t i=0; i<4; i++) {
2065        if (TMath::Abs(z-7.3*i)<0.60) {
2066           d+=dd;
2067           if (TMath::Abs(y-0.00)>3.30) d+=dd; 
2068           break;
2069        }
2070        if (TMath::Abs(z+7.3*i)<0.60) {
2071           d+=dd; 
2072           if (TMath::Abs(y-0.00)>3.30) d+=dd; 
2073           break;
2074        }
2075      }
2076   } else
2077   if (6<fR&&fR<8) {   //SPD2
2078      Double_t dd=0.0063; x0=21.5;
2079      d=dd;
2080      if (TMath::Abs(y-3.08)>0.5) d+=dd;
2081      if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2082   } else
2083   if (3<fR&&fR<5) {   //SPD1
2084      Double_t dd=0.0063; x0=21.5;
2085      d=dd;
2086      if (TMath::Abs(y+0.21)>0.6) d+=dd;
2087      if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2088   }
2089
2090   return d;
2091 }
2092 //------------------------------------------------------------------------
2093 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2094 fR(det.fR),
2095 fRmisal(det.fRmisal),
2096 fPhi(det.fPhi),
2097 fSinPhi(det.fSinPhi),
2098 fCosPhi(det.fCosPhi),
2099 fYmin(det.fYmin),
2100 fYmax(det.fYmax),
2101 fZmin(det.fZmin),
2102 fZmax(det.fZmax),
2103 fIsBad(det.fIsBad),
2104 fNChips(det.fNChips),
2105 fChipIsBad(det.fChipIsBad)
2106 {
2107   //Copy constructor
2108 }
2109 //------------------------------------------------------------------------
2110 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2111                                                const AliITSDetTypeRec *detTypeRec)
2112 {
2113   //--------------------------------------------------------------------
2114   // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2115   //--------------------------------------------------------------------
2116
2117   // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2118   // while in the tracker they start from 0 for each layer
2119   for(Int_t il=0; il<ilayer; il++) 
2120     idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2121
2122   Int_t detType;
2123   if (ilayer==0 || ilayer==1) {        // ----------  SPD
2124     detType = 0;
2125   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
2126     detType = 1;
2127   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
2128     detType = 2;
2129   } else {
2130     printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2131     return;
2132   }
2133
2134   // Get calibration from AliITSDetTypeRec
2135   AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2136   calib->SetModuleIndex(idet);
2137   AliITSCalibration *calibSPDdead = 0;
2138   if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2139   if (calib->IsBad() ||
2140       (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2141     {
2142       SetBad();
2143       //      printf("lay %d bad %d\n",ilayer,idet);
2144     }
2145
2146   // Get segmentation from AliITSDetTypeRec
2147   AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2148
2149   // Read info about bad chips
2150   fNChips = segm->GetMaximumChipIndex()+1;
2151   //printf("ilayer %d  detType %d idet %d fNChips %d %d  GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2152   if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2153   fChipIsBad = new Bool_t[fNChips];
2154   for (Int_t iCh=0;iCh<fNChips;iCh++) {
2155     fChipIsBad[iCh] = calib->IsChipBad(iCh);
2156     if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2157     //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2158   }
2159
2160   return;
2161 }
2162 //------------------------------------------------------------------------
2163 Double_t AliITStrackerMI::GetEffectiveThickness()
2164 {
2165   //--------------------------------------------------------------------
2166   // Returns the thickness between the current layer and the vertex (units X0)
2167   //--------------------------------------------------------------------
2168
2169   if(fUseTGeo!=0) {
2170     if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2171     if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2172     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2173   }
2174
2175   // beam pipe
2176   Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2177   Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2178
2179   // layers
2180   Double_t x0=0;
2181   Double_t xn=fgLayers[fI].GetR();
2182   for (Int_t i=0; i<fI; i++) {
2183     Double_t xi=fgLayers[i].GetR();
2184     Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2185     d+=dLayer*xi*xi;
2186   }
2187
2188   // shields
2189   if (fI>1) {
2190     Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2191     d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2192   }
2193   if (fI>3) {
2194     Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2195     d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2196   }
2197   return d/(xn*xn);
2198 }
2199 //------------------------------------------------------------------------
2200 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2201   //-------------------------------------------------------------------
2202   // This function returns number of clusters within the "window" 
2203   //--------------------------------------------------------------------
2204   Int_t ncl=0;
2205   for (Int_t i=fI; i<fN; i++) {
2206     const AliITSRecPoint *c=fClusters[i];
2207     if (c->GetZ() > fZmax) break;
2208     if (c->IsUsed()) continue;
2209     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());    
2210     Double_t y=fR*det.GetPhi() + c->GetY();
2211
2212     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2213     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2214
2215     if (y<fYmin) continue;
2216     if (y>fYmax) continue;
2217     ncl++;
2218   }
2219   return ncl;
2220 }
2221 //------------------------------------------------------------------------
2222 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2223                                 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff) 
2224 {
2225   //--------------------------------------------------------------------
2226   // This function refits the track "track" at the position "x" using
2227   // the clusters from "clusters"
2228   // If "extra"==kTRUE, 
2229   //    the clusters from overlapped modules get attached to "track" 
2230   // If "planeff"==kTRUE,
2231   //    special approach for plane efficiency evaluation is applyed
2232   //--------------------------------------------------------------------
2233
2234   Int_t index[AliITSgeomTGeo::kNLayers];
2235   Int_t k;
2236   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2237   Int_t nc=clusters->GetNumberOfClusters();
2238   for (k=0; k<nc; k++) { 
2239     Int_t idx=clusters->GetClusterIndex(k);
2240     Int_t ilayer=(idx&0xf0000000)>>28;
2241     index[ilayer]=idx; 
2242   }
2243
2244   return RefitAt(xx,track,index,extra,planeeff); // call the method below
2245 }
2246 //------------------------------------------------------------------------
2247 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2248                                 const Int_t *clusters,Bool_t extra, Bool_t planeeff) 
2249 {
2250   //--------------------------------------------------------------------
2251   // This function refits the track "track" at the position "x" using
2252   // the clusters from array
2253   // If "extra"==kTRUE, 
2254   //    the clusters from overlapped modules get attached to "track" 
2255   // If "planeff"==kTRUE,
2256   //    special approach for plane efficiency evaluation is applyed
2257   //--------------------------------------------------------------------
2258   Int_t index[AliITSgeomTGeo::kNLayers];
2259   Int_t k;
2260   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2261   //
2262   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) { 
2263     index[k]=clusters[k]; 
2264   }
2265
2266   // special for cosmics and TPC prolonged tracks: 
2267   // propagate to the innermost of:
2268   // - innermost layer crossed by the track
2269   // - innermost layer where a cluster was associated to the track
2270   static AliITSRecoParam *repa = NULL;
2271   if(!repa){
2272     repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
2273     if(!repa){
2274       repa = AliITSRecoParam::GetHighFluxParam();
2275       AliWarning("Using default AliITSRecoParam class");
2276     }
2277   }
2278   Int_t evsp=repa->GetEventSpecie();
2279   ULong_t trStatus=0;
2280   if(track->GetESDtrack()) trStatus=track->GetStatus();
2281   Int_t innermostlayer=0;
2282   if((evsp&AliRecoParam::kCosmic) || (trStatus&AliESDtrack::kTPCin))  {
2283     innermostlayer=5;
2284     Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2285     for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2286       if( (drphi < (fgLayers[innermostlayer].GetR()+1.)) ||
2287           index[innermostlayer] >= 0 ) break;
2288     }
2289
2290     AliDebug(2,Form(" drphi  %f  innermost %d",drphi,innermostlayer));
2291   }
2292
2293   Int_t modstatus=1; // found
2294   Float_t xloc,zloc;
2295   Int_t from, to, step;
2296   if (xx > track->GetX()) {
2297       from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2298       step=+1;
2299   } else {
2300       from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2301       step=-1;
2302   }
2303   TString dir = (step>0 ? "outward" : "inward");
2304
2305   for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2306      AliITSlayer &layer=fgLayers[ilayer];
2307      Double_t r=layer.GetR();
2308
2309      if (step<0 && xx>r) break;
2310
2311      // material between SSD and SDD, SDD and SPD
2312      Double_t hI=ilayer-0.5*step; 
2313      if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2314        if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2315      if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2316        if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2317
2318
2319      Double_t oldGlobXYZ[3];
2320      if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2321
2322      // continue if we are already beyond this layer
2323      Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2324      if(step>0 && oldGlobR > r) continue; // going outward
2325      if(step<0 && oldGlobR < r) continue; // going inward
2326
2327      Double_t phi,z;
2328      if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2329
2330      Int_t idet=layer.FindDetectorIndex(phi,z);
2331
2332      // check if we allow a prolongation without point for large-eta tracks
2333      Int_t skip = CheckSkipLayer(track,ilayer,idet);
2334      if (skip==2) {
2335        modstatus = 4; // out in z
2336        if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2337          track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2338        }
2339        // cross layer
2340        // apply correction for material of the current layer
2341        // add time if going outward
2342        CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2343        continue;
2344      }
2345
2346      if (idet<0) return kFALSE;
2347
2348
2349      const AliITSdetector &det=layer.GetDetector(idet);
2350      // only for ITS-SA tracks refit
2351      if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2352      // 
2353      if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2354
2355      track->SetDetectorIndex(idet);
2356      if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2357
2358      Double_t dz,zmin,zmax,dy,ymin,ymax;
2359
2360      const AliITSRecPoint *clAcc=0;
2361      Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2362
2363      Int_t idx=index[ilayer];
2364      if (idx>=0) { // cluster in this layer
2365        modstatus = 6; // no refit
2366        const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx); 
2367        if (cl) {
2368          if (idet != cl->GetDetectorIndex()) {
2369            idet=cl->GetDetectorIndex();
2370            const AliITSdetector &detc=layer.GetDetector(idet);
2371            if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2372            track->SetDetectorIndex(idet);
2373            if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2374          }
2375          Int_t cllayer = (idx & 0xf0000000) >> 28;;
2376          Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2377          if (chi2<maxchi2) { 
2378            clAcc=cl; 
2379            maxchi2=chi2; 
2380            modstatus = 1; // found
2381          } else {
2382             return kFALSE; //
2383          }
2384        }
2385      } else { // no cluster in this layer
2386        if (skip==1) {
2387          modstatus = 3; // skipped
2388          // Plane Eff determination:
2389          if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2390            if (IsOKForPlaneEff(track,clusters,ilayer))  // only adequate track for plane eff. evaluation
2391               UseTrackForPlaneEff(track,ilayer);
2392          }
2393        } else {
2394          modstatus = 5; // no cls in road
2395          // check dead
2396          if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2397          dz = 0.5*(zmax-zmin);
2398          dy = 0.5*(ymax-ymin);
2399          Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2400          if (dead==1) modstatus = 7; // holes in z in SPD
2401          if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2402        }
2403      }
2404      
2405      if (clAcc) {
2406        if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2407        track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2408      }
2409      track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2410
2411
2412      if (extra) { // search for extra clusters in overlapped modules
2413        AliITStrackV2 tmp(*track);
2414        if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2415        layer.SelectClusters(zmin,zmax,ymin,ymax);
2416        
2417        const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2418        Int_t idetExtra=-1;  
2419        maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2420        Double_t tolerance=0.1;
2421        while ((clExtra=layer.GetNextCluster(ci))!=0) {
2422          // only clusters in another module! (overlaps)
2423          idetExtra = clExtra->GetDetectorIndex();
2424          if (idet == idetExtra) continue;
2425          
2426          const AliITSdetector &detx=layer.GetDetector(idetExtra);
2427          
2428          if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2429          if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2430          if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2431          if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2432          
2433          Double_t chi2=tmp.GetPredictedChi2(clExtra);
2434          if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2435        }
2436        if (cci>=0) {
2437          track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2438          track->SetExtraModule(ilayer,idetExtra);
2439        }
2440      } // end search for extra clusters in overlapped modules
2441      
2442      // Correct for material of the current layer
2443      // cross material
2444      // add time if going outward
2445      if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2446      track->SetCheckInvariant(kTRUE);
2447   } // end loop on layers
2448
2449   if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2450
2451   return kTRUE;
2452 }
2453 //------------------------------------------------------------------------
2454 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2455 {
2456   //
2457   // calculate normalized chi2
2458   //  return NormalizedChi2(track,0);
2459   Float_t chi2 = 0;
2460   Float_t sum=0;
2461   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2462   //  track->fdEdxMismatch=0;
2463   Float_t dedxmismatch =0;
2464   Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack); 
2465   if (mode<100){
2466     for (Int_t i = 0;i<6;i++){
2467       if (track->GetClIndex(i)>=0){
2468         Float_t cerry, cerrz;
2469         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2470         else 
2471           { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2472         cerry*=cerry;
2473         cerrz*=cerrz;   
2474         Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2475         if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2476           Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2477           if (ratio<0.5) {
2478             cchi2+=(0.5-ratio)*10.;
2479             //track->fdEdxMismatch+=(0.5-ratio)*10.;
2480             dedxmismatch+=(0.5-ratio)*10.;          
2481           }
2482         }
2483         if (i<2 ||i>3){
2484           AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));  
2485           Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2486           if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.); 
2487           if (i<2) chi2+=2*cl->GetDeltaProbability();
2488         }
2489         chi2+=cchi2;
2490         sum++;
2491       }
2492     }
2493     if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2494       track->SetdEdxMismatch(dedxmismatch);
2495     }
2496   }
2497   else{
2498     for (Int_t i = 0;i<4;i++){
2499       if (track->GetClIndex(i)>=0){
2500         Float_t cerry, cerrz;
2501         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2502         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2503         cerry*=cerry;
2504         cerrz*=cerrz;
2505         chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2506         chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);      
2507         sum++;
2508       }
2509     }
2510     for (Int_t i = 4;i<6;i++){
2511       if (track->GetClIndex(i)>=0){     
2512         Float_t cerry, cerrz;
2513         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2514         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2515         cerry*=cerry;
2516         cerrz*=cerrz;   
2517         Float_t cerryb, cerrzb;
2518         if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2519         else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2520         cerryb*=cerryb;
2521         cerrzb*=cerrzb;
2522         chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2523         chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);      
2524         sum++;
2525       }
2526     }
2527   }
2528   if (track->GetESDtrack()->GetTPCsignal()>85){
2529     Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2530     if (ratio<0.5) {
2531       chi2+=(0.5-ratio)*5.;      
2532     }
2533     if (ratio>2){
2534       chi2+=(ratio-2.0)*3; 
2535     }
2536   }
2537   //
2538   Double_t match = TMath::Sqrt(track->GetChi22());
2539   if (track->GetConstrain())  match/=track->GetNumberOfClusters();
2540   if (!track->GetConstrain()) { 
2541     if (track->GetNumberOfClusters()>2) {
2542       match/=track->GetNumberOfClusters()-2.;
2543     } else {
2544       match=0;
2545     }
2546   }
2547   if (match<0) match=0;
2548
2549   // penalty factor for missing points (NDeadZone>0), but no penalty
2550   // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2551   // or there is a dead from OCDB)
2552   Float_t deadzonefactor = 0.; 
2553   if (track->GetNDeadZone()>0.) {    
2554     Int_t sumDeadZoneProbability=0; 
2555     for(Int_t ilay=0;ilay<6;ilay++) {
2556       if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2557     }
2558     Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2559     if(nDeadZoneWithProbNot1>0) {
2560       Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
2561       AliDebug(2,Form("nDeadZone %f sumDZProbability %d nDZWithProbNot1 %d deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
2562       deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2563       Float_t one = 1.;
2564       deadZoneProbability = TMath::Min(deadZoneProbability,one);
2565       deadzonefactor = 3.*(1.1-deadZoneProbability);  
2566     }
2567   }  
2568
2569   Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2570     (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2571                                 1./(1.+track->GetNSkipped()));     
2572   AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped %f\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
2573   AliDebug(2,Form("NormChi2 %f  cls %d\n",normchi2,track->GetNumberOfClusters()));
2574   return normchi2;
2575 }
2576 //------------------------------------------------------------------------
2577 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2578 {
2579   //
2580   // return matching chi2 between two tracks
2581   Double_t largeChi2=1000.;
2582
2583   AliITStrackMI track3(*track2);
2584   if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2585   TMatrixD vec(5,1);
2586   vec(0,0)=track1->GetY()   - track3.GetY();
2587   vec(1,0)=track1->GetZ()   - track3.GetZ();
2588   vec(2,0)=track1->GetSnp() - track3.GetSnp();
2589   vec(3,0)=track1->GetTgl() - track3.GetTgl();
2590   vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2591   //
2592   TMatrixD cov(5,5);
2593   cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2594   cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2595   cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2596   cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2597   cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2598   
2599   cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2600   cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2601   cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2602   cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2603   //
2604   cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2605   cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2606   cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2607   //
2608   cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2609   cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2610   //
2611   cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2612   
2613   cov.Invert();
2614   TMatrixD vec2(cov,TMatrixD::kMult,vec);
2615   TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2616   return chi2(0,0);
2617 }
2618 //------------------------------------------------------------------------
2619 Double_t  AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2620 {
2621   //
2622   //  return probability that given point (characterized by z position and error) 
2623   //  is in SPD dead zone
2624   //     This method assumes that fSPDdetzcentre is ordered from -z to +z
2625   //
2626   Double_t probability = 0.;
2627   Double_t nearestz = 0.,distz=0.;
2628   Int_t    nearestzone = -1;
2629   Double_t mindistz = 1000.;
2630
2631   // find closest dead zone
2632   for (Int_t i=0; i<3; i++) {
2633     distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2634     if (distz<mindistz) {
2635       nearestzone=i;
2636       nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2637       mindistz=distz;
2638     }
2639   }
2640
2641   // too far from dead zone
2642   if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2643
2644
2645   Double_t zmin, zmax;   
2646   if (nearestzone==0) { // dead zone at z = -7
2647     zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2648     zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2649   } else if (nearestzone==1) { // dead zone at z = 0
2650     zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2651     zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2652   } else if (nearestzone==2) { // dead zone at z = +7
2653     zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2654     zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2655   } else {
2656     zmin = 0.;
2657     zmax = 0.;
2658   }
2659   // probability that the true z is in the range [zmin,zmax] (i.e. inside 
2660   // dead zone)
2661   probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) - 
2662                       AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2663   AliDebug(2,Form("zpos %f +- %f nearestzone %d  zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
2664   return probability;
2665 }
2666 //------------------------------------------------------------------------
2667 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2668 {
2669   //
2670   // calculate normalized chi2
2671   Float_t chi2[6];
2672   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2673   Float_t ncl = 0;
2674   for (Int_t i = 0;i<6;i++){
2675     if (TMath::Abs(track->GetDy(i))>0){      
2676       chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2677       chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2678       ncl++;
2679     }
2680     else{chi2[i]=10000;}
2681   }
2682   Int_t index[6];
2683   TMath::Sort(6,chi2,index,kFALSE);
2684   Float_t max = float(ncl)*fac-1.;
2685   Float_t sumchi=0, sumweight=0; 
2686   for (Int_t i=0;i<max+1;i++){
2687     Float_t weight = (i<max)?1.:(max+1.-i);
2688     sumchi+=weight*chi2[index[i]];
2689     sumweight+=weight;
2690   }
2691   Double_t normchi2 = sumchi/sumweight;
2692   return normchi2;
2693 }
2694 //------------------------------------------------------------------------
2695 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2696 {
2697   //
2698   // calculate normalized chi2
2699   //  if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2700   Int_t npoints = 0;
2701   Double_t res =0;
2702   for (Int_t i=0;i<6;i++){
2703     if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2704     Double_t sy1 = forwardtrack->GetSigmaY(i);
2705     Double_t sz1 = forwardtrack->GetSigmaZ(i);
2706     Double_t sy2 = backtrack->GetSigmaY(i);
2707     Double_t sz2 = backtrack->GetSigmaZ(i);
2708     if (i<2){ sy2=1000.;sz2=1000;}
2709     //    
2710     Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2711     Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2712     // 
2713     Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2714     Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2715     //
2716     res+= nz0*nz0+ny0*ny0;
2717     npoints++;
2718   }
2719   if (npoints>1) return 
2720                    TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2721                    //2*forwardtrack->fNUsed+
2722                    res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2723                                   1./(1.+forwardtrack->GetNSkipped()));
2724   return 1000;
2725 }
2726 //------------------------------------------------------------------------
2727 Float_t  *AliITStrackerMI::GetWeight(Int_t index) {
2728   //--------------------------------------------------------------------
2729   //       Return pointer to a given cluster
2730   //--------------------------------------------------------------------
2731   Int_t l=(index & 0xf0000000) >> 28;
2732   Int_t c=(index & 0x0fffffff) >> 00;
2733   return fgLayers[l].GetWeight(c);
2734 }
2735 //------------------------------------------------------------------------
2736 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2737 {
2738   //---------------------------------------------
2739   // register track to the list
2740   //
2741   if (track->GetESDtrack()->GetKinkIndex(0)!=0) return;  //don't register kink tracks
2742   //
2743   //
2744   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2745     Int_t index = track->GetClusterIndex(icluster);
2746     Int_t l=(index & 0xf0000000) >> 28;
2747     Int_t c=(index & 0x0fffffff) >> 00;
2748     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2749     for (Int_t itrack=0;itrack<4;itrack++){
2750       if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2751         fgLayers[l].SetClusterTracks(itrack,c,id);
2752         break;
2753       }
2754     }
2755   }
2756 }
2757 //------------------------------------------------------------------------
2758 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2759 {
2760   //---------------------------------------------
2761   // unregister track from the list
2762   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2763     Int_t index = track->GetClusterIndex(icluster);
2764     Int_t l=(index & 0xf0000000) >> 28;
2765     Int_t c=(index & 0x0fffffff) >> 00;
2766     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2767     for (Int_t itrack=0;itrack<4;itrack++){
2768       if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2769         fgLayers[l].SetClusterTracks(itrack,c,-1);
2770       }
2771     }
2772   }
2773 }
2774 //------------------------------------------------------------------------
2775 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2776 {
2777   //-------------------------------------------------------------
2778   //get number of shared clusters
2779   //-------------------------------------------------------------
2780   Float_t shared=0;
2781   for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2782   // mean  number of clusters
2783   Float_t *ny = GetNy(id), *nz = GetNz(id); 
2784
2785  
2786   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2787     Int_t index = track->GetClusterIndex(icluster);
2788     Int_t l=(index & 0xf0000000) >> 28;
2789     Int_t c=(index & 0x0fffffff) >> 00;
2790     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2791     // if (ny[l]<1.e-13){
2792     //   printf("problem\n");
2793     // }
2794     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2795     Float_t weight=1;
2796     //
2797     Float_t deltan = 0;
2798     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2799     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2800       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2801     if (l<2 || l>3){      
2802       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2803     }
2804     else{
2805       deltan = (cl->GetNz()-nz[l]);
2806     }
2807     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2808     weight = 2./TMath::Max(3.+deltan,2.);
2809     //
2810     for (Int_t itrack=0;itrack<4;itrack++){
2811       if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2812         list[l]=index;
2813         clist[l] = (AliITSRecPoint*)GetCluster(index);
2814         shared+=weight; 
2815         break;
2816       }
2817     }
2818   }
2819   track->SetNUsed(shared);
2820   return shared;
2821 }
2822 //------------------------------------------------------------------------
2823 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2824 {
2825   //
2826   // find first shared track 
2827   //
2828   // mean  number of clusters
2829   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID); 
2830   //
2831   for (Int_t i=0;i<6;i++) overlist[i]=-1;
2832   Int_t sharedtrack=100000;
2833   Int_t tracks[24],trackindex=0;
2834   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2835   //
2836   for (Int_t icluster=0;icluster<6;icluster++){
2837     if (clusterlist[icluster]<0) continue;
2838     Int_t index = clusterlist[icluster];
2839     Int_t l=(index & 0xf0000000) >> 28;
2840     Int_t c=(index & 0x0fffffff) >> 00;
2841     // if (ny[l]<1.e-13){
2842     //   printf("problem\n");
2843     // }
2844     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2845     //if (l>3) continue;
2846     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2847     //
2848     Float_t deltan = 0;
2849     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2850     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2851       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2852     if (l<2 || l>3){      
2853       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2854     }
2855     else{
2856       deltan = (cl->GetNz()-nz[l]);
2857     }
2858     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2859     //
2860     for (Int_t itrack=3;itrack>=0;itrack--){
2861       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2862       if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2863        tracks[trackindex]  = fgLayers[l].GetClusterTracks(itrack,c);
2864        trackindex++;
2865       }
2866     }
2867   }
2868   if (trackindex==0) return -1;
2869   if (trackindex==1){    
2870     sharedtrack = tracks[0];
2871   }else{
2872     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2873     else{
2874       //
2875       Int_t tracks2[24], cluster[24];
2876       for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2877       Int_t index =0;
2878       //
2879       for (Int_t i=0;i<trackindex;i++){
2880         if (tracks[i]<0) continue;
2881         tracks2[index] = tracks[i];
2882         cluster[index]++;       
2883         for (Int_t j=i+1;j<trackindex;j++){
2884           if (tracks[j]<0) continue;
2885           if (tracks[j]==tracks[i]){
2886             cluster[index]++;
2887             tracks[j]=-1;
2888           }
2889         }
2890         index++;
2891       }
2892       Int_t max=0;
2893       for (Int_t i=0;i<index;i++){
2894         if (cluster[index]>max) {
2895           sharedtrack=tracks2[index];
2896           max=cluster[index];
2897         }
2898       }
2899     }
2900   }
2901   
2902   if (sharedtrack>=100000) return -1;
2903   //
2904   // make list of overlaps
2905   shared =0;
2906   for (Int_t icluster=0;icluster<6;icluster++){
2907     if (clusterlist[icluster]<0) continue;
2908     Int_t index = clusterlist[icluster];
2909     Int_t l=(index & 0xf0000000) >> 28;
2910     Int_t c=(index & 0x0fffffff) >> 00;
2911     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2912     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2913     if (l==0 || l==1){
2914       if (cl->GetNy()>2) continue;
2915       if (cl->GetNz()>2) continue;
2916     }
2917     if (l==4 || l==5){
2918       if (cl->GetNy()>3) continue;
2919       if (cl->GetNz()>3) continue;
2920     }
2921     //
2922     for (Int_t itrack=3;itrack>=0;itrack--){
2923       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2924       if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2925         overlist[l]=index;
2926         shared++;      
2927       }
2928     }
2929   }
2930   return sharedtrack;
2931 }
2932 //------------------------------------------------------------------------
2933 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2934   //
2935   // try to find track hypothesys without conflicts
2936   // with minimal chi2;
2937   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2938   Int_t entries1 = arr1->GetEntriesFast();
2939   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2940   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2941   Int_t entries2 = arr2->GetEntriesFast();
2942   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2943   //
2944   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2945   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2946   if (track10->Pt()>0.5+track20->Pt()) return track10;
2947
2948   for (Int_t itrack=0;itrack<entries1;itrack++){
2949     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2950     UnRegisterClusterTracks(track,trackID1);
2951   }
2952   //
2953   for (Int_t itrack=0;itrack<entries2;itrack++){
2954     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2955     UnRegisterClusterTracks(track,trackID2);
2956   }
2957   Int_t index1=0;
2958   Int_t index2=0;
2959   Float_t maxconflicts=6;
2960   Double_t maxchi2 =1000.;
2961   //
2962   // get weight of hypothesys - using best hypothesy
2963   Double_t w1,w2;
2964  
2965   Int_t list1[6],list2[6];
2966   AliITSRecPoint *clist1[6], *clist2[6] ;
2967   RegisterClusterTracks(track10,trackID1);
2968   RegisterClusterTracks(track20,trackID2);
2969   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2970   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2971   UnRegisterClusterTracks(track10,trackID1);
2972   UnRegisterClusterTracks(track20,trackID2);
2973   //
2974   // normalized chi2
2975   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2976   Float_t nerry[6],nerrz[6];
2977   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2978   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2979   for (Int_t i=0;i<6;i++){
2980      if ( (erry1[i]>0) && (erry2[i]>0)) {
2981        nerry[i] = TMath::Min(erry1[i],erry2[i]);
2982        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2983      }else{
2984        nerry[i] = TMath::Max(erry1[i],erry2[i]);
2985        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2986      }
2987      if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2988        chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2989        chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2990        ncl1++;
2991      }
2992      if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2993        chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2994        chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2995        ncl2++;
2996      }
2997   }
2998   chi21/=ncl1;
2999   chi22/=ncl2;
3000   //
3001   // 
3002   Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
3003   Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
3004   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
3005   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
3006   //
3007   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
3008         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
3009         +1.*track10->Pt()/(track10->Pt()+track20->Pt())
3010         );
3011   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
3012         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
3013         +1.*track20->Pt()/(track10->Pt()+track20->Pt())
3014         );
3015
3016   Double_t sumw = w1+w2;
3017   w1/=sumw;
3018   w2/=sumw;
3019   if (w1<w2*0.5) {
3020     w1 = (d2+0.5)/(d1+d2+1.);
3021     w2 = (d1+0.5)/(d1+d2+1.);
3022   }
3023   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
3024   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
3025   //
3026   // get pair of "best" hypothesys
3027   //  
3028   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1); 
3029   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2); 
3030
3031   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3032     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3033     //if (track1->fFakeRatio>0) continue;
3034     RegisterClusterTracks(track1,trackID1);
3035     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3036       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3037
3038       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3039       //if (track2->fFakeRatio>0) continue;
3040       Float_t nskipped=0;            
3041       RegisterClusterTracks(track2,trackID2);
3042       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3043       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3044       UnRegisterClusterTracks(track2,trackID2);
3045       //
3046       if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3047       if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3048       if (nskipped>0.5) continue;
3049       //
3050       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3051       if (conflict1+1<cconflict1) continue;
3052       if (conflict2+1<cconflict2) continue;
3053       Float_t conflict=0;
3054       Float_t sumchi2=0;
3055       Float_t sum=0;
3056       for (Int_t i=0;i<6;i++){
3057         //
3058         Float_t c1 =0.; // conflict coeficients
3059         Float_t c2 =0.; 
3060         if (clist1[i]&&clist2[i]){
3061           Float_t deltan = 0;
3062           if (i<2 || i>3){      
3063             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3064           }
3065           else{
3066             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3067           }
3068           c1  = 2./TMath::Max(3.+deltan,2.);      
3069           c2  = 2./TMath::Max(3.+deltan,2.);      
3070         }
3071         else{
3072           if (clist1[i]){
3073             Float_t deltan = 0;
3074             if (i<2 || i>3){      
3075               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3076             }
3077             else{
3078               deltan = (clist1[i]->GetNz()-nz1[i]);
3079             }
3080             c1  = 2./TMath::Max(3.+deltan,2.);    
3081             c2  = 0;
3082           }
3083
3084           if (clist2[i]){
3085             Float_t deltan = 0;
3086             if (i<2 || i>3){      
3087               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3088             }
3089             else{
3090               deltan = (clist2[i]->GetNz()-nz2[i]);
3091             }
3092             c2  = 2./TMath::Max(3.+deltan,2.);    
3093             c1  = 0;
3094           }       
3095         }
3096         //
3097         chi21=0;chi22=0;
3098         if (TMath::Abs(track1->GetDy(i))>0.) {
3099           chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3100             (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3101           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3102           //  (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3103         }else{
3104           if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3105         }
3106         //
3107         if (TMath::Abs(track2->GetDy(i))>0.) {
3108           chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3109             (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3110           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3111           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3112         }
3113         else{
3114           if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3115         }
3116         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3117         if (chi21>0) sum+=w1;
3118         if (chi22>0) sum+=w2;
3119         conflict+=(c1+c2);
3120       }
3121       Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3122       if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3123       Double_t normchi2 = 2*conflict+sumchi2/sum;
3124       if ( normchi2 <maxchi2 ){   
3125         index1 = itrack1;
3126         index2 = itrack2;
3127         maxconflicts = conflict;
3128         maxchi2 = normchi2;
3129       }      
3130     }
3131     UnRegisterClusterTracks(track1,trackID1);
3132   }
3133   //
3134   //  if (maxconflicts<4 && maxchi2<th0){   
3135   if (maxchi2<th0*2.){   
3136     Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3137     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3138     track1->SetChi2MIP(5,maxconflicts);
3139     track1->SetChi2MIP(6,maxchi2);
3140     track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3141     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
3142     track1->SetChi2MIP(8,index1);
3143     fBestTrackIndex[trackID1] =index1;
3144     UpdateESDtrack(track1, AliESDtrack::kITSin);
3145   }  
3146   else if (track10->GetChi2MIP(0)<th1){
3147     track10->SetChi2MIP(5,maxconflicts);
3148     track10->SetChi2MIP(6,maxchi2);    
3149     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
3150     UpdateESDtrack(track10,AliESDtrack::kITSin);
3151   }   
3152   
3153   for (Int_t itrack=0;itrack<entries1;itrack++){
3154     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3155     UnRegisterClusterTracks(track,trackID1);
3156   }
3157   //
3158   for (Int_t itrack=0;itrack<entries2;itrack++){
3159     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3160     UnRegisterClusterTracks(track,trackID2);
3161   }
3162
3163   if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3164       &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){ 
3165     //  if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3166   //    &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){ 
3167     RegisterClusterTracks(track10,trackID1);
3168   }
3169   if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3170       &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){ 
3171     //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3172     //  &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){ 
3173     RegisterClusterTracks(track20,trackID2);  
3174   }
3175   return track10; 
3176  
3177 }
3178 //------------------------------------------------------------------------
3179 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3180   //--------------------------------------------------------------------
3181   // This function marks clusters assigned to the track
3182   //--------------------------------------------------------------------
3183   AliTracker::UseClusters(t,from);
3184
3185   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3186   //if (c->GetQ()>2) c->Use();
3187   if (c->GetSigmaZ2()>0.1) c->Use();
3188   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3189   //if (c->GetQ()>2) c->Use();
3190   if (c->GetSigmaZ2()>0.1) c->Use();
3191
3192 }
3193 //------------------------------------------------------------------------
3194 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3195 {
3196   //------------------------------------------------------------------
3197   // add track to the list of hypothesys
3198   //------------------------------------------------------------------
3199
3200   if (esdindex>=fTrackHypothesys.GetEntriesFast()) 
3201     fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3202   //
3203   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3204   if (!array) {
3205     array = new TObjArray(10);
3206     fTrackHypothesys.AddAt(array,esdindex);
3207   }
3208   array->AddLast(track);
3209 }
3210 //------------------------------------------------------------------------
3211 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3212 {
3213   //-------------------------------------------------------------------
3214   // compress array of track hypothesys
3215   // keep only maxsize best hypothesys
3216   //-------------------------------------------------------------------
3217   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3218   if (! (fTrackHypothesys.At(esdindex)) ) return;
3219   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3220   Int_t entries = array->GetEntriesFast();
3221   //
3222   //- find preliminary besttrack as a reference
3223   Float_t minchi2=10000;
3224   Int_t maxn=0;
3225   AliITStrackMI * besttrack=0;
3226   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3227     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3228     if (!track) continue;
3229     Float_t chi2 = NormalizedChi2(track,0);
3230     //
3231     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3232     track->SetLabel(tpcLabel);
3233     CookdEdx(track);
3234     track->SetFakeRatio(1.);
3235     CookLabel(track,0.); //For comparison only
3236     //
3237     //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3238     if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3239       if (track->GetNumberOfClusters()<maxn) continue;
3240       maxn = track->GetNumberOfClusters();
3241       if (chi2<minchi2){
3242         minchi2=chi2;
3243         besttrack=track;
3244       }
3245     }
3246     else{
3247       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3248         delete array->RemoveAt(itrack);
3249       }  
3250     }
3251   }
3252   if (!besttrack) return;
3253   //
3254   //
3255   //take errors of best track as a reference
3256   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3257   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3258   for (Int_t j=0;j<6;j++) {
3259     if (besttrack->GetClIndex(j)>=0){
3260       erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3261       errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3262       ny[j]   = besttrack->GetNy(j);
3263       nz[j]   = besttrack->GetNz(j);
3264     }
3265   }
3266   //
3267   // calculate normalized chi2
3268   //
3269   Float_t * chi2        = new Float_t[entries];
3270   Int_t * index         = new Int_t[entries];  
3271   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3272   for (Int_t itrack=0;itrack<entries;itrack++){
3273     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3274     if (track){
3275       AliDebug(2,Form("track %d  ncls %d\n",itrack,track->GetNumberOfClusters()));
3276       track->SetChi2MIP(0,GetNormalizedChi2(track, mode));            
3277       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) 
3278         chi2[itrack] = track->GetChi2MIP(0);
3279       else{
3280         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3281           delete array->RemoveAt(itrack);            
3282         }
3283       }
3284     }
3285   }
3286   //
3287   TMath::Sort(entries,chi2,index,kFALSE);
3288   besttrack = (AliITStrackMI*)array->At(index[0]);
3289   if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3290   if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3291     for (Int_t j=0;j<6;j++){
3292       if (besttrack->GetClIndex(j)>=0){
3293         erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3294         errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3295         ny[j]   = besttrack->GetNy(j);
3296         nz[j]   = besttrack->GetNz(j);
3297       }
3298     }
3299   }
3300   //
3301   // calculate one more time with updated normalized errors
3302   for (Int_t i=0;i<entries;i++) chi2[i] =10000;  
3303   for (Int_t itrack=0;itrack<entries;itrack++){
3304     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3305     if (track){      
3306       track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3307       AliDebug(2,Form("track %d  ncls %d\n",itrack,track->GetNumberOfClusters()));            
3308       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) 
3309         chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone()); 
3310       else
3311         {
3312           if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3313             delete array->RemoveAt(itrack);     
3314           }
3315         }
3316     }   
3317   }
3318   entries = array->GetEntriesFast();  
3319   //
3320   //
3321   if (entries>0){
3322     TObjArray * newarray = new TObjArray();  
3323     TMath::Sort(entries,chi2,index,kFALSE);
3324     besttrack = (AliITStrackMI*)array->At(index[0]);
3325     if (besttrack){
3326       AliDebug(2,Form("ncls best track %d     %f   %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3327       //
3328       for (Int_t j=0;j<6;j++){
3329         if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3330           erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3331           errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3332           ny[j]   = besttrack->GetNy(j);
3333           nz[j]   = besttrack->GetNz(j);
3334         }
3335       }
3336       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3337       minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3338       Float_t minn = besttrack->GetNumberOfClusters()-3;
3339       Int_t accepted=0;
3340       for (Int_t i=0;i<entries;i++){
3341         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);    
3342         if (!track) continue;
3343         if (accepted>maxcut) break;
3344         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3345         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3346           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3347             delete array->RemoveAt(index[i]);
3348             continue;
3349           }
3350         }
3351         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3352         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3353           if (!shortbest) accepted++;
3354           //
3355           newarray->AddLast(array->RemoveAt(index[i]));      
3356           for (Int_t j=0;j<6;j++){
3357             if (nz[j]==0){
3358               erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3359               errz[j] = track->GetSigmaZ(j); errz[j]   = track->GetSigmaZ(j+6);
3360               ny[j]   = track->GetNy(j);
3361               nz[j]   = track->GetNz(j);
3362             }
3363           }
3364         }
3365         else{
3366           delete array->RemoveAt(index[i]);
3367         }
3368       }
3369       array->Delete();
3370       delete fTrackHypothesys.RemoveAt(esdindex);
3371       fTrackHypothesys.AddAt(newarray,esdindex);
3372     }
3373     else{
3374       array->Delete();
3375       delete fTrackHypothesys.RemoveAt(esdindex);
3376     }
3377   }
3378   delete [] chi2;
3379   delete [] index;
3380 }
3381 //------------------------------------------------------------------------
3382 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3383 {
3384   //-------------------------------------------------------------
3385   // try to find best hypothesy
3386   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3387   //-------------------------------------------------------------
3388   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3389   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3390   if (!array) return 0;
3391   Int_t entries = array->GetEntriesFast();
3392   if (!entries) return 0;  
3393   Float_t minchi2 = 100000;
3394   AliITStrackMI * besttrack=0;
3395   //
3396   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
3397   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3398   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};     
3399   Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3400   //
3401   for (Int_t i=0;i<entries;i++){    
3402     AliITStrackMI * track = (AliITStrackMI*)array->At(i);    
3403     if (!track) continue;
3404     Float_t sigmarfi,sigmaz;
3405     GetDCASigma(track,sigmarfi,sigmaz);
3406     track->SetDnorm(0,sigmarfi);
3407     track->SetDnorm(1,sigmaz);
3408     //
3409     track->SetChi2MIP(1,1000000);
3410     track->SetChi2MIP(2,1000000);
3411     track->SetChi2MIP(3,1000000);
3412     //
3413     // backtrack
3414     backtrack = new(backtrack) AliITStrackMI(*track); 
3415     if (track->GetConstrain()) {
3416       if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3417       if (!backtrack->Improve(0,xyzVtx,ersVtx))         continue;     
3418       backtrack->ResetCovariance(10.);      
3419     }else{
3420       backtrack->ResetCovariance(10.);
3421     }
3422     backtrack->ResetClusters();
3423
3424     Double_t x = original->GetX();
3425     if (!RefitAt(x,backtrack,track)) continue;
3426     //
3427     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3428     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3429     if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.)  continue;
3430     track->SetChi22(GetMatchingChi2(backtrack,original));
3431
3432     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
3433     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
3434     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
3435
3436
3437     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1))  continue;
3438     //
3439     //forward track - without constraint
3440     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3441     forwardtrack->ResetClusters();
3442     x = track->GetX();
3443     RefitAt(x,forwardtrack,track);
3444     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));    
3445     if  (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0)  continue;
3446     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2))  continue;
3447     
3448     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3449     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3450     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
3451     forwardtrack->SetD(0,track->GetD(0));
3452     forwardtrack->SetD(1,track->GetD(1));    
3453     {
3454       Int_t list[6];
3455       AliITSRecPoint* clist[6];
3456       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));      
3457       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3458     }
3459     
3460     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3461     if  ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;    
3462     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3463       track->SetChi2MIP(3,1000);
3464       continue; 
3465     }
3466     Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();    
3467     //
3468     for (Int_t ichi=0;ichi<5;ichi++){
3469       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3470     }
3471     if (chi2 < minchi2){
3472       //besttrack = new AliITStrackMI(*forwardtrack);
3473       besttrack = track;
3474       besttrack->SetLabel(track->GetLabel());
3475       besttrack->SetFakeRatio(track->GetFakeRatio());
3476       minchi2   = chi2;
3477       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3478       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3479       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
3480     }    
3481   }
3482   delete backtrack;
3483   delete forwardtrack;
3484   Int_t accepted=0;
3485   for (Int_t i=0;i<entries;i++){    
3486     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3487    
3488     if (!track) continue;
3489     
3490     if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. || 
3491         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3492         track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3493       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3494         delete array->RemoveAt(i);    
3495         continue;
3496       }
3497     }
3498     else{
3499       accepted++;
3500     }
3501   }
3502   //
3503   array->Compress();
3504   SortTrackHypothesys(esdindex,checkmax,1);
3505
3506   array = (TObjArray*) fTrackHypothesys.At(esdindex);
3507   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3508   besttrack = (AliITStrackMI*)array->At(0);  
3509   if (!besttrack)  return 0;
3510   besttrack->SetChi2MIP(8,0);
3511   fBestTrackIndex[esdindex]=0;
3512   entries = array->GetEntriesFast();
3513   AliITStrackMI *longtrack =0;
3514   minchi2 =1000;
3515   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3516   for (Int_t itrack=entries-1;itrack>0;itrack--) {