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