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