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