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