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