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;
1448   if (fZOffset<0)            // old geometry
1449     dphi = -(phi-fPhiOffset);
1450   else                       // new geometry
1451     dphi = phi-fPhiOffset;
1452
1453   if      (dphi <  0) dphi += 2*TMath::Pi();
1454   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1455   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1456   if (np>=fNladders) np-=fNladders;
1457   if (np<0)          np+=fNladders;
1458
1459   Double_t dz=fZOffset-z;
1460   Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
1461   if (nz>=fNdetectors) return -1;
1462   if (nz<0)            return -1;
1463
1464   return np*fNdetectors + nz;
1465 }
1466
1467
1468 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci){
1469   //--------------------------------------------------------------------
1470   // This function returns clusters within the "window"
1471   //--------------------------------------------------------------------
1472
1473   if (fCurrentSlice<0){
1474     Double_t rpi2 = 2.*fR*TMath::Pi();
1475     for (Int_t i=fI; i<fImax; i++) {
1476       Double_t y = fY[i];
1477       if (fYmax<y) y -= rpi2;
1478       if (fYmin>y) y += rpi2;
1479       if (y<fYmin) continue;
1480       if (y>fYmax) continue;
1481       if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1482       ci=i;
1483       fI=i+1;
1484       return fClusters[i];
1485     }
1486   }
1487   else{
1488     for (Int_t i=fI; i<fImax; i++) {
1489       if (fYcs[i]<fYmin) continue;
1490       if (fYcs[i]>fYmax) continue;
1491       if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1492       ci=fClusterIndexCs[i];
1493       fI=i+1;
1494       return fClustersCs[i];
1495     }
1496   }
1497   return 0;
1498 }
1499
1500
1501
1502 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1503 const {
1504   //--------------------------------------------------------------------
1505   //This function returns the layer thickness at this point (units X0)
1506   //--------------------------------------------------------------------
1507   Double_t d=0.0085;
1508   x0=21.82;
1509   if (43<fR&&fR<45) { //SSD2
1510      Double_t dd=0.0034;
1511      d=dd;
1512      if (TMath::Abs(y-0.00)>3.40) d+=dd;
1513      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1514      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1515      for (Int_t i=0; i<12; i++) {
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.9*(i+0.5))<0.15) {
1522           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1523           d+=0.0034;
1524           break;
1525        }
1526        if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1527        if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1528      }
1529   } else
1530   if (37<fR&&fR<41) { //SSD1
1531      Double_t dd=0.0034;
1532      d=dd;
1533      if (TMath::Abs(y-0.00)>3.40) d+=dd;
1534      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1535      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1536      for (Int_t i=0; i<11; i++) {
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+3.9*i)<0.15) {
1543           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1544           d+=dd;
1545           break;
1546        }
1547        if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1548        if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1549      }
1550   } else
1551   if (13<fR&&fR<26) { //SDD
1552      Double_t dd=0.0033;
1553      d=dd;
1554      if (TMath::Abs(y-0.00)>3.30) d+=dd;
1555
1556      if (TMath::Abs(y-1.80)<0.55) {
1557         d+=0.016;
1558         for (Int_t j=0; j<20; j++) {
1559           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1560           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1561         }
1562      }
1563      if (TMath::Abs(y+1.80)<0.55) {
1564         d+=0.016;
1565         for (Int_t j=0; j<20; j++) {
1566           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1567           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1568         }
1569      }
1570
1571      for (Int_t i=0; i<4; i++) {
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        if (TMath::Abs(z+7.3*i)<0.60) {
1578           d+=dd;
1579           if (TMath::Abs(y-0.00)>3.30) d+=dd;
1580           break;
1581        }
1582      }
1583   } else
1584   if (6<fR&&fR<8) {   //SPD2
1585      Double_t dd=0.0063; x0=21.5;
1586      d=dd;
1587      if (TMath::Abs(y-3.08)>0.5) d+=dd;
1588      //if (TMath::Abs(y-3.08)>0.45) d+=dd;
1589      if (TMath::Abs(y-3.03)<0.10) {d+=0.014;}
1590   } else
1591   if (3<fR&&fR<5) {   //SPD1
1592      Double_t dd=0.0063; x0=21.5;
1593      d=dd;
1594      if (TMath::Abs(y+0.21)>0.6) d+=dd;
1595      //if (TMath::Abs(y+0.21)>0.45) d+=dd;
1596      if (TMath::Abs(y+0.10)<0.10) {d+=0.014;}
1597   }
1598
1599   return d;
1600 }
1601
1602 Double_t AliITStrackerMI::GetEffectiveThickness(Double_t y,Double_t z) const
1603 {
1604   //--------------------------------------------------------------------
1605   //Returns the thickness between the current layer and the vertex (units X0)
1606   //--------------------------------------------------------------------
1607   Double_t d=0.0028*3*3; //beam pipe
1608   Double_t x0=0;
1609
1610   Double_t xn=fgLayers[fI].GetR();
1611   for (Int_t i=0; i<fI; i++) {
1612     Double_t xi=fgLayers[i].GetR();
1613     d+=fgLayers[i].GetThickness(y,z,x0)*xi*xi;
1614   }
1615
1616   if (fI>1) {
1617     Double_t xi=9.;
1618     d+=0.0097*xi*xi;
1619   }
1620
1621   if (fI>3) {
1622     Double_t xi=0.5*(fgLayers[3].GetR()+fgLayers[4].GetR());
1623     d+=0.0034*xi*xi;
1624   }
1625
1626   return d/(xn*xn);
1627 }
1628
1629 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
1630   //--------------------------------------------------------------------
1631   // This function returns number of clusters within the "window"
1632   //--------------------------------------------------------------------
1633   Int_t ncl=0;
1634   for (Int_t i=fI; i<fN; i++) {
1635     const AliITSRecPoint *c=fClusters[i];
1636     if (c->GetZ() > fZmax) break;
1637     if (c->IsUsed()) continue;
1638     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
1639     Double_t y=fR*det.GetPhi() + c->GetY();
1640
1641     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
1642     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
1643
1644     if (y<fYmin) continue;
1645     if (y>fYmax) continue;
1646     ncl++;
1647   }
1648   return ncl;
1649 }
1650
1651 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *t,
1652                                 const AliITStrackMI *c, Bool_t extra) {
1653   //--------------------------------------------------------------------
1654   // This function refits the track "t" at the position "x" using
1655   // the clusters from "c"
1656   // If "extra"==kTRUE,
1657   //    the clusters from overlapped modules get attached to "t"
1658   //--------------------------------------------------------------------
1659   Int_t index[kMaxLayer];
1660   Int_t k;
1661   for (k=0; k<kMaxLayer; k++) index[k]=-1;
1662   Int_t nc=c->GetNumberOfClusters();
1663   for (k=0; k<nc; k++) {
1664     Int_t idx=c->GetClusterIndex(k),nl=(idx&0xf0000000)>>28;
1665     index[nl]=idx;
1666   }
1667
1668   Int_t from, to, step;
1669   if (xx > t->GetX()) {
1670       from=0; to=kMaxLayer;
1671       step=+1;
1672   } else {
1673       from=kMaxLayer-1; to=-1;
1674       step=-1;
1675   }
1676
1677   for (Int_t i=from; i != to; i += step) {
1678      AliITSlayer &layer=fgLayers[i];
1679      Double_t r=layer.GetR();
1680
1681      {
1682      Double_t hI=i-0.5*step;
1683      if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) {
1684         Double_t rs=0.5*(fgLayers[i-step].GetR() + r);
1685         Double_t d=0.0034, x0=38.6;
1686         if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;}
1687         if (!t->PropagateTo(rs,-step*d,x0)) {
1688           return kFALSE;
1689         }
1690      }
1691      }
1692
1693      // remember old position [SR, GSI 18.02.2003]
1694      Double_t oldX=0., oldY=0., oldZ=0.;
1695      if (t->IsStartedTimeIntegral() && step==1) {
1696         t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ);
1697      }
1698      //
1699
1700      Double_t x,y,z;
1701      if (!t->GetGlobalXYZat(r,x,y,z)) {
1702        return kFALSE;
1703      }
1704      Double_t phi=TMath::ATan2(y,x);
1705      Int_t idet=layer.FindDetectorIndex(phi,z);
1706      if (idet<0) {
1707        return kFALSE;
1708      }
1709      const AliITSdetector &det=layer.GetDetector(idet);
1710      phi=det.GetPhi();
1711      if (!t->Propagate(phi,det.GetR())) {
1712        return kFALSE;
1713      }
1714      t->SetDetectorIndex(idet);
1715
1716      const AliITSRecPoint *cl=0;
1717      Double_t maxchi2=1000.*kMaxChi2;
1718
1719      Int_t idx=index[i];
1720      if (idx>0) {
1721         const AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(idx);
1722         if (c){
1723           if (idet != c->GetDetectorIndex()) {
1724             idet=c->GetDetectorIndex();
1725             const AliITSdetector &det=layer.GetDetector(idet);
1726             if (!t->Propagate(det.GetPhi(),det.GetR())) {
1727               return kFALSE;
1728             }
1729             t->SetDetectorIndex(idet);
1730           }
1731           //Double_t chi2=t->GetPredictedChi2(c);
1732           Int_t layer = (idx & 0xf0000000) >> 28;;
1733           Double_t chi2=GetPredictedChi2MI(t,c,layer);
1734           if (chi2<maxchi2) {
1735             cl=c;
1736             maxchi2=chi2;
1737           } else {
1738             return kFALSE;
1739           }
1740         }
1741      }
1742
1743      if (cl) {
1744        //if (!t->Update(cl,maxchi2,idx)) {
1745        if (!UpdateMI(t,cl,maxchi2,idx)) {
1746           return kFALSE;
1747        }
1748        t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1);
1749      }
1750
1751      {
1752      Double_t x0;
1753      Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0);
1754      t->CorrectForMaterial(-step*d,x0);
1755      }
1756
1757      if (extra) { //search for extra clusters
1758         AliITStrackV2 tmp(*t);
1759         Double_t dz=4*TMath::Sqrt(tmp.GetSigmaZ2()+kSigmaZ2[i]);
1760         if (dz < 0.5*TMath::Abs(tmp.GetTgl())) dz=0.5*TMath::Abs(tmp.GetTgl());
1761         Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
1762         if (dy < 0.5*TMath::Abs(tmp.GetSnp())) dy=0.5*TMath::Abs(tmp.GetSnp());
1763         Double_t zmin=t->GetZ() - dz;
1764         Double_t zmax=t->GetZ() + dz;
1765         Double_t ymin=t->GetY() + phi*r - dy;
1766         Double_t ymax=t->GetY() + phi*r + dy;
1767         layer.SelectClusters(zmin,zmax,ymin,ymax);
1768
1769         const AliITSRecPoint *c=0; Int_t ci=-1,cci=-1;
1770         Double_t maxchi2=1000.*kMaxChi2, tolerance=0.1;
1771         while ((c=layer.GetNextCluster(ci))!=0) {
1772            if (idet == c->GetDetectorIndex()) continue;
1773
1774            const AliITSdetector &det=layer.GetDetector(c->GetDetectorIndex());
1775
1776            if (!tmp.Propagate(det.GetPhi(),det.GetR())) continue;
1777
1778            if (TMath::Abs(tmp.GetZ() - c->GetZ()) > tolerance) continue;
1779            if (TMath::Abs(tmp.GetY() - c->GetY()) > tolerance) continue;
1780
1781            Double_t chi2=tmp.GetPredictedChi2(c);
1782            if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
1783         }
1784         if (cci>=0) t->SetExtraCluster(i,(i<<28)+cci);
1785      }
1786
1787      // track time update [SR, GSI 17.02.2003]
1788      if (t->IsStartedTimeIntegral() && step==1) {
1789         Double_t newX, newY, newZ;
1790         t->GetGlobalXYZat(t->GetX(),newX,newY,newZ);
1791         Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
1792                        (oldZ-newZ)*(oldZ-newZ);
1793         t->AddTimeStep(TMath::Sqrt(dL2));
1794      }
1795      //
1796
1797   }
1798
1799   if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
1800   return kTRUE;
1801 }
1802
1803 Bool_t
1804 AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *t,const Int_t *clindex) {
1805   //--------------------------------------------------------------------
1806   // This function refits the track "t" at the position "x" using
1807   // the clusters from array
1808   //--------------------------------------------------------------------
1809   Int_t index[kMaxLayer];
1810   Int_t k;
1811   for (k=0; k<kMaxLayer; k++) index[k]=-1;
1812   //
1813   for (k=0; k<kMaxLayer; k++) {
1814     index[k]=clindex[k];
1815   }
1816
1817   Int_t from, to, step;
1818   if (xx > t->GetX()) {
1819       from=0; to=kMaxLayer;
1820       step=+1;
1821   } else {
1822       from=kMaxLayer-1; to=-1;
1823       step=-1;
1824   }
1825
1826   for (Int_t i=from; i != to; i += step) {
1827      AliITSlayer &layer=fgLayers[i];
1828      Double_t r=layer.GetR();
1829      if (step<0 && xx>r) break;  //
1830      {
1831      Double_t hI=i-0.5*step;
1832      if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) {
1833         Double_t rs=0.5*(fgLayers[i-step].GetR() + r);
1834         Double_t d=0.0034, x0=38.6;
1835         if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;}
1836         if (!t->PropagateTo(rs,-step*d,x0)) {
1837           return kFALSE;
1838         }
1839      }
1840      }
1841
1842      // remember old position [SR, GSI 18.02.2003]
1843      Double_t oldX=0., oldY=0., oldZ=0.;
1844      if (t->IsStartedTimeIntegral() && step==1) {
1845         t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ);
1846      }
1847      //
1848
1849      Double_t x,y,z;
1850      if (!t->GetGlobalXYZat(r,x,y,z)) {
1851        return kFALSE;
1852      }
1853      Double_t phi=TMath::ATan2(y,x);
1854      Int_t idet=layer.FindDetectorIndex(phi,z);
1855      if (idet<0) {
1856        return kFALSE;
1857      }
1858      const AliITSdetector &det=layer.GetDetector(idet);
1859      phi=det.GetPhi();
1860      if (!t->Propagate(phi,det.GetR())) {
1861        return kFALSE;
1862      }
1863      t->SetDetectorIndex(idet);
1864
1865      const AliITSRecPoint *cl=0;
1866      Double_t maxchi2=1000.*kMaxChi2;
1867
1868      Int_t idx=index[i];
1869      if (idx>0) {
1870         const AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(idx);
1871         if (c){
1872           if (idet != c->GetDetectorIndex()) {
1873             idet=c->GetDetectorIndex();
1874             const AliITSdetector &det=layer.GetDetector(idet);
1875             if (!t->Propagate(det.GetPhi(),det.GetR())) {
1876               return kFALSE;
1877             }
1878             t->SetDetectorIndex(idet);
1879           }
1880           //Double_t chi2=t->GetPredictedChi2(c);
1881           Int_t layer = (idx & 0xf0000000) >> 28;;
1882           Double_t chi2=GetPredictedChi2MI(t,c,layer);
1883           if (chi2<maxchi2) {
1884             cl=c;
1885             maxchi2=chi2;
1886           } else {
1887             return kFALSE;
1888           }
1889         }
1890      }
1891      /*
1892      if (cl==0)
1893      if (t->GetNumberOfClusters()>2) {
1894         Double_t dz=4*TMath::Sqrt(t->GetSigmaZ2()+kSigmaZ2[i]);
1895         Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
1896         Double_t zmin=t->GetZ() - dz;
1897         Double_t zmax=t->GetZ() + dz;
1898         Double_t ymin=t->GetY() + phi*r - dy;
1899         Double_t ymax=t->GetY() + phi*r + dy;
1900         layer.SelectClusters(zmin,zmax,ymin,ymax);
1901
1902         const AliITSRecPoint *c=0; Int_t ci=-1;
1903         while ((c=layer.GetNextCluster(ci))!=0) {
1904            if (idet != c->GetDetectorIndex()) continue;
1905            Double_t chi2=t->GetPredictedChi2(c);
1906            if (chi2<maxchi2) { cl=c; maxchi2=chi2; idx=ci; }
1907         }
1908      }
1909      */
1910      if (cl) {
1911        //if (!t->Update(cl,maxchi2,idx)) {
1912        if (!UpdateMI(t,cl,maxchi2,idx)) {
1913           return kFALSE;
1914        }
1915        t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1);
1916      }
1917
1918      {
1919      Double_t x0;
1920      Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0);
1921      t->CorrectForMaterial(-step*d,x0);
1922      }
1923
1924      // track time update [SR, GSI 17.02.2003]
1925      if (t->IsStartedTimeIntegral() && step==1) {
1926         Double_t newX, newY, newZ;
1927         t->GetGlobalXYZat(t->GetX(),newX,newY,newZ);
1928         Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
1929                        (oldZ-newZ)*(oldZ-newZ);
1930         t->AddTimeStep(TMath::Sqrt(dL2));
1931      }
1932      //
1933
1934   }
1935
1936   if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
1937   return kTRUE;
1938 }
1939
1940 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
1941 {
1942   //
1943   // calculate normalized chi2
1944   //  return NormalizedChi2(track,0);
1945   Float_t chi2 = 0;
1946   Float_t sum=0;
1947   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
1948   //  track->fdEdxMismatch=0;
1949   Float_t dedxmismatch =0;
1950   Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
1951   if (mode<100){
1952     for (Int_t i = 0;i<6;i++){
1953       if (track->GetClIndex(i)>0){
1954         Float_t cerry, cerrz;
1955         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
1956         else
1957           { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
1958         cerry*=cerry;
1959         cerrz*=cerrz;
1960         Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
1961         if (i>1){
1962           Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
1963           if (ratio<0.5) {
1964             cchi2+=(0.5-ratio)*10.;
1965             //track->fdEdxMismatch+=(0.5-ratio)*10.;
1966             dedxmismatch+=(0.5-ratio)*10.;
1967           }
1968         }
1969         if (i<2 ||i>3){
1970           AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
1971           Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
1972           if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
1973           if (i<2) chi2+=2*cl->GetDeltaProbability();
1974         }
1975         chi2+=cchi2;
1976         sum++;
1977       }
1978     }
1979     if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
1980       track->SetdEdxMismatch(dedxmismatch);
1981     }
1982   }
1983   else{
1984     for (Int_t i = 0;i<4;i++){
1985       if (track->GetClIndex(i)>0){
1986         Float_t cerry, cerrz;
1987         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
1988         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
1989         cerry*=cerry;
1990         cerrz*=cerrz;
1991         chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
1992         chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
1993         sum++;
1994       }
1995     }
1996     for (Int_t i = 4;i<6;i++){
1997       if (track->GetClIndex(i)>0){
1998         Float_t cerry, cerrz;
1999         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2000         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2001         cerry*=cerry;
2002         cerrz*=cerrz;
2003         Float_t cerryb, cerrzb;
2004         if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2005         else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2006         cerryb*=cerryb;
2007         cerrzb*=cerrzb;
2008         chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2009         chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2010         sum++;
2011       }
2012     }
2013   }
2014   if (track->GetESDtrack()->GetTPCsignal()>85){
2015     Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2016     if (ratio<0.5) {
2017       chi2+=(0.5-ratio)*5.;
2018     }
2019     if (ratio>2){
2020       chi2+=(ratio-2.0)*3;
2021     }
2022   }
2023   //
2024   Double_t match = TMath::Sqrt(track->GetChi22());
2025   if (track->GetConstrain())  match/=track->GetNumberOfClusters();
2026   if (!track->GetConstrain()) match/=track->GetNumberOfClusters()-2.;
2027   if (match<0) match=0;
2028   Float_t deadzonefactor = (track->GetNDeadZone()>0) ? 3*(1.1-track->GetDeadZoneProbability()):0.;
2029   Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2030     (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2031                                 1./(1.+track->GetNSkipped()));
2032
2033  return normchi2;
2034 }
2035
2036
2037 Double_t AliITStrackerMI::GetMatchingChi2(AliITStrackMI * track1, AliITStrackMI * track2)
2038 {
2039   //
2040   // return matching chi2 between two tracks
2041   AliITStrackMI track3(*track2);
2042   track3.Propagate(track1->GetAlpha(),track1->GetX());
2043   TMatrixD vec(5,1);
2044   vec(0,0)=track1->GetY()   - track3.GetY();
2045   vec(1,0)=track1->GetZ()   - track3.GetZ();
2046   vec(2,0)=track1->GetSnp() - track3.GetSnp();
2047   vec(3,0)=track1->GetTgl() - track3.GetTgl();
2048   vec(4,0)=track1->Get1Pt() - track3.Get1Pt();
2049   //
2050   TMatrixD cov(5,5);
2051   cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2052   cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2053   cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2054   cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2055   cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2056
2057   cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2058   cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2059   cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2060   cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2061   //
2062   cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2063   cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2064   cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2065   //
2066   cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2067   cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2068   //
2069   cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2070
2071   cov.Invert();
2072   TMatrixD vec2(cov,TMatrixD::kMult,vec);
2073   TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2074   return chi2(0,0);
2075 }
2076
2077 Double_t  AliITStrackerMI::GetDeadZoneProbability(Double_t zpos, Double_t zerr)
2078 {
2079   //
2080   //  return probability that given point - characterized by z position and error  is in dead zone
2081   //
2082   Double_t probability =0;
2083   Double_t absz = TMath::Abs(zpos);
2084   Double_t nearestz = (absz<2)? 0.:7.1;
2085   if (TMath::Abs(absz-nearestz)>0.25+3*zerr) return 0;
2086   Double_t zmin=0, zmax=0;
2087   if (zpos<-6.){
2088     zmin = -7.25; zmax = -6.95;
2089   }
2090   if (zpos>6){
2091     zmin = 7.0; zmax =7.3;
2092   }
2093   if (absz<2){
2094     zmin = -0.75; zmax = 1.5;
2095   }
2096   probability = (TMath::Erf((zpos-zmin)/zerr) - TMath::Erf((zpos-zmax)/zerr))*0.5;
2097   return probability;
2098 }
2099
2100
2101 Double_t AliITStrackerMI::GetTruncatedChi2(AliITStrackMI * track, Float_t fac)
2102 {
2103   //
2104   // calculate normalized chi2
2105   Float_t chi2[6];
2106   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2107   Float_t ncl = 0;
2108   for (Int_t i = 0;i<6;i++){
2109     if (TMath::Abs(track->GetDy(i))>0){
2110       chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2111       chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2112       ncl++;
2113     }
2114     else{chi2[i]=10000;}
2115   }
2116   Int_t index[6];
2117   TMath::Sort(6,chi2,index,kFALSE);
2118   Float_t max = float(ncl)*fac-1.;
2119   Float_t sumchi=0, sumweight=0;
2120   for (Int_t i=0;i<max+1;i++){
2121     Float_t weight = (i<max)?1.:(max+1.-i);
2122     sumchi+=weight*chi2[index[i]];
2123     sumweight+=weight;
2124   }
2125   Double_t normchi2 = sumchi/sumweight;
2126   return normchi2;
2127 }
2128
2129
2130 Double_t AliITStrackerMI::GetInterpolatedChi2(AliITStrackMI * forwardtrack, AliITStrackMI * backtrack)
2131 {
2132   //
2133   // calculate normalized chi2
2134   //  if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2135   Int_t npoints = 0;
2136   Double_t res =0;
2137   for (Int_t i=0;i<6;i++){
2138     if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2139     Double_t sy1 = forwardtrack->GetSigmaY(i);
2140     Double_t sz1 = forwardtrack->GetSigmaZ(i);
2141     Double_t sy2 = backtrack->GetSigmaY(i);
2142     Double_t sz2 = backtrack->GetSigmaZ(i);
2143     if (i<2){ sy2=1000.;sz2=1000;}
2144     //
2145     Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2146     Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2147     //
2148     Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2149     Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2150     //
2151     res+= nz0*nz0+ny0*ny0;
2152     npoints++;
2153   }
2154   if (npoints>1) return
2155                    TMath::Max(TMath::Abs(0.3*forwardtrack->Get1Pt())-0.5,0.)+
2156                    //2*forwardtrack->fNUsed+
2157                    res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2158                                   1./(1.+forwardtrack->GetNSkipped()));
2159   return 1000;
2160 }
2161
2162
2163
2164
2165
2166 Float_t  *AliITStrackerMI::GetWeight(Int_t index) {
2167   //--------------------------------------------------------------------
2168   //       Return pointer to a given cluster
2169   //--------------------------------------------------------------------
2170   Int_t l=(index & 0xf0000000) >> 28;
2171   Int_t c=(index & 0x0fffffff) >> 00;
2172   return fgLayers[l].GetWeight(c);
2173 }
2174
2175 void AliITStrackerMI::RegisterClusterTracks(AliITStrackMI* track,Int_t id)
2176 {
2177   //---------------------------------------------
2178   // register track to the list
2179   //
2180   if (track->GetESDtrack()->GetKinkIndex(0)!=0) return;  //don't register kink tracks
2181   //
2182   //
2183   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2184     Int_t index = track->GetClusterIndex(icluster);
2185     Int_t l=(index & 0xf0000000) >> 28;
2186     Int_t c=(index & 0x0fffffff) >> 00;
2187     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2188     for (Int_t itrack=0;itrack<4;itrack++){
2189       if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2190         fgLayers[l].SetClusterTracks(itrack,c,id);
2191         break;
2192       }
2193     }
2194   }
2195 }
2196 void AliITStrackerMI::UnRegisterClusterTracks(AliITStrackMI* track, Int_t id)
2197 {
2198   //---------------------------------------------
2199   // unregister track from the list
2200   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2201     Int_t index = track->GetClusterIndex(icluster);
2202     Int_t l=(index & 0xf0000000) >> 28;
2203     Int_t c=(index & 0x0fffffff) >> 00;
2204     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2205     for (Int_t itrack=0;itrack<4;itrack++){
2206       if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2207         fgLayers[l].SetClusterTracks(itrack,c,-1);
2208       }
2209     }
2210   }
2211 }
2212 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2213 {
2214   //-------------------------------------------------------------
2215   //get number of shared clusters
2216   //-------------------------------------------------------------
2217   Float_t shared=0;
2218   for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2219   // mean  number of clusters
2220   Float_t *ny = GetNy(id), *nz = GetNz(id);
2221
2222
2223   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2224     Int_t index = track->GetClusterIndex(icluster);
2225     Int_t l=(index & 0xf0000000) >> 28;
2226     Int_t c=(index & 0x0fffffff) >> 00;
2227     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2228     if (ny[l]==0){
2229       printf("problem\n");
2230     }
2231     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2232     Float_t weight=1;
2233     //
2234     Float_t deltan = 0;
2235     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2236     if (l>2&&track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2237     if (l<2 || l>3){
2238       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2239     }
2240     else{
2241       deltan = (cl->GetNz()-nz[l]);
2242     }
2243     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2244     weight = 2./TMath::Max(3.+deltan,2.);
2245     //
2246     for (Int_t itrack=0;itrack<4;itrack++){
2247       if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2248         list[l]=index;
2249         clist[l] = (AliITSRecPoint*)GetCluster(index);
2250         shared+=weight;
2251         break;
2252       }
2253     }
2254   }
2255   track->SetNUsed(shared);
2256   return shared;
2257 }
2258
2259 Int_t AliITStrackerMI::GetOverlapTrack(AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2260 {
2261   //
2262   // find first shared track
2263   //
2264   // mean  number of clusters
2265   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2266   //
2267   for (Int_t i=0;i<6;i++) overlist[i]=-1;
2268   Int_t sharedtrack=100000;
2269   Int_t tracks[24],trackindex=0;
2270   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2271   //
2272   for (Int_t icluster=0;icluster<6;icluster++){
2273     if (clusterlist[icluster]<0) continue;
2274     Int_t index = clusterlist[icluster];
2275     Int_t l=(index & 0xf0000000) >> 28;
2276     Int_t c=(index & 0x0fffffff) >> 00;
2277     if (ny[l]==0){
2278       printf("problem\n");
2279     }
2280     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2281     //if (l>3) continue;
2282     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2283     //
2284     Float_t deltan = 0;
2285     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2286     if (l>2&&track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2287     if (l<2 || l>3){
2288       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2289     }
2290     else{
2291       deltan = (cl->GetNz()-nz[l]);
2292     }
2293     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2294     //
2295     for (Int_t itrack=3;itrack>=0;itrack--){
2296       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2297       if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2298        tracks[trackindex]  = fgLayers[l].GetClusterTracks(itrack,c);
2299        trackindex++;
2300       }
2301     }
2302   }
2303   if (trackindex==0) return -1;
2304   if (trackindex==1){
2305     sharedtrack = tracks[0];
2306   }else{
2307     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2308     else{
2309       //
2310       Int_t track[24], cluster[24];
2311       for (Int_t i=0;i<trackindex;i++){ track[i]=-1; cluster[i]=0;}
2312       Int_t index =0;
2313       //
2314       for (Int_t i=0;i<trackindex;i++){
2315         if (tracks[i]<0) continue;
2316         track[index] = tracks[i];
2317         cluster[index]++;
2318         for (Int_t j=i+1;j<trackindex;j++){
2319           if (tracks[j]<0) continue;
2320           if (tracks[j]==tracks[i]){
2321             cluster[index]++;
2322             tracks[j]=-1;
2323           }
2324         }
2325         index++;
2326       }
2327       Int_t max=0;
2328       for (Int_t i=0;i<index;i++){
2329         if (cluster[index]>max) {
2330           sharedtrack=track[index];
2331           max=cluster[index];
2332         }
2333       }
2334     }
2335   }
2336
2337   if (sharedtrack>=100000) return -1;
2338   //
2339   // make list of overlaps
2340   shared =0;
2341   for (Int_t icluster=0;icluster<6;icluster++){
2342     if (clusterlist[icluster]<0) continue;
2343     Int_t index = clusterlist[icluster];
2344     Int_t l=(index & 0xf0000000) >> 28;
2345     Int_t c=(index & 0x0fffffff) >> 00;
2346     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2347     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2348     if (l==0 || l==1){
2349       if (cl->GetNy()>2) continue;
2350       if (cl->GetNz()>2) continue;
2351     }
2352     if (l==4 || l==5){
2353       if (cl->GetNy()>3) continue;
2354       if (cl->GetNz()>3) continue;
2355     }
2356     //
2357     for (Int_t itrack=3;itrack>=0;itrack--){
2358       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2359       if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2360         overlist[l]=index;
2361         shared++;
2362       }
2363     }
2364   }
2365   return sharedtrack;
2366 }
2367
2368
2369 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2370   //
2371   // try to find track hypothesys without conflicts
2372   // with minimal chi2;
2373   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2374   Int_t entries1 = arr1->GetEntriesFast();
2375   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2376   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2377   Int_t entries2 = arr2->GetEntriesFast();
2378   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2379   //
2380   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2381   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2382   if (TMath::Abs(1./track10->Get1Pt())>0.5+TMath::Abs(1/track20->Get1Pt())) return track10;
2383
2384   for (Int_t itrack=0;itrack<entries1;itrack++){
2385     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2386     UnRegisterClusterTracks(track,trackID1);
2387   }
2388   //
2389   for (Int_t itrack=0;itrack<entries2;itrack++){
2390     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2391     UnRegisterClusterTracks(track,trackID2);
2392   }
2393   Int_t index1=0;
2394   Int_t index2=0;
2395   Float_t maxconflicts=6;
2396   Double_t maxchi2 =1000.;
2397   //
2398   // get weight of hypothesys - using best hypothesy
2399   Double_t w1,w2;
2400
2401   Int_t list1[6],list2[6];
2402   AliITSRecPoint *clist1[6], *clist2[6] ;
2403   RegisterClusterTracks(track10,trackID1);
2404   RegisterClusterTracks(track20,trackID2);
2405   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2406   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2407   UnRegisterClusterTracks(track10,trackID1);
2408   UnRegisterClusterTracks(track20,trackID2);
2409   //
2410   // normalized chi2
2411   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2412   Float_t nerry[6],nerrz[6];
2413   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2414   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2415   for (Int_t i=0;i<6;i++){
2416      if ( (erry1[i]>0) && (erry2[i]>0)) {
2417        nerry[i] = TMath::Min(erry1[i],erry2[i]);
2418        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2419      }else{
2420        nerry[i] = TMath::Max(erry1[i],erry2[i]);
2421        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2422      }
2423      if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2424        chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2425        chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2426        ncl1++;
2427      }
2428      if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2429        chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2430        chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2431        ncl2++;
2432      }
2433   }
2434   chi21/=ncl1;
2435   chi22/=ncl2;
2436   //
2437   //
2438   Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2439   Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2440   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2441   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2442   //
2443   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2444         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2445         +1.*TMath::Abs(1./track10->Get1Pt())/(TMath::Abs(1./track10->Get1Pt())+TMath::Abs(1./track20->Get1Pt()))
2446         );
2447   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2448         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2449         +1.*TMath::Abs(1./track20->Get1Pt())/(TMath::Abs(1./track10->Get1Pt())+TMath::Abs(1./track20->Get1Pt()))
2450         );
2451
2452   Double_t sumw = w1+w2;
2453   w1/=sumw;
2454   w2/=sumw;
2455   if (w1<w2*0.5) {
2456     w1 = (d2+0.5)/(d1+d2+1.);
2457     w2 = (d1+0.5)/(d1+d2+1.);
2458   }
2459   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2460   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2461   //
2462   // get pair of "best" hypothesys
2463   //
2464   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2465   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2466
2467   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2468     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2469     //if (track1->fFakeRatio>0) continue;
2470     RegisterClusterTracks(track1,trackID1);
2471     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2472       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2473
2474       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2475       //if (track2->fFakeRatio>0) continue;
2476       Float_t nskipped=0;
2477       RegisterClusterTracks(track2,trackID2);
2478       Int_t list1[6],list2[6];
2479       AliITSRecPoint *clist1[6], *clist2[6] ;
2480       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2481       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2482       UnRegisterClusterTracks(track2,trackID2);
2483       //
2484       if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2485       if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2486       if (nskipped>0.5) continue;
2487       //
2488       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2489       if (conflict1+1<cconflict1) continue;
2490       if (conflict2+1<cconflict2) continue;
2491       Float_t conflict=0;
2492       Float_t sumchi2=0;
2493       Float_t sum=0;
2494       for (Int_t i=0;i<6;i++){
2495         //
2496         Float_t c1 =0.; // conflict coeficients
2497         Float_t c2 =0.;
2498         if (clist1[i]&&clist2[i]){
2499           Float_t deltan = 0;
2500           if (i<2 || i>3){
2501             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2502           }
2503           else{
2504             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2505           }
2506           c1  = 2./TMath::Max(3.+deltan,2.);
2507           c2  = 2./TMath::Max(3.+deltan,2.);
2508         }
2509         else{
2510           if (clist1[i]){
2511             Float_t deltan = 0;
2512             if (i<2 || i>3){
2513               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
2514             }
2515             else{
2516               deltan = (clist1[i]->GetNz()-nz1[i]);
2517             }
2518             c1  = 2./TMath::Max(3.+deltan,2.);
2519             c2  = 0;
2520           }
2521
2522           if (clist2[i]){
2523             Float_t deltan = 0;
2524             if (i<2 || i>3){
2525               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
2526             }
2527             else{
2528               deltan = (clist2[i]->GetNz()-nz2[i]);
2529             }
2530             c2  = 2./TMath::Max(3.+deltan,2.);
2531             c1  = 0;
2532           }
2533         }
2534         //
2535         Double_t chi21=0,chi22=0;
2536         if (TMath::Abs(track1->GetDy(i))>0.) {
2537           chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
2538             (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
2539           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
2540           //  (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
2541         }else{
2542           if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
2543         }
2544         //
2545         if (TMath::Abs(track2->GetDy(i))>0.) {
2546           chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
2547             (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
2548           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
2549           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
2550         }
2551         else{
2552           if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
2553         }
2554         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
2555         if (chi21>0) sum+=w1;
2556         if (chi22>0) sum+=w2;
2557         conflict+=(c1+c2);
2558       }
2559       Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
2560       if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
2561       Double_t normchi2 = 2*conflict+sumchi2/sum;
2562       if ( normchi2 <maxchi2 ){
2563         index1 = itrack1;
2564         index2 = itrack2;
2565         maxconflicts = conflict;
2566         maxchi2 = normchi2;
2567       }
2568     }
2569     UnRegisterClusterTracks(track1,trackID1);
2570   }
2571   //
2572   //  if (maxconflicts<4 && maxchi2<th0){
2573   if (maxchi2<th0*2.){
2574     Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
2575     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
2576     track1->SetChi2MIP(5,maxconflicts);
2577     track1->SetChi2MIP(6,maxchi2);
2578     track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
2579     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
2580     track1->SetChi2MIP(8,index1);
2581     fBestTrackIndex[trackID1] =index1;
2582     UpdateESDtrack(track1, AliESDtrack::kITSin);
2583   }
2584   else if (track10->GetChi2MIP(0)<th1){
2585     track10->SetChi2MIP(5,maxconflicts);
2586     track10->SetChi2MIP(6,maxchi2);
2587     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
2588     UpdateESDtrack(track10,AliESDtrack::kITSin);
2589   }
2590
2591   for (Int_t itrack=0;itrack<entries1;itrack++){
2592     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2593     UnRegisterClusterTracks(track,trackID1);
2594   }
2595   //
2596   for (Int_t itrack=0;itrack<entries2;itrack++){
2597     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2598     UnRegisterClusterTracks(track,trackID2);
2599   }
2600
2601   if (track10->GetConstrain()&&track10->GetChi2MIP(0)<kMaxChi2PerCluster[0]&&track10->GetChi2MIP(1)<kMaxChi2PerCluster[1]
2602       &&track10->GetChi2MIP(2)<kMaxChi2PerCluster[2]&&track10->GetChi2MIP(3)<kMaxChi2PerCluster[3]){
2603     //  if (track10->fChi2MIP[0]<kMaxChi2PerCluster[0]&&track10->fChi2MIP[1]<kMaxChi2PerCluster[1]
2604   //    &&track10->fChi2MIP[2]<kMaxChi2PerCluster[2]&&track10->fChi2MIP[3]<kMaxChi2PerCluster[3]){
2605     RegisterClusterTracks(track10,trackID1);
2606   }
2607   if (track20->GetConstrain()&&track20->GetChi2MIP(0)<kMaxChi2PerCluster[0]&&track20->GetChi2MIP(1)<kMaxChi2PerCluster[1]
2608       &&track20->GetChi2MIP(2)<kMaxChi2PerCluster[2]&&track20->GetChi2MIP(3)<kMaxChi2PerCluster[3]){
2609     //if (track20->fChi2MIP[0]<kMaxChi2PerCluster[0]&&track20->fChi2MIP[1]<kMaxChi2PerCluster[1]
2610     //  &&track20->fChi2MIP[2]<kMaxChi2PerCluster[2]&&track20->fChi2MIP[3]<kMaxChi2PerCluster[3]){
2611     RegisterClusterTracks(track20,trackID2);
2612   }
2613   return track10;
2614
2615 }
2616
2617 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
2618   //--------------------------------------------------------------------
2619   // This function marks clusters assigned to the track
2620   //--------------------------------------------------------------------
2621   AliTracker::UseClusters(t,from);
2622
2623   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
2624   //if (c->GetQ()>2) c->Use();
2625   if (c->GetSigmaZ2()>0.1) c->Use();
2626   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
2627   //if (c->GetQ()>2) c->Use();
2628   if (c->GetSigmaZ2()>0.1) c->Use();
2629
2630 }
2631
2632
2633 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
2634 {
2635   //------------------------------------------------------------------
2636   // add track to the list of hypothesys
2637   //------------------------------------------------------------------
2638
2639   if (esdindex>=fTrackHypothesys.GetEntriesFast()) fTrackHypothesys.Expand(esdindex*2+10);
2640   //
2641   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2642   if (!array) {
2643     array = new TObjArray(10);
2644     fTrackHypothesys.AddAt(array,esdindex);
2645   }
2646   array->AddLast(track);
2647 }
2648
2649 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
2650 {
2651   //-------------------------------------------------------------------
2652   // compress array of track hypothesys
2653   // keep only maxsize best hypothesys
2654   //-------------------------------------------------------------------
2655   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
2656   if (! (fTrackHypothesys.At(esdindex)) ) return;
2657   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2658   Int_t entries = array->GetEntriesFast();
2659   //
2660   //- find preliminary besttrack as a reference
2661   Float_t minchi2=10000;
2662   Int_t maxn=0;
2663   AliITStrackMI * besttrack=0;
2664   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
2665     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2666     if (!track) continue;
2667     Float_t chi2 = NormalizedChi2(track,0);
2668     //
2669     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
2670     track->SetLabel(tpcLabel);
2671     CookdEdx(track);
2672     track->SetFakeRatio(1.);
2673     CookLabel(track,0.); //For comparison only
2674     //
2675     //if (chi2<kMaxChi2PerCluster[0]&&track->fFakeRatio==0){
2676     if (chi2<kMaxChi2PerCluster[0]){
2677       if (track->GetNumberOfClusters()<maxn) continue;
2678       maxn = track->GetNumberOfClusters();
2679       if (chi2<minchi2){
2680         minchi2=chi2;
2681         besttrack=track;
2682       }
2683     }
2684     else{
2685       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
2686         delete array->RemoveAt(itrack);
2687       }
2688     }
2689   }
2690   if (!besttrack) return;
2691   //
2692   //
2693   //take errors of best track as a reference
2694   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
2695   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
2696   for (Int_t i=0;i<6;i++) {
2697     if (besttrack->GetClIndex(i)>0){
2698       erry[i] = besttrack->GetSigmaY(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2699       errz[i] = besttrack->GetSigmaZ(i); errz[i+6] = besttrack->GetSigmaZ(i+6);
2700       ny[i]   = besttrack->GetNy(i);
2701       nz[i]   = besttrack->GetNz(i);
2702     }
2703   }
2704   //
2705   // calculate normalized chi2
2706   //
2707   Float_t * chi2        = new Float_t[entries];
2708   Int_t * index         = new Int_t[entries];
2709   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
2710   for (Int_t itrack=0;itrack<entries;itrack++){
2711     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2712     if (track){
2713       track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
2714       if (track->GetChi2MIP(0)<kMaxChi2PerCluster[0])
2715         chi2[itrack] = track->GetChi2MIP(0);
2716       else{
2717         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
2718           delete array->RemoveAt(itrack);
2719         }
2720       }
2721     }
2722   }
2723   //
2724   TMath::Sort(entries,chi2,index,kFALSE);
2725   besttrack = (AliITStrackMI*)array->At(index[0]);
2726   if (besttrack&&besttrack->GetChi2MIP(0)<kMaxChi2PerCluster[0]){
2727     for (Int_t i=0;i<6;i++){
2728       if (besttrack->GetClIndex(i)>0){
2729         erry[i] = besttrack->GetSigmaY(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2730         errz[i] = besttrack->GetSigmaZ(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2731         ny[i]   = besttrack->GetNy(i);
2732         nz[i]   = besttrack->GetNz(i);
2733       }
2734     }
2735   }
2736   //
2737   // calculate one more time with updated normalized errors
2738   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
2739   for (Int_t itrack=0;itrack<entries;itrack++){
2740     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2741     if (track){
2742       track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
2743       if (track->GetChi2MIP(0)<kMaxChi2PerCluster[0])
2744         chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
2745       else
2746         {
2747           if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
2748             delete array->RemoveAt(itrack);
2749           }
2750         }
2751     }
2752   }
2753   entries = array->GetEntriesFast();
2754   //
2755   //
2756   if (entries>0){
2757     TObjArray * newarray = new TObjArray();
2758     TMath::Sort(entries,chi2,index,kFALSE);
2759     besttrack = (AliITStrackMI*)array->At(index[0]);
2760     if (besttrack){
2761       //
2762       for (Int_t i=0;i<6;i++){
2763         if (besttrack->GetNz(i)>0&&besttrack->GetNy(i)>0){
2764           erry[i] = besttrack->GetSigmaY(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2765           errz[i] = besttrack->GetSigmaZ(i); errz[i+6] = besttrack->GetSigmaZ(i+6);
2766           ny[i]   = besttrack->GetNy(i);
2767           nz[i]   = besttrack->GetNz(i);
2768         }
2769       }
2770       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
2771       Float_t minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(kMaxChi2PerCluster[0]));
2772       Float_t minn = besttrack->GetNumberOfClusters()-3;
2773       Int_t accepted=0;
2774       for (Int_t i=0;i<entries;i++){
2775         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
2776         if (!track) continue;
2777         if (accepted>maxcut) break;
2778         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
2779         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
2780           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
2781             delete array->RemoveAt(index[i]);
2782             continue;
2783           }
2784         }
2785         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
2786         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
2787           if (!shortbest) accepted++;
2788           //
2789           newarray->AddLast(array->RemoveAt(index[i]));
2790           for (Int_t i=0;i<6;i++){
2791             if (nz[i]==0){
2792               erry[i] = track->GetSigmaY(i); erry[i+6] = track->GetSigmaY(i+6);
2793               errz[i] = track->GetSigmaZ(i); errz[i]   = track->GetSigmaZ(i+6);
2794               ny[i]   = track->GetNy(i);
2795               nz[i]   = track->GetNz(i);
2796             }
2797           }
2798         }
2799         else{
2800           delete array->RemoveAt(index[i]);
2801         }
2802       }
2803       array->Delete();
2804       delete fTrackHypothesys.RemoveAt(esdindex);
2805       fTrackHypothesys.AddAt(newarray,esdindex);
2806     }
2807     else{
2808       array->Delete();
2809       delete fTrackHypothesys.RemoveAt(esdindex);
2810     }
2811   }
2812   delete [] chi2;
2813   delete [] index;
2814 }
2815
2816
2817
2818 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
2819 {
2820   //-------------------------------------------------------------
2821   // try to find best hypothesy
2822   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
2823   //-------------------------------------------------------------
2824   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
2825   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2826   if (!array) return 0;
2827   Int_t entries = array->GetEntriesFast();
2828   if (!entries) return 0;
2829   Float_t minchi2 = 100000;
2830   AliITStrackMI * besttrack=0;
2831   //
2832   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
2833   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
2834   Double_t xyzv[]={GetX(),GetY(),GetZ()};
2835   Double_t ersv[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
2836   //
2837   for (Int_t i=0;i<entries;i++){
2838     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
2839     if (!track) continue;
2840     Float_t sigmarfi,sigmaz;
2841     GetDCASigma(track,sigmarfi,sigmaz);
2842     track->SetDnorm(0,sigmarfi);
2843     track->SetDnorm(1,sigmaz);
2844     //
2845     track->SetChi2MIP(1,1000000);
2846     track->SetChi2MIP(2,1000000);
2847     track->SetChi2MIP(3,1000000);
2848     //
2849     // backtrack
2850     backtrack = new(backtrack) AliITStrackMI(*track);
2851     if (track->GetConstrain()){
2852       if (!backtrack->PropagateTo(3.,0.0028,65.19)) continue;
2853       if (!backtrack->Improve(0,xyzv,ersv))         continue;
2854       //if (!backtrack->PropagateTo(2.,0.0028,0))     continue; // This
2855       //if (!backtrack->Improve(0,xyzv,ersv))         continue; // is
2856       //if (!backtrack->PropagateTo(1.,0.0028,0))     continue; // an over-kill
2857       //if (!backtrack->Improve(0,xyzv,ersv))         continue; //   (I.B.)
2858       //if (!backtrack->PropagateToVertex())          continue; //
2859       backtrack->ResetCovariance(10.);
2860       //if (!backtrack->Improve(0,xyzv,ersv))         continue;
2861     }else{
2862       backtrack->ResetCovariance(10.);
2863     }
2864     backtrack->ResetClusters();
2865
2866     Double_t x = original->GetX();
2867     if (!RefitAt(x,backtrack,track)) continue;
2868     //
2869     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
2870     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
2871     if (track->GetChi2MIP(1)>kMaxChi2PerCluster[1]*6.)  continue;
2872     track->SetChi22(GetMatchingChi2(backtrack,original));
2873
2874     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
2875     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
2876     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
2877
2878
2879     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>kMaxChi2PerCluster[1])  continue;
2880     Bool_t isOK=kTRUE;
2881     if(!isOK) continue;
2882     //
2883     //forward track - without constraint
2884     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
2885     forwardtrack->ResetClusters();
2886     x = track->GetX();
2887     RefitAt(x,forwardtrack,track);
2888     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
2889     if  (track->GetChi2MIP(2)>kMaxChi2PerCluster[2]*6.0)  continue;
2890     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>kMaxChi2PerCluster[2])  continue;
2891
2892     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
2893     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
2894     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
2895     forwardtrack->SetD(0,track->GetD(0));
2896     forwardtrack->SetD(1,track->GetD(1));
2897     {
2898       Int_t list[6];
2899       AliITSRecPoint* clist[6];
2900       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
2901       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
2902     }
2903
2904     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
2905     if  ( (track->GetChi2MIP(3)>6.*kMaxChi2PerCluster[3])) continue;
2906     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*kMaxChi2PerCluster[3])) {
2907       track->SetChi2MIP(3,1000);
2908       continue;
2909     }
2910     Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
2911     //
2912     for (Int_t ichi=0;ichi<5;ichi++){
2913       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
2914     }
2915     if (chi2 < minchi2){
2916       //besttrack = new AliITStrackMI(*forwardtrack);
2917       besttrack = track;
2918       besttrack->SetLabel(track->GetLabel());
2919       besttrack->SetFakeRatio(track->GetFakeRatio());
2920       minchi2   = chi2;
2921       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
2922       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
2923       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
2924     }
2925   }
2926   delete backtrack;
2927   delete forwardtrack;
2928   Int_t accepted=0;
2929   for (Int_t i=0;i<entries;i++){
2930     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
2931     if (!track) continue;
2932
2933     if (accepted>checkmax || track->GetChi2MIP(3)>kMaxChi2PerCluster[3]*6. ||
2934         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
2935         track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
2936       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
2937         delete array->RemoveAt(i);
2938         continue;
2939       }
2940     }
2941     else{
2942       accepted++;
2943     }
2944   }
2945   //
2946   array->Compress();
2947   SortTrackHypothesys(esdindex,checkmax,1);
2948   array = (TObjArray*) fTrackHypothesys.At(esdindex);
2949   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
2950   besttrack = (AliITStrackMI*)array->At(0);
2951   if (!besttrack)  return 0;
2952   besttrack->SetChi2MIP(8,0);
2953   fBestTrackIndex[esdindex]=0;
2954   entries = array->GetEntriesFast();
2955   AliITStrackMI *longtrack =0;
2956   minchi2 =1000;
2957   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
2958   for (Int_t itrack=entries-1;itrack>0;itrack--){
2959     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2960     if (!track->GetConstrain()) continue;
2961     if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
2962     if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
2963     if (track->GetChi2MIP(0)>4.) continue;
2964     minn = track->GetNumberOfClusters()+track->GetNDeadZone();
2965     longtrack =track;
2966   }
2967   //if (longtrack) besttrack=longtrack;
2968
2969   Int_t list[6];
2970   AliITSRecPoint * clist[6];
2971   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
2972   if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<kMaxChi2PerCluster[0]&&besttrack->GetChi2MIP(1)<kMaxChi2PerCluster[1]
2973       &&besttrack->GetChi2MIP(2)<kMaxChi2PerCluster[2]&&besttrack->GetChi2MIP(3)<kMaxChi2PerCluster[3]){
2974     RegisterClusterTracks(besttrack,esdindex);
2975   }
2976   //
2977   //
2978   if (shared>0.0){
2979     Int_t nshared;
2980     Int_t overlist[6];
2981     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
2982     if (sharedtrack>=0){
2983       //
2984       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
2985       if (besttrack){
2986         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
2987       }
2988       else return 0;
2989     }
2990   }
2991
2992   if (shared>2.5) return 0;
2993   if (shared>1.0) return besttrack;
2994   //
2995   // Don't sign clusters if not gold track
2996   //
2997   if (!besttrack->IsGoldPrimary()) return besttrack;
2998   if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
2999   //
3000   if (fConstraint[fPass]){
3001     //
3002     // sign clusters
3003     //
3004     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3005     for (Int_t i=0;i<6;i++){
3006       Int_t index = besttrack->GetClIndex(i);
3007       if (index<=0) continue;
3008       Int_t ilayer =  (index & 0xf0000000) >> 28;
3009       if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3010       AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3011       if (!c) continue;
3012       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3013       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3014       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
3015       if ( ilayer>2&& besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3016       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
3017
3018       Bool_t cansign = kTRUE;
3019       for (Int_t itrack=0;itrack<entries; itrack++){
3020         AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3021         if (!track) continue;
3022         if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3023         if ( (track->GetClIndex(ilayer)>0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3024           cansign = kFALSE;
3025           break;
3026         }
3027       }
3028       if (cansign){
3029         if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3030         if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3031         if (!c->IsUsed()) c->Use();
3032       }
3033     }
3034   }
3035   return besttrack;
3036 }
3037
3038
3039
3040 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3041 {
3042   //
3043   // get "best" hypothesys
3044   //
3045
3046   Int_t nentries = itsTracks.GetEntriesFast();
3047   for (Int_t i=0;i<nentries;i++){
3048     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3049     if (!track) continue;
3050     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3051     if (!array) continue;
3052     if (array->GetEntriesFast()<=0) continue;
3053     //
3054     AliITStrackMI* longtrack=0;
3055     Float_t minn=0;
3056     Float_t maxchi2=1000;
3057     for (Int_t j=0;j<array->GetEntriesFast();j++){
3058       AliITStrackMI* track = (AliITStrackMI*)array->At(j);
3059       if (!track) continue;
3060       if (track->GetGoldV0()) {
3061         longtrack = track;   //gold V0 track taken
3062         break;
3063       }
3064       if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3065       Float_t chi2 = track->GetChi2MIP(0);
3066       if (fAfterV0){
3067         if (!track->GetGoldV0()&&track->GetConstrain()==kFALSE) chi2+=5;
3068       }
3069       if (track->GetNumberOfClusters()+track->GetNDeadZone()>minn) maxchi2 = track->GetChi2MIP(0);
3070       //
3071       if (chi2 > maxchi2) continue;
3072       minn= track->GetNumberOfClusters()+track->GetNDeadZone();
3073       maxchi2 = chi2;
3074       longtrack=track;
3075     }
3076     //
3077     //
3078     //
3079     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3080     if (!longtrack) {longtrack = besttrack;}
3081     else besttrack= longtrack;
3082     //
3083     if (besttrack){
3084       Int_t list[6];
3085       AliITSRecPoint * clist[6];
3086       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3087       //
3088       track->SetNUsed(shared);
3089       track->SetNSkipped(besttrack->GetNSkipped());
3090       track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3091       if (shared>0){
3092         Int_t nshared;
3093         Int_t overlist[6];
3094         //
3095         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3096         //if (sharedtrack==-1) sharedtrack=0;
3097         if (sharedtrack>=0){
3098           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3099         }
3100       }
3101       if (besttrack&&fAfterV0){
3102         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3103       }
3104       if (besttrack&&fConstraint[fPass])
3105         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3106       //if (besttrack&&besttrack->fConstrain)
3107       //        UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3108       if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5){
3109         if ( (TMath::Abs(besttrack->GetD(0))>0.1) && fConstraint[fPass]) {
3110           track->SetReconstructed(kFALSE);
3111         }
3112         if ( (TMath::Abs(besttrack->GetD(1))>0.1) && fConstraint[fPass]){
3113           track->SetReconstructed(kFALSE);
3114         }
3115       }
3116
3117     }
3118   }
3119 }
3120
3121
3122 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3123   //--------------------------------------------------------------------
3124   //This function "cooks" a track label. If label<0, this track is fake.
3125   //--------------------------------------------------------------------
3126   Int_t tpcLabel=-1;
3127
3128   if ( track->GetESDtrack())   tpcLabel =  TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3129
3130    track->SetChi2MIP(9,0);
3131    Int_t nwrong=0;
3132    for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3133      Int_t cindex = track->GetClusterIndex(i);
3134      Int_t l=(cindex & 0xf0000000) >> 28;
3135      AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3136      Int_t isWrong=1;
3137      for (Int_t ind=0;ind<3;ind++){
3138        if (tpcLabel>0)
3139          if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
3140      }
3141      track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3142      nwrong+=isWrong;
3143    }
3144    track->SetFakeRatio(double(nwrong)/double(track->GetNumberOfClusters()));
3145    if (tpcLabel>0){
3146      if (track->GetFakeRatio()>wrong) track->SetLabel(-tpcLabel);
3147      else
3148        track->SetLabel(tpcLabel);
3149    }
3150
3151 }
3152
3153
3154
3155 void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
3156 {
3157   //
3158   //
3159   //  Int_t list[6];
3160   //AliITSRecPoint * clist[6];
3161   //  Int_t shared = GetNumberOfSharedClusters(track,index,list,clist);
3162   Float_t dedx[4];
3163   Int_t accepted=0;
3164   track->SetChi2MIP(9,0);
3165   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3166     Int_t cindex = track->GetClusterIndex(i);
3167     Int_t l=(cindex & 0xf0000000) >> 28;
3168     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3169     Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3170     Int_t isWrong=1;
3171     for (Int_t ind=0;ind<3;ind++){
3172       if (cl->GetLabel(ind)==lab) isWrong=0;
3173     }
3174     track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3175     if (l<2) continue;
3176     //if (l>3 && (cl->GetNy()>4) || (cl->GetNz()>4)) continue;  //shared track
3177     //if (l>3&& !(cl->GetType()==1||cl->GetType()==10)) continue;
3178     //if (l<4&& !(cl->GetType()==1)) continue;
3179     dedx[accepted]= track->GetSampledEdx(i);
3180     //dedx[accepted]= track->fNormQ[l];
3181     accepted++;
3182   }
3183   if (accepted<1) {
3184     track->SetdEdx(0);
3185     return;
3186   }
3187   Int_t indexes[4];
3188   TMath::Sort(accepted,dedx,indexes,kFALSE);
3189   Double_t low=0.;
3190   Double_t up=0.51;
3191   Double_t nl=low*accepted, nu =up*accepted;
3192   Float_t sumamp = 0;
3193   Float_t sumweight =0;
3194   for (Int_t i=0; i<accepted; i++) {
3195     Float_t weight =1;
3196     if (i<nl+0.1)     weight = TMath::Max(1.-(nl-i),0.);
3197     if (i>nu-1)     weight = TMath::Max(nu-i,0.);
3198     sumamp+= dedx[indexes[i]]*weight;
3199     sumweight+=weight;
3200   }
3201   track->SetdEdx(sumamp/sumweight);
3202 }
3203
3204
3205 void  AliITStrackerMI::MakeCoeficients(Int_t ntracks){
3206   //
3207   //
3208   if (fCoeficients) delete []fCoeficients;
3209   fCoeficients = new Float_t[ntracks*48];
3210   for (Int_t i=0;i<ntracks*48;i++) fCoeficients[i]=-1.;
3211 }
3212
3213
3214 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3215 {
3216   //
3217   //
3218   //
3219   Float_t erry,errz;
3220   Float_t theta = track->GetTgl();
3221   Float_t phi   = track->GetSnp();
3222   phi = TMath::Sqrt(phi*phi/(1.-phi*phi));
3223   GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz);
3224   Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz);
3225   Float_t ny,nz;
3226   GetNTeor(layer,cluster, theta,phi,ny,nz);
3227   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3228   if (delta>1){
3229     chi2+=0.5*TMath::Min(delta/2,2.);
3230     chi2+=2.*cluster->GetDeltaProbability();
3231   }
3232   //
3233   track->SetNy(layer,ny);
3234   track->SetNz(layer,nz);
3235   track->SetSigmaY(layer,erry);
3236   track->SetSigmaZ(layer, errz);
3237   //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3238   track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/(1.- track->GetSnp()*track->GetSnp())));
3239   return chi2;
3240
3241 }
3242
3243 Int_t    AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3244 {
3245   //
3246   //
3247   //
3248   Int_t layer = (index & 0xf0000000) >> 28;
3249   track->SetClIndex(layer, index);
3250   if ( (layer>1) &&track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3251     chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3252     track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3253   }
3254
3255   if (cl->GetQ()<=0) return 0;  // ingore the "virtual" clusters
3256
3257   // Take into account the mis-alignment
3258   Double_t x=track->GetX()+cl->GetX();
3259   if (!track->PropagateTo(x,0.,0.)) return 0;
3260
3261   return track->UpdateMI(cl->GetY(),cl->GetZ(),track->GetSigmaY(layer),track->GetSigmaZ(layer),chi2,index);
3262 }
3263
3264 void AliITStrackerMI::GetNTeor(Int_t layer, const AliITSRecPoint* /*cl*/, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz)
3265 {
3266   //
3267   //get "mean shape"
3268   //
3269   if (layer==0){
3270     ny = 1.+TMath::Abs(phi)*3.2;
3271     nz = 1.+TMath::Abs(theta)*0.34;
3272     return;
3273   }
3274   if (layer==1){
3275     ny = 1.+TMath::Abs(phi)*3.2;
3276     nz = 1.+TMath::Abs(theta)*0.28;
3277     return;
3278   }
3279
3280   if (layer>3){
3281     ny = 2.02+TMath::Abs(phi)*1.95;
3282     nz = 2.02+TMath::Abs(phi)*2.35;
3283     return;
3284   }
3285   ny  = 6.6-2.7*TMath::Abs(phi);
3286   nz  = 2.8-3.11*TMath::Abs(phi)+0.45*TMath::Abs(theta);
3287 }
3288
3289
3290
3291 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)
3292 {
3293   //calculate cluster position error
3294   //
3295   Float_t nz,ny;
3296   GetNTeor(layer, cl,theta,phi,ny,nz);
3297   erry   = TMath::Sqrt(cl->GetSigmaY2());
3298   errz   = TMath::Sqrt(cl->GetSigmaZ2());
3299   //
3300   // PIXELS
3301   if (layer<2){
3302
3303     if (TMath::Abs(ny-cl->GetNy())>0.6)  {
3304       if (ny<cl->GetNy()){
3305         erry*=0.4+TMath::Abs(ny-cl->GetNy());
3306         errz*=0.4+TMath::Abs(ny-cl->GetNy());
3307       }else{
3308         erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
3309         errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
3310       }
3311     }
3312     if (TMath::Abs(nz-cl->GetNz())>1.)  {
3313       erry*=TMath::Abs(nz-cl->GetNz());
3314       errz*=TMath::Abs(nz-cl->GetNz());
3315     }
3316     erry*=0.85;
3317     errz*=0.85;
3318     erry= TMath::Min(erry,float(0.005));
3319     errz= TMath::Min(errz,float(0.03));
3320     return 10;
3321   }
3322
3323 //STRIPS
3324   if (layer>3){
3325     //factor 1.8 appears in new simulation
3326     //
3327     Float_t scale=1.8;
3328     if (cl->GetNy()==100||cl->GetNz()==100){
3329       erry = 0.004*scale;
3330       errz = 0.2*scale;
3331       return 100;
3332     }
3333     if (cl->GetNy()+cl->GetNz()>12){
3334       erry = 0.06*scale;
3335       errz = 0.57*scale;
3336       return 100;
3337     }
3338     Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
3339     Float_t chargematch = TMath::Max(double(normq/expQ),2.);
3340     //
3341     if (cl->GetType()==1 || cl->GetType()==10 ){
3342       if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
3343         errz = 0.043*scale;
3344         erry = 0.00094*scale;
3345         return 101;
3346       }
3347       if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
3348         errz = 0.06*scale;
3349         erry =0.0013*scale;
3350         return 102;
3351       }
3352       erry = 0.0027*scale;
3353       errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15)*scale;
3354       return 103;
3355     }
3356     if (cl->GetType()==2 || cl->GetType()==11 ){
3357       erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05)*scale;
3358       errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5)*scale;
3359       return 104;
3360     }
3361
3362     if (cl->GetType()>100 ){
3363       if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
3364         errz = 0.05*scale;
3365         erry = 0.00096*scale;
3366         return 105;
3367       }
3368       if (cl->GetNy()+cl->GetNz()-nz-ny<1){
3369         errz = 0.10*scale;
3370         erry = 0.0025*scale;
3371         return 106;
3372       }
3373
3374       errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4)*scale;
3375       erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05)*scale;
3376       return 107;
3377     }
3378     Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
3379     if (diff<1) diff=1;
3380     if (diff>4) diff=4;
3381
3382     if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){
3383       errz = 0.14*diff;
3384       erry = 0.003*diff;
3385       return 108;
3386     }
3387     erry = 0.04*diff;
3388     errz = 0.06*diff;
3389     return 109;
3390   }
3391   //DRIFTS
3392   Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
3393   Float_t chargematch = normq/expQ;
3394   Float_t factorz=1;
3395   Int_t   cnz = cl->GetNz()%10;
3396   //charge match
3397   if (cl->GetType()==1){
3398     if (chargematch<1.25){
3399       erry =  0.0028*(1.+6./cl->GetQ());  // gold clusters
3400     }
3401     else{
3402       erry = 0.003*chargematch;
3403       if (cl->GetNz()==3) erry*=1.5;
3404     }
3405     if (chargematch<1.0){
3406       errz =  0.0011*(1.+6./cl->GetQ());
3407     }
3408     else{
3409       errz = 0.002*(1+2*(chargematch-1.));
3410     }
3411     if (cnz>nz+0.6) {
3412       erry*=(cnz-nz+0.5);
3413       errz*=1.4*(cnz-nz+0.5);
3414     }
3415   }
3416   if (cl->GetType()>1){
3417     if (chargematch<1){
3418       erry =  0.00385*(1.+6./cl->GetQ());  // gold clusters
3419       errz =  0.0016*(1.+6./cl->GetQ());
3420     }
3421     else{
3422       errz = 0.0014*(1+3*(chargematch-1.));
3423       erry = 0.003*(1+3*(chargematch-1.));
3424     }
3425     if (cnz>nz+0.6) {
3426       erry*=(cnz-nz+0.5);
3427       errz*=1.4*(cnz-nz+0.5);
3428     }
3429   }
3430
3431   if (TMath::Abs(cl->GetY())>2.5){
3432     factorz*=1+2*(TMath::Abs(cl->GetY())-2.5);
3433   }
3434   if (TMath::Abs(cl->GetY())<1){
3435     factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.);
3436   }
3437   factorz= TMath::Min(factorz,float(4.));
3438   errz*=factorz;
3439
3440   erry= TMath::Min(erry,float(0.05));
3441   errz= TMath::Min(errz,float(0.05));
3442   return 200;
3443 }
3444
3445
3446
3447 void   AliITStrackerMI::GetDCASigma(AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3448 {
3449   //
3450   //DCA sigmas parameterization
3451   //to be paramterized using external parameters in future
3452   //
3453   //
3454   sigmarfi = 0.004+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3455   sigmaz   = 0.011+4.37*TMath::Abs(track->GetC());
3456 }
3457
3458
3459 void AliITStrackerMI::SignDeltas( TObjArray *ClusterArray, Float_t vz)
3460 {
3461   //
3462   //
3463   Int_t entries = ClusterArray->GetEntriesFast();
3464   if (entries<4) return;
3465   AliITSRecPoint* cluster = (AliITSRecPoint*)ClusterArray->At(0);
3466   Int_t layer = cluster->GetLayer();
3467   if (layer>1) return;
3468   Int_t index[10000];
3469   Int_t ncandidates=0;
3470   Float_t r = (layer>0)? 7:4;
3471   //
3472   for (Int_t i=0;i<entries;i++){
3473     AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(i);
3474     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3475     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3476     index[ncandidates] = i;  //candidate to belong to delta electron track
3477     ncandidates++;
3478     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3479       cl0->SetDeltaProbability(1);
3480     }
3481   }
3482   //
3483   //
3484   //
3485   for (Int_t i=0;i<ncandidates;i++){
3486     AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(index[i]);
3487     if (cl0->GetDeltaProbability()>0.8) continue;
3488     //
3489     Int_t ncl = 0;
3490     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3491     sumy=sumz=sumy2=sumyz=sumw=0.0;
3492     for (Int_t j=0;j<ncandidates;j++){
3493       if (i==j) continue;
3494       AliITSRecPoint* cl1 = (AliITSRecPoint*)ClusterArray->At(index[j]);
3495       //
3496       Float_t dz = cl0->GetZ()-cl1->GetZ();
3497       Float_t dy = cl0->GetY()-cl1->GetY();
3498       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3499         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3500         y[ncl] = cl1->GetY();
3501         z[ncl] = cl1->GetZ();
3502         sumy+= y[ncl]*weight;
3503         sumz+= z[ncl]*weight;
3504         sumy2+=y[ncl]*y[ncl]*weight;
3505         sumyz+=y[ncl]*z[ncl]*weight;
3506         sumw+=weight;
3507         ncl++;
3508       }
3509     }
3510     if (ncl<4) continue;
3511     Float_t det = sumw*sumy2  - sumy*sumy;
3512     Float_t delta=1000;
3513     if (TMath::Abs(det)>0.01){
3514       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
3515       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
3516       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3517     }
3518     else{
3519       Float_t z0  = sumyz/sumy;
3520       delta = TMath::Abs(cl0->GetZ()-z0);
3521     }
3522     if ( delta<0.05) {
3523       cl0->SetDeltaProbability(1-20.*delta);
3524     }
3525   }
3526 }
3527
3528
3529 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3530 {
3531   //
3532   //
3533   track->UpdateESDtrack(flags);
3534   AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3535   if (oldtrack) delete oldtrack;
3536   track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3537   if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3538     printf("Problem\n");
3539   }
3540 }
3541
3542
3543
3544 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3545   //
3546   //Get nearest upper layer close to the point xr.
3547   // rough approximation
3548   //
3549   const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3550   Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3551   Int_t res =6;
3552   for (Int_t i=0;i<6;i++){
3553     if (radius<kRadiuses[i]){
3554       res =i;
3555       break;
3556     }
3557   }
3558   return res;
3559 }
3560
3561
3562 void AliITStrackerMI::UpdateTPCV0(AliESD *event){
3563   //
3564   //try to update, or reject TPC  V0s
3565   //
3566   Int_t nv0s = event->GetNumberOfV0s();
3567   Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3568
3569   for (Int_t i=0;i<nv0s;i++){
3570     AliESDv0 * vertex = event->GetV0(i);
3571     Int_t ip = vertex->GetIndex(0);
3572     Int_t im = vertex->GetIndex(1);
3573     //
3574     TObjArray * arrayp = (ip<nitstracks) ? (TObjArray*)fTrackHypothesys.At(ip):0;
3575     TObjArray * arraym = (im<nitstracks) ? (TObjArray*)fTrackHypothesys.At(im):0;
3576     AliITStrackMI * trackp = (arrayp!=0) ? (AliITStrackMI*)arrayp->At(0):0;
3577     AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0;
3578     //
3579     //
3580     if (trackp){
3581       if (trackp->GetNumberOfClusters()+trackp->GetNDeadZone()>5.5){
3582         if (trackp->GetConstrain()&&trackp->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3583         if (!trackp->GetConstrain()&&trackp->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3584       }
3585     }
3586
3587     if (trackm){
3588       if (trackm->GetNumberOfClusters()+trackm->GetNDeadZone()>5.5){
3589         if (trackm->GetConstrain()&&trackm->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3590         if (!trackm->GetConstrain()&&trackm->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3591       }
3592     }
3593     if (vertex->GetStatus()==-100) continue;
3594     //
3595     Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
3596     Int_t clayer = GetNearestLayer(xrp);                    //I.B.
3597     vertex->SetNBefore(clayer);        //
3598     vertex->SetChi2Before(9*clayer);   //
3599     vertex->SetNAfter(6-clayer);       //
3600     vertex->SetChi2After(0);           //
3601     //
3602     if (clayer >1 ){ // calculate chi2 before vertex
3603       Float_t chi2p = 0, chi2m=0;
3604       //
3605       if (trackp){
3606         for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3607           if (trackp->GetClIndex(ilayer)>0){
3608             chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3609               trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3610           }
3611           else{
3612             chi2p+=9;
3613           }
3614         }
3615       }else{
3616         chi2p = 9*clayer;
3617       }
3618       //
3619       if (trackm){
3620         for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3621           if (trackm->GetClIndex(ilayer)>0){
3622             chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3623               trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3624           }
3625           else{
3626             chi2m+=9;
3627           }
3628         }
3629       }else{
3630         chi2m = 9*clayer;
3631       }
3632       vertex->SetChi2Before(TMath::Min(chi2p,chi2m));
3633       if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10);  // track exist before vertex
3634     }
3635
3636     if (clayer < 5 ){ // calculate chi2 after vertex
3637       Float_t chi2p = 0, chi2m=0;
3638       //
3639       if (trackp&&TMath::Abs(trackp->GetTgl())<1.){
3640         for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3641           if (trackp->GetClIndex(ilayer)>0){
3642             chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3643               trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3644           }
3645           else{
3646             chi2p+=9;
3647           }
3648         }
3649       }else{
3650         chi2p = 0;
3651       }
3652       //
3653       if (trackm&&TMath::Abs(trackm->GetTgl())<1.){
3654         for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3655           if (trackm->GetClIndex(ilayer)>0){
3656             chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3657               trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3658           }
3659           else{
3660             chi2m+=9;
3661           }
3662         }
3663       }else{
3664         chi2m = 0;
3665       }
3666       vertex->SetChi2After(TMath::Max(chi2p,chi2m));
3667       if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20);  // track not found in ITS
3668     }
3669   }
3670   //
3671 }
3672
3673
3674
3675 void  AliITStrackerMI::FindV02(AliESD *event)
3676 {
3677   //
3678   // V0 finder
3679   //
3680   //  Cuts on DCA -  R dependent
3681   //          max distance DCA between 2 tracks cut
3682   //          maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist);
3683   //
3684   const Float_t kMaxDist0      = 0.1;
3685   const Float_t kMaxDist1      = 0.1;
3686   const Float_t kMaxDist       = 1;
3687   const Float_t kMinPointAngle  = 0.85;
3688   const Float_t kMinPointAngle2 = 0.99;
3689   const Float_t kMinR           = 0.5;
3690   const Float_t kMaxR           = 220;
3691   //const Float_t kCausality0Cut   = 0.19;
3692   //const Float_t kLikelihood01Cut = 0.25;
3693   //const Float_t kPointAngleCut   = 0.9996;
3694   const Float_t kCausality0Cut   = 0.19;
3695   const Float_t kLikelihood01Cut = 0.45;
3696   const Float_t kLikelihood1Cut  = 0.5;
3697   const Float_t kCombinedCut     = 0.55;
3698
3699   //
3700   //
3701   TTreeSRedirector &cstream = *fDebugStreamer;
3702   Int_t ntracks    = event->GetNumberOfTracks();
3703   Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3704   fOriginal.Expand(ntracks);
3705   fTrackHypothesys.Expand(ntracks);
3706   fBestHypothesys.Expand(ntracks);
3707   //
3708   AliHelix * helixes   = new AliHelix[ntracks+2];
3709   TObjArray trackarray(ntracks+2);     //array with tracks - with vertex constrain
3710   TObjArray trackarrayc(ntracks+2);    //array of "best    tracks" - without vertex constrain
3711   TObjArray trackarrayl(ntracks+2);    //array of "longest tracks" - without vertex constrain
3712   Bool_t * forbidden   = new Bool_t [ntracks+2];
3713   Int_t   *itsmap      = new Int_t  [ntracks+2];
3714   Float_t *dist        = new Float_t[ntracks+2];
3715   Float_t *normdist0   = new Float_t[ntracks+2];
3716   Float_t *normdist1   = new Float_t[ntracks+2];
3717   Float_t *normdist    = new Float_t[ntracks+2];
3718   Float_t *norm        = new Float_t[ntracks+2];
3719   Float_t *maxr        = new Float_t[ntracks+2];
3720   Float_t *minr        = new Float_t[ntracks+2];
3721   Float_t *minPointAngle= new Float_t[ntracks+2];
3722   //
3723   AliV0 *pvertex      = new AliV0;
3724   AliITStrackMI * dummy= new AliITStrackMI;
3725   dummy->SetLabel(0);
3726   AliITStrackMI  trackat0;    //temporary track for DCA calculation
3727   //
3728   Float_t primvertex[3]={GetX(),GetY(),GetZ()};
3729   //
3730   // make its -  esd map
3731   //
3732   for (Int_t itrack=0;itrack<ntracks+2;itrack++) {
3733     itsmap[itrack]        = -1;
3734     forbidden[itrack]     = kFALSE;
3735     maxr[itrack]          = kMaxR;
3736     minr[itrack]          = kMinR;
3737     minPointAngle[itrack] = kMinPointAngle;
3738   }
3739   for (Int_t itrack=0;itrack<nitstracks;itrack++){
3740     AliITStrackMI * original =   (AliITStrackMI*)(fOriginal.At(itrack));
3741     Int_t           esdindex =   original->GetESDtrack()->GetID();
3742     itsmap[esdindex]         =   itrack;
3743   }
3744   //
3745   // create its tracks from esd tracks if not done before
3746   //
3747   for (Int_t itrack=0;itrack<ntracks;itrack++){
3748     if (itsmap[itrack]>=0) continue;
3749     AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack)));
3750     //tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY());
3751     //tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ();
3752     tpctrack->GetDZ(GetX(),GetY(),GetZ(),tpctrack->GetDP());   //I.B.
3753     if (tpctrack->GetD(0)<20 && tpctrack->GetD(1)<20){
3754       // tracks which can reach inner part of ITS
3755       // propagate track to outer its volume - with correction for material
3756       CorrectForDeadZoneMaterial(tpctrack);
3757     }
3758     itsmap[itrack] = nitstracks;
3759     fOriginal.AddAt(tpctrack,nitstracks);
3760     nitstracks++;
3761   }
3762   //
3763   // fill temporary arrays
3764   //
3765   for (Int_t itrack=0;itrack<ntracks;itrack++){
3766     AliESDtrack *  esdtrack = event->GetTrack(itrack);
3767     Int_t          itsindex = itsmap[itrack];
3768     AliITStrackMI *original = (AliITStrackMI*)fOriginal.At(itsmap[itrack]);
3769     if (!original) continue;
3770     AliITStrackMI *bestConst  = 0;
3771     AliITStrackMI *bestLong   = 0;
3772     AliITStrackMI *best       = 0;
3773     //
3774     //
3775     TObjArray * array    = (TObjArray*)  fTrackHypothesys.At(itsindex);
3776     Int_t       hentries = (array==0) ?  0 : array->GetEntriesFast();
3777     // Get best track with vertex constrain
3778     for (Int_t ih=0;ih<hentries;ih++){
3779       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
3780       if (!trackh->GetConstrain()) continue;
3781       if (!bestConst) bestConst = trackh;
3782       if (trackh->GetNumberOfClusters()>5.0){
3783         bestConst  = trackh;                         // full track -  with minimal chi2
3784         break;
3785       }
3786       if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone())  continue;
3787       bestConst = trackh;
3788       break;
3789     }
3790     // Get best long track without vertex constrain and best track without vertex constrain
3791     for (Int_t ih=0;ih<hentries;ih++){
3792       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
3793       if (trackh->GetConstrain()) continue;
3794       if (!best)     best     = trackh;
3795       if (!bestLong) bestLong = trackh;
3796       if (trackh->GetNumberOfClusters()>5.0){
3797         bestLong  = trackh;                         // full track -  with minimal chi2
3798         break;
3799       }
3800       if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone())  continue;
3801       bestLong = trackh;
3802     }
3803     if (!best) {
3804       best     = original;
3805       bestLong = original;
3806     }
3807     //I.B. trackat0 = *bestLong;
3808     new (&trackat0) AliITStrackMI(*bestLong);
3809     Double_t xx,yy,zz,alpha;
3810     bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz);
3811     alpha = TMath::ATan2(yy,xx);
3812     trackat0.Propagate(alpha,0);
3813     // calculate normalized distances to the vertex
3814     //
3815     Float_t ptfac  = (1.+100.*TMath::Abs(trackat0.GetC()));
3816     if ( bestLong->GetNumberOfClusters()>3 ){
3817       dist[itsindex]      = trackat0.GetY();
3818       norm[itsindex]      = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
3819       normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
3820       normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
3821       normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
3822       if (!bestConst){
3823         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<6) normdist[itsindex]*=2.;
3824         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<5) normdist[itsindex]*=2.;
3825         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<4) normdist[itsindex]*=2.;
3826       }else{
3827         if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<6) normdist[itsindex]*=1.5;
3828         if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<5) normdist[itsindex]*=1.5;
3829       }
3830     }
3831     else{
3832       if (bestConst&&bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5){
3833         dist[itsindex] = bestConst->GetD(0);
3834         norm[itsindex] = bestConst->GetDnorm(0);
3835         normdist0[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
3836         normdist1[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
3837         normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
3838       }else{
3839         dist[itsindex]      = trackat0.GetY();
3840         norm[itsindex]      = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
3841         normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
3842         normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
3843         normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
3844         if (TMath::Abs(trackat0.GetTgl())>1.05){
3845           if (normdist[itsindex]<3) forbidden[itsindex]=kTRUE;
3846           if (normdist[itsindex]>3) {
3847             minr[itsindex] = TMath::Max(Float_t(40.),minr[itsindex]);
3848           }
3849         }
3850       }
3851     }
3852     //
3853     //-----------------------------------------------------------
3854     //Forbid primary track candidates -
3855     //
3856     //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5");
3857     //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5");
3858     //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0");
3859     //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0");
3860     //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2");
3861     //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1");
3862     //-----------------------------------------------------------
3863     if (bestConst){
3864       if (bestLong->GetNumberOfClusters()<4       && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5)               forbidden[itsindex]=kTRUE;
3865       if (normdist[itsindex]<3 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>5.5)               forbidden[itsindex]=kTRUE;
3866       if (normdist[itsindex]<2 && bestConst->GetClIndex(0)>0 && bestConst->GetClIndex(1)>0 ) forbidden[itsindex]=kTRUE;
3867       if (normdist[itsindex]<1 && bestConst->GetClIndex(0)>0)                              forbidden[itsindex]=kTRUE;
3868       if (normdist[itsindex]<4 && bestConst->GetNormChi2(0)<2)                             forbidden[itsindex]=kTRUE;
3869       if (normdist[itsindex]<5 && bestConst->GetNormChi2(0)<1)                             forbidden[itsindex]=kTRUE;
3870       if (bestConst->GetNormChi2(0)<2.5) {
3871         minPointAngle[itsindex]= 0.9999;
3872         maxr[itsindex]         = 10;
3873       }
3874     }
3875     //
3876     //forbid daughter kink candidates
3877     //
3878     if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
3879     Bool_t isElectron = kTRUE;
3880     Bool_t isProton   = kTRUE;
3881     Double_t pid[5];
3882     esdtrack->GetESDpid(pid);
3883     for (Int_t i=1;i<5;i++){
3884       if (pid[0]<pid[i]) isElectron= kFALSE;
3885       if (pid[4]<pid[i]) isProton= kFALSE;
3886     }
3887     if (isElectron){
3888       forbidden[itsindex]=kFALSE;
3889       normdist[itsindex]*=-1;
3890     }
3891     if (isProton){
3892       if (normdist[itsindex]>2) forbidden[itsindex]=kFALSE;
3893       normdist[itsindex]*=-1;
3894     }
3895
3896     //
3897     // Causality cuts in TPC volume
3898     //
3899     if (esdtrack->GetTPCdensity(0,10) >0.6)  maxr[itsindex] = TMath::Min(Float_t(110),maxr[itsindex]);
3900     if (esdtrack->GetTPCdensity(10,30)>0.6)  maxr[itsindex] = TMath::Min(Float_t(120),maxr[itsindex]);
3901     if (esdtrack->GetTPCdensity(20,40)>0.6)  maxr[itsindex] = TMath::Min(Float_t(130),maxr[itsindex]);
3902     if (esdtrack->GetTPCdensity(30,50)>0.6)  maxr[itsindex] = TMath::Min(Float_t(140),maxr[itsindex]);
3903     //
3904     if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->GetNumberOfClusters()<3) minr[itsindex]=100;
3905     //
3906     //
3907     if (kFALSE){
3908       cstream<<"Track"<<
3909         "Tr0.="<<best<<
3910         "Tr1.="<<((bestConst)? bestConst:dummy)<<
3911         "Tr2.="<<bestLong<<
3912         "Tr3.="<<&trackat0<<
3913         "Esd.="<<esdtrack<<
3914         "Dist="<<dist[itsindex]<<
3915         "ND0="<<normdist0[itsindex]<<
3916         "ND1="<<normdist1[itsindex]<<
3917         "ND="<<normdist[itsindex]<<
3918         "Pz="<<primvertex[2]<<
3919         "Forbid="<<forbidden[itsindex]<<
3920         "\n";
3921       //
3922     }
3923     trackarray.AddAt(best,itsindex);
3924     trackarrayc.AddAt(bestConst,itsindex);
3925     trackarrayl.AddAt(bestLong,itsindex);
3926     new (&helixes[itsindex]) AliHelix(*best);
3927   }
3928   //
3929   //
3930   //
3931   // first iterration of V0 finder
3932   //
3933   for (Int_t iesd0=0;iesd0<ntracks;iesd0++){
3934     Int_t itrack0 = itsmap[iesd0];
3935     if (forbidden[itrack0]) continue;
3936     AliITStrackMI * btrack0 = (AliITStrackMI*)trackarray.At(itrack0);
3937     if (!btrack0) continue;
3938     if (btrack0->GetSign()>0) continue;
3939     AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
3940     //
3941     for (Int_t iesd1=0;iesd1<ntracks;iesd1++){
3942       Int_t itrack1 = itsmap[iesd1];
3943       if (forbidden[itrack1]) continue;
3944
3945       AliITStrackMI * btrack1 = (AliITStrackMI*)trackarray.At(itrack1);
3946       if (!btrack1) continue;
3947       if (btrack1->GetSign()<0) continue;
3948       Bool_t isGold = kFALSE;
3949       if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){
3950         isGold = kTRUE;
3951       }
3952       AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
3953       AliHelix &h1 = helixes[itrack0];
3954       AliHelix &h2 = helixes[itrack1];
3955       //
3956       // find linear distance
3957       Double_t rmin =0;
3958       //
3959       //
3960       //
3961       Double_t phase[2][2],radius[2];
3962       Int_t  points = h1.GetRPHIintersections(h2, phase, radius);
3963       if    (points==0)  continue;
3964       Double_t delta[2]={1000000,1000000};
3965       rmin = radius[0];
3966       h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]);
3967       if (points==2){
3968         if (radius[1]<rmin) rmin = radius[1];
3969         h1.ParabolicDCA(h2,phase[1][0],phase[1][1],radius[1],delta[1]);
3970       }
3971       rmin = TMath::Sqrt(rmin);
3972       Double_t distance = 0;
3973       Double_t radiusC  = 0;
3974       Int_t    iphase   = 0;
3975       if (points==1 || delta[0]<delta[1]){
3976         distance = TMath::Sqrt(delta[0]);
3977         radiusC  = TMath::Sqrt(radius[0]);
3978       }else{
3979         distance = TMath::Sqrt(delta[1]);
3980         radiusC  = TMath::Sqrt(radius[1]);
3981         iphase=1;
3982       }
3983       if (radiusC<TMath::Max(minr[itrack0],minr[itrack1]))    continue;
3984       if (radiusC>TMath::Min(maxr[itrack0],maxr[itrack1]))     continue;
3985       Float_t maxDist  = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1));
3986       if (distance>maxDist) continue;
3987       Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex);
3988       if (pointAngle<TMath::Max(minPointAngle[itrack0],minPointAngle[itrack1])) continue;
3989       //
3990       //
3991       //      Double_t distance = TestV0(h1,h2,pvertex,rmin);
3992       //
3993       //       if (distance>maxDist)           continue;
3994       //       if (pvertex->GetRr()<kMinR)     continue;
3995       //       if (pvertex->GetRr()>kMaxR)     continue;
3996       AliITStrackMI * track0=btrack0;
3997       AliITStrackMI * track1=btrack1;
3998       //      if (pvertex->GetRr()<3.5){
3999       if (radiusC<3.5){
4000         //use longest tracks inside the pipe
4001         track0 = (AliITStrackMI*)trackarrayl.At(itrack0);
4002         track1 = (AliITStrackMI*)trackarrayl.At(itrack1);
4003       }
4004       //
4005       //
4006       pvertex->SetParamN(*track0);
4007       pvertex->SetParamP(*track1);
4008       pvertex->Update(primvertex);
4009       pvertex->SetClusters(track0->ClIndex(),track1->ClIndex());  // register clusters
4010
4011       if (pvertex->GetRr()<kMinR) continue;
4012       if (pvertex->GetRr()>kMaxR) continue;
4013       if (pvertex->GetV0CosineOfPointingAngle()<kMinPointAngle) continue;
4014 //Bo:      if (pvertex->GetDist2()>maxDist) continue;
4015       if (pvertex->GetDcaV0Daughters()>maxDist) continue;
4016 //Bo:        pvertex->SetLab(0,track0->GetLabel());
4017 //Bo:        pvertex->SetLab(1,track1->GetLabel());
4018       pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4019       pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4020       //
4021       AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
4022       AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
4023
4024       //
4025       //
4026       TObjArray * array0b     = (TObjArray*)fBestHypothesys.At(itrack0);
4027       if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->GetTgl())<1.1)
4028         FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack0),itrack0, kFALSE);
4029       TObjArray * array1b    = (TObjArray*)fBestHypothesys.At(itrack1);
4030       if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->GetTgl())<1.1)
4031         FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack1),itrack1, kFALSE);
4032       //
4033       AliITStrackMI * track0b = (AliITStrackMI*)fOriginal.At(itrack0);
4034       AliITStrackMI * track1b = (AliITStrackMI*)fOriginal.At(itrack1);
4035       AliITStrackMI * track0l = (AliITStrackMI*)fOriginal.At(itrack0);
4036       AliITStrackMI * track1l = (AliITStrackMI*)fOriginal.At(itrack1);
4037
4038       Float_t minchi2before0=16;
4039       Float_t minchi2before1=16;
4040       Float_t minchi2after0 =16;
4041       Float_t minchi2after1 =16;
4042       Double_t xrp[3]; pvertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
4043       Int_t maxLayer = GetNearestLayer(xrp);                   //I.B.
4044
4045       if (array0b) for (Int_t i=0;i<5;i++){
4046         // best track after vertex
4047         AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i);
4048         if (!btrack) continue;
4049         if (btrack->GetNumberOfClusters()>track0l->GetNumberOfClusters()) track0l = btrack;
4050         //      if (btrack->fX<pvertex->GetRr()-2.-0.5/(0.1+pvertex->GetAnglep()[2])) {
4051         if (btrack->GetX()<pvertex->GetRr()-2.) {
4052           if ( (maxLayer>i+2|| (i==0)) && btrack->GetNumberOfClusters()==(6-i)&&i<3){
4053             Float_t sumchi2= 0;
4054             Float_t sumn   = 0;
4055             if (maxLayer<3){   // take prim vertex as additional measurement
4056               if (normdist[itrack0]>0 && htrackc0){
4057                 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack0]*normdist[itrack0],16.);
4058               }else{
4059                 sumchi2 +=  TMath::Min((3.-maxLayer)*(3*normdist[itrack0]*normdist[itrack0]+3.),16.);
4060               }
4061               sumn    +=  3-maxLayer;
4062             }
4063             for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4064               sumn+=1.;
4065               if (!btrack->GetClIndex(ilayer)){
4066                 sumchi2+=25;
4067                 continue;
4068               }else{
4069                 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4070                 for (Int_t itrack=0;itrack<4;itrack++){
4071                   if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack0){
4072                     sumchi2+=18.;  //shared cluster
4073                     break;
4074                   }
4075                 }
4076                 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4077                 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4078               }
4079             }
4080             sumchi2/=sumn;
4081             if (sumchi2<minchi2before0) minchi2before0=sumchi2;
4082           }
4083           continue;   //safety space - Geo manager will give exact layer
4084         }
4085         track0b       = btrack;
4086         minchi2after0 = btrack->GetNormChi2(i);
4087         break;
4088       }
4089       if (array1b) for (Int_t i=0;i<5;i++){
4090         // best track after vertex
4091         AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i);
4092         if (!btrack) continue;
4093         if (btrack->GetNumberOfClusters()>track1l->GetNumberOfClusters()) track1l = btrack;
4094         //      if (btrack->fX<pvertex->GetRr()-2-0.5/(0.1+pvertex->GetAnglep()[2])){
4095         if (btrack->GetX()<pvertex->GetRr()-2){
4096           if ((maxLayer>i+2 || (i==0))&&btrack->GetNumberOfClusters()==(6-i)&&(i<3)){
4097             Float_t sumchi2= 0;
4098             Float_t sumn   = 0;
4099             if (maxLayer<3){   // take prim vertex as additional measurement
4100               if (normdist[itrack1]>0 && htrackc1){
4101                 sumchi2 +=  TMath::Min((3.-maxLayer)*normdist[itrack1]*normdist[itrack1],16.);
4102               }else{
4103                 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack1]*normdist[itrack1]+3.),16.);
4104               }
4105               sumn    +=  3-maxLayer;
4106             }
4107             for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4108               sumn+=1.;
4109               if (!btrack->GetClIndex(ilayer)){
4110                 sumchi2+=30;
4111                 continue;
4112               }else{
4113                 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4114                 for (Int_t itrack=0;itrack<4;itrack++){
4115                   if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack1){
4116                     sumchi2+=18.;  //shared cluster
4117                     break;
4118                   }
4119                 }
4120                 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4121                 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4122               }
4123             }
4124             sumchi2/=sumn;
4125             if (sumchi2<minchi2before1) minchi2before1=sumchi2;
4126           }
4127           continue;   //safety space - Geo manager will give exact layer
4128         }
4129         track1b = btrack;
4130         minchi2after1 = btrack->GetNormChi2(i);
4131         break;
4132       }
4133       //
4134       // position resolution - used for DCA cut
4135       Float_t sigmad = track0b->GetSigmaY2()+track0b->GetSigmaZ2()+track1b->GetSigmaY2()+track1b->GetSigmaZ2()+
4136         (track0b->GetX()-pvertex->GetRr())*(track0b->GetX()-pvertex->GetRr())*(track0b->GetSigmaSnp2()+track0b->GetSigmaTgl2())+
4137         (track1b->GetX()-pvertex->GetRr())*(track1b->GetX()-pvertex->GetRr())*(track1b->GetSigmaSnp2()+track1b->GetSigmaTgl2());
4138       sigmad =TMath::Sqrt(sigmad)+0.04;
4139       if (pvertex->GetRr()>50){
4140         Double_t cov0[15],cov1[15];
4141         track0b->GetESDtrack()->GetInnerExternalCovariance(cov0);
4142         track1b->GetESDtrack()->GetInnerExternalCovariance(cov1);
4143         sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+
4144           (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+
4145           (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]);
4146         sigmad =TMath::Sqrt(sigmad)+0.05;
4147       }
4148       //
4149       AliV0 vertex2;
4150       vertex2.SetParamN(*track0b);
4151       vertex2.SetParamP(*track1b);
4152       vertex2.Update(primvertex);
4153       //Bo:      if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4154       if (vertex2.GetDcaV0Daughters()<=pvertex->GetDcaV0Daughters()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4155         pvertex->SetParamN(*track0b);
4156         pvertex->SetParamP(*track1b);
4157         pvertex->Update(primvertex);
4158         pvertex->SetClusters(track0b->ClIndex(),track1b->ClIndex());  // register clusters
4159         pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4160         pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4161       }
4162       pvertex->SetDistSigma(sigmad);
4163       //Bo:      pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);
4164       pvertex->SetNormDCAPrim(normdist[itrack0],normdist[itrack1]);
4165       //
4166       // define likelihhod and causalities
4167       //
4168       Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;
4169       if (maxLayer<1){
4170         Float_t fnorm0 = normdist[itrack0];
4171         if (fnorm0<0) fnorm0*=-3;
4172         Float_t fnorm1 = normdist[itrack1];
4173         if (fnorm1<0) fnorm1*=-3;
4174         if (pvertex->GetAnglep()[2]>0.1 ||  (pvertex->GetRr()<10.5)&& pvertex->GetAnglep()[2]>0.05 || pvertex->GetRr()<3){
4175           pb0    =  TMath::Exp(-TMath::Min(fnorm0,Float_t(16.))/12.);
4176           pb1    =  TMath::Exp(-TMath::Min(fnorm1,Float_t(16.))/12.);
4177         }
4178         pvertex->SetChi2Before(normdist[itrack0]);
4179         pvertex->SetChi2After(normdist[itrack1]);
4180         pvertex->SetNAfter(0);
4181         pvertex->SetNBefore(0);
4182       }else{
4183         pvertex->SetChi2Before(minchi2before0);
4184         pvertex->SetChi2After(minchi2before1);
4185          if (pvertex->GetAnglep()[2]>0.1 || ( pvertex->GetRr()<10.5 && pvertex->GetAnglep()[2]>0.05) || pvertex->GetRr()<3){
4186            pb0    =  TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.);
4187            pb1    =  TMath::Exp(-TMath::Min(minchi2before1,Float_t(16))/12.);
4188          }
4189          pvertex->SetNAfter(maxLayer);
4190          pvertex->SetNBefore(maxLayer);
4191       }
4192       if (pvertex->GetRr()<90){
4193         pa0  *= TMath::Min(track0->GetESDtrack()->GetTPCdensity(0,60),Float_t(1.));
4194         pa1  *= TMath::Min(track1->GetESDtrack()->GetTPCdensity(0,60),Float_t(1.));
4195       }
4196       if (pvertex->GetRr()<20){
4197         pa0  *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2;
4198         pa1  *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2;
4199       }
4200       //
4201       pvertex->SetCausality(pb0,pb1,pa0,pa1);
4202       //
4203       //  Likelihood calculations  - apply cuts
4204       //
4205       Bool_t v0OK = kTRUE;
4206       Float_t p12 = pvertex->GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4];
4207       p12        += pvertex->GetParamN()->GetParameter()[4]*pvertex->GetParamN()->GetParameter()[4];
4208       p12         = TMath::Sqrt(p12);                             // "mean" momenta
4209       Float_t    sigmap0   = 0.0001+0.001/(0.1+pvertex->GetRr());
4210       Float_t    sigmap    = 0.5*sigmap0*(0.6+0.4*p12);           // "resolution: of point angle - as a function of radius and momenta
4211
4212       Float_t causalityA  = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]);
4213       Float_t causalityB  = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Float_t(0.7))*
4214                                         TMath::Min(pvertex->GetCausalityP()[3],Float_t(0.7)));
4215       //
4216       //Bo:      Float_t likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5);
4217       Float_t lDistNorm = pvertex->GetDcaV0Daughters()/pvertex->GetDistSigma();
4218       Float_t likelihood0 = (TMath::Exp(-lDistNorm)+0.1) *(pvertex->GetDcaV0Daughters()<0.5)*(lDistNorm<5);
4219
4220       Float_t likelihood1 = TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/sigmap)+
4221         0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(4.*sigmap))+
4222         0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(8.*sigmap))+
4223         0.1*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/0.01);
4224       //
4225       if (causalityA<kCausality0Cut)                                          v0OK = kFALSE;
4226       if (TMath::Sqrt(likelihood0*likelihood1)<kLikelihood01Cut)              v0OK = kFALSE;
4227       if (likelihood1<kLikelihood1Cut)                                        v0OK = kFALSE;
4228       if (TMath::Power(likelihood0*likelihood1*causalityB,0.33)<kCombinedCut) v0OK = kFALSE;
4229
4230       //
4231       //
4232       if (kFALSE){
4233         Bool_t gold = TMath::Abs(TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel()))==1;
4234         cstream<<"It0"<<
4235           "Tr0.="<<track0<<                       //best without constrain
4236           "Tr1.="<<track1<<                       //best without constrain
4237           "Tr0B.="<<track0b<<                     //best without constrain  after vertex
4238           "Tr1B.="<<track1b<<                     //best without constrain  after vertex
4239           "Tr0C.="<<htrackc0<<                    //best with constrain     if exist
4240           "Tr1C.="<<htrackc1<<                    //best with constrain     if exist
4241           "Tr0L.="<<track0l<<                     //longest best
4242           "Tr1L.="<<track1l<<                     //longest best
4243           "Esd0.="<<track0->GetESDtrack()<<           // esd track0 params
4244           "Esd1.="<<track1->GetESDtrack()<<           // esd track1 params
4245           "V0.="<<pvertex<<                       //vertex properties
4246           "V0b.="<<&vertex2<<                       //vertex properties at "best" track
4247           "ND0="<<normdist[itrack0]<<             //normalize distance for track0
4248           "ND1="<<normdist[itrack1]<<             //normalize distance for track1
4249           "Gold="<<gold<<                         //
4250           //      "RejectBase="<<rejectBase<<             //rejection in First itteration
4251           "OK="<<v0OK<<
4252           "rmin="<<rmin<<
4253           "sigmad="<<sigmad<<
4254           "\n";
4255       }
4256       //if (rejectBase) continue;
4257       //
4258       pvertex->SetStatus(0);
4259       //      if (rejectBase) {
4260       //        pvertex->SetStatus(-100);
4261       //}
4262       if (pvertex->GetV0CosineOfPointingAngle()>kMinPointAngle2) {
4263         //Bo:     pvertex->SetESDindexes(track0->GetESDtrack()->GetID(),track1->GetESDtrack()->GetID());
4264         pvertex->SetIndex(0,track0->GetESDtrack()->GetID());//Bo: consistency 0 for neg
4265         pvertex->SetIndex(1,track1->GetESDtrack()->GetID());//Bo: consistency 1 for pos
4266         if (v0OK){
4267           //      AliV0vertex vertexjuri(*track0,*track1);
4268           //      vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
4269           //      event->AddV0(&vertexjuri);
4270           pvertex->SetStatus(100);
4271         }
4272         pvertex->SetOnFlyStatus(kTRUE);
4273         pvertex->ChangeMassHypothesis(kK0Short);
4274         event->AddV0(pvertex);
4275       }
4276     }
4277   }
4278   //
4279   //
4280   // delete temporary arrays
4281   //
4282   delete[] minPointAngle;
4283   delete[] maxr;
4284   delete[] minr;
4285   delete[] norm;
4286   delete[] normdist;
4287   delete[] normdist1;
4288   delete[] normdist0;
4289   delete[] dist;
4290   delete[] itsmap;
4291   delete[] helixes;
4292   delete   pvertex;
4293 }
4294
4295
4296 void AliITStrackerMI::RefitV02(AliESD *event)
4297 {
4298   //
4299   //try to refit  V0s in the third path of the reconstruction
4300   //
4301   TTreeSRedirector &cstream = *fDebugStreamer;
4302   //
4303   Int_t  nv0s = event->GetNumberOfV0s();
4304   Float_t primvertex[3]={GetX(),GetY(),GetZ()};
4305   AliV0 v0temp;
4306   for (Int_t iv0 = 0; iv0<nv0s;iv0++){
4307     AliV0 * v0mi = (AliV0*)event->GetV0(iv0);
4308     if (!v0mi) continue;
4309     Int_t     itrack0   = v0mi->GetIndex(0);
4310     Int_t     itrack1   = v0mi->GetIndex(1);
4311     AliESDtrack *esd0   = event->GetTrack(itrack0);
4312     AliESDtrack *esd1   = event->GetTrack(itrack1);
4313     if (!esd0||!esd1) continue;
4314     AliITStrackMI tpc0(*esd0);
4315     AliITStrackMI tpc1(*esd1);
4316     Double_t x,y,z; v0mi->GetXYZ(x,y,z); //I.B.
4317     Double_t alpha =TMath::ATan2(y,x);   //I.B.
4318     if (v0mi->GetRr()>85){
4319       if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4320         v0temp.SetParamN(tpc0);
4321         v0temp.SetParamP(tpc1);
4322         v0temp.Update(primvertex);
4323         if (kFALSE) cstream<<"Refit"<<
4324           "V0.="<<v0mi<<
4325           "V0refit.="<<&v0temp<<
4326           "Tr0.="<<&tpc0<<
4327           "Tr1.="<<&tpc1<<
4328           "\n";
4329         //Bo:   if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4330         if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4331           v0mi->SetParamN(tpc0);
4332           v0mi->SetParamP(tpc1);
4333           v0mi->Update(primvertex);
4334         }
4335       }
4336       continue;
4337     }
4338     if (v0mi->GetRr()>35){
4339       CorrectForDeadZoneMaterial(&tpc0);
4340       CorrectForDeadZoneMaterial(&tpc1);
4341       if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4342         v0temp.SetParamN(tpc0);
4343         v0temp.SetParamP(tpc1);
4344         v0temp.Update(primvertex);
4345         if (kFALSE) cstream<<"Refit"<<
4346           "V0.="<<v0mi<<
4347           "V0refit.="<<&v0temp<<
4348           "Tr0.="<<&tpc0<<
4349           "Tr1.="<<&tpc1<<
4350           "\n";
4351         //Bo:   if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4352         if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4353           v0mi->SetParamN(tpc0);
4354           v0mi->SetParamP(tpc1);
4355           v0mi->Update(primvertex);
4356         }
4357       }
4358       continue;
4359     }
4360     CorrectForDeadZoneMaterial(&tpc0);
4361     CorrectForDeadZoneMaterial(&tpc1);
4362     //    if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4363     if (RefitAt(v0mi->GetRr(),&tpc0, v0mi->GetClusters(0)) && RefitAt(v0mi->GetRr(),&tpc1, v0mi->GetClusters(1))){
4364       v0temp.SetParamN(tpc0);
4365       v0temp.SetParamP(tpc1);
4366       v0temp.Update(primvertex);
4367       if (kFALSE) cstream<<"Refit"<<
4368         "V0.="<<v0mi<<
4369         "V0refit.="<<&v0temp<<
4370         "Tr0.="<<&tpc0<<
4371         "Tr1.="<<&tpc1<<
4372         "\n";
4373       //Bo:      if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4374       if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4375         v0mi->SetParamN(tpc0);
4376         v0mi->SetParamP(tpc1);
4377         v0mi->Update(primvertex);
4378       }
4379     }
4380   }
4381 }
4382
4383
4384
4385
4386
4387
4388