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