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