ITS/AliITStrackerMI.cxx

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