ITS/AliITStrackerMI.cxx

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