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