ITS/AliITStrackerMI.cxx

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