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