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