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