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 /* $Id$ */
  17
  18 //-------------------------------------------------------------------------
  19 //               Implementation of the ITS tracker class
  20 //    It reads AliITSRecPoint clusters and creates AliITStrackMI tracks
  21 //                   and fills with them the ESD
  22 //          Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
  23 //     dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
  24 //
  25 //-------------------------------------------------------------------------
  26
  27 #include <TMatrixD.h>
  28 #include <TTree.h>
  29 #include <TTreeStream.h>
  30 #include <TTree.h>
  31
  32 #include "AliESD.h"
  33 #include "AliV0.h"
  34 #include "AliHelix.h"
  35 #include "AliITSRecPoint.h"
  36 #include "AliITSgeomTGeo.h"
  37 #include "AliITStrackerMI.h"
  38 #include "AliTrackPointArray.h"
  39 #include "AliAlignObj.h"
  40
  41 ClassImp(AliITStrackerMI)
  42
  43
  44
  45 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[kMaxLayer]; // ITS layers
  46 AliITStrackerMI::AliITStrackerMI():AliTracker(),
  47 fI(0),
  48 fBestTrack(),
  49 fTrackToFollow(),
  50 fTrackHypothesys(),
  51 fBestHypothesys(),
  52 fOriginal(),
  53 fCurrentEsdTrack(),
  54 fPass(0),
  55 fAfterV0(kFALSE),
  56 fLastLayerToTrackTo(0),
  57 fCoeficients(0),
  58 fEsd(0),
  59 fDebugStreamer(0){
  60   //Default constructor
  61 }
  62
  63
  64 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
  65 fI(kMaxLayer),
  66 fBestTrack(),
  67 fTrackToFollow(),
  68 fTrackHypothesys(),
  69 fBestHypothesys(),
  70 fOriginal(),
  71 fCurrentEsdTrack(),
  72 fPass(0),
  73 fAfterV0(kFALSE),
  74 fLastLayerToTrackTo(kLastLayerToTrackTo),
  75 fCoeficients(0),
  76 fEsd(0),
  77 fDebugStreamer(0){
  78   //--------------------------------------------------------------------
  79   //This is the AliITStrackerMI constructor
  80   //--------------------------------------------------------------------
  81   if (geom) {
  82     AliWarning("\"geom\" is actually a dummy argument !");
  83   }
  84
  85   fCoeficients = 0;
  86   fAfterV0     = kFALSE;
  87
  88   for (Int_t i=1; i<kMaxLayer+1; i++) {
  89     Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
  90     Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
  91
  92     Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
  93     AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
  94     Double_t poff=TMath::ATan2(y,x);
  95     Double_t zoff=z;
  96     Double_t r=TMath::Sqrt(x*x + y*y);
  97
  98     AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
  99     r += TMath::Sqrt(x*x + y*y);
 100     AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
 101     r += TMath::Sqrt(x*x + y*y);
 102     AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
 103     r += TMath::Sqrt(x*x + y*y);
 104     r*=0.25;
 105
 106     new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
 107
 108     for (Int_t j=1; j<nlad+1; j++) {
 109       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
 110         TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
 111         const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
 112         m.Multiply(tm);
 113         Double_t txyz[3]={0.}, xyz[3]={0.};
 114         m.LocalToMaster(txyz,xyz);
 115         Double_t r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 116         Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
 117
 118         if (phi<0) phi+=TMath::TwoPi();
 119         else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
 120
 121         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 122         new(&det) AliITSdetector(r,phi);
 123       }
 124     }
 125
 126   }
 127
 128   fI=kMaxLayer;
 129
 130   fPass=0;
 131   fConstraint[0]=1; fConstraint[1]=0;
 132
 133   Double_t xyz[]={kXV,kYV,kZV}, ers[]={kSigmaXV,kSigmaYV,kSigmaZV};
 134   SetVertex(xyz,ers);
 135
 136   for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=kLayersNotToSkip[i];
 137   fLastLayerToTrackTo=kLastLayerToTrackTo;
 138   for (Int_t i=0;i<100000;i++){
 139     fBestTrackIndex[i]=0;
 140   }
 141   //
 142   fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
 143
 144 }
 145
 146 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
 147 fI(tracker.fI),
 148 fBestTrack(tracker.fBestTrack),
 149 fTrackToFollow(tracker.fTrackToFollow),
 150 fTrackHypothesys(tracker.fTrackHypothesys),
 151 fBestHypothesys(tracker.fBestHypothesys),
 152 fOriginal(tracker.fOriginal),
 153 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
 154 fPass(tracker.fPass),
 155 fAfterV0(tracker.fAfterV0),
 156 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
 157 fCoeficients(tracker.fCoeficients),
 158 fEsd(tracker.fEsd),
 159 fDebugStreamer(tracker.fDebugStreamer){
 160   //Copy constructor
 161 }
 162
 163 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
 164   //Assignment operator
 165   this->~AliITStrackerMI();
 166   new(this) AliITStrackerMI(tracker);
 167   return *this;
 168 }
 169
 170
 171 AliITStrackerMI::~AliITStrackerMI()
 172 {
 173   //
 174   //destructor
 175   //
 176   if (fCoeficients) delete []fCoeficients;
 177   if (fDebugStreamer) {
 178     //fDebugStreamer->Close();
 179     delete fDebugStreamer;
 180   }
 181 }
 182
 183 void AliITStrackerMI::SetLayersNotToSkip(Int_t *l) {
 184   //--------------------------------------------------------------------
 185   //This function set masks of the layers which must be not skipped
 186   //--------------------------------------------------------------------
 187   for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=l[i];
 188 }
 189
 190 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
 191   //--------------------------------------------------------------------
 192   //This function loads ITS clusters
 193   //--------------------------------------------------------------------
 194   TBranch *branch=cTree->GetBranch("ITSRecPoints");
 195   if (!branch) {
 196     Error("LoadClusters"," can't get the branch !\n");
 197     return 1;
 198   }
 199
 200   TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
 201   branch->SetAddress(&clusters);
 202
 203   Int_t j=0;
 204   Int_t detector=0;
 205   for (Int_t i=0; i<kMaxLayer; i++) {
 206     Int_t ndet=fgLayers[i].GetNdetectors();
 207     Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
 208     for (; j<jmax; j++) {
 209       if (!cTree->GetEvent(j)) continue;
 210       Int_t ncl=clusters->GetEntriesFast();
 211       SignDeltas(clusters,GetZ());
 212
 213       //*** The Delta transformation
 214       TGeoHMatrix delta;
 215       AliITSgeomTGeo::GetTrackingMatrix(j,delta);
 216       TGeoHMatrix m;
 217       AliITSgeomTGeo::GetOrigMatrix(j,m);
 218       delta.MultiplyLeft(&(m.Inverse()));
 219       const TGeoHMatrix *tm = AliITSgeomTGeo::GetTracking2LocalMatrix(j);
 220       delta.MultiplyLeft(&(tm->Inverse()));
 221       //***
 222
 223       while (ncl--) {
 224         AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
 225
 226         //*** Shift the cluster to the misaligned position
 227         Double_t  xyz[3]={c->GetX(),c->GetY(),c->GetZ()};
 228         Double_t mxyz[3]={0.};
 229
 230         delta.LocalToMaster(xyz,mxyz);
 231
 232         c->SetX(mxyz[0]);
 233         c->AliCluster::SetY(mxyz[1]);
 234         c->AliCluster::SetZ(mxyz[2]);
 235         //***
 236
 237         fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
 238       }
 239       clusters->Delete();
 240       //add dead zone virtual "cluster"
 241       if (i<2){
 242         for (Float_t ydead = 0; ydead < 1.31 ; ydead+=(i+1.)*0.018){
 243           Int_t lab[4] = {0,0,0,detector};
 244           Int_t info[3] = {0,0,0};
 245           Float_t hit[5]={0,0,0.004/12.,0.001/12.,0};
 246           if (i==0) hit[0] =ydead-0.4;
 247           if (i==1) hit[0]=ydead-3.75;
 248           hit[1] =-0.04;
 249           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
 250             fgLayers[i].InsertCluster(new AliITSRecPoint(lab, hit, info));
 251           hit[1]=-7.05;
 252           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
 253             fgLayers[i].InsertCluster(new AliITSRecPoint(lab, hit, info));
 254           hit[1]=-7.15;
 255           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
 256             fgLayers[i].InsertCluster(new AliITSRecPoint(lab, hit, info));
 257           hit[1] =0.06;
 258           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
 259             fgLayers[i].InsertCluster(new AliITSRecPoint(lab, hit, info));
 260           hit[1]=7.05;
 261           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
 262             fgLayers[i].InsertCluster(new AliITSRecPoint(lab, hit, info));
 263           hit[1]=7.25;
 264           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
 265             fgLayers[i].InsertCluster(new AliITSRecPoint(lab, hit, info));
 266         }
 267       }
 268
 269     }
 270     //
 271     fgLayers[i].ResetRoad(); //road defined by the cluster density
 272     fgLayers[i].SortClusters();
 273   }
 274
 275   return 0;
 276 }
 277
 278 void AliITStrackerMI::UnloadClusters() {
 279   //--------------------------------------------------------------------
 280   //This function unloads ITS clusters
 281   //--------------------------------------------------------------------
 282   for (Int_t i=0; i<kMaxLayer; i++) fgLayers[i].ResetClusters();
 283 }
 284
 285 static Int_t CorrectForDeadZoneMaterial(AliITStrackMI *t) {
 286   //--------------------------------------------------------------------
 287   // Correction for the material between the TPC and the ITS
 288   // (should it belong to the TPC code ?)
 289   //--------------------------------------------------------------------
 290   Double_t riw=80., diw=0.0053, x0iw=30; // TPC inner wall ?
 291   Double_t rcd=61., dcd=0.0053, x0cd=30; // TPC "central drum" ?
 292   Double_t yr=12.8, dr=0.03; // rods ?
 293   Double_t zm=0.2, dm=0.40;  // membrane
 294   //Double_t rr=52., dr=0.19, x0r=24., yyr=7.77; //rails
 295   Double_t rs=50., ds=0.001; // something belonging to the ITS (screen ?)
 296
 297   if (t->GetX() > riw) {
 298      if (!t->PropagateTo(riw,diw,x0iw)) return 1;
 299      if (TMath::Abs(t->GetY())>yr) t->CorrectForMaterial(dr);
 300      if (TMath::Abs(t->GetZ())<zm) t->CorrectForMaterial(dm);
 301      if (!t->PropagateTo(rcd,dcd,x0cd)) return 1;
 302      //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
 303      //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,dr,x0r);
 304      if (!t->PropagateTo(rs,ds)) return 1;
 305   } else if (t->GetX() < rs) {
 306      if (!t->PropagateTo(rs,-ds)) return 1;
 307      //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
 308      //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,-dr,x0r);
 309      if (!t->PropagateTo(rcd,-dcd,x0cd)) return 1;
 310      if (!t->PropagateTo(riw+0.001,-diw,x0iw)) return 1;
 311   } else {
 312   ::Error("CorrectForDeadZoneMaterial","track is already in the dead zone !");
 313     return 1;
 314   }
 315
 316   return 0;
 317 }
 318
 319 Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) {
 320   //--------------------------------------------------------------------
 321   // This functions reconstructs ITS tracks
 322   // The clusters must be already loaded !
 323   //--------------------------------------------------------------------
 324   TObjArray itsTracks(15000);
 325   fOriginal.Clear();
 326   fEsd = event;         // store pointer to the esd
 327   {/* Read ESD tracks */
 328     Int_t nentr=event->GetNumberOfTracks();
 329     Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
 330     while (nentr--) {
 331       AliESDtrack *esd=event->GetTrack(nentr);
 332
 333       if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
 334       if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
 335       if (esd->GetStatus()&AliESDtrack::kITSin) continue;
 336       if (esd->GetKinkIndex(0)>0) continue;   //kink daughter
 337       AliITStrackMI *t=0;
 338       try {
 339         t=new AliITStrackMI(*esd);
 340       } catch (const Char_t *msg) {
 341         //Warning("Clusters2Tracks",msg);
 342         delete t;
 343         continue;
 344       }
 345       //t->fD[0] = t->GetD(GetX(),GetY());
 346       //t->fD[1] = t->GetZat(GetX())-GetZ();
 347       t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP());              //I.B.
 348       Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
 349       if (t->GetMass()<0.13) t->SetMass(0.13957); // MI look to the esd - mass hypothesys  !!!!!!!!!!!
 350       // write expected q
 351       t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 352
 353       if (esd->GetV0Index(0)>0 && t->GetD(0)<30){
 354         //track - can be  V0 according to TPC
 355       }
 356       else{
 357         if (TMath::Abs(t->GetD(0))>10) {
 358           delete t;
 359           continue;
 360         }
 361
 362         if (TMath::Abs(vdist)>20) {
 363           delete t;
 364           continue;
 365         }
 366         if (TMath::Abs(1/t->Get1Pt())<0.120) {
 367           delete t;
 368           continue;
 369         }
 370
 371         if (CorrectForDeadZoneMaterial(t)!=0) {
 372           //Warning("Clusters2Tracks",
 373           //        "failed to correct for the material in the dead zone !\n");
 374           delete t;
 375           continue;
 376         }
 377       }
 378       t->SetReconstructed(kFALSE);
 379       itsTracks.AddLast(t);
 380       fOriginal.AddLast(t);
 381     }
 382   } /* End Read ESD tracks */
 383
 384   itsTracks.Sort();
 385   fOriginal.Sort();
 386   Int_t nentr=itsTracks.GetEntriesFast();
 387   fTrackHypothesys.Expand(nentr);
 388   fBestHypothesys.Expand(nentr);
 389   MakeCoeficients(nentr);
 390   Int_t ntrk=0;
 391   for (fPass=0; fPass<2; fPass++) {
 392      Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
 393      for (Int_t i=0; i<nentr; i++) {
 394 //       cerr<<fPass<<"    "<<i<<'\r';
 395        fCurrentEsdTrack = i;
 396        AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(i);
 397        if (t==0) continue;              //this track has been already tracked
 398        if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue;  //this track was  already  "succesfully" reconstructed
 399        //if ( (TMath::Abs(t->GetD(GetX(),GetY()))  >3.) && fConstraint[fPass]) continue;
 400        //if ( (TMath::Abs(t->GetZat(GetX())-GetZ())>3.) && fConstraint[fPass]) continue;
 401        Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz);              //I.B.
 402        if ( (TMath::Abs(dz[0])>3.) && fConstraint[fPass]) continue;
 403        if ( (TMath::Abs(dz[1])>3.) && fConstraint[fPass]) continue;
 404
 405        Int_t tpcLabel=t->GetLabel(); //save the TPC track label
 406        fI = 6;
 407        ResetTrackToFollow(*t);
 408        ResetBestTrack();
 409        FollowProlongationTree(t,i,fConstraint[fPass]);
 410
 411        SortTrackHypothesys(fCurrentEsdTrack,20,0);  //MI change
 412        //
 413        AliITStrackMI * besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
 414        if (!besttrack) continue;
 415        besttrack->SetLabel(tpcLabel);
 416        //       besttrack->CookdEdx();
 417        CookdEdx(besttrack);
 418        besttrack->SetFakeRatio(1.);
 419        CookLabel(besttrack,0.); //For comparison only
 420        UpdateESDtrack(besttrack,AliESDtrack::kITSin);
 421
 422        /*
 423        if ( besttrack->GetNumberOfClusters()<6 && fConstraint[fPass]) {
 424          continue;
 425        }
 426        if (besttrack->fChi2MIP[0]+besttrack->fNUsed>3.5) continue;
 427        if ( (TMath::Abs(besttrack->fD[0]*besttrack->fD[0]+besttrack->fD[1]*besttrack->fD[1])>0.1) && fConstraint[fPass])  continue;
 428        //delete itsTracks.RemoveAt(i);
 429        */
 430        if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue;  //to be tracked also without vertex constrain
 431
 432
 433        t->SetReconstructed(kTRUE);
 434        ntrk++;
 435      }
 436      GetBestHypothesysMIP(itsTracks);
 437   }
 438
 439   //GetBestHypothesysMIP(itsTracks);
 440   UpdateTPCV0(event);
 441   FindV02(event);
 442   fAfterV0 = kTRUE;
 443   //GetBestHypothesysMIP(itsTracks);
 444   //
 445   itsTracks.Delete();
 446   //
 447   Int_t entries = fTrackHypothesys.GetEntriesFast();
 448   for (Int_t ientry=0;ientry<entries;ientry++){
 449     TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
 450     if (array) array->Delete();
 451     delete fTrackHypothesys.RemoveAt(ientry);
 452   }
 453
 454   fTrackHypothesys.Delete();
 455   fBestHypothesys.Delete();
 456   fOriginal.Clear();
 457   delete []fCoeficients;
 458   fCoeficients=0;
 459   Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
 460
 461   return 0;
 462 }
 463
 464
 465 Int_t AliITStrackerMI::PropagateBack(AliESD *event) {
 466   //--------------------------------------------------------------------
 467   // This functions propagates reconstructed ITS tracks back
 468   // The clusters must be loaded !
 469   //--------------------------------------------------------------------
 470   Int_t nentr=event->GetNumberOfTracks();
 471   Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
 472
 473   Int_t ntrk=0;
 474   for (Int_t i=0; i<nentr; i++) {
 475      AliESDtrack *esd=event->GetTrack(i);
 476
 477      if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
 478      if (esd->GetStatus()&AliESDtrack::kITSout) continue;
 479
 480      AliITStrackMI *t=0;
 481      try {
 482         t=new AliITStrackMI(*esd);
 483      } catch (const Char_t *msg) {
 484        //Warning("PropagateBack",msg);
 485         delete t;
 486         continue;
 487      }
 488      t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 489
 490      ResetTrackToFollow(*t);
 491
 492      // propagete to vertex [SR, GSI 17.02.2003]
 493      // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
 494      if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {
 495        if (fTrackToFollow.PropagateToVertex(event->GetVertex())) {
 496           fTrackToFollow.StartTimeIntegral();
 497        }
 498        fTrackToFollow.PropagateTo(3.,-0.0028,65.19);
 499      }
 500
 501      fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
 502      if (RefitAt(49.,&fTrackToFollow,t)) {
 503         if (CorrectForDeadZoneMaterial(&fTrackToFollow)!=0) {
 504           //Warning("PropagateBack",
 505           //        "failed to correct for the material in the dead zone !\n");
 506           delete t;
 507           continue;
 508         }
 509         fTrackToFollow.SetLabel(t->GetLabel());
 510         //fTrackToFollow.CookdEdx();
 511         CookLabel(&fTrackToFollow,0.); //For comparison only
 512         fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
 513         //UseClusters(&fTrackToFollow);
 514         ntrk++;
 515      }
 516      delete t;
 517   }
 518
 519   Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
 520
 521   return 0;
 522 }
 523
 524 Int_t AliITStrackerMI::RefitInward(AliESD *event) {
 525   //--------------------------------------------------------------------
 526   // This functions refits ITS tracks using the
 527   // "inward propagated" TPC tracks
 528   // The clusters must be loaded !
 529   //--------------------------------------------------------------------
 530   RefitV02(event);
 531   Int_t nentr=event->GetNumberOfTracks();
 532   Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
 533
 534   Int_t ntrk=0;
 535   for (Int_t i=0; i<nentr; i++) {
 536     AliESDtrack *esd=event->GetTrack(i);
 537
 538     if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
 539     if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
 540     if (esd->GetStatus()&AliESDtrack::kTPCout)
 541       if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
 542
 543     AliITStrackMI *t=0;
 544     try {
 545         t=new AliITStrackMI(*esd);
 546     } catch (const Char_t *msg) {
 547       //Warning("RefitInward",msg);
 548         delete t;
 549         continue;
 550     }
 551     t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 552     if (CorrectForDeadZoneMaterial(t)!=0) {
 553       //Warning("RefitInward",
 554       //         "failed to correct for the material in the dead zone !\n");
 555        delete t;
 556        continue;
 557     }
 558
 559     ResetTrackToFollow(*t);
 560     fTrackToFollow.ResetClusters();
 561
 562     if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
 563       fTrackToFollow.ResetCovariance(10.);
 564
 565     //Refitting...
 566     if (RefitAt(3.7, &fTrackToFollow, t,kTRUE)) {
 567        fTrackToFollow.SetLabel(t->GetLabel());
 568        //       fTrackToFollow.CookdEdx();
 569        CookdEdx(&fTrackToFollow);
 570
 571        CookLabel(&fTrackToFollow,0.0); //For comparison only
 572
 573        if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {//The beam pipe
 574          AliESDtrack  *esdTrack =fTrackToFollow.GetESDtrack();
 575          esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit);
 576          Float_t r[3]={0.,0.,0.};
 577          Double_t maxD=3.;
 578          esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
 579          ntrk++;
 580        }
 581     }
 582     delete t;
 583   }
 584
 585   Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
 586
 587   return 0;
 588 }
 589
 590 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
 591   //--------------------------------------------------------------------
 592   //       Return pointer to a given cluster
 593   //--------------------------------------------------------------------
 594   Int_t l=(index & 0xf0000000) >> 28;
 595   Int_t c=(index & 0x0fffffff) >> 00;
 596   return fgLayers[l].GetCluster(c);
 597 }
 598
 599 #include "TGeoManager.h"
 600 #include "TGeoMatrix.h"
 601 #include "TGeoPhysicalNode.h"
 602
 603 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
 604   //
 605   // Get track space point with index i
 606   //
 607   Int_t l=(index & 0xf0000000) >> 28;
 608   Int_t c=(index & 0x0fffffff) >> 00;
 609   AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 610   Int_t idet = cl->GetDetectorIndex();
 611
 612   const char* name = AliAlignObj::SymName((AliAlignObj::ELayerID)
 613                                           (l+AliAlignObj::kFirstLayer), idet);
 614   TGeoPNEntry *mapPN = gGeoManager->GetAlignableEntry(name);
 615
 616   if (!mapPN) return kFALSE;
 617   TGeoPhysicalNode *node = mapPN->GetPhysicalNode();
 618   if (!node) {
 619     gGeoManager->MakeAlignablePN(name);
 620     node = mapPN->GetPhysicalNode();
 621   }
 622   if (!node) return kFALSE;
 623   TGeoHMatrix* matrix = node->GetMatrix();
 624   if (!matrix) return kFALSE;
 625
 626   //
 627   // Calculate the global coordinates
 628   //
 629   Double_t localCoord[]  = {cl->GetDetLocalX(), 0, cl->GetDetLocalZ()};
 630   // LG  AliAlignObj makes life simple but keep in mind that alignable
 631   // LG  volume doesn't mean sensitive volume. There might be a shift
 632   // LG  between the 2, has it is here for the SPD :
 633   if (l<2) localCoord[1] = 0.01;
 634   // LG   !!!   Check for this when the new geometry comes   !!!
 635   Double_t globalCoord[3] = {0};
 636   matrix->LocalToMaster(localCoord, globalCoord);
 637   Float_t xyz[3]= {globalCoord[0], globalCoord[1], globalCoord[2]};
 638
 639   //
 640   // Calculate the cov matrix
 641   //
 642   TGeoRotation rotMatrix(*matrix);
 643   TGeoRotation rotMatrixTr(rotMatrix.Inverse());
 644   Double_t sigmaArray[] = {cl->GetSigmaY2(),0,0, 0,0,0, 0,0,cl->GetSigmaZ2()};
 645   TGeoRotation sigmMatrix;
 646   sigmMatrix.SetMatrix( sigmaArray );
 647   sigmMatrix.MultiplyBy(&rotMatrixTr, kFALSE);
 648   sigmMatrix.MultiplyBy(&rotMatrix, kTRUE);
 649   const Double_t *globalSigma =  sigmMatrix.GetRotationMatrix();
 650   Float_t cov[6]= { globalSigma[0], globalSigma[1], globalSigma[2],
 651                                     globalSigma[4], globalSigma[5],
 652                                                     globalSigma[8] };
 653
 654   p.SetXYZ(xyz[0],xyz[1],xyz[2],cov);
 655   AliAlignObj::ELayerID iLayer = AliAlignObj::kInvalidLayer;
 656   switch (l) {
 657   case 0:
 658     iLayer = AliAlignObj::kSPD1;
 659     break;
 660   case 1:
 661     iLayer = AliAlignObj::kSPD2;
 662     break;
 663   case 2:
 664     iLayer = AliAlignObj::kSDD1;
 665     break;
 666   case 3:
 667     iLayer = AliAlignObj::kSDD2;
 668     break;
 669   case 4:
 670     iLayer = AliAlignObj::kSSD1;
 671     break;
 672   case 5:
 673     iLayer = AliAlignObj::kSSD2;
 674     break;
 675   default:
 676     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 677     break;
 678   };
 679   UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,idet);
 680   p.SetVolumeID((UShort_t)volid);
 681   return kTRUE;
 682 }
 683
 684 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
 685 {
 686   //--------------------------------------------------------------------
 687   // Follow prolongation tree
 688   //--------------------------------------------------------------------
 689   //
 690   AliESDtrack * esd = otrack->GetESDtrack();
 691   if (esd->GetV0Index(0)>0){
 692     //
 693     // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
 694     //                      mapping of esd track is different as its track in Containers
 695     //                      Need something more stable
 696     //                      Indexes are set back againg to the ESD track indexes in UpdateTPCV0
 697     for (Int_t i=0;i<3;i++){
 698       Int_t  index = esd->GetV0Index(i);
 699       if (index==0) break;
 700       AliESDv0 * vertex = fEsd->GetV0(index);
 701       if (vertex->GetStatus()<0) continue;     // rejected V0
 702       //
 703             if (esd->GetSign()>0) {
 704               vertex->SetIndex(0,esdindex);
 705             }
 706             else{
 707               vertex->SetIndex(1,esdindex);
 708             }
 709     }
 710   }
 711   TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
 712   if (!bestarray){
 713     bestarray = new TObjArray(5);
 714     fBestHypothesys.AddAt(bestarray,esdindex);
 715   }
 716
 717   //
 718   //setup tree of the prolongations
 719   //
 720   static AliITStrackMI tracks[7][100];
 721   AliITStrackMI *currenttrack;
 722   static AliITStrackMI currenttrack1;
 723   static AliITStrackMI currenttrack2;
 724   static AliITStrackMI backuptrack;
 725   Int_t ntracks[7];
 726   Int_t nindexes[7][100];
 727   Float_t normalizedchi2[100];
 728   for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
 729   otrack->SetNSkipped(0);
 730   new (&(tracks[6][0])) AliITStrackMI(*otrack);
 731   ntracks[6]=1;
 732   for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
 733   //
 734   //
 735   // follow prolongations
 736   for (Int_t ilayer=5;ilayer>=0;ilayer--){
 737     //
 738     AliITSlayer &layer=fgLayers[ilayer];
 739     Double_t r=layer.GetR();
 740     ntracks[ilayer]=0;
 741     //
 742     //
 743    Int_t nskipped=0;
 744     Float_t nused =0;
 745     for (Int_t itrack =0;itrack<ntracks[ilayer+1];itrack++){
 746       //set current track
 747       if (ntracks[ilayer]>=100) break;
 748       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
 749       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
 750       if (ntracks[ilayer]>15+ilayer){
 751         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
 752         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
 753       }
 754
 755       new(&currenttrack1)  AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
 756       if (ilayer==3 || ilayer==1) {
 757         Double_t rs=0.5*(fgLayers[ilayer+1].GetR() + r);
 758         Double_t d=0.0034, x0=38.6;
 759         if (ilayer==1) {rs=9.; d=0.0097; x0=42;}
 760         if (!currenttrack1.PropagateTo(rs,d,x0)) {
 761           continue;
 762         }
 763       }
 764       //
 765       //find intersection with layer
 766       Double_t x,y,z;
 767       if (!currenttrack1.GetGlobalXYZat(r,x,y,z)) {
 768         continue;
 769       }
 770       Double_t phi=TMath::ATan2(y,x);
 771       Int_t idet=layer.FindDetectorIndex(phi,z);
 772       if (idet<0) {
 773         continue;
 774       }
 775       //propagate to the intersection
 776       const AliITSdetector &det=layer.GetDetector(idet);
 777       phi=det.GetPhi();
 778       new(&currenttrack2)  AliITStrackMI(currenttrack1);
 779       if (!currenttrack1.Propagate(phi,det.GetR())) {
 780         continue;
 781       }
 782       currenttrack2.Propagate(phi,det.GetR());  //
 783       currenttrack1.SetDetectorIndex(idet);
 784       currenttrack2.SetDetectorIndex(idet);
 785
 786       //
 787       //
 788       Double_t dz=7.5*TMath::Sqrt(currenttrack1.GetSigmaZ2() + 16.*kSigmaZ2[ilayer]);
 789       Double_t dy=7.5*TMath::Sqrt(currenttrack1.GetSigmaY2() + 16.*kSigmaY2[ilayer]);
 790       //
 791       Bool_t isBoundary=kFALSE;
 792       if (currenttrack1.GetY()-dy< det.GetYmin()+0.2) isBoundary = kTRUE;
 793       if (currenttrack1.GetY()+dy> det.GetYmax()-0.2) isBoundary = kTRUE;
 794       if (currenttrack1.GetZ()-dz< det.GetZmin()+0.2) isBoundary = kTRUE;
 795       if (currenttrack1.GetZ()+dz> det.GetZmax()-0.2) isBoundary = kTRUE;
 796
 797       if (isBoundary){ // track at boundary between detectors
 798         Float_t maxtgl = TMath::Abs(currenttrack1.GetTgl());
 799         if (maxtgl>1) maxtgl=1;
 800         dz = TMath::Sqrt(dz*dz+0.25*maxtgl*maxtgl);
 801         //
 802         Float_t maxsnp = TMath::Abs(currenttrack1.GetSnp());
 803         if (maxsnp>0.95) continue;
 804         //if (maxsnp>0.5) maxsnp=0.5;
 805         dy=TMath::Sqrt(dy*dy+0.25*maxsnp*maxsnp);
 806       }
 807
 808       Double_t zmin=currenttrack1.GetZ() - dz;
 809       Double_t zmax=currenttrack1.GetZ() + dz;
 810       Double_t ymin=currenttrack1.GetY() + r*phi - dy;
 811       Double_t ymax=currenttrack1.GetY() + r*phi + dy;
 812       layer.SelectClusters(zmin,zmax,ymin,ymax);
 813       //
 814       // loop over all possible prolongations
 815       //
 816       Double_t msz=1./((currenttrack1.GetSigmaZ2() + 16.*kSigmaZ2[ilayer]));
 817       Double_t msy=1./((currenttrack1.GetSigmaY2() + 16.*kSigmaY2[ilayer]));
 818       if (constrain){
 819         msy/=60; msz/=60.;
 820       }
 821       else{
 822         msy/=50; msz/=50.;
 823       }
 824       //
 825       const AliITSRecPoint *c=0; Int_t ci=-1;
 826       Double_t chi2=12345.;
 827       Int_t deadzone=0;
 828       currenttrack = &currenttrack1;
 829       while ((c=layer.GetNextCluster(ci))!=0) {
 830         if (ntracks[ilayer]>95) break; //space for skipped clusters
 831         Bool_t change =kFALSE;
 832         if (c->GetQ()==0 && (deadzone==1)) continue;
 833         Int_t idet=c->GetDetectorIndex();
 834         if (currenttrack->GetDetectorIndex()!=idet) {
 835           const AliITSdetector &det=layer.GetDetector(idet);
 836           Double_t y,z;
 837           if (!currenttrack2.GetProlongationFast(det.GetPhi(),det.GetR(),y,z)) continue;
 838           Float_t pz = (z - c->GetZ()) , py=(y - c->GetY());
 839           if (pz*pz*msz+py*py*msy>1.) continue;
 840           //
 841           new (&backuptrack) AliITStrackMI(currenttrack2);
 842           change = kTRUE;
 843           currenttrack =&currenttrack2;
 844           if (!currenttrack->Propagate(det.GetPhi(),det.GetR())) {
 845             new (currenttrack) AliITStrackMI(backuptrack);
 846             change = kFALSE;
 847             continue;
 848           }
 849           currenttrack->SetDetectorIndex(idet);
 850         }
 851         else{
 852           Float_t pz = (currenttrack->GetZ() - c->GetZ()) , py=(currenttrack->GetY() - c->GetY());
 853           if (pz*pz*msz+py*py*msy>1.) continue;
 854         }
 855
 856         chi2=GetPredictedChi2MI(currenttrack,c,ilayer);
 857         if (chi2<kMaxChi2s[ilayer]){
 858           if (c->GetQ()==0) deadzone=1;     // take dead zone only once
 859           if (ntracks[ilayer]>=100) continue;
 860           AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
 861           updatetrack->SetClIndex(ilayer,0);
 862           if (change){
 863             new (&currenttrack2) AliITStrackMI(backuptrack);
 864           }
 865           if (c->GetQ()!=0){
 866             if (!UpdateMI(updatetrack,c,chi2,(ilayer<<28)+ci)) continue;
 867             updatetrack->SetSampledEdx(c->GetQ(),updatetrack->GetNumberOfClusters()-1); //b.b.
 868           }
 869           else {
 870             updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
 871             updatetrack->SetDeadZoneProbability(GetDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
 872           }
 873           if (c->IsUsed()){
 874             updatetrack->IncrementNUsed();
 875           }
 876           Double_t x0;
 877           Double_t d=layer.GetThickness(updatetrack->GetY(),updatetrack->GetZ(),x0);
 878           updatetrack->CorrectForMaterial(d,x0);
 879           if (constrain) {
 880             updatetrack->SetConstrain(constrain);
 881             fI = ilayer;
 882             Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
 883             Double_t xyz[]={GetX(),GetY(),GetZ()};
 884             Double_t ptfactor = 1;
 885             Double_t ers[]={GetSigmaX()*ptfactor,GetSigmaY()*ptfactor,GetSigmaZ()};
 886             Bool_t isPrim = kTRUE;
 887             if (ilayer<4){
 888               //updatetrack->fD[0] = updatetrack->GetD(GetX(),GetY());
 889               //updatetrack->fD[1] = updatetrack->GetZat(GetX())-GetZ();
 890               updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
 891               if ( TMath::Abs(updatetrack->GetD(0)/(1.+ilayer))>0.4 ||  TMath::Abs(updatetrack->GetD(1)/(1.+ilayer))>0.4) isPrim=kFALSE;
 892             }
 893             if (isPrim) updatetrack->Improve(d,xyz,ers);
 894           } //apply vertex constrain
 895           ntracks[ilayer]++;
 896         }  // create new hypothesy
 897       } // loop over possible cluster prolongation
 898       //      if (constrain&&itrack<2&&currenttrack1.fNSkipped==0 && deadzone==0){
 899       if (constrain&&itrack<2&&currenttrack1.GetNSkipped()==0 && deadzone==0&&ntracks[ilayer]<100){
 900         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
 901         vtrack->SetClIndex(ilayer,0);
 902         fI = ilayer;
 903         Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
 904         Double_t xyz[]={GetX(),GetY(),GetZ()};
 905         Double_t ers[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
 906         vtrack->Improve(d,xyz,ers);
 907         vtrack->IncrementNSkipped();
 908         ntracks[ilayer]++;
 909       }
 910
 911       if (constrain&&itrack<1&&TMath::Abs(currenttrack1.GetTgl())>1.1){  //big theta -- for low mult. runs
 912         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
 913         vtrack->SetClIndex(ilayer,0);
 914         fI = ilayer;
 915         Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
 916         Double_t xyz[]={GetX(),GetY(),GetZ()};
 917         Double_t ers[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
 918         vtrack->Improve(d,xyz,ers);
 919         vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
 920         ntracks[ilayer]++;
 921       }
 922
 923
 924     } //loop over track candidates
 925     //
 926     //
 927     Int_t accepted=0;
 928
 929     Int_t golds=0;
 930     for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
 931       normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
 932       if ( normalizedchi2[itrack]<3+0.5*ilayer) golds++;
 933       if (ilayer>4) accepted++;
 934       else{
 935         if ( constrain && normalizedchi2[itrack]<kMaxNormChi2C[ilayer]+1) accepted++;
 936         if (!constrain && normalizedchi2[itrack]<kMaxNormChi2NonC[ilayer]+1) accepted++;
 937       }
 938     }
 939     TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
 940     ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
 941     if (ntracks[ilayer]<golds+2+ilayer) ntracks[ilayer]=TMath::Min(golds+2+ilayer,accepted);
 942     if (ntracks[ilayer]>90) ntracks[ilayer]=90;
 943   } //loop over layers
 944   //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]);
 945   Int_t max = constrain? 20: 5;
 946
 947   for (Int_t i=0;i<TMath::Min(max,ntracks[0]);i++) {
 948     AliITStrackMI & track= tracks[0][nindexes[0][i]];
 949     if (track.GetNumberOfClusters()<2) continue;
 950     if (!constrain&&track.GetNormChi2(0)>7.)continue;
 951     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
 952   }
 953   for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
 954     AliITStrackMI & track= tracks[1][nindexes[1][i]];
 955     if (track.GetNumberOfClusters()<4) continue;
 956     if (!constrain&&track.GetNormChi2(1)>7.)continue;
 957     if (constrain) track.IncrementNSkipped();
 958     if (!constrain) {
 959       track.SetD(0,track.GetD(GetX(),GetY()));
 960       track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
 961       if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
 962         track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
 963       }
 964     }
 965     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
 966   }
 967   //}
 968
 969   if (!constrain){
 970     for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
 971       AliITStrackMI & track= tracks[2][nindexes[2][i]];
 972       if (track.GetNumberOfClusters()<3) continue;
 973       if (!constrain&&track.GetNormChi2(2)>7.)continue;
 974       if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
 975       if (!constrain){
 976         track.SetD(0,track.GetD(GetX(),GetY()));
 977         track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
 978         if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
 979           track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
 980         }
 981       }
 982       AddTrackHypothesys(new AliITStrackMI(track), esdindex);
 983     }
 984   }
 985
 986   if (!constrain){
 987     //
 988     // register best tracks - important for V0 finder
 989     //
 990     for (Int_t ilayer=0;ilayer<5;ilayer++){
 991       if (ntracks[ilayer]==0) continue;
 992       AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
 993       if (track.GetNumberOfClusters()<1) continue;
 994       CookLabel(&track,0);
 995       bestarray->AddAt(new AliITStrackMI(track),ilayer);
 996     }
 997   }
 998   //
 999   // update TPC V0 information
