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