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