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