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