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