1000   //
1001   if (otrack->GetESDtrack()->GetV0Index(0)>0){
1002     Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1003     for (Int_t i=0;i<3;i++){
1004       Int_t  index = otrack->GetESDtrack()->GetV0Index(i);
1005       if (index==0) break;
1006       AliV0 * vertex = (AliV0*)fEsd->GetV0(index);
1007       if (vertex->GetStatus()<0) continue;     // rejected V0
1008       //
1009       if (otrack->GetSign()>0) {
1010         vertex->SetIndex(0,esdindex);
1011       }
1012       else{
1013         vertex->SetIndex(1,esdindex);
1014       }
1015       //find nearest layer with track info
1016       Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
1017       Int_t nearestold  = GetNearestLayer(xrp);               //I.B.
1018       Int_t nearest     = nearestold;
1019       for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1020         if (ntracks[nearest]==0){
1021           nearest = ilayer;
1022         }
1023       }
1024       //
1025       AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1026       if (nearestold<5&&nearest<5){
1027         Bool_t accept = track.GetNormChi2(nearest)<10;
1028         if (accept){
1029           if (track.GetSign()>0) {
1030             vertex->SetParamP(track);
1031             vertex->Update(fprimvertex);
1032             //      vertex->SetIndex(0,track.fESDtrack->GetID());
1033             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1034           }else{
1035             vertex->SetParamN(track);
1036             vertex->Update(fprimvertex);
1037             //vertex->SetIndex(1,track.fESDtrack->GetID());
1038             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1039           }
1040           vertex->SetStatus(vertex->GetStatus()+1);
1041         }else{
1042           //  vertex->SetStatus(-2);  // reject V0  - not enough clusters
1043         }
1044       }
1045       // if (nearestold>3){
1046 //      Int_t indexlayer = (ntracks[0]>0)? 0:1;
1047 //      if (ntracks[indexlayer]>0){
1048 //        AliITStrackMI & track= tracks[indexlayer][nindexes[indexlayer][0]];
1049 //        if (track.GetNumberOfClusters()>4&&track.fNormChi2[indexlayer]<4){
1050 //          vertex->SetStatus(-1);  // reject V0 - clusters before
1051 //        }
1052 //      }
1053 //      }
1054     }
1055   }
1056 }
1057
1058
1059 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1060 {
1061   //--------------------------------------------------------------------
1062   //
1063   //
1064   return fgLayers[layer];
1065 }
1066
1067 AliITStrackerMI::AliITSlayer::AliITSlayer():
1068 fR(0),
1069 fPhiOffset(0),
1070 fNladders(0),
1071 fZOffset(0),
1072 fNdetectors(0),
1073 fDetectors(0),
1074 fN(0),
1075 fDy5(0),
1076 fDy10(0),
1077 fDy20(0),
1078 fClustersCs(0),
1079 fClusterIndexCs(0),
1080 fYcs(0),
1081 fZcs(0),
1082 fNcs(0),
1083 fCurrentSlice(-1),
1084 fZmax(0),
1085 fYmin(0),
1086 fYmax(0),
1087 fI(0),
1088 fImax(0),
1089 fSkip(0),
1090 fAccepted(0),
1091 fRoad(0){
1092   //--------------------------------------------------------------------
1093   //default AliITSlayer constructor
1094   //--------------------------------------------------------------------
1095   for (Int_t i=0; i<kMaxClusterPerLayer;i++) {
1096     fClusterWeight[i]=0;
1097     fClusterTracks[0][i]=-1;
1098     fClusterTracks[1][i]=-1;
1099     fClusterTracks[2][i]=-1;
1100     fClusterTracks[3][i]=-1;
1101   }
1102 }
1103
1104 AliITStrackerMI::AliITSlayer::
1105 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1106 fR(r),
1107 fPhiOffset(p),
1108 fNladders(nl),
1109 fZOffset(z),
1110 fNdetectors(nd),
1111 fDetectors(0),
1112 fN(0),
1113 fDy5(0),
1114 fDy10(0),
1115 fDy20(0),
1116 fClustersCs(0),
1117 fClusterIndexCs(0),
1118 fYcs(0),
1119 fZcs(0),
1120 fNcs(0),
1121 fCurrentSlice(-1),
1122 fZmax(0),
1123 fYmin(0),
1124 fYmax(0),
1125 fI(0),
1126 fImax(0),
1127 fSkip(0),
1128 fAccepted(0),
1129 fRoad(0) {
1130   //--------------------------------------------------------------------
1131   //main AliITSlayer constructor
1132   //--------------------------------------------------------------------
1133   fDetectors=new AliITSdetector[fNladders*fNdetectors];
1134   fRoad=2*fR*TMath::Sqrt(3.14/1.);//assuming that there's only one cluster
1135 }
1136
1137 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1138 fR(layer.fR),
1139 fPhiOffset(layer.fPhiOffset),
1140 fNladders(layer.fNladders),
1141 fZOffset(layer.fZOffset),
1142 fNdetectors(layer.fNdetectors),
1143 fDetectors(layer.fDetectors),
1144 fN(layer.fN),
1145 fDy5(layer.fDy5),
1146 fDy10(layer.fDy10),
1147 fDy20(layer.fDy20),
1148 fClustersCs(layer.fClustersCs),
1149 fClusterIndexCs(layer.fClusterIndexCs),
1150 fYcs(layer.fYcs),
1151 fZcs(layer.fZcs),
1152 fNcs(layer.fNcs),
1153 fCurrentSlice(layer.fCurrentSlice),
1154 fZmax(layer.fZmax),
1155 fYmin(layer.fYmin),
1156 fYmax(layer.fYmax),
1157 fI(layer.fI),
1158 fImax(layer.fImax),
1159 fSkip(layer.fSkip),
1160 fAccepted(layer.fAccepted),
1161 fRoad(layer.fRoad){
1162   //Copy constructor
1163 }
1164
1165
1166 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1167   //--------------------------------------------------------------------
1168   // AliITSlayer destructor
1169   //--------------------------------------------------------------------
1170   delete[] fDetectors;
1171   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1172   for (Int_t i=0; i<kMaxClusterPerLayer;i++) {
1173     fClusterWeight[i]=0;
1174     fClusterTracks[0][i]=-1;
1175     fClusterTracks[1][i]=-1;
1176     fClusterTracks[2][i]=-1;
1177     fClusterTracks[3][i]=-1;
1178   }
1179 }
1180
1181 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1182   //--------------------------------------------------------------------
1183   // This function removes loaded clusters
1184   //--------------------------------------------------------------------
1185   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1186   for (Int_t i=0; i<kMaxClusterPerLayer;i++){
1187     fClusterWeight[i]=0;
1188     fClusterTracks[0][i]=-1;
1189     fClusterTracks[1][i]=-1;
1190     fClusterTracks[2][i]=-1;
1191     fClusterTracks[3][i]=-1;
1192   }
1193
1194   fN=0;
1195   fI=0;
1196 }
1197
1198 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1199   //--------------------------------------------------------------------
1200   // This function reset weights of the clusters
1201   //--------------------------------------------------------------------
1202   for (Int_t i=0; i<kMaxClusterPerLayer;i++) {
1203     fClusterWeight[i]=0;
1204     fClusterTracks[0][i]=-1;
1205     fClusterTracks[1][i]=-1;
1206     fClusterTracks[2][i]=-1;
1207     fClusterTracks[3][i]=-1;
1208   }
1209   for (Int_t i=0; i<fN;i++) {
1210     AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1211     if (cl&&cl->IsUsed()) cl->Use();
1212   }
1213
1214 }
1215
1216 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1217   //--------------------------------------------------------------------
1218   // This function calculates the road defined by the cluster density
1219   //--------------------------------------------------------------------
1220   Int_t n=0;
1221   for (Int_t i=0; i<fN; i++) {
1222      if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1223   }
1224   //if (n>1) fRoad=2*fR*TMath::Sqrt(3.14/n);
1225   if (n>1) fRoad=2*fR*TMath::Sqrt(3.14/n);
1226 }
1227
1228
1229 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *c) {
1230   //--------------------------------------------------------------------
1231   //This function adds a cluster to this layer
1232   //--------------------------------------------------------------------
1233   if (fN==kMaxClusterPerLayer) {
1234     ::Error("InsertCluster","Too many clusters !\n");
1235     return 1;
1236   }
1237   fCurrentSlice=-1;
1238   fClusters[fN]=c;
1239   fN++;
1240   AliITSdetector &det=GetDetector(c->GetDetectorIndex());
1241   if (c->GetY()<det.GetYmin()) det.SetYmin(c->GetY());
1242   if (c->GetY()>det.GetYmax()) det.SetYmax(c->GetY());
1243   if (c->GetZ()<det.GetZmin()) det.SetZmin(c->GetZ());
1244   if (c->GetZ()>det.GetZmax()) det.SetZmax(c->GetZ());
1245
1246   return 0;
1247 }
1248
1249 void  AliITStrackerMI::AliITSlayer::SortClusters()
1250 {
1251   //
1252   //sort clusters
1253   //
1254   AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1255   Float_t *z                = new Float_t[fN];
1256   Int_t   * index           = new Int_t[fN];
1257   //
1258   for (Int_t i=0;i<fN;i++){
1259     z[i] = fClusters[i]->GetZ();
1260   }
1261   TMath::Sort(fN,z,index,kFALSE);
1262   for (Int_t i=0;i<fN;i++){
1263     clusters[i] = fClusters[index[i]];
1264   }
1265   //
1266   for (Int_t i=0;i<fN;i++){
1267     fClusters[i] = clusters[i];
1268     fZ[i]        = fClusters[i]->GetZ();
1269     AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1270     Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1271     if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1272     fY[i] = y;
1273   }
1274   delete[] index;
1275   delete[] z;
1276   delete[] clusters;
1277   //
1278
1279   fYB[0]=10000000;
1280   fYB[1]=-10000000;
1281   for (Int_t i=0;i<fN;i++){
1282     if (fY[i]<fYB[0]) fYB[0]=fY[i];
1283     if (fY[i]>fYB[1]) fYB[1]=fY[i];
1284     fClusterIndex[i] = i;
1285   }
1286   //
1287   // fill slices
1288   fDy5 = (fYB[1]-fYB[0])/5.;
1289   fDy10 = (fYB[1]-fYB[0])/10.;
1290   fDy20 = (fYB[1]-fYB[0])/20.;
1291   for (Int_t i=0;i<6;i++)  fN5[i] =0;
1292   for (Int_t i=0;i<11;i++) fN10[i]=0;
1293   for (Int_t i=0;i<21;i++) fN20[i]=0;
1294   //
1295   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;}
1296   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;}
1297   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;}
1298   //
1299   //
1300   for (Int_t i=0;i<fN;i++)
1301     for (Int_t irot=-1;irot<=1;irot++){
1302       Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1303       // slice 5
1304       for (Int_t slice=0; slice<6;slice++){
1305         if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<kMaxClusterPerLayer5){
1306           fClusters5[slice][fN5[slice]] = fClusters[i];
1307           fY5[slice][fN5[slice]] = curY;
1308           fZ5[slice][fN5[slice]] = fZ[i];
1309           fClusterIndex5[slice][fN5[slice]]=i;
1310           fN5[slice]++;
1311         }
1312       }
1313       // slice 10
1314       for (Int_t slice=0; slice<11;slice++){
1315         if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<kMaxClusterPerLayer10){
1316           fClusters10[slice][fN10[slice]] = fClusters[i];
1317           fY10[slice][fN10[slice]] = curY;
1318           fZ10[slice][fN10[slice]] = fZ[i];
1319           fClusterIndex10[slice][fN10[slice]]=i;
1320           fN10[slice]++;
1321         }
1322       }
1323       // slice 20
1324       for (Int_t slice=0; slice<21;slice++){
1325         if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<kMaxClusterPerLayer20){
1326           fClusters20[slice][fN20[slice]] = fClusters[i];
1327           fY20[slice][fN20[slice]] = curY;
1328           fZ20[slice][fN20[slice]] = fZ[i];
1329           fClusterIndex20[slice][fN20[slice]]=i;
1330           fN20[slice]++;
1331         }
1332       }
1333     }
1334
1335   //
1336   // consistency check
1337   //
1338   for (Int_t i=0;i<fN-1;i++){
1339     if (fZ[i]>fZ[i+1]){
1340       printf("Bugg\n");
1341     }
1342   }
1343   //
1344   for (Int_t slice=0;slice<21;slice++)
1345   for (Int_t i=0;i<fN20[slice]-1;i++){
1346     if (fZ20[slice][i]>fZ20[slice][i+1]){
1347       printf("Bugg\n");
1348     }
1349   }
1350
1351
1352 }
1353
1354
1355 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1356   //--------------------------------------------------------------------
1357   // This function returns the index of the nearest cluster
1358   //--------------------------------------------------------------------
1359   Int_t ncl=0;
1360   const Float_t *zcl;
1361   if (fCurrentSlice<0) {
1362     ncl = fN;
1363     zcl   = fZ;
1364   }
1365   else{
1366     ncl   = fNcs;
1367     zcl   = fZcs;;
1368   }
1369
1370   if (ncl==0) return 0;
1371   Int_t b=0, e=ncl-1, m=(b+e)/2;
1372   for (; b<e; m=(b+e)/2) {
1373     //    if (z > fClusters[m]->GetZ()) b=m+1;
1374     if (z > zcl[m]) b=m+1;
1375     else e=m;
1376   }
1377   return m;
1378 }
1379
1380
1381 void AliITStrackerMI::AliITSlayer::
1382 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1383   //--------------------------------------------------------------------
1384   // This function sets the "window"
1385   //--------------------------------------------------------------------
1386
1387   Double_t circle=2*TMath::Pi()*fR;
1388   fYmin = ymin; fYmax =ymax;
1389   Float_t ymiddle = (fYmin+fYmax)*0.5;
1390   if (ymiddle<fYB[0]) {fYmin+=circle; fYmax+=circle;ymiddle+=circle;}
1391   else{
1392     if (ymiddle>fYB[1]) {fYmin-=circle; fYmax-=circle;ymiddle-=circle;}
1393   }
1394   //
1395   fCurrentSlice =-1;
1396   // defualt take all
1397   fClustersCs = fClusters;
1398   fClusterIndexCs = fClusterIndex;
1399   fYcs  = fY;
1400   fZcs  = fZ;
1401   fNcs  = fN;
1402   //
1403   //is in 20 slice?
1404   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1405     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1406     if (slice<0) slice=0;
1407     if (slice>20) slice=20;
1408     Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1409     if (isOK) {
1410       fCurrentSlice=slice;
1411       fClustersCs = fClusters20[fCurrentSlice];
1412       fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1413       fYcs  = fY20[fCurrentSlice];
1414       fZcs  = fZ20[fCurrentSlice];
1415       fNcs  = fN20[fCurrentSlice];
1416     }
1417   }
1418   //
1419   //is in 10 slice?
1420   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1421     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1422     if (slice<0) slice=0;
1423     if (slice>10) slice=10;
1424     Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1425     if (isOK) {
1426       fCurrentSlice=slice;
1427       fClustersCs = fClusters10[fCurrentSlice];
1428       fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1429       fYcs  = fY10[fCurrentSlice];
1430       fZcs  = fZ10[fCurrentSlice];
1431       fNcs  = fN10[fCurrentSlice];
1432     }
1433   }
1434   //
1435   //is in 5 slice?
1436   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1437     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1438     if (slice<0) slice=0;
1439     if (slice>5) slice=5;
1440     Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1441     if ( isOK){
1442       fCurrentSlice=slice;
1443       fClustersCs = fClusters5[fCurrentSlice];
1444       fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1445       fYcs  = fY5[fCurrentSlice];
1446       fZcs  = fZ5[fCurrentSlice];
1447       fNcs  = fN5[fCurrentSlice];
1448     }
1449   }
1450   //
1451   fI=FindClusterIndex(zmin); fZmax=zmax;
1452   fImax = TMath::Min(FindClusterIndex(zmax)+1,fNcs);
1453   fSkip = 0;
1454   fAccepted =0;
1455 }
1456
1457
1458
1459
1460 Int_t AliITStrackerMI::AliITSlayer::
1461 FindDetectorIndex(Double_t phi, Double_t z) const {
1462   //--------------------------------------------------------------------
1463   //This function finds the detector crossed by the track
1464   //--------------------------------------------------------------------
1465   Double_t dphi=-(phi-fPhiOffset);
1466   if      (dphi <  0) dphi += 2*TMath::Pi();
1467   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1468   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1469   if (np>=fNladders) np-=fNladders;
1470   if (np<0)          np+=fNladders;
1471
1472   Double_t dz=fZOffset-z;
1473   Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
1474   if (nz>=fNdetectors) return -1;
1475   if (nz<0)            return -1;
1476
1477   return np*fNdetectors + nz;
1478 }
1479
1480
1481 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci){
1482   //--------------------------------------------------------------------
1483   // This function returns clusters within the "window"
1484   //--------------------------------------------------------------------
1485
1486   if (fCurrentSlice<0){
1487     Double_t rpi2 = 2.*fR*TMath::Pi();
1488     for (Int_t i=fI; i<fImax; i++) {
1489       Double_t y = fY[i];
1490       if (fYmax<y) y -= rpi2;
1491       if (fYmin>y) y += rpi2;
1492       if (y<fYmin) continue;
1493       if (y>fYmax) continue;
1494       if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1495       ci=i;
1496       fI=i+1;
1497       return fClusters[i];
1498     }
1499   }
1500   else{
1501     for (Int_t i=fI; i<fImax; i++) {
1502       if (fYcs[i]<fYmin) continue;
1503       if (fYcs[i]>fYmax) continue;
1504       if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1505       ci=fClusterIndexCs[i];
1506       fI=i+1;
1507       return fClustersCs[i];
1508     }
1509   }
1510   return 0;
1511 }
1512
1513
1514
1515 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1516 const {
1517   //--------------------------------------------------------------------
1518   //This function returns the layer thickness at this point (units X0)
1519   //--------------------------------------------------------------------
1520   Double_t d=0.0085;
1521   x0=21.82;
1522   if (43<fR&&fR<45) { //SSD2
1523      Double_t dd=0.0034;
1524      d=dd;
1525      if (TMath::Abs(y-0.00)>3.40) d+=dd;
1526      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1527      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1528      for (Int_t i=0; i<12; i++) {
1529        if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1530           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1531           d+=0.0034;
1532           break;
1533        }
1534        if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
1535           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1536           d+=0.0034;
1537           break;
1538        }
1539        if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1540        if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1541      }
1542   } else
1543   if (37<fR&&fR<41) { //SSD1
1544      Double_t dd=0.0034;
1545      d=dd;
1546      if (TMath::Abs(y-0.00)>3.40) d+=dd;
1547      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1548      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1549      for (Int_t i=0; i<11; i++) {
1550        if (TMath::Abs(z-3.9*i)<0.15) {
1551           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1552           d+=dd;
1553           break;
1554        }
1555        if (TMath::Abs(z+3.9*i)<0.15) {
1556           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1557           d+=dd;
1558           break;
1559        }
1560        if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1561        if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1562      }
1563   } else
1564   if (13<fR&&fR<26) { //SDD
1565      Double_t dd=0.0033;
1566      d=dd;
1567      if (TMath::Abs(y-0.00)>3.30) d+=dd;
1568
1569      if (TMath::Abs(y-1.80)<0.55) {
1570         d+=0.016;
1571         for (Int_t j=0; j<20; j++) {
1572           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1573           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1574         }
1575      }
1576      if (TMath::Abs(y+1.80)<0.55) {
1577         d+=0.016;
1578         for (Int_t j=0; j<20; j++) {
1579           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1580           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1581         }
1582      }
1583
1584      for (Int_t i=0; i<4; i++) {
1585        if (TMath::Abs(z-7.3*i)<0.60) {
1586           d+=dd;
1587           if (TMath::Abs(y-0.00)>3.30) d+=dd;
1588           break;
1589        }
1590        if (TMath::Abs(z+7.3*i)<0.60) {
1591           d+=dd;
1592           if (TMath::Abs(y-0.00)>3.30) d+=dd;
1593           break;
1594        }
1595      }
1596   } else
1597   if (6<fR&&fR<8) {   //SPD2
1598      Double_t dd=0.0063; x0=21.5;
1599      d=dd;
1600      if (TMath::Abs(y-3.08)>0.5) d+=dd;
1601      //if (TMath::Abs(y-3.08)>0.45) d+=dd;
1602      if (TMath::Abs(y-3.03)<0.10) {d+=0.014;}
1603   } else
1604   if (3<fR&&fR<5) {   //SPD1
1605      Double_t dd=0.0063; x0=21.5;
1606      d=dd;
1607      if (TMath::Abs(y+0.21)>0.6) d+=dd;
1608      //if (TMath::Abs(y+0.21)>0.45) d+=dd;
1609      if (TMath::Abs(y+0.10)<0.10) {d+=0.014;}
1610   }
1611
1612   return d;
1613 }
1614
1615 Double_t AliITStrackerMI::GetEffectiveThickness(Double_t y,Double_t z) const
1616 {
1617   //--------------------------------------------------------------------
1618   //Returns the thickness between the current layer and the vertex (units X0)
1619   //--------------------------------------------------------------------
1620   Double_t d=0.0028*3*3; //beam pipe
1621   Double_t x0=0;
1622
1623   Double_t xn=fgLayers[fI].GetR();
1624   for (Int_t i=0; i<fI; i++) {
1625     Double_t xi=fgLayers[i].GetR();
1626     d+=fgLayers[i].GetThickness(y,z,x0)*xi*xi;
1627   }
1628
1629   if (fI>1) {
1630     Double_t xi=9.;
1631     d+=0.0097*xi*xi;
1632   }
1633
1634   if (fI>3) {
1635     Double_t xi=0.5*(fgLayers[3].GetR()+fgLayers[4].GetR());
1636     d+=0.0034*xi*xi;
1637   }
1638
1639   return d/(xn*xn);
1640 }
1641
1642 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
1643   //--------------------------------------------------------------------
1644   // This function returns number of clusters within the "window"
1645   //--------------------------------------------------------------------
1646   Int_t ncl=0;
1647   for (Int_t i=fI; i<fN; i++) {
1648     const AliITSRecPoint *c=fClusters[i];
1649     if (c->GetZ() > fZmax) break;
1650     if (c->IsUsed()) continue;
1651     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
1652     Double_t y=fR*det.GetPhi() + c->GetY();
1653
1654     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
1655     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
1656
1657     if (y<fYmin) continue;
1658     if (y>fYmax) continue;
1659     ncl++;
1660   }
1661   return ncl;
1662 }
1663
1664 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *t,
1665                                 const AliITStrackMI *c, Bool_t extra) {
1666   //--------------------------------------------------------------------
1667   // This function refits the track "t" at the position "x" using
1668   // the clusters from "c"
1669   // If "extra"==kTRUE,
1670   //    the clusters from overlapped modules get attached to "t"
1671   //--------------------------------------------------------------------
1672   Int_t index[kMaxLayer];
1673   Int_t k;
1674   for (k=0; k<kMaxLayer; k++) index[k]=-1;
1675   Int_t nc=c->GetNumberOfClusters();
1676   for (k=0; k<nc; k++) {
1677     Int_t idx=c->GetClusterIndex(k),nl=(idx&0xf0000000)>>28;
1678     index[nl]=idx;
1679   }
1680
1681   Int_t from, to, step;
1682   if (xx > t->GetX()) {
1683       from=0; to=kMaxLayer;
1684       step=+1;
1685   } else {
1686       from=kMaxLayer-1; to=-1;
1687       step=-1;
1688   }
1689
1690   for (Int_t i=from; i != to; i += step) {
1691      AliITSlayer &layer=fgLayers[i];
1692      Double_t r=layer.GetR();
1693
1694      {
1695      Double_t hI=i-0.5*step;
1696      if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) {
1697         Double_t rs=0.5*(fgLayers[i-step].GetR() + r);
1698         Double_t d=0.0034, x0=38.6;
1699         if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;}
1700         if (!t->PropagateTo(rs,-step*d,x0)) {
1701           return kFALSE;
1702         }
1703      }
1704      }
1705
1706      // remember old position [SR, GSI 18.02.2003]
1707      Double_t oldX=0., oldY=0., oldZ=0.;
1708      if (t->IsStartedTimeIntegral() && step==1) {
1709         t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ);
1710      }
1711      //
1712
1713      Double_t x,y,z;
1714      if (!t->GetGlobalXYZat(r,x,y,z)) {
1715        return kFALSE;
1716      }
1717      Double_t phi=TMath::ATan2(y,x);
1718      Int_t idet=layer.FindDetectorIndex(phi,z);
1719      if (idet<0) {
1720        return kFALSE;
1721      }
1722      const AliITSdetector &det=layer.GetDetector(idet);
1723      phi=det.GetPhi();
1724      if (!t->Propagate(phi,det.GetR())) {
1725        return kFALSE;
1726      }
1727      t->SetDetectorIndex(idet);
1728
1729      const AliITSRecPoint *cl=0;
1730      Double_t maxchi2=1000.*kMaxChi2;
1731
1732      Int_t idx=index[i];
1733      if (idx>0) {
1734         const AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(idx);
1735         if (c){
1736           if (idet != c->GetDetectorIndex()) {
1737             idet=c->GetDetectorIndex();
1738             const AliITSdetector &det=layer.GetDetector(idet);
1739             if (!t->Propagate(det.GetPhi(),det.GetR())) {
1740               return kFALSE;
1741             }
1742             t->SetDetectorIndex(idet);
1743           }
1744           //Double_t chi2=t->GetPredictedChi2(c);
1745           Int_t layer = (idx & 0xf0000000) >> 28;;
1746           Double_t chi2=GetPredictedChi2MI(t,c,layer);
1747           if (chi2<maxchi2) {
1748             cl=c;
1749             maxchi2=chi2;
1750           } else {
1751             return kFALSE;
1752           }
1753         }
1754      }
1755
1756      if (cl) {
1757        //if (!t->Update(cl,maxchi2,idx)) {
1758        if (!UpdateMI(t,cl,maxchi2,idx)) {
1759           return kFALSE;
1760        }
1761        t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1);
1762      }
1763
1764      {
1765      Double_t x0;
1766      Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0);
1767      t->CorrectForMaterial(-step*d,x0);
1768      }
1769
1770      if (extra) { //search for extra clusters
1771         AliITStrackV2 tmp(*t);
1772         Double_t dz=4*TMath::Sqrt(tmp.GetSigmaZ2()+kSigmaZ2[i]);
1773         if (dz < 0.5*TMath::Abs(tmp.GetTgl())) dz=0.5*TMath::Abs(tmp.GetTgl());
1774         Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
1775         if (dy < 0.5*TMath::Abs(tmp.GetSnp())) dy=0.5*TMath::Abs(tmp.GetSnp());
1776         Double_t zmin=t->GetZ() - dz;
1777         Double_t zmax=t->GetZ() + dz;
1778         Double_t ymin=t->GetY() + phi*r - dy;
1779         Double_t ymax=t->GetY() + phi*r + dy;
1780         layer.SelectClusters(zmin,zmax,ymin,ymax);
1781
1782         const AliITSRecPoint *c=0; Int_t ci=-1,cci=-1;
1783         Double_t maxchi2=1000.*kMaxChi2, tolerance=0.1;
1784         while ((c=layer.GetNextCluster(ci))!=0) {
1785            if (idet == c->GetDetectorIndex()) continue;
1786
1787            const AliITSdetector &det=layer.GetDetector(c->GetDetectorIndex());
1788
1789            if (!tmp.Propagate(det.GetPhi(),det.GetR())) continue;
1790
1791            if (TMath::Abs(tmp.GetZ() - c->GetZ()) > tolerance) continue;
1792            if (TMath::Abs(tmp.GetY() - c->GetY()) > tolerance) continue;
1793
1794            Double_t chi2=tmp.GetPredictedChi2(c);
1795            if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
1796         }
1797         if (cci>=0) t->SetExtraCluster(i,(i<<28)+cci);
1798      }
1799
1800      // track time update [SR, GSI 17.02.2003]
1801      if (t->IsStartedTimeIntegral() && step==1) {
1802         Double_t newX, newY, newZ;
1803         t->GetGlobalXYZat(t->GetX(),newX,newY,newZ);
1804         Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
1805                        (oldZ-newZ)*(oldZ-newZ);
1806         t->AddTimeStep(TMath::Sqrt(dL2));
1807      }
1808      //
1809
1810   }
1811
1812   if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
1813   return kTRUE;
1814 }
1815
1816 Bool_t
1817 AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *t,const Int_t *clindex) {
1818   //--------------------------------------------------------------------
1819   // This function refits the track "t" at the position "x" using
1820   // the clusters from array
1821   //--------------------------------------------------------------------
1822   Int_t index[kMaxLayer];
1823   Int_t k;
1824   for (k=0; k<kMaxLayer; k++) index[k]=-1;
1825   //
1826   for (k=0; k<kMaxLayer; k++) {
1827     index[k]=clindex[k];
1828   }
1829
1830   Int_t from, to, step;
1831   if (xx > t->GetX()) {
1832       from=0; to=kMaxLayer;
1833       step=+1;
1834   } else {
1835       from=kMaxLayer-1; to=-1;
1836       step=-1;
1837   }
1838
1839   for (Int_t i=from; i != to; i += step) {
1840      AliITSlayer &layer=fgLayers[i];
1841      Double_t r=layer.GetR();
1842      if (step<0 && xx>r) break;  //
1843      {
1844      Double_t hI=i-0.5*step;
1845      if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) {
1846         Double_t rs=0.5*(fgLayers[i-step].GetR() + r);
1847         Double_t d=0.0034, x0=38.6;
1848         if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;}
1849         if (!t->PropagateTo(rs,-step*d,x0)) {
1850           return kFALSE;
1851         }
1852      }
1853      }
1854
1855      // remember old position [SR, GSI 18.02.2003]
1856      Double_t oldX=0., oldY=0., oldZ=0.;
1857      if (t->IsStartedTimeIntegral() && step==1) {
1858         t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ);
1859      }
1860      //
1861
1862      Double_t x,y,z;
1863      if (!t->GetGlobalXYZat(r,x,y,z)) {
1864        return kFALSE;
1865      }
1866      Double_t phi=TMath::ATan2(y,x);
1867      Int_t idet=layer.FindDetectorIndex(phi,z);
1868      if (idet<0) {
1869        return kFALSE;
1870      }
1871      const AliITSdetector &det=layer.GetDetector(idet);
1872      phi=det.GetPhi();
1873      if (!t->Propagate(phi,det.GetR())) {
1874        return kFALSE;
1875      }
1876      t->SetDetectorIndex(idet);
1877
1878      const AliITSRecPoint *cl=0;
1879      Double_t maxchi2=1000.*kMaxChi2;
1880
1881      Int_t idx=index[i];
1882      if (idx>0) {
1883         const AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(idx);
1884         if (c){
1885           if (idet != c->GetDetectorIndex()) {
1886             idet=c->GetDetectorIndex();
1887             const AliITSdetector &det=layer.GetDetector(idet);
1888             if (!t->Propagate(det.GetPhi(),det.GetR())) {
1889               return kFALSE;
1890             }
1891             t->SetDetectorIndex(idet);
1892           }
1893           //Double_t chi2=t->GetPredictedChi2(c);
1894           Int_t layer = (idx & 0xf0000000) >> 28;;
1895           Double_t chi2=GetPredictedChi2MI(t,c,layer);
1896           if (chi2<maxchi2) {
1897             cl=c;
1898             maxchi2=chi2;
1899           } else {
1900             return kFALSE;
1901           }
1902         }
1903      }
1904      /*
1905      if (cl==0)
1906      if (t->GetNumberOfClusters()>2) {
1907         Double_t dz=4*TMath::Sqrt(t->GetSigmaZ2()+kSigmaZ2[i]);
1908         Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
1909         Double_t zmin=t->GetZ() - dz;
1910         Double_t zmax=t->GetZ() + dz;
1911         Double_t ymin=t->GetY() + phi*r - dy;
1912         Double_t ymax=t->GetY() + phi*r + dy;
1913         layer.SelectClusters(zmin,zmax,ymin,ymax);
1914
1915         const AliITSRecPoint *c=0; Int_t ci=-1;
1916         while ((c=layer.GetNextCluster(ci))!=0) {
1917            if (idet != c->GetDetectorIndex()) continue;
1918            Double_t chi2=t->GetPredictedChi2(c);
1919            if (chi2<maxchi2) { cl=c; maxchi2=chi2; idx=ci; }
1920         }
1921      }
1922      */
1923      if (cl) {
1924        //if (!t->Update(cl,maxchi2,idx)) {
1925        if (!UpdateMI(t,cl,maxchi2,idx)) {
1926           return kFALSE;
1927        }
1928        t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1);
1929      }
1930
1931      {
1932      Double_t x0;
1933      Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0);
1934      t->CorrectForMaterial(-step*d,x0);
1935      }
1936
1937      // track time update [SR, GSI 17.02.2003]
1938      if (t->IsStartedTimeIntegral() && step==1) {
1939         Double_t newX, newY, newZ;
1940         t->GetGlobalXYZat(t->GetX(),newX,newY,newZ);
1941         Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
1942                        (oldZ-newZ)*(oldZ-newZ);
1943         t->AddTimeStep(TMath::Sqrt(dL2));
1944      }
1945      //
1946
1947   }
1948
1949   if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
1950   return kTRUE;
1951 }
1952
1953 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
1954 {
1955   //
1956   // calculate normalized chi2
1957   //  return NormalizedChi2(track,0);
1958   Float_t chi2 = 0;
1959   Float_t sum=0;
1960   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
1961   //  track->fdEdxMismatch=0;
1962   Float_t dedxmismatch =0;
1963   Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
1964   if (mode<100){
1965     for (Int_t i = 0;i<6;i++){
1966       if (track->GetClIndex(i)>0){
1967         Float_t cerry, cerrz;
1968         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
1969         else
1970           { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
1971         cerry*=cerry;
1972         cerrz*=cerrz;
1973         Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
1974         if (i>1){
1975           Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
1976           if (ratio<0.5) {
1977             cchi2+=(0.5-ratio)*10.;
1978             //track->fdEdxMismatch+=(0.5-ratio)*10.;
1979             dedxmismatch+=(0.5-ratio)*10.;
1980           }
1981         }
1982         if (i<2 ||i>3){
1983           AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
1984           Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
1985           if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
1986           if (i<2) chi2+=2*cl->GetDeltaProbability();
1987         }
1988         chi2+=cchi2;
1989         sum++;
1990       }
1991     }
1992     if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
1993       track->SetdEdxMismatch(dedxmismatch);
1994     }
1995   }
1996   else{
1997     for (Int_t i = 0;i<4;i++){
1998       if (track->GetClIndex(i)>0){
1999         Float_t cerry, cerrz;
2000         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2001         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2002         cerry*=cerry;
2003         cerrz*=cerrz;
2004         chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2005         chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2006         sum++;
2007       }
2008     }
2009     for (Int_t i = 4;i<6;i++){
2010       if (track->GetClIndex(i)>0){
2011         Float_t cerry, cerrz;
2012         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2013         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2014         cerry*=cerry;
2015         cerrz*=cerrz;
2016         Float_t cerryb, cerrzb;
2017         if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2018         else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2019         cerryb*=cerryb;
2020         cerrzb*=cerrzb;
2021         chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2022         chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2023         sum++;
2024       }
2025     }
2026   }
2027   if (track->GetESDtrack()->GetTPCsignal()>85){
2028     Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2029     if (ratio<0.5) {
2030       chi2+=(0.5-ratio)*5.;
2031     }
2032     if (ratio>2){
2033       chi2+=(ratio-2.0)*3;
2034     }
2035   }
2036   //
2037   Double_t match = TMath::Sqrt(track->GetChi22());
2038   if (track->GetConstrain())  match/=track->GetNumberOfClusters();
2039   if (!track->GetConstrain()) match/=track->GetNumberOfClusters()-2.;
2040   if (match<0) match=0;
2041   Float_t deadzonefactor = (track->GetNDeadZone()>0) ? 3*(1.1-track->GetDeadZoneProbability()):0.;
2042   Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2043     (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2044                                 1./(1.+track->GetNSkipped()));
2045
2046  return normchi2;
2047 }
2048
2049
2050 Double_t AliITStrackerMI::GetMatchingChi2(AliITStrackMI * track1, AliITStrackMI * track2)
2051 {
2052   //
2053   // return matching chi2 between two tracks
2054   AliITStrackMI track3(*track2);
2055   track3.Propagate(track1->GetAlpha(),track1->GetX());
2056   TMatrixD vec(5,1);
2057   vec(0,0)=track1->GetY()   - track3.GetY();
2058   vec(1,0)=track1->GetZ()   - track3.GetZ();
2059   vec(2,0)=track1->GetSnp() - track3.GetSnp();
2060   vec(3,0)=track1->GetTgl() - track3.GetTgl();
2061   vec(4,0)=track1->Get1Pt() - track3.Get1Pt();
2062   //
2063   TMatrixD cov(5,5);
2064   cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2065   cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2066   cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2067   cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2068   cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2069
2070   cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2071   cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2072   cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2073   cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2074   //
2075   cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2076   cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2077   cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2078   //
2079   cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2080   cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2081   //
2082   cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2083
2084   cov.Invert();
2085   TMatrixD vec2(cov,TMatrixD::kMult,vec);
2086   TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2087   return chi2(0,0);
2088 }
2089
2090 Double_t  AliITStrackerMI::GetDeadZoneProbability(Double_t zpos, Double_t zerr)
2091 {
2092   //
2093   //  return probability that given point - characterized by z position and error  is in dead zone
2094   //
2095   Double_t probability =0;
2096   Double_t absz = TMath::Abs(zpos);
2097   Double_t nearestz = (absz<2)? 0.:7.1;
2098   if (TMath::Abs(absz-nearestz)>0.25+3*zerr) return 0;
2099   Double_t zmin=0, zmax=0;
2100   if (zpos<-6.){
2101     zmin = -7.25; zmax = -6.95;
2102   }
2103   if (zpos>6){
2104     zmin = 7.0; zmax =7.3;
2105   }
2106   if (absz<2){
2107     zmin = -0.75; zmax = 1.5;
2108   }
2109   probability = (TMath::Erf((zpos-zmin)/zerr) - TMath::Erf((zpos-zmax)/zerr))*0.5;
2110   return probability;
2111 }
2112
2113
2114 Double_t AliITStrackerMI::GetTruncatedChi2(AliITStrackMI * track, Float_t fac)
2115 {
2116   //
2117   // calculate normalized chi2
2118   Float_t chi2[6];
2119   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2120   Float_t ncl = 0;
2121   for (Int_t i = 0;i<6;i++){
2122     if (TMath::Abs(track->GetDy(i))>0){
2123       chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2124       chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2125       ncl++;
2126     }
2127     else{chi2[i]=10000;}
2128   }
2129   Int_t index[6];
2130   TMath::Sort(6,chi2,index,kFALSE);
2131   Float_t max = float(ncl)*fac-1.;
2132   Float_t sumchi=0, sumweight=0;
2133   for (Int_t i=0;i<max+1;i++){
2134     Float_t weight = (i<max)?1.:(max+1.-i);
2135     sumchi+=weight*chi2[index[i]];
2136     sumweight+=weight;
2137   }
2138   Double_t normchi2 = sumchi/sumweight;
2139   return normchi2;
2140 }
2141
2142
2143 Double_t AliITStrackerMI::GetInterpolatedChi2(AliITStrackMI * forwardtrack, AliITStrackMI * backtrack)
2144 {
2145   //
2146   // calculate normalized chi2
2147   //  if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2148   Int_t npoints = 0;
2149   Double_t res =0;
2150   for (Int_t i=0;i<6;i++){
2151     if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2152     Double_t sy1 = forwardtrack->GetSigmaY(i);
2153     Double_t sz1 = forwardtrack->GetSigmaZ(i);
2154     Double_t sy2 = backtrack->GetSigmaY(i);
2155     Double_t sz2 = backtrack->GetSigmaZ(i);
2156     if (i<2){ sy2=1000.;sz2=1000;}
2157     //
2158     Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2159     Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2160     //
2161     Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2162     Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2163     //
2164     res+= nz0*nz0+ny0*ny0;
2165     npoints++;
2166   }
2167   if (npoints>1) return
2168                    TMath::Max(TMath::Abs(0.3*forwardtrack->Get1Pt())-0.5,0.)+
2169                    //2*forwardtrack->fNUsed+
2170                    res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2171                                   1./(1.+forwardtrack->GetNSkipped()));
2172   return 1000;
2173 }
2174
2175
2176
2177
2178
2179 Float_t  *AliITStrackerMI::GetWeight(Int_t index) {
2180   //--------------------------------------------------------------------
2181   //       Return pointer to a given cluster
2182   //--------------------------------------------------------------------
2183   Int_t l=(index & 0xf0000000) >> 28;
2184   Int_t c=(index & 0x0fffffff) >> 00;
2185   return fgLayers[l].GetWeight(c);
2186 }
2187
2188 void AliITStrackerMI::RegisterClusterTracks(AliITStrackMI* track,Int_t id)
2189 {
2190   //---------------------------------------------
2191   // register track to the list
2192   //
2193   if (track->GetESDtrack()->GetKinkIndex(0)!=0) return;  //don't register kink tracks
2194   //
2195   //
2196   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2197     Int_t index = track->GetClusterIndex(icluster);
2198     Int_t l=(index & 0xf0000000) >> 28;
2199     Int_t c=(index & 0x0fffffff) >> 00;
2200     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2201     for (Int_t itrack=0;itrack<4;itrack++){
2202       if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2203         fgLayers[l].SetClusterTracks(itrack,c,id);
2204         break;
2205       }
2206     }
2207   }
2208 }
2209 void AliITStrackerMI::UnRegisterClusterTracks(AliITStrackMI* track, Int_t id)
2210 {
2211   //---------------------------------------------
2212   // unregister track from the list
2213   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2214     Int_t index = track->GetClusterIndex(icluster);
2215     Int_t l=(index & 0xf0000000) >> 28;
2216     Int_t c=(index & 0x0fffffff) >> 00;
2217     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2218     for (Int_t itrack=0;itrack<4;itrack++){
2219       if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2220         fgLayers[l].SetClusterTracks(itrack,c,-1);
2221       }
2222     }
2223   }
2224 }
2225 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2226 {
2227   //-------------------------------------------------------------
2228   //get number of shared clusters
2229   //-------------------------------------------------------------
2230   Float_t shared=0;
2231   for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2232   // mean  number of clusters
2233   Float_t *ny = GetNy(id), *nz = GetNz(id);
2234
2235
2236   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2237     Int_t index = track->GetClusterIndex(icluster);
2238     Int_t l=(index & 0xf0000000) >> 28;
2239     Int_t c=(index & 0x0fffffff) >> 00;
2240     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2241     if (ny[l]==0){
2242       printf("problem\n");
2243     }
2244     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2245     Float_t weight=1;
2246     //
2247     Float_t deltan = 0;
2248     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2249     if (l>2&&track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2250     if (l<2 || l>3){
2251       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2252     }
2253     else{
2254       deltan = (cl->GetNz()-nz[l]);
2255     }
2256     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2257     weight = 2./TMath::Max(3.+deltan,2.);
2258     //
2259     for (Int_t itrack=0;itrack<4;itrack++){
2260       if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2261         list[l]=index;
2262         clist[l] = (AliITSRecPoint*)GetCluster(index);
2263         shared+=weight;
2264         break;
2265       }
2266     }
2267   }
2268   track->SetNUsed(shared);
2269   return shared;
2270 }
2271
2272 Int_t AliITStrackerMI::GetOverlapTrack(AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2273 {
2274   //
2275   // find first shared track
2276   //
2277   // mean  number of clusters
2278   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2279   //
2280   for (Int_t i=0;i<6;i++) overlist[i]=-1;
2281   Int_t sharedtrack=100000;
2282   Int_t tracks[24],trackindex=0;
2283   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2284   //
2285   for (Int_t icluster=0;icluster<6;icluster++){
2286     if (clusterlist[icluster]<0) continue;
2287     Int_t index = clusterlist[icluster];
2288     Int_t l=(index & 0xf0000000) >> 28;
2289     Int_t c=(index & 0x0fffffff) >> 00;
2290     if (ny[l]==0){
2291       printf("problem\n");
2292     }
2293     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2294     //if (l>3) continue;
2295     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2296     //
2297     Float_t deltan = 0;
2298     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2299     if (l>2&&track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2300     if (l<2 || l>3){
2301       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2302     }
2303     else{
2304       deltan = (cl->GetNz()-nz[l]);
2305     }
2306     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2307     //
2308     for (Int_t itrack=3;itrack>=0;itrack--){
2309       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2310       if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2311        tracks[trackindex]  = fgLayers[l].GetClusterTracks(itrack,c);
2312        trackindex++;
2313       }
2314     }
2315   }
2316   if (trackindex==0) return -1;
2317   if (trackindex==1){
2318     sharedtrack = tracks[0];
2319   }else{
2320     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2321     else{
2322       //
2323       Int_t track[24], cluster[24];
2324       for (Int_t i=0;i<trackindex;i++){ track[i]=-1; cluster[i]=0;}
2325       Int_t index =0;
2326       //
2327       for (Int_t i=0;i<trackindex;i++){
2328         if (tracks[i]<0) continue;
2329         track[index] = tracks[i];
2330         cluster[index]++;
2331         for (Int_t j=i+1;j<trackindex;j++){
2332           if (tracks[j]<0) continue;
2333           if (tracks[j]==tracks[i]){
2334             cluster[index]++;
2335             tracks[j]=-1;
2336           }
2337         }
2338         index++;
2339       }
2340       Int_t max=0;
2341       for (Int_t i=0;i<index;i++){
2342         if (cluster[index]>max) {
2343           sharedtrack=track[index];
2344           max=cluster[index];
2345         }
2346       }
2347     }
2348   }
2349
2350   if (sharedtrack>=100000) return -1;
2351   //
2352   // make list of overlaps
2353   shared =0;
2354   for (Int_t icluster=0;icluster<6;icluster++){
2355     if (clusterlist[icluster]<0) continue;
2356     Int_t index = clusterlist[icluster];
2357     Int_t l=(index & 0xf0000000) >> 28;
2358     Int_t c=(index & 0x0fffffff) >> 00;
2359     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2360     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2361     if (l==0 || l==1){
2362       if (cl->GetNy()>2) continue;
2363       if (cl->GetNz()>2) continue;
2364     }
2365     if (l==4 || l==5){
2366       if (cl->GetNy()>3) continue;
2367       if (cl->GetNz()>3) continue;
2368     }
2369     //
2370     for (Int_t itrack=3;itrack>=0;itrack--){
2371       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2372       if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2373         overlist[l]=index;
2374         shared++;
2375       }
2376     }
2377   }
2378   return sharedtrack;
2379 }
2380
2381
2382 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2383   //
2384   // try to find track hypothesys without conflicts
2385   // with minimal chi2;
2386   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2387   Int_t entries1 = arr1->GetEntriesFast();
2388   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2389   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2390   Int_t entries2 = arr2->GetEntriesFast();
2391   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2392   //
2393   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2394   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2395   if (TMath::Abs(1./track10->Get1Pt())>0.5+TMath::Abs(1/track20->Get1Pt())) return track10;
2396
2397   for (Int_t itrack=0;itrack<entries1;itrack++){
2398     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2399     UnRegisterClusterTracks(track,trackID1);
2400   }
2401   //
2402   for (Int_t itrack=0;itrack<entries2;itrack++){
2403     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2404     UnRegisterClusterTracks(track,trackID2);
2405   }
2406   Int_t index1=0;
2407   Int_t index2=0;
2408   Float_t maxconflicts=6;
2409   Double_t maxchi2 =1000.;
2410   //
2411   // get weight of hypothesys - using best hypothesy
2412   Double_t w1,w2;
2413
2414   Int_t list1[6],list2[6];
2415   AliITSRecPoint *clist1[6], *clist2[6] ;
2416   RegisterClusterTracks(track10,trackID1);
2417   RegisterClusterTracks(track20,trackID2);
2418   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2419   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2420   UnRegisterClusterTracks(track10,trackID1);
2421   UnRegisterClusterTracks(track20,trackID2);
2422   //
2423   // normalized chi2
2424   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2425   Float_t nerry[6],nerrz[6];
2426   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2427   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2428   for (Int_t i=0;i<6;i++){
2429      if ( (erry1[i]>0) && (erry2[i]>0)) {
2430        nerry[i] = TMath::Min(erry1[i],erry2[i]);
2431        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2432      }else{
2433        nerry[i] = TMath::Max(erry1[i],erry2[i]);
2434        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2435      }
2436      if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2437        chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2438        chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2439        ncl1++;
2440      }
2441      if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2442        chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2443        chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2444        ncl2++;
2445      }
2446   }
2447   chi21/=ncl1;
2448   chi22/=ncl2;
2449   //
2450   //
2451   Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2452   Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2453   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2454   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2455   //
2456   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2457         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2458         +1.*TMath::Abs(1./track10->Get1Pt())/(TMath::Abs(1./track10->Get1Pt())+TMath::Abs(1./track20->Get1Pt()))
2459         );
2460   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2461         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2462         +1.*TMath::Abs(1./track20->Get1Pt())/(TMath::Abs(1./track10->Get1Pt())+TMath::Abs(1./track20->Get1Pt()))
2463         );
2464
2465   Double_t sumw = w1+w2;
2466   w1/=sumw;
2467   w2/=sumw;
2468   if (w1<w2*0.5) {
2469     w1 = (d2+0.5)/(d1+d2+1.);
2470     w2 = (d1+0.5)/(d1+d2+1.);
2471   }
2472   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2473   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2474   //
2475   // get pair of "best" hypothesys
2476   //
2477   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2478   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2479
2480   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2481     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2482     //if (track1->fFakeRatio>0) continue;
2483     RegisterClusterTracks(track1,trackID1);
2484     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2485       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2486
2487       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2488       //if (track2->fFakeRatio>0) continue;
2489       Float_t nskipped=0;
2490       RegisterClusterTracks(track2,trackID2);
2491       Int_t list1[6],list2[6];
2492       AliITSRecPoint *clist1[6], *clist2[6] ;
2493       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2494       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2495       UnRegisterClusterTracks(track2,trackID2);
2496       //
2497       if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2498       if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2499       if (nskipped>0.5) continue;
2500       //
2501       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2502       if (conflict1+1<cconflict1) continue;
2503       if (conflict2+1<cconflict2) continue;
2504       Float_t conflict=0;
2505       Float_t sumchi2=0;
2506       Float_t sum=0;
2507       for (Int_t i=0;i<6;i++){
2508         //
2509         Float_t c1 =0.; // conflict coeficients
2510         Float_t c2 =0.;
2511         if (clist1[i]&&clist2[i]){
2512           Float_t deltan = 0;
2513           if (i<2 || i>3){
2514             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2515           }
2516           else{
2517             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2518           }
2519           c1  = 2./TMath::Max(3.+deltan,2.);
2520           c2  = 2./TMath::Max(3.+deltan,2.);
2521         }
2522         else{
2523           if (clist1[i]){
2524             Float_t deltan = 0;
2525             if (i<2 || i>3){
2526               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
2527             }
2528             else{
2529               deltan = (clist1[i]->GetNz()-nz1[i]);
2530             }
2531             c1  = 2./TMath::Max(3.+deltan,2.);
2532             c2  = 0;
2533           }
2534
2535           if (clist2[i]){
2536             Float_t deltan = 0;
2537             if (i<2 || i>3){
2538               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
2539             }
2540             else{
2541               deltan = (clist2[i]->GetNz()-nz2[i]);
2542             }
2543             c2  = 2./TMath::Max(3.+deltan,2.);
2544             c1  = 0;
2545           }
2546         }
2547         //
2548         Double_t chi21=0,chi22=0;
2549         if (TMath::Abs(track1->GetDy(i))>0.) {
2550           chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
2551             (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
2552           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
2553           //  (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
2554         }else{
2555           if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
2556         }
2557         //
2558         if (TMath::Abs(track2->GetDy(i))>0.) {
2559           chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
2560             (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
2561           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
2562           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
2563         }
2564         else{
2565           if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
2566         }
2567         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
2568         if (chi21>0) sum+=w1;
2569         if (chi22>0) sum+=w2;
2570         conflict+=(c1+c2);
2571       }
2572       Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
2573       if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
2574       Double_t normchi2 = 2*conflict+sumchi2/sum;
2575       if ( normchi2 <maxchi2 ){
2576         index1 = itrack1;
2577         index2 = itrack2;
2578         maxconflicts = conflict;
2579         maxchi2 = normchi2;
2580       }
2581     }
2582     UnRegisterClusterTracks(track1,trackID1);
2583   }
2584   //
2585   //  if (maxconflicts<4 && maxchi2<th0){
2586   if (maxchi2<th0*2.){
2587     Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
2588     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
2589     track1->SetChi2MIP(5,maxconflicts);
2590     track1->SetChi2MIP(6,maxchi2);
2591     track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
2592     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
2593     track1->SetChi2MIP(8,index1);
2594     fBestTrackIndex[trackID1] =index1;
2595     UpdateESDtrack(track1, AliESDtrack::kITSin);
2596   }
2597   else if (track10->GetChi2MIP(0)<th1){
2598     track10->SetChi2MIP(5,maxconflicts);
2599     track10->SetChi2MIP(6,maxchi2);
2600     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
2601     UpdateESDtrack(track10,AliESDtrack::kITSin);
2602   }
2603
2604   for (Int_t itrack=0;itrack<entries1;itrack++){
2605     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2606     UnRegisterClusterTracks(track,trackID1);
2607   }
2608   //
2609   for (Int_t itrack=0;itrack<entries2;itrack++){
2610     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2611     UnRegisterClusterTracks(track,trackID2);
2612   }
2613
2614   if (track10->GetConstrain()&&track10->GetChi2MIP(0)<kMaxChi2PerCluster[0]&&track10->GetChi2MIP(1)<kMaxChi2PerCluster[1]
2615       &&track10->GetChi2MIP(2)<kMaxChi2PerCluster[2]&&track10->GetChi2MIP(3)<kMaxChi2PerCluster[3]){
2616     //  if (track10->fChi2MIP[0]<kMaxChi2PerCluster[0]&&track10->fChi2MIP[1]<kMaxChi2PerCluster[1]
2617   //    &&track10->fChi2MIP[2]<kMaxChi2PerCluster[2]&&track10->fChi2MIP[3]<kMaxChi2PerCluster[3]){
2618     RegisterClusterTracks(track10,trackID1);
2619   }
2620   if (track20->GetConstrain()&&track20->GetChi2MIP(0)<kMaxChi2PerCluster[0]&&track20->GetChi2MIP(1)<kMaxChi2PerCluster[1]
2621       &&track20->GetChi2MIP(2)<kMaxChi2PerCluster[2]&&track20->GetChi2MIP(3)<kMaxChi2PerCluster[3]){
2622     //if (track20->fChi2MIP[0]<kMaxChi2PerCluster[0]&&track20->fChi2MIP[1]<kMaxChi2PerCluster[1]
2623     //  &&track20->fChi2MIP[2]<kMaxChi2PerCluster[2]&&track20->fChi2MIP[3]<kMaxChi2PerCluster[3]){
2624     RegisterClusterTracks(track20,trackID2);
2625   }
2626   return track10;
2627
2628 }
2629
2630 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
2631   //--------------------------------------------------------------------
2632   // This function marks clusters assigned to the track
2633   //--------------------------------------------------------------------
2634   AliTracker::UseClusters(t,from);
2635
2636   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
2637   //if (c->GetQ()>2) c->Use();
2638   if (c->GetSigmaZ2()>0.1) c->Use();
2639   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
2640   //if (c->GetQ()>2) c->Use();
2641   if (c->GetSigmaZ2()>0.1) c->Use();
2642
2643 }
2644
2645
2646 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
2647 {
2648   //------------------------------------------------------------------
2649   // add track to the list of hypothesys
2650   //------------------------------------------------------------------
2651
2652   if (esdindex>=fTrackHypothesys.GetEntriesFast()) fTrackHypothesys.Expand(esdindex*2+10);
2653   //
2654   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2655   if (!array) {
2656     array = new TObjArray(10);
2657     fTrackHypothesys.AddAt(array,esdindex);
2658   }
2659   array->AddLast(track);
2660 }
2661
2662 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
2663 {
2664   //-------------------------------------------------------------------
2665   // compress array of track hypothesys
2666   // keep only maxsize best hypothesys
2667   //-------------------------------------------------------------------
2668   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
2669   if (! (fTrackHypothesys.At(esdindex)) ) return;
2670   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2671   Int_t entries = array->GetEntriesFast();
2672   //
2673   //- find preliminary besttrack as a reference
2674   Float_t minchi2=10000;
2675   Int_t maxn=0;
2676   AliITStrackMI * besttrack=0;
2677   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
2678     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2679     if (!track) continue;
2680     Float_t chi2 = NormalizedChi2(track,0);
2681     //
2682     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
2683     track->SetLabel(tpcLabel);
2684     CookdEdx(track);
2685     track->SetFakeRatio(1.);
2686     CookLabel(track,0.); //For comparison only
2687     //
2688     //if (chi2<kMaxChi2PerCluster[0]&&track->fFakeRatio==0){
2689     if (chi2<kMaxChi2PerCluster[0]){
2690       if (track->GetNumberOfClusters()<maxn) continue;
2691       maxn = track->GetNumberOfClusters();
2692       if (chi2<minchi2){
2693         minchi2=chi2;
2694         besttrack=track;
2695       }
2696     }
2697     else{
2698       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
2699         delete array->RemoveAt(itrack);
2700       }
2701     }
2702   }
2703   if (!besttrack) return;
2704   //
2705   //
2706   //take errors of best track as a reference
2707   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
2708   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
2709   for (Int_t i=0;i<6;i++) {
2710     if (besttrack->GetClIndex(i)>0){
2711       erry[i] = besttrack->GetSigmaY(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2712       errz[i] = besttrack->GetSigmaZ(i); errz[i+6] = besttrack->GetSigmaZ(i+6);
2713       ny[i]   = besttrack->GetNy(i);
2714       nz[i]   = besttrack->GetNz(i);
2715     }
2716   }
2717   //
2718   // calculate normalized chi2
2719   //
2720   Float_t * chi2        = new Float_t[entries];
2721   Int_t * index         = new Int_t[entries];
2722   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
2723   for (Int_t itrack=0;itrack<entries;itrack++){
2724     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2725     if (track){
2726       track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
2727       if (track->GetChi2MIP(0)<kMaxChi2PerCluster[0])
2728         chi2[itrack] = track->GetChi2MIP(0);
2729       else{
2730         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
2731           delete array->RemoveAt(itrack);
2732         }
2733       }
2734     }
2735   }
2736   //
2737   TMath::Sort(entries,chi2,index,kFALSE);
2738   besttrack = (AliITStrackMI*)array->At(index[0]);
2739   if (besttrack&&besttrack->GetChi2MIP(0)<kMaxChi2PerCluster[0]){
2740     for (Int_t i=0;i<6;i++){
2741       if (besttrack->GetClIndex(i)>0){
2742         erry[i] = besttrack->GetSigmaY(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2743         errz[i] = besttrack->GetSigmaZ(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2744         ny[i]   = besttrack->GetNy(i);
2745         nz[i]   = besttrack->GetNz(i);
2746       }
2747     }
2748   }
2749   //
2750   // calculate one more time with updated normalized errors
2751   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
2752   for (Int_t itrack=0;itrack<entries;itrack++){
2753     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2754     if (track){
2755       track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
2756       if (track->GetChi2MIP(0)<kMaxChi2PerCluster[0])
2757         chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
2758       else
2759         {
2760           if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
2761             delete array->RemoveAt(itrack);
2762           }
2763         }
2764     }
2765   }
2766   entries = array->GetEntriesFast();
2767   //
2768   //
2769   if (entries>0){
2770     TObjArray * newarray = new TObjArray();
2771     TMath::Sort(entries,chi2,index,kFALSE);
2772     besttrack = (AliITStrackMI*)array->At(index[0]);
2773     if (besttrack){
2774       //
2775       for (Int_t i=0;i<6;i++){
2776         if (besttrack->GetNz(i)>0&&besttrack->GetNy(i)>0){
2777           erry[i] = besttrack->GetSigmaY(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2778           errz[i] = besttrack->GetSigmaZ(i); errz[i+6] = besttrack->GetSigmaZ(i+6);
2779           ny[i]   = besttrack->GetNy(i);
2780           nz[i]   = besttrack->GetNz(i);
2781         }
2782       }
2783       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
2784       Float_t minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(kMaxChi2PerCluster[0]));
2785       Float_t minn = besttrack->GetNumberOfClusters()-3;
2786       Int_t accepted=0;
2787       for (Int_t i=0;i<entries;i++){
2788         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
2789         if (!track) continue;
2790         if (accepted>maxcut) break;
2791         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
2792         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
2793           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
2794             delete array->RemoveAt(index[i]);
2795             continue;
2796           }
2797         }
2798         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
2799         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
2800           if (!shortbest) accepted++;
2801           //
2802           newarray->AddLast(array->RemoveAt(index[i]));
2803           for (Int_t i=0;i<6;i++){
2804             if (nz[i]==0){
2805               erry[i] = track->GetSigmaY(i); erry[i+6] = track->GetSigmaY(i+6);
2806               errz[i] = track->GetSigmaZ(i); errz[i]   = track->GetSigmaZ(i+6);
2807               ny[i]   = track->GetNy(i);
2808               nz[i]   = track->GetNz(i);
2809             }
2810           }
2811         }
2812         else{
2813           delete array->RemoveAt(index[i]);
2814         }
2815       }
2816       array->Delete();
2817       delete fTrackHypothesys.RemoveAt(esdindex);
2818       fTrackHypothesys.AddAt(newarray,esdindex);
2819     }
2820     else{
2821       array->Delete();
2822       delete fTrackHypothesys.RemoveAt(esdindex);
2823     }
2824   }
2825   delete [] chi2;
2826   delete [] index;
2827 }
2828
2829
2830
2831 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
2832 {
2833   //-------------------------------------------------------------
2834   // try to find best hypothesy
2835   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
2836   //-------------------------------------------------------------
2837   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
2838   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2839   if (!array) return 0;
2840   Int_t entries = array->GetEntriesFast();
2841   if (!entries) return 0;
2842   Float_t minchi2 = 100000;
2843   AliITStrackMI * besttrack=0;
2844   //
2845   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
2846   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
2847   Double_t xyzv[]={GetX(),GetY(),GetZ()};
2848   Double_t ersv[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
2849   //
2850   for (Int_t i=0;i<entries;i++){
2851     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
2852     if (!track) continue;
2853     Float_t sigmarfi,sigmaz;
2854     GetDCASigma(track,sigmarfi,sigmaz);
2855     track->SetDnorm(0,sigmarfi);
2856     track->SetDnorm(1,sigmaz);
2857     //
2858     track->SetChi2MIP(1,1000000);
2859     track->SetChi2MIP(2,1000000);
2860     track->SetChi2MIP(3,1000000);
2861     //
2862     // backtrack
2863     backtrack = new(backtrack) AliITStrackMI(*track);
2864     if (track->GetConstrain()){
2865       if (!backtrack->PropagateTo(3.,0.0028,65.19)) continue;
2866       if (!backtrack->Improve(0,xyzv,ersv))         continue;
2867       //if (!backtrack->PropagateTo(2.,0.0028,0))     continue; // This
2868       //if (!backtrack->Improve(0,xyzv,ersv))         continue; // is
2869       //if (!backtrack->PropagateTo(1.,0.0028,0))     continue; // an over-kill
2870       //if (!backtrack->Improve(0,xyzv,ersv))         continue; //   (I.B.)
2871       //if (!backtrack->PropagateToVertex())          continue; //
2872       backtrack->ResetCovariance(10.);
2873       //if (!backtrack->Improve(0,xyzv,ersv))         continue;
2874     }else{
2875       backtrack->ResetCovariance(10.);
2876     }
2877     backtrack->ResetClusters();
2878
2879     Double_t x = original->GetX();
2880     if (!RefitAt(x,backtrack,track)) continue;
2881     //
2882     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
2883     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
2884     if (track->GetChi2MIP(1)>kMaxChi2PerCluster[1]*6.)  continue;
2885     track->SetChi22(GetMatchingChi2(backtrack,original));
2886
2887     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
2888     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
2889     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
2890
2891
2892     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>kMaxChi2PerCluster[1])  continue;
2893     Bool_t isOK=kTRUE;
2894     if(!isOK) continue;
2895     //
2896     //forward track - without constraint
2897     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
2898     forwardtrack->ResetClusters();
2899     x = track->GetX();
2900     RefitAt(x,forwardtrack,track);
2901     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
2902     if  (track->GetChi2MIP(2)>kMaxChi2PerCluster[2]*6.0)  continue;
2903     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>kMaxChi2PerCluster[2])  continue;
2904
2905     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
2906     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
2907     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
2908     forwardtrack->SetD(0,track->GetD(0));
2909     forwardtrack->SetD(1,track->GetD(1));
2910     {
2911       Int_t list[6];
2912       AliITSRecPoint* clist[6];
2913       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
2914       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
2915     }
2916
2917     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
2918     if  ( (track->GetChi2MIP(3)>6.*kMaxChi2PerCluster[3])) continue;
2919     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*kMaxChi2PerCluster[3])) {
2920       track->SetChi2MIP(3,1000);
2921       continue;
2922     }
2923     Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
2924     //
2925     for (Int_t ichi=0;ichi<5;ichi++){
2926       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
2927     }
2928     if (chi2 < minchi2){
2929       //besttrack = new AliITStrackMI(*forwardtrack);
2930       besttrack = track;
2931       besttrack->SetLabel(track->GetLabel());
2932       besttrack->SetFakeRatio(track->GetFakeRatio());
2933       minchi2   = chi2;
2934       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
2935       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
2936       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
2937     }
2938   }
2939   delete backtrack;
2940   delete forwardtrack;
2941   Int_t accepted=0;
2942   for (Int_t i=0;i<entries;i++){
2943     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
2944     if (!track) continue;
2945
2946     if (accepted>checkmax || track->GetChi2MIP(3)>kMaxChi2PerCluster[3]*6. ||
2947         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
2948         track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
2949       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
2950         delete array->RemoveAt(i);
2951         continue;
2952       }
2953     }
2954     else{
2955       accepted++;
2956     }
2957   }
2958   //
2959   array->Compress();
2960   SortTrackHypothesys(esdindex,checkmax,1);
2961   array = (TObjArray*) fTrackHypothesys.At(esdindex);
2962   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
2963   besttrack = (AliITStrackMI*)array->At(0);
2964   if (!besttrack)  return 0;
2965   besttrack->SetChi2MIP(8,0);
2966   fBestTrackIndex[esdindex]=0;
2967   entries = array->GetEntriesFast();
2968   AliITStrackMI *longtrack =0;
2969   minchi2 =1000;
2970   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
2971   for (Int_t itrack=entries-1;itrack>0;itrack--){
2972     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2973     if (!track->GetConstrain()) continue;
2974     if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
2975     if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
2976     if (track->GetChi2MIP(0)>4.) continue;
2977     minn = track->GetNumberOfClusters()+track->GetNDeadZone();
2978     longtrack =track;
2979   }
2980   //if (longtrack) besttrack=longtrack;
2981
2982   Int_t list[6];
2983   AliITSRecPoint * clist[6];
2984   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
2985   if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<kMaxChi2PerCluster[0]&&besttrack->GetChi2MIP(1)<kMaxChi2PerCluster[1]
2986       &&besttrack->GetChi2MIP(2)<kMaxChi2PerCluster[2]&&besttrack->GetChi2MIP(3)<kMaxChi2PerCluster[3]){
2987     RegisterClusterTracks(besttrack,esdindex);
2988   }
2989   //
2990   //
2991   if (shared>0.0){
2992     Int_t nshared;
2993     Int_t overlist[6];
2994     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
2995     if (sharedtrack>=0){
2996       //
2997       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
2998       if (besttrack){
2999         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3000       }
3001       else return 0;
3002     }
3003   }
3004
3005   if (shared>2.5) return 0;
3006   if (shared>1.0) return besttrack;
3007   //
3008   // Don't sign clusters if not gold track
3009   //
3010   if (!besttrack->IsGoldPrimary()) return besttrack;
3011   if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
3012   //
3013   if (fConstraint[fPass]){
3014     //
3015     // sign clusters
3016     //
3017     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3018     for (Int_t i=0;i<6;i++){
3019       Int_t index = besttrack->GetClIndex(i);
3020       if (index<=0) continue;
3021       Int_t ilayer =  (index & 0xf0000000) >> 28;
3022       if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3023       AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3024       if (!c) continue;
3025       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3026       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3027       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
3028       if ( ilayer>2&& besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3029       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
3030
3031       Bool_t cansign = kTRUE;
3032       for (Int_t itrack=0;itrack<entries; itrack++){
3033         AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3034         if (!track) continue;
3035         if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3036         if ( (track->GetClIndex(ilayer)>0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3037           cansign = kFALSE;
3038           break;
3039         }
3040       }
3041       if (cansign){
3042         if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3043         if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3044         if (!c->IsUsed()) c->Use();
3045       }
3046     }
3047   }
3048   return besttrack;
3049 }
3050
3051
3052
3053 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3054 {
3055   //
3056   // get "best" hypothesys
3057   //
3058
3059   Int_t nentries = itsTracks.GetEntriesFast();
3060   for (Int_t i=0;i<nentries;i++){
3061     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3062     if (!track) continue;
3063     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3064     if (!array) continue;
3065     if (array->GetEntriesFast()<=0) continue;
3066     //
3067     AliITStrackMI* longtrack=0;
3068     Float_t minn=0;
3069     Float_t maxchi2=1000;
3070     for (Int_t j=0;j<array->GetEntriesFast();j++){
3071       AliITStrackMI* track = (AliITStrackMI*)array->At(j);
3072       if (!track) continue;
3073       if (track->GetGoldV0()) {
3074         longtrack = track;   //gold V0 track taken
3075         break;
3076       }
3077       if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3078       Float_t chi2 = track->GetChi2MIP(0);
3079       if (fAfterV0){
3080         if (!track->GetGoldV0()&&track->GetConstrain()==kFALSE) chi2+=5;
3081       }
3082       if (track->GetNumberOfClusters()+track->GetNDeadZone()>minn) maxchi2 = track->GetChi2MIP(0);
3083       //
3084       if (chi2 > maxchi2) continue;
3085       minn= track->GetNumberOfClusters()+track->GetNDeadZone();
3086       maxchi2 = chi2;
3087       longtrack=track;
3088     }
3089     //
3090     //
3091     //
3092     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3093     if (!longtrack) {longtrack = besttrack;}
3094     else besttrack= longtrack;
3095     //
3096     if (besttrack){
3097       Int_t list[6];
3098       AliITSRecPoint * clist[6];
3099       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3100       //
3101       track->SetNUsed(shared);
3102       track->SetNSkipped(besttrack->GetNSkipped());
3103       track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3104       if (shared>0){
3105         Int_t nshared;
3106         Int_t overlist[6];
3107         //
3108         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3109         //if (sharedtrack==-1) sharedtrack=0;
3110         if (sharedtrack>=0){
3111           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3112         }
3113       }
3114       if (besttrack&&fAfterV0){
3115         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3116       }
3117       if (besttrack&&fConstraint[fPass])
3118         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3119       //if (besttrack&&besttrack->fConstrain)
3120       //        UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3121       if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5){
3122         if ( (TMath::Abs(besttrack->GetD(0))>0.1) && fConstraint[fPass]) {
3123           track->SetReconstructed(kFALSE);
3124         }
3125         if ( (TMath::Abs(besttrack->GetD(1))>0.1) && fConstraint[fPass]){
3126           track->SetReconstructed(kFALSE);
3127         }
3128       }
3129
3130     }
3131   }
3132 }
3133
3134
3135 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3136   //--------------------------------------------------------------------
3137   //This function "cooks" a track label. If label<0, this track is fake.
3138   //--------------------------------------------------------------------
3139   Int_t tpcLabel=-1;
3140
3141   if ( track->GetESDtrack())   tpcLabel =  TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3142
3143    track->SetChi2MIP(9,0);
3144    Int_t nwrong=0;
3145    for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3146      Int_t cindex = track->GetClusterIndex(i);
3147      Int_t l=(cindex & 0xf0000000) >> 28;
3148      AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3149      Int_t isWrong=1;
3150      for (Int_t ind=0;ind<3;ind++){
3151        if (tpcLabel>0)
3152          if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
3153      }
3154      track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3155      nwrong+=isWrong;
3156    }
3157    track->SetFakeRatio(double(nwrong)/double(track->GetNumberOfClusters()));
3158    if (tpcLabel>0){
3159      if (track->GetFakeRatio()>wrong) track->SetLabel(-tpcLabel);
3160      else
3161        track->SetLabel(tpcLabel);
3162    }
3163
3164 }
3165
3166
3167
3168 void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
3169 {
3170   //
3171   //
3172   //  Int_t list[6];
3173   //AliITSRecPoint * clist[6];
3174   //  Int_t shared = GetNumberOfSharedClusters(track,index,list,clist);
3175   Float_t dedx[4];
3176   Int_t accepted=0;
3177   track->SetChi2MIP(9,0);
3178   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3179     Int_t cindex = track->GetClusterIndex(i);
3180     Int_t l=(cindex & 0xf0000000) >> 28;
3181     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3182     Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3183     Int_t isWrong=1;
3184     for (Int_t ind=0;ind<3;ind++){
3185       if (cl->GetLabel(ind)==lab) isWrong=0;
3186     }
3187     track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3188     if (l<2) continue;
3189     //if (l>3 && (cl->GetNy()>4) || (cl->GetNz()>4)) continue;  //shared track
3190     //if (l>3&& !(cl->GetType()==1||cl->GetType()==10)) continue;
3191     //if (l<4&& !(cl->GetType()==1)) continue;
3192     dedx[accepted]= track->GetSampledEdx(i);
3193     //dedx[accepted]= track->fNormQ[l];
3194     accepted++;
3195   }
3196   if (accepted<1) {
3197     track->SetdEdx(0);
3198     return;
3199   }
3200   Int_t indexes[4];
3201   TMath::Sort(accepted,dedx,indexes,kFALSE);
3202   Double_t low=0.;
3203   Double_t up=0.51;
3204   Double_t nl=low*accepted, nu =up*accepted;
3205   Float_t sumamp = 0;
3206   Float_t sumweight =0;
3207   for (Int_t i=0; i<accepted; i++) {
3208     Float_t weight =1;
3209     if (i<nl+0.1)     weight = TMath::Max(1.-(nl-i),0.);
3210     if (i>nu-1)     weight = TMath::Max(nu-i,0.);
3211     sumamp+= dedx[indexes[i]]*weight;
3212     sumweight+=weight;
3213   }
3214   track->SetdEdx(sumamp/sumweight);
3215 }
3216
3217
3218 void  AliITStrackerMI::MakeCoeficients(Int_t ntracks){
3219   //
3220   //
3221   if (fCoeficients) delete []fCoeficients;
3222   fCoeficients = new Float_t[ntracks*48];
3223   for (Int_t i=0;i<ntracks*48;i++) fCoeficients[i]=-1.;
3224 }
3225
3226
3227 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3228 {
3229   //
3230   //
3231   //
3232   Float_t erry,errz;
3233   Float_t theta = track->GetTgl();
3234   Float_t phi   = track->GetSnp();
3235   phi = TMath::Sqrt(phi*phi/(1.-phi*phi));
3236   GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz);
3237   Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz);
3238   Float_t ny,nz;
3239   GetNTeor(layer,cluster, theta,phi,ny,nz);
3240   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3241   if (delta>1){
3242     chi2+=0.5*TMath::Min(delta/2,2.);
3243     chi2+=2.*cluster->GetDeltaProbability();
3244   }
3245   //
3246   track->SetNy(layer,ny);
3247   track->SetNz(layer,nz);
3248   track->SetSigmaY(layer,erry);
3249   track->SetSigmaZ(layer, errz);
3250   //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3251   track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/(1.- track->GetSnp()*track->GetSnp())));
3252   return chi2;
3253
3254 }
3255
3256 Int_t    AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3257 {
3258   //
3259   //
3260   //
3261   Int_t layer = (index & 0xf0000000) >> 28;
3262   track->SetClIndex(layer, index);
3263   if ( (layer>1) &&track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3264     chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3265     track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3266   }
3267
3268   if (cl->GetQ()<=0) return 0;  // ingore the "virtual" clusters
3269
3270   // Take into account the mis-alignment
3271   Double_t x=track->GetX()+cl->GetX();
3272   if (!track->PropagateTo(x,0.,0.)) return 0;
3273
3274   return track->UpdateMI(cl->GetY(),cl->GetZ(),track->GetSigmaY(layer),track->GetSigmaZ(layer),chi2,index);
3275 }
3276
3277 void AliITStrackerMI::GetNTeor(Int_t layer, const AliITSRecPoint* /*cl*/, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz)
3278 {
3279   //
3280   //get "mean shape"
3281   //
3282   if (layer==0){
3283     ny = 1.+TMath::Abs(phi)*3.2;
3284     nz = 1.+TMath::Abs(theta)*0.34;
3285     return;
3286   }
3287   if (layer==1){
3288     ny = 1.+TMath::Abs(phi)*3.2;
3289     nz = 1.+TMath::Abs(theta)*0.28;
3290     return;
3291   }
3292
3293   if (layer>3){
3294     ny = 2.02+TMath::Abs(phi)*1.95;
3295     nz = 2.02+TMath::Abs(phi)*2.35;
3296     return;
3297   }
3298   ny  = 6.6-2.7*TMath::Abs(phi);
3299   nz  = 2.8-3.11*TMath::Abs(phi)+0.45*TMath::Abs(theta);
3300 }
3301
3302
3303
3304 Int_t AliITStrackerMI::GetError(Int_t layer, const AliITSRecPoint*cl, Float_t theta, Float_t phi,Float_t expQ, Float_t &erry, Float_t &errz)
3305 {
3306   //calculate cluster position error
3307   //
3308   Float_t nz,ny;
3309   GetNTeor(layer, cl,theta,phi,ny,nz);
3310   erry   = TMath::Sqrt(cl->GetSigmaY2());
3311   errz   = TMath::Sqrt(cl->GetSigmaZ2());
3312   //
3313   // PIXELS
3314   if (layer<2){
3315
3316     if (TMath::Abs(ny-cl->GetNy())>0.6)  {
3317       if (ny<cl->GetNy()){
3318         erry*=0.4+TMath::Abs(ny-cl->GetNy());
3319         errz*=0.4+TMath::Abs(ny-cl->GetNy());
3320       }else{
3321         erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
3322         errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
3323       }
3324     }
3325     if (TMath::Abs(nz-cl->GetNz())>1.)  {
3326       erry*=TMath::Abs(nz-cl->GetNz());
3327       errz*=TMath::Abs(nz-cl->GetNz());
3328     }
3329     erry*=0.85;
3330     errz*=0.85;
3331     erry= TMath::Min(erry,float(0.005));
3332     errz= TMath::Min(errz,float(0.03));
3333     return 10;
3334   }
3335
3336 //STRIPS
3337   if (layer>3){
3338     //factor 1.8 appears in new simulation
3339     //
3340     Float_t scale=1.8;
3341     if (cl->GetNy()==100||cl->GetNz()==100){
3342       erry = 0.004*scale;
3343       errz = 0.2*scale;
3344       return 100;
3345     }
3346     if (cl->GetNy()+cl->GetNz()>12){
3347       erry = 0.06*scale;
3348       errz = 0.57*scale;
3349       return 100;
3350     }
3351     Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
3352     Float_t chargematch = TMath::Max(double(normq/expQ),2.);
3353     //
3354     if (cl->GetType()==1 || cl->GetType()==10 ){
3355       if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
3356         errz = 0.043*scale;
3357         erry = 0.00094*scale;
3358         return 101;
3359       }
3360       if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
3361         errz = 0.06*scale;
3362         erry =0.0013*scale;
3363         return 102;
3364       }
3365       erry = 0.0027*scale;
3366       errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15)*scale;
3367       return 103;
3368     }
3369     if (cl->GetType()==2 || cl->GetType()==11 ){
3370       erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05)*scale;
3371       errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5)*scale;
3372       return 104;
3373     }
3374
3375     if (cl->GetType()>100 ){
3376       if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
3377         errz = 0.05*scale;
3378         erry = 0.00096*scale;
3379         return 105;
3380       }
3381       if (cl->GetNy()+cl->GetNz()-nz-ny<1){
3382         errz = 0.10*scale;
3383         erry = 0.0025*scale;
3384         return 106;
3385       }
3386
3387       errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4)*scale;
3388       erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05)*scale;
3389       return 107;
3390     }
3391     Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
3392     if (diff<1) diff=1;
3393     if (diff>4) diff=4;
3394
3395     if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){
3396       errz = 0.14*diff;
3397       erry = 0.003*diff;
3398       return 108;
3399     }
3400     erry = 0.04*diff;
3401     errz = 0.06*diff;
3402     return 109;
3403   }
3404   //DRIFTS
3405   Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
3406   Float_t chargematch = normq/expQ;
3407   Float_t factorz=1;
3408   Int_t   cnz = cl->GetNz()%10;
3409   //charge match
3410   if (cl->GetType()==1){
3411     if (chargematch<1.25){
3412       erry =  0.0028*(1.+6./cl->GetQ());  // gold clusters
3413     }
3414     else{
3415       erry = 0.003*chargematch;
3416       if (cl->GetNz()==3) erry*=1.5;
3417     }
3418     if (chargematch<1.0){
3419       errz =  0.0011*(1.+6./cl->GetQ());
3420     }
3421     else{
3422       errz = 0.002*(1+2*(chargematch-1.));
3423     }
3424     if (cnz>nz+0.6) {
3425       erry*=(cnz-nz+0.5);
3426       errz*=1.4*(cnz-nz+0.5);
3427     }
3428   }
3429   if (cl->GetType()>1){
3430     if (chargematch<1){
3431       erry =  0.00385*(1.+6./cl->GetQ());  // gold clusters
3432       errz =  0.0016*(1.+6./cl->GetQ());
3433     }
3434     else{
3435       errz = 0.0014*(1+3*(chargematch-1.));
3436       erry = 0.003*(1+3*(chargematch-1.));
3437     }
3438     if (cnz>nz+0.6) {
3439       erry*=(cnz-nz+0.5);
3440       errz*=1.4*(cnz-nz+0.5);
3441     }
3442   }
3443
3444   if (TMath::Abs(cl->GetY())>2.5){
3445     factorz*=1+2*(TMath::Abs(cl->GetY())-2.5);
3446   }
3447   if (TMath::Abs(cl->GetY())<1){
3448     factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.);
3449   }
3450   factorz= TMath::Min(factorz,float(4.));
3451   errz*=factorz;
3452
3453   erry= TMath::Min(erry,float(0.05));
3454   errz= TMath::Min(errz,float(0.05));
3455   return 200;
3456 }
3457
3458
3459
3460 void   AliITStrackerMI::GetDCASigma(AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3461 {
3462   //
3463   //DCA sigmas parameterization
3464   //to be paramterized using external parameters in future
3465   //
3466   //
3467   sigmarfi = 0.004+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3468   sigmaz   = 0.011+4.37*TMath::Abs(track->GetC());
3469 }
3470
3471
3472 void AliITStrackerMI::SignDeltas( TObjArray *ClusterArray, Float_t vz)
3473 {
3474   //
3475   //
3476   Int_t entries = ClusterArray->GetEntriesFast();
3477   if (entries<4) return;
3478   AliITSRecPoint* cluster = (AliITSRecPoint*)ClusterArray->At(0);
3479   Int_t layer = cluster->GetLayer();
3480   if (layer>1) return;
3481   Int_t index[10000];
3482   Int_t ncandidates=0;
3483   Float_t r = (layer>0)? 7:4;
3484   //
3485   for (Int_t i=0;i<entries;i++){
3486     AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(i);
3487     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3488     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3489     index[ncandidates] = i;  //candidate to belong to delta electron track
3490     ncandidates++;
3491     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3492       cl0->SetDeltaProbability(1);
3493     }
3494   }
3495   //
3496   //
3497   //
3498   for (Int_t i=0;i<ncandidates;i++){
3499     AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(index[i]);
3500     if (cl0->GetDeltaProbability()>0.8) continue;
3501     //
3502     Int_t ncl = 0;
3503     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3504     sumy=sumz=sumy2=sumyz=sumw=0.0;
3505     for (Int_t j=0;j<ncandidates;j++){
3506       if (i==j) continue;
3507       AliITSRecPoint* cl1 = (AliITSRecPoint*)ClusterArray->At(index[j]);
3508       //
3509       Float_t dz = cl0->GetZ()-cl1->GetZ();
3510       Float_t dy = cl0->GetY()-cl1->GetY();
3511       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3512         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3513         y[ncl] = cl1->GetY();
3514         z[ncl] = cl1->GetZ();
3515         sumy+= y[ncl]*weight;
3516         sumz+= z[ncl]*weight;
3517         sumy2+=y[ncl]*y[ncl]*weight;
3518         sumyz+=y[ncl]*z[ncl]*weight;
3519         sumw+=weight;
3520         ncl++;
3521       }
3522     }
3523     if (ncl<4) continue;
3524     Float_t det = sumw*sumy2  - sumy*sumy;
3525     Float_t delta=1000;
3526     if (TMath::Abs(det)>0.01){
3527       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
3528       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
3529       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3530     }
3531     else{
3532       Float_t z0  = sumyz/sumy;
3533       delta = TMath::Abs(cl0->GetZ()-z0);
3534     }
3535     if ( delta<0.05) {
3536       cl0->SetDeltaProbability(1-20.*delta);
3537     }
3538   }
3539 }
3540
3541
3542 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3543 {
3544   //
3545   //
3546   track->UpdateESDtrack(flags);
3547   AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3548   if (oldtrack) delete oldtrack;
3549   track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3550   if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3551     printf("Problem\n");
3552   }
3553 }
3554
3555
3556
3557 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3558   //
3559   //Get nearest upper layer close to the point xr.
3560   // rough approximation
3561   //
3562   const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3563   Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3564   Int_t res =6;
3565   for (Int_t i=0;i<6;i++){
3566     if (radius<kRadiuses[i]){
3567       res =i;
3568       break;
3569     }
3570   }
3571   return res;
3572 }
3573
3574
3575 void AliITStrackerMI::UpdateTPCV0(AliESD *event){
3576   //
3577   //try to update, or reject TPC  V0s
3578   //
3579   Int_t nv0s = event->GetNumberOfV0s();
3580   Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3581
3582   for (Int_t i=0;i<nv0s;i++){
3583     AliESDv0 * vertex = event->GetV0(i);
3584     Int_t ip = vertex->GetIndex(0);
3585     Int_t im = vertex->GetIndex(1);
3586     //
3587     TObjArray * arrayp = (ip<nitstracks) ? (TObjArray*)fTrackHypothesys.At(ip):0;
3588     TObjArray * arraym = (im<nitstracks) ? (TObjArray*)fTrackHypothesys.At(im):0;
3589     AliITStrackMI * trackp = (arrayp!=0) ? (AliITStrackMI*)arrayp->At(0):0;
3590     AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0;
3591     //
3592     //
3593     if (trackp){
3594       if (trackp->GetNumberOfClusters()+trackp->GetNDeadZone()>5.5){
3595         if (trackp->GetConstrain()&&trackp->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3596         if (!trackp->GetConstrain()&&trackp->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3597       }
3598     }
3599
3600     if (trackm){
3601       if (trackm->GetNumberOfClusters()+trackm->GetNDeadZone()>5.5){
3602         if (trackm->GetConstrain()&&trackm->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3603         if (!trackm->GetConstrain()&&trackm->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3604       }
3605     }
3606     if (vertex->GetStatus()==-100) continue;
3607     //
3608     Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
3609     Int_t clayer = GetNearestLayer(xrp);                    //I.B.
3610     vertex->SetNBefore(clayer);        //
3611     vertex->SetChi2Before(9*clayer);   //
3612     vertex->SetNAfter(6-clayer);       //
3613     vertex->SetChi2After(0);           //
3614     //
3615     if (clayer >1 ){ // calculate chi2 before vertex
3616       Float_t chi2p = 0, chi2m=0;
3617       //
3618       if (trackp){
3619         for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3620           if (trackp->GetClIndex(ilayer)>0){
3621             chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3622               trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3623           }
3624           else{
3625             chi2p+=9;
3626           }
3627         }
3628       }else{
3629         chi2p = 9*clayer;
3630       }
3631       //
3632       if (trackm){
3633         for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3634           if (trackm->GetClIndex(ilayer)>0){
3635             chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3636               trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3637           }
3638           else{
3639             chi2m+=9;
3640           }
3641         }
3642       }else{
3643         chi2m = 9*clayer;
3644       }
3645       vertex->SetChi2Before(TMath::Min(chi2p,chi2m));
3646       if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10);  // track exist before vertex
3647     }
3648
3649     if (clayer < 5 ){ // calculate chi2 after vertex
3650       Float_t chi2p = 0, chi2m=0;
3651       //
3652       if (trackp&&TMath::Abs(trackp->GetTgl())<1.){
3653         for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3654           if (trackp->GetClIndex(ilayer)>0){
3655             chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3656               trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3657           }
3658           else{
3659             chi2p+=9;
3660           }
3661         }
3662       }else{
3663         chi2p = 0;
3664       }
3665       //
3666       if (trackm&&TMath::Abs(trackm->GetTgl())<1.){
3667         for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3668           if (trackm->GetClIndex(ilayer)>0){
3669             chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3670               trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3671           }
3672           else{
3673             chi2m+=9;
3674           }
3675         }
3676       }else{
3677         chi2m = 0;
3678       }
3679       vertex->SetChi2After(TMath::Max(chi2p,chi2m));
3680       if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20);  // track not found in ITS
3681     }
3682   }
3683   //
3684 }
3685
3686
3687
3688 void  AliITStrackerMI::FindV02(AliESD *event)
3689 {
3690   //
3691   // V0 finder
3692   //
3693   //  Cuts on DCA -  R dependent
3694   //          max distance DCA between 2 tracks cut
3695   //          maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist);
3696   //
3697   const Float_t kMaxDist0      = 0.1;
3698   const Float_t kMaxDist1      = 0.1;
3699   const Float_t kMaxDist       = 1;
3700   const Float_t kMinPointAngle  = 0.85;
3701   const Float_t kMinPointAngle2 = 0.99;
3702   const Float_t kMinR           = 0.5;
3703   const Float_t kMaxR           = 220;
3704   //const Float_t kCausality0Cut   = 0.19;
3705   //const Float_t kLikelihood01Cut = 0.25;
3706   //const Float_t kPointAngleCut   = 0.9996;
3707   const Float_t kCausality0Cut   = 0.19;
3708   const Float_t kLikelihood01Cut = 0.45;
3709   const Float_t kLikelihood1Cut  = 0.5;
3710   const Float_t kCombinedCut     = 0.55;
3711
3712   //
3713   //
3714   TTreeSRedirector &cstream = *fDebugStreamer;
3715   Int_t ntracks    = event->GetNumberOfTracks();
3716   Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3717   fOriginal.Expand(ntracks);
3718   fTrackHypothesys.Expand(ntracks);
3719   fBestHypothesys.Expand(ntracks);
3720   //
3721   AliHelix * helixes   = new AliHelix[ntracks+2];
3722   TObjArray trackarray(ntracks+2);     //array with tracks - with vertex constrain
3723   TObjArray trackarrayc(ntracks+2);    //array of "best    tracks" - without vertex constrain
3724   TObjArray trackarrayl(ntracks+2);    //array of "longest tracks" - without vertex constrain
3725   Bool_t * forbidden   = new Bool_t [ntracks+2];
3726   Int_t   *itsmap      = new Int_t  [ntracks+2];
3727   Float_t *dist        = new Float_t[ntracks+2];
3728   Float_t *normdist0   = new Float_t[ntracks+2];
3729   Float_t *normdist1   = new Float_t[ntracks+2];
3730   Float_t *normdist    = new Float_t[ntracks+2];
3731   Float_t *norm        = new Float_t[ntracks+2];
3732   Float_t *maxr        = new Float_t[ntracks+2];
3733   Float_t *minr        = new Float_t[ntracks+2];
3734   Float_t *minPointAngle= new Float_t[ntracks+2];
3735   //
3736   AliV0 *pvertex      = new AliV0;
3737   AliITStrackMI * dummy= new AliITStrackMI;
3738   dummy->SetLabel(0);
3739   AliITStrackMI  trackat0;    //temporary track for DCA calculation
3740   //
3741   Float_t primvertex[3]={GetX(),GetY(),GetZ()};
3742   //
3743   // make its -  esd map
3744   //
3745   for (Int_t itrack=0;itrack<ntracks+2;itrack++) {
3746     itsmap[itrack]        = -1;
3747     forbidden[itrack]     = kFALSE;
3748     maxr[itrack]          = kMaxR;
3749     minr[itrack]          = kMinR;
3750     minPointAngle[itrack] = kMinPointAngle;
3751   }
3752   for (Int_t itrack=0;itrack<nitstracks;itrack++){
3753     AliITStrackMI * original =   (AliITStrackMI*)(fOriginal.At(itrack));
3754     Int_t           esdindex =   original->GetESDtrack()->GetID();
3755     itsmap[esdindex]         =   itrack;
3756   }
3757   //
3758   // create its tracks from esd tracks if not done before
3759   //
3760   for (Int_t itrack=0;itrack<ntracks;itrack++){
3761     if (itsmap[itrack]>=0) continue;
3762     AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack)));
3763     //tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY());
3764     //tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ();
3765     tpctrack->GetDZ(GetX(),GetY(),GetZ(),tpctrack->GetDP());   //I.B.
3766     if (tpctrack->GetD(0)<20 && tpctrack->GetD(1)<20){
3767       // tracks which can reach inner part of ITS
3768       // propagate track to outer its volume - with correction for material
3769       CorrectForDeadZoneMaterial(tpctrack);
3770     }
3771     itsmap[itrack] = nitstracks;
3772     fOriginal.AddAt(tpctrack,nitstracks);
3773     nitstracks++;
3774   }
3775   //
3776   // fill temporary arrays
3777   //
3778   for (Int_t itrack=0;itrack<ntracks;itrack++){
3779     AliESDtrack *  esdtrack = event->GetTrack(itrack);
3780     Int_t          itsindex = itsmap[itrack];
3781     AliITStrackMI *original = (AliITStrackMI*)fOriginal.At(itsmap[itrack]);
3782     if (!original) continue;
3783     AliITStrackMI *bestConst  = 0;
3784     AliITStrackMI *bestLong   = 0;
3785     AliITStrackMI *best       = 0;
3786     //
3787     //
3788     TObjArray * array    = (TObjArray*)  fTrackHypothesys.At(itsindex);
3789     Int_t       hentries = (array==0) ?  0 : array->GetEntriesFast();
3790     // Get best track with vertex constrain
3791     for (Int_t ih=0;ih<hentries;ih++){
3792       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
3793       if (!trackh->GetConstrain()) continue;
3794       if (!bestConst) bestConst = trackh;
3795       if (trackh->GetNumberOfClusters()>5.0){
3796         bestConst  = trackh;                         // full track -  with minimal chi2
3797         break;
3798       }
3799       if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone())  continue;
3800       bestConst = trackh;
3801       break;
3802     }
3803     // Get best long track without vertex constrain and best track without vertex constrain
3804     for (Int_t ih=0;ih<hentries;ih++){
3805       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
3806       if (trackh->GetConstrain()) continue;
3807       if (!best)     best     = trackh;
3808       if (!bestLong) bestLong = trackh;
3809       if (trackh->GetNumberOfClusters()>5.0){
3810         bestLong  = trackh;                         // full track -  with minimal chi2
3811         break;
3812       }
3813       if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone())  continue;
3814       bestLong = trackh;
3815     }
3816     if (!best) {
3817       best     = original;
3818       bestLong = original;
3819     }
3820     //I.B. trackat0 = *bestLong;
3821     new (&trackat0) AliITStrackMI(*bestLong);
3822     Double_t xx,yy,zz,alpha;
3823     bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz);
3824     alpha = TMath::ATan2(yy,xx);
3825     trackat0.Propagate(alpha,0);
3826     // calculate normalized distances to the vertex
3827     //
3828     Float_t ptfac  = (1.+100.*TMath::Abs(trackat0.GetC()));
3829     if ( bestLong->GetNumberOfClusters()>3 ){
3830       dist[itsindex]      = trackat0.GetY();
3831       norm[itsindex]      = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
3832       normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
3833       normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
3834       normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
3835       if (!bestConst){
3836         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<6) normdist[itsindex]*=2.;
3837         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<5) normdist[itsindex]*=2.;
3838         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<4) normdist[itsindex]*=2.;
3839       }else{
3840         if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<6) normdist[itsindex]*=1.5;
3841         if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<5) normdist[itsindex]*=1.5;
3842       }
3843     }
3844     else{
3845       if (bestConst&&bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5){
3846         dist[itsindex] = bestConst->GetD(0);
3847         norm[itsindex] = bestConst->GetDnorm(0);
3848         normdist0[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
3849         normdist1[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
3850         normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
3851       }else{
3852         dist[itsindex]      = trackat0.GetY();
3853         norm[itsindex]      = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
3854         normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
3855         normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
3856         normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
3857         if (TMath::Abs(trackat0.GetTgl())>1.05){
3858           if (normdist[itsindex]<3) forbidden[itsindex]=kTRUE;
3859           if (normdist[itsindex]>3) {
3860             minr[itsindex] = TMath::Max(Float_t(40.),minr[itsindex]);
3861           }
3862         }
3863       }
3864     }
3865     //
3866     //-----------------------------------------------------------
3867     //Forbid primary track candidates -
3868     //
3869     //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5");
3870     //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5");
3871     //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0");
3872     //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0");
3873     //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2");
3874     //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1");
3875     //-----------------------------------------------------------
3876     if (bestConst){
3877       if (bestLong->GetNumberOfClusters()<4       && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5)               forbidden[itsindex]=kTRUE;
3878       if (normdist[itsindex]<3 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>5.5)               forbidden[itsindex]=kTRUE;
3879       if (normdist[itsindex]<2 && bestConst->GetClIndex(0)>0 && bestConst->GetClIndex(1)>0 ) forbidden[itsindex]=kTRUE;
3880       if (normdist[itsindex]<1 && bestConst->GetClIndex(0)>0)                              forbidden[itsindex]=kTRUE;
3881       if (normdist[itsindex]<4 && bestConst->GetNormChi2(0)<2)                             forbidden[itsindex]=kTRUE;
3882       if (normdist[itsindex]<5 && bestConst->GetNormChi2(0)<1)                             forbidden[itsindex]=kTRUE;
3883       if (bestConst->GetNormChi2(0)<2.5) {
3884         minPointAngle[itsindex]= 0.9999;
3885         maxr[itsindex]         = 10;
3886       }
3887     }
3888     //
3889     //forbid daughter kink candidates
3890     //
3891     if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
3892     Bool_t isElectron = kTRUE;
3893     Bool_t isProton   = kTRUE;
3894     Double_t pid[5];
3895     esdtrack->GetESDpid(pid);
3896     for (Int_t i=1;i<5;i++){
3897       if (pid[0]<pid[i]) isElectron= kFALSE;
3898       if (pid[4]<pid[i]) isProton= kFALSE;
3899     }
3900     if (isElectron){
3901       forbidden[itsindex]=kFALSE;
3902       normdist[itsindex]*=-1;
3903     }
3904     if (isProton){
3905       if (normdist[itsindex]>2) forbidden[itsindex]=kFALSE;
3906       normdist[itsindex]*=-1;
3907     }
3908
3909     //
3910     // Causality cuts in TPC volume
3911     //
3912     if (esdtrack->GetTPCdensity(0,10) >0.6)  maxr[itsindex] = TMath::Min(Float_t(110),maxr[itsindex]);
3913     if (esdtrack->GetTPCdensity(10,30)>0.6)  maxr[itsindex] = TMath::Min(Float_t(120),maxr[itsindex]);
3914     if (esdtrack->GetTPCdensity(20,40)>0.6)  maxr[itsindex] = TMath::Min(Float_t(130),maxr[itsindex]);
3915     if (esdtrack->GetTPCdensity(30,50)>0.6)  maxr[itsindex] = TMath::Min(Float_t(140),maxr[itsindex]);
3916     //
3917     if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->GetNumberOfClusters()<3) minr[itsindex]=100;
3918     //
3919     //
3920     if (kFALSE){
3921       cstream<<"Track"<<
3922         "Tr0.="<<best<<
3923         "Tr1.="<<((bestConst)? bestConst:dummy)<<
3924         "Tr2.="<<bestLong<<
3925         "Tr3.="<<&trackat0<<
3926         "Esd.="<<esdtrack<<
3927         "Dist="<<dist[itsindex]<<
3928         "ND0="<<normdist0[itsindex]<<
3929         "ND1="<<normdist1[itsindex]<<
3930         "ND="<<normdist[itsindex]<<
3931         "Pz="<<primvertex[2]<<
3932         "Forbid="<<forbidden[itsindex]<<
3933         "\n";
3934       //
3935     }
3936     trackarray.AddAt(best,itsindex);
3937     trackarrayc.AddAt(bestConst,itsindex);
3938     trackarrayl.AddAt(bestLong,itsindex);
3939     new (&helixes[itsindex]) AliHelix(*best);
3940   }
3941   //
3942   //
3943   //
3944   // first iterration of V0 finder
3945   //
3946   for (Int_t iesd0=0;iesd0<ntracks;iesd0++){
3947     Int_t itrack0 = itsmap[iesd0];
3948     if (forbidden[itrack0]) continue;
3949     AliITStrackMI * btrack0 = (AliITStrackMI*)trackarray.At(itrack0);
3950     if (!btrack0) continue;
3951     if (btrack0->GetSign()>0) continue;
3952     AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
3953     //
3954     for (Int_t iesd1=0;iesd1<ntracks;iesd1++){
3955       Int_t itrack1 = itsmap[iesd1];
3956       if (forbidden[itrack1]) continue;
3957
3958       AliITStrackMI * btrack1 = (AliITStrackMI*)trackarray.At(itrack1);
3959       if (!btrack1) continue;
3960       if (btrack1->GetSign()<0) continue;
3961       Bool_t isGold = kFALSE;
3962       if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){
3963         isGold = kTRUE;
3964       }
3965       AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
3966       AliHelix &h1 = helixes[itrack0];
3967       AliHelix &h2 = helixes[itrack1];
3968       //
3969       // find linear distance
3970       Double_t rmin =0;
3971       //
3972       //
3973       //
3974       Double_t phase[2][2],radius[2];
3975       Int_t  points = h1.GetRPHIintersections(h2, phase, radius);
3976       if    (points==0)  continue;
3977       Double_t delta[2]={1000000,1000000};
3978       rmin = radius[0];
3979       h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]);
3980       if (points==2){
3981         if (radius[1]<rmin) rmin = radius[1];
3982         h1.ParabolicDCA(h2,phase[1][0],phase[1][1],radius[1],delta[1]);
3983       }
3984       rmin = TMath::Sqrt(rmin);
3985       Double_t distance = 0;
3986       Double_t radiusC  = 0;
3987       Int_t    iphase   = 0;
3988       if (points==1 || delta[0]<delta[1]){
3989         distance = TMath::Sqrt(delta[0]);
3990         radiusC  = TMath::Sqrt(radius[0]);
3991       }else{
3992         distance = TMath::Sqrt(delta[1]);
3993         radiusC  = TMath::Sqrt(radius[1]);
3994         iphase=1;
3995       }
3996       if (radiusC<TMath::Max(minr[itrack0],minr[itrack1]))    continue;
3997       if (radiusC>TMath::Min(maxr[itrack0],maxr[itrack1]))     continue;
3998       Float_t maxDist  = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1));
3999       if (distance>maxDist) continue;
4000       Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex);
4001       if (pointAngle<TMath::Max(minPointAngle[itrack0],minPointAngle[itrack1])) continue;
4002       //
4003       //
4004       //      Double_t distance = TestV0(h1,h2,pvertex,rmin);
4005       //
4006       //       if (distance>maxDist)           continue;
4007       //       if (pvertex->GetRr()<kMinR)     continue;
4008       //       if (pvertex->GetRr()>kMaxR)     continue;
4009       AliITStrackMI * track0=btrack0;
4010       AliITStrackMI * track1=btrack1;
4011       //      if (pvertex->GetRr()<3.5){
4012       if (radiusC<3.5){
4013         //use longest tracks inside the pipe
4014         track0 = (AliITStrackMI*)trackarrayl.At(itrack0);
4015         track1 = (AliITStrackMI*)trackarrayl.At(itrack1);
4016       }
4017       //
4018       //
4019       pvertex->SetParamN(*track0);
4020       pvertex->SetParamP(*track1);
4021       pvertex->Update(primvertex);
4022       pvertex->SetClusters(track0->ClIndex(),track1->ClIndex());  // register clusters
4023
4024       if (pvertex->GetRr()<kMinR) continue;
4025       if (pvertex->GetRr()>kMaxR) continue;
4026       if (pvertex->GetV0CosineOfPointingAngle()<kMinPointAngle) continue;
4027 //Bo:      if (pvertex->GetDist2()>maxDist) continue;
4028       if (pvertex->GetDcaV0Daughters()>maxDist) continue;
4029 //Bo:        pvertex->SetLab(0,track0->GetLabel());
4030 //Bo:        pvertex->SetLab(1,track1->GetLabel());
4031       pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4032       pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4033       //
4034       AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
4035       AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
4036
4037       //
4038       //
4039       TObjArray * array0b     = (TObjArray*)fBestHypothesys.At(itrack0);
4040       if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->GetTgl())<1.1)
4041         FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack0),itrack0, kFALSE);
4042       TObjArray * array1b    = (TObjArray*)fBestHypothesys.At(itrack1);
4043       if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->GetTgl())<1.1)
4044         FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack1),itrack1, kFALSE);
4045       //
4046       AliITStrackMI * track0b = (AliITStrackMI*)fOriginal.At(itrack0);
4047       AliITStrackMI * track1b = (AliITStrackMI*)fOriginal.At(itrack1);
4048       AliITStrackMI * track0l = (AliITStrackMI*)fOriginal.At(itrack0);
4049       AliITStrackMI * track1l = (AliITStrackMI*)fOriginal.At(itrack1);
4050
4051       Float_t minchi2before0=16;
4052       Float_t minchi2before1=16;
4053       Float_t minchi2after0 =16;
4054       Float_t minchi2after1 =16;
4055       Double_t xrp[3]; pvertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
4056       Int_t maxLayer = GetNearestLayer(xrp);                   //I.B.
4057
4058       if (array0b) for (Int_t i=0;i<5;i++){
4059         // best track after vertex
4060         AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i);
4061         if (!btrack) continue;
4062         if (btrack->GetNumberOfClusters()>track0l->GetNumberOfClusters()) track0l = btrack;
4063         //      if (btrack->fX<pvertex->GetRr()-2.-0.5/(0.1+pvertex->GetAnglep()[2])) {
4064         if (btrack->GetX()<pvertex->GetRr()-2.) {
4065           if ( (maxLayer>i+2|| (i==0)) && btrack->GetNumberOfClusters()==(6-i)&&i<3){
4066             Float_t sumchi2= 0;
4067             Float_t sumn   = 0;
4068             if (maxLayer<3){   // take prim vertex as additional measurement
4069               if (normdist[itrack0]>0 && htrackc0){
4070                 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack0]*normdist[itrack0],16.);
4071               }else{
4072                 sumchi2 +=  TMath::Min((3.-maxLayer)*(3*normdist[itrack0]*normdist[itrack0]+3.),16.);
4073               }
4074               sumn    +=  3-maxLayer;
4075             }
4076             for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4077               sumn+=1.;
4078               if (!btrack->GetClIndex(ilayer)){
4079                 sumchi2+=25;
4080                 continue;
4081               }else{
4082                 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4083                 for (Int_t itrack=0;itrack<4;itrack++){
4084                   if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack0){
4085                     sumchi2+=18.;  //shared cluster
4086                     break;
4087                   }
4088                 }
4089                 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4090                 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4091               }
4092             }
4093             sumchi2/=sumn;
4094             if (sumchi2<minchi2before0) minchi2before0=sumchi2;
4095           }
4096           continue;   //safety space - Geo manager will give exact layer
4097         }
4098         track0b       = btrack;
4099         minchi2after0 = btrack->GetNormChi2(i);
4100         break;
4101       }
4102       if (array1b) for (Int_t i=0;i<5;i++){
4103         // best track after vertex
4104         AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i);
4105         if (!btrack) continue;
4106         if (btrack->GetNumberOfClusters()>track1l->GetNumberOfClusters()) track1l = btrack;
4107         //      if (btrack->fX<pvertex->GetRr()-2-0.5/(0.1+pvertex->GetAnglep()[2])){
4108         if (btrack->GetX()<pvertex->GetRr()-2){
4109           if ((maxLayer>i+2 || (i==0))&&btrack->GetNumberOfClusters()==(6-i)&&(i<3)){
4110             Float_t sumchi2= 0;
4111             Float_t sumn   = 0;
4112             if (maxLayer<3){   // take prim vertex as additional measurement
4113               if (normdist[itrack1]>0 && htrackc1){
4114                 sumchi2 +=  TMath::Min((3.-maxLayer)*normdist[itrack1]*normdist[itrack1],16.);
4115               }else{
4116                 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack1]*normdist[itrack1]+3.),16.);
4117               }
4118               sumn    +=  3-maxLayer;
4119             }
4120             for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4121               sumn+=1.;
4122               if (!btrack->GetClIndex(ilayer)){
4123                 sumchi2+=30;
4124                 continue;
4125               }else{
4126                 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4127                 for (Int_t itrack=0;itrack<4;itrack++){
4128                   if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack1){
4129                     sumchi2+=18.;  //shared cluster
4130                     break;
4131                   }
4132                 }
4133                 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4134                 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4135               }
4136             }
4137             sumchi2/=sumn;
4138             if (sumchi2<minchi2before1) minchi2before1=sumchi2;
4139           }
4140           continue;   //safety space - Geo manager will give exact layer
4141         }
4142         track1b = btrack;
4143         minchi2after1 = btrack->GetNormChi2(i);
4144         break;
4145       }
4146       //
4147       // position resolution - used for DCA cut
4148       Float_t sigmad = track0b->GetSigmaY2()+track0b->GetSigmaZ2()+track1b->GetSigmaY2()+track1b->GetSigmaZ2()+
4149         (track0b->GetX()-pvertex->GetRr())*(track0b->GetX()-pvertex->GetRr())*(track0b->GetSigmaSnp2()+track0b->GetSigmaTgl2())+
4150         (track1b->GetX()-pvertex->GetRr())*(track1b->GetX()-pvertex->GetRr())*(track1b->GetSigmaSnp2()+track1b->GetSigmaTgl2());
4151       sigmad =TMath::Sqrt(sigmad)+0.04;
4152       if (pvertex->GetRr()>50){
4153         Double_t cov0[15],cov1[15];
4154         track0b->GetESDtrack()->GetInnerExternalCovariance(cov0);
4155         track1b->GetESDtrack()->GetInnerExternalCovariance(cov1);
4156         sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+
4157           (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+
4158           (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]);
4159         sigmad =TMath::Sqrt(sigmad)+0.05;
4160       }
4161       //
4162       AliV0 vertex2;
4163       vertex2.SetParamN(*track0b);
4164       vertex2.SetParamP(*track1b);
4165       vertex2.Update(primvertex);
4166       //Bo:      if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4167       if (vertex2.GetDcaV0Daughters()<=pvertex->GetDcaV0Daughters()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4168         pvertex->SetParamN(*track0b);
4169         pvertex->SetParamP(*track1b);
4170         pvertex->Update(primvertex);
4171         pvertex->SetClusters(track0b->ClIndex(),track1b->ClIndex());  // register clusters
4172         pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4173         pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4174       }
4175       pvertex->SetDistSigma(sigmad);
4176       //Bo:      pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);
4177       pvertex->SetNormDCAPrim(normdist[itrack0],normdist[itrack1]);
4178       //
4179       // define likelihhod and causalities
4180       //
4181       Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;
4182       if (maxLayer<1){
4183         Float_t fnorm0 = normdist[itrack0];
4184         if (fnorm0<0) fnorm0*=-3;
4185         Float_t fnorm1 = normdist[itrack1];
4186         if (fnorm1<0) fnorm1*=-3;
4187         if (pvertex->GetAnglep()[2]>0.1 ||  (pvertex->GetRr()<10.5)&& pvertex->GetAnglep()[2]>0.05 || pvertex->GetRr()<3){
4188           pb0    =  TMath::Exp(-TMath::Min(fnorm0,Float_t(16.))/12.);
4189           pb1    =  TMath::Exp(-TMath::Min(fnorm1,Float_t(16.))/12.);
4190         }
4191         pvertex->SetChi2Before(normdist[itrack0]);
4192         pvertex->SetChi2After(normdist[itrack1]);
4193         pvertex->SetNAfter(0);
4194         pvertex->SetNBefore(0);
4195       }else{
4196         pvertex->SetChi2Before(minchi2before0);
4197         pvertex->SetChi2After(minchi2before1);
4198          if (pvertex->GetAnglep()[2]>0.1 || ( pvertex->GetRr()<10.5 && pvertex->GetAnglep()[2]>0.05) || pvertex->GetRr()<3){
4199            pb0    =  TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.);
4200            pb1    =  TMath::Exp(-TMath::Min(minchi2before1,Float_t(16))/12.);
4201          }
4202          pvertex->SetNAfter(maxLayer);
4203          pvertex->SetNBefore(maxLayer);
4204       }
4205       if (pvertex->GetRr()<90){
4206         pa0  *= TMath::Min(track0->GetESDtrack()->GetTPCdensity(0,60),Float_t(1.));
4207         pa1  *= TMath::Min(track1->GetESDtrack()->GetTPCdensity(0,60),Float_t(1.));
4208       }
4209       if (pvertex->GetRr()<20){
4210         pa0  *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2;
4211         pa1  *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2;
4212       }
4213       //
4214       pvertex->SetCausality(pb0,pb1,pa0,pa1);
4215       //
4216       //  Likelihood calculations  - apply cuts
4217       //
4218       Bool_t v0OK = kTRUE;
4219       Float_t p12 = pvertex->GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4];
4220       p12        += pvertex->GetParamN()->GetParameter()[4]*pvertex->GetParamN()->GetParameter()[4];
4221       p12         = TMath::Sqrt(p12);                             // "mean" momenta
4222       Float_t    sigmap0   = 0.0001+0.001/(0.1+pvertex->GetRr());
4223       Float_t    sigmap    = 0.5*sigmap0*(0.6+0.4*p12);           // "resolution: of point angle - as a function of radius and momenta
4224
4225       Float_t causalityA  = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]);
4226       Float_t causalityB  = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Float_t(0.7))*
4227                                         TMath::Min(pvertex->GetCausalityP()[3],Float_t(0.7)));
4228       //
4229       //Bo:      Float_t likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5);
4230       Float_t lDistNorm = pvertex->GetDcaV0Daughters()/pvertex->GetDistSigma();
4231       Float_t likelihood0 = (TMath::Exp(-lDistNorm)+0.1) *(pvertex->GetDcaV0Daughters()<0.5)*(lDistNorm<5);
4232
4233       Float_t likelihood1 = TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/sigmap)+
4234         0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(4.*sigmap))+
4235         0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(8.*sigmap))+
4236         0.1*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/0.01);
4237       //
4238       if (causalityA<kCausality0Cut)                                          v0OK = kFALSE;
4239       if (TMath::Sqrt(likelihood0*likelihood1)<kLikelihood01Cut)              v0OK = kFALSE;
4240       if (likelihood1<kLikelihood1Cut)                                        v0OK = kFALSE;
4241       if (TMath::Power(likelihood0*likelihood1*causalityB,0.33)<kCombinedCut) v0OK = kFALSE;
4242
4243       //
4244       //
4245       if (kFALSE){
4246         Bool_t gold = TMath::Abs(TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel()))==1;
4247         cstream<<"It0"<<
4248           "Tr0.="<<track0<<                       //best without constrain
4249           "Tr1.="<<track1<<                       //best without constrain
4250           "Tr0B.="<<track0b<<                     //best without constrain  after vertex
4251           "Tr1B.="<<track1b<<                     //best without constrain  after vertex
4252           "Tr0C.="<<htrackc0<<                    //best with constrain     if exist
4253           "Tr1C.="<<htrackc1<<                    //best with constrain     if exist
4254           "Tr0L.="<<track0l<<                     //longest best
4255           "Tr1L.="<<track1l<<                     //longest best
4256           "Esd0.="<<track0->GetESDtrack()<<           // esd track0 params
4257           "Esd1.="<<track1->GetESDtrack()<<           // esd track1 params
4258           "V0.="<<pvertex<<                       //vertex properties
4259           "V0b.="<<&vertex2<<                       //vertex properties at "best" track
4260           "ND0="<<normdist[itrack0]<<             //normalize distance for track0
4261           "ND1="<<normdist[itrack1]<<             //normalize distance for track1
4262           "Gold="<<gold<<                         //
4263           //      "RejectBase="<<rejectBase<<             //rejection in First itteration
4264           "OK="<<v0OK<<
4265           "rmin="<<rmin<<
4266           "sigmad="<<sigmad<<
4267           "\n";
4268       }
4269       //if (rejectBase) continue;
4270       //
4271       pvertex->SetStatus(0);
4272       //      if (rejectBase) {
4273       //        pvertex->SetStatus(-100);
4274       //}
4275       if (pvertex->GetV0CosineOfPointingAngle()>kMinPointAngle2) {
4276         //Bo:     pvertex->SetESDindexes(track0->GetESDtrack()->GetID(),track1->GetESDtrack()->GetID());
4277         pvertex->SetIndex(0,track0->GetESDtrack()->GetID());//Bo: consistency 0 for neg
4278         pvertex->SetIndex(1,track1->GetESDtrack()->GetID());//Bo: consistency 1 for pos
4279         if (v0OK){
4280           //      AliV0vertex vertexjuri(*track0,*track1);
4281           //      vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
4282           //      event->AddV0(&vertexjuri);
4283           pvertex->SetStatus(100);
4284         }
4285         pvertex->SetOnFlyStatus(kTRUE);
4286         pvertex->ChangeMassHypothesis(kK0Short);
4287         event->AddV0(pvertex);
4288       }
4289     }
4290   }
4291   //
4292   //
4293   // delete temporary arrays
4294   //
4295   delete[] minPointAngle;
4296   delete[] maxr;
4297   delete[] minr;
4298   delete[] norm;
4299   delete[] normdist;
4300   delete[] normdist1;
4301   delete[] normdist0;
4302   delete[] dist;
4303   delete[] itsmap;
4304   delete[] helixes;
4305   delete   pvertex;
4306 }
4307
4308
4309 void AliITStrackerMI::RefitV02(AliESD *event)
4310 {
4311   //
4312   //try to refit  V0s in the third path of the reconstruction
4313   //
4314   TTreeSRedirector &cstream = *fDebugStreamer;
4315   //
4316   Int_t  nv0s = event->GetNumberOfV0s();
4317   Float_t primvertex[3]={GetX(),GetY(),GetZ()};
4318   AliV0 v0temp;
4319   for (Int_t iv0 = 0; iv0<nv0s;iv0++){
4320     AliV0 * v0mi = (AliV0*)event->GetV0(iv0);
4321     if (!v0mi) continue;
4322     Int_t     itrack0   = v0mi->GetIndex(0);
4323     Int_t     itrack1   = v0mi->GetIndex(1);
4324     AliESDtrack *esd0   = event->GetTrack(itrack0);
4325     AliESDtrack *esd1   = event->GetTrack(itrack1);
4326     if (!esd0||!esd1) continue;
4327     AliITStrackMI tpc0(*esd0);
4328     AliITStrackMI tpc1(*esd1);
4329     Double_t x,y,z; v0mi->GetXYZ(x,y,z); //I.B.
4330     Double_t alpha =TMath::ATan2(y,x);   //I.B.
4331     if (v0mi->GetRr()>85){
4332       if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4333         v0temp.SetParamN(tpc0);
4334         v0temp.SetParamP(tpc1);
4335         v0temp.Update(primvertex);
4336         if (kFALSE) cstream<<"Refit"<<
4337           "V0.="<<v0mi<<
4338           "V0refit.="<<&v0temp<<
4339           "Tr0.="<<&tpc0<<
4340           "Tr1.="<<&tpc1<<
4341           "\n";
4342         //Bo:   if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4343         if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4344           v0mi->SetParamN(tpc0);
4345           v0mi->SetParamP(tpc1);
4346           v0mi->Update(primvertex);
4347         }
4348       }
4349       continue;
4350     }
4351     if (v0mi->GetRr()>35){
4352       CorrectForDeadZoneMaterial(&tpc0);
4353       CorrectForDeadZoneMaterial(&tpc1);
4354       if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4355         v0temp.SetParamN(tpc0);
4356         v0temp.SetParamP(tpc1);
4357         v0temp.Update(primvertex);
4358         if (kFALSE) cstream<<"Refit"<<
4359           "V0.="<<v0mi<<
4360           "V0refit.="<<&v0temp<<
4361           "Tr0.="<<&tpc0<<
4362           "Tr1.="<<&tpc1<<
4363           "\n";
4364         //Bo:   if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4365         if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4366           v0mi->SetParamN(tpc0);
4367           v0mi->SetParamP(tpc1);
4368           v0mi->Update(primvertex);
4369         }
4370       }
4371       continue;
4372     }
4373     CorrectForDeadZoneMaterial(&tpc0);
4374     CorrectForDeadZoneMaterial(&tpc1);
4375     //    if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4376     if (RefitAt(v0mi->GetRr(),&tpc0, v0mi->GetClusters(0)) && RefitAt(v0mi->GetRr(),&tpc1, v0mi->GetClusters(1))){
4377       v0temp.SetParamN(tpc0);
4378       v0temp.SetParamP(tpc1);
4379       v0temp.Update(primvertex);
4380       if (kFALSE) cstream<<"Refit"<<
4381         "V0.="<<v0mi<<
4382         "V0refit.="<<&v0temp<<
4383         "Tr0.="<<&tpc0<<
4384         "Tr1.="<<&tpc1<<
4385         "\n";
4386       //Bo:      if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4387       if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4388         v0mi->SetParamN(tpc0);
4389         v0mi->SetParamP(tpc1);
4390         v0mi->Update(primvertex);
4391       }
4392     }
4393   }
4394 }
4395
4396
4397
4398
4399
4400
4401