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