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