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