ITS/AliITStrackerMI.cxx

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