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   AliCluster c(*cl);
3439   c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3440   c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3441
3442   return track->UpdateMI(&c,chi2,index);
3443 }
3444
3445 //------------------------------------------------------------------------
3446 void AliITStrackerMI::GetDCASigma(AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3447 {
3448   //
3449   //DCA sigmas parameterization
3450   //to be paramterized using external parameters in future
3451   //
3452   //
3453   sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3454   sigmaz   = 0.0110+4.37*TMath::Abs(track->GetC());
3455 }
3456 //------------------------------------------------------------------------
3457 void AliITStrackerMI::SignDeltas( TObjArray *ClusterArray, Float_t vz)
3458 {
3459   //
3460   //
3461   Int_t entries = ClusterArray->GetEntriesFast();
3462   if (entries<4) return;
3463   AliITSRecPoint* cluster = (AliITSRecPoint*)ClusterArray->At(0);
3464   Int_t layer = cluster->GetLayer();
3465   if (layer>1) return;
3466   Int_t index[10000];
3467   Int_t ncandidates=0;
3468   Float_t r = (layer>0)? 7:4;
3469   //
3470   for (Int_t i=0;i<entries;i++){
3471     AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(i);
3472     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3473     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3474     index[ncandidates] = i;  //candidate to belong to delta electron track
3475     ncandidates++;
3476     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3477       cl0->SetDeltaProbability(1);
3478     }
3479   }
3480   //
3481   //
3482   //
3483   for (Int_t i=0;i<ncandidates;i++){
3484     AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(index[i]);
3485     if (cl0->GetDeltaProbability()>0.8) continue;
3486     //
3487     Int_t ncl = 0;
3488     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3489     sumy=sumz=sumy2=sumyz=sumw=0.0;
3490     for (Int_t j=0;j<ncandidates;j++){
3491       if (i==j) continue;
3492       AliITSRecPoint* cl1 = (AliITSRecPoint*)ClusterArray->At(index[j]);
3493       //
3494       Float_t dz = cl0->GetZ()-cl1->GetZ();
3495       Float_t dy = cl0->GetY()-cl1->GetY();
3496       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3497         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3498         y[ncl] = cl1->GetY();
3499         z[ncl] = cl1->GetZ();
3500         sumy+= y[ncl]*weight;
3501         sumz+= z[ncl]*weight;
3502         sumy2+=y[ncl]*y[ncl]*weight;
3503         sumyz+=y[ncl]*z[ncl]*weight;
3504         sumw+=weight;
3505         ncl++;
3506       }
3507     }
3508     if (ncl<4) continue;
3509     Float_t det = sumw*sumy2  - sumy*sumy;
3510     Float_t delta=1000;
3511     if (TMath::Abs(det)>0.01){
3512       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
3513       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
3514       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3515     }
3516     else{
3517       Float_t z0  = sumyz/sumy;
3518       delta = TMath::Abs(cl0->GetZ()-z0);
3519     }
3520     if ( delta<0.05) {
3521       cl0->SetDeltaProbability(1-20.*delta);
3522     }
3523   }
3524 }
3525 //------------------------------------------------------------------------
3526 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3527 {
3528   //
3529   //
3530   track->UpdateESDtrack(flags);
3531   AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3532   if (oldtrack) delete oldtrack;
3533   track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3534   if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3535     printf("Problem\n");
3536   }
3537 }
3538 //------------------------------------------------------------------------
3539 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3540   //
3541   // Get nearest upper layer close to the point xr.
3542   // rough approximation
3543   //
3544   const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3545   Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3546   Int_t res =6;
3547   for (Int_t i=0;i<6;i++){
3548     if (radius<kRadiuses[i]){
3549       res =i;
3550       break;
3551     }
3552   }
3553   return res;
3554 }
3555 //------------------------------------------------------------------------
3556 void AliITStrackerMI::UpdateTPCV0(AliESDEvent *event){
3557   //
3558   //try to update, or reject TPC  V0s
3559   //
3560   Int_t nv0s = event->GetNumberOfV0s();
3561   Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3562
3563   for (Int_t i=0;i<nv0s;i++){
3564     AliESDv0 * vertex = event->GetV0(i);
3565     Int_t ip = vertex->GetIndex(0);
3566     Int_t im = vertex->GetIndex(1);
3567     //
3568     TObjArray * arrayp = (ip<nitstracks) ? (TObjArray*)fTrackHypothesys.At(ip):0;
3569     TObjArray * arraym = (im<nitstracks) ? (TObjArray*)fTrackHypothesys.At(im):0;
3570     AliITStrackMI * trackp = (arrayp!=0) ? (AliITStrackMI*)arrayp->At(0):0;
3571     AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0;
3572     //
3573     //
3574     if (trackp){
3575       if (trackp->GetNumberOfClusters()+trackp->GetNDeadZone()>5.5){
3576         if (trackp->GetConstrain()&&trackp->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3577         if (!trackp->GetConstrain()&&trackp->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3578       }
3579     }
3580
3581     if (trackm){
3582       if (trackm->GetNumberOfClusters()+trackm->GetNDeadZone()>5.5){
3583         if (trackm->GetConstrain()&&trackm->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3584         if (!trackm->GetConstrain()&&trackm->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3585       }
3586     }
3587     if (vertex->GetStatus()==-100) continue;
3588     //
3589     Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
3590     Int_t clayer = GetNearestLayer(xrp);                    //I.B.
3591     vertex->SetNBefore(clayer);        //
3592     vertex->SetChi2Before(9*clayer);   //
3593     vertex->SetNAfter(6-clayer);       //
3594     vertex->SetChi2After(0);           //
3595     //
3596     if (clayer >1 ){ // calculate chi2 before vertex
3597       Float_t chi2p = 0, chi2m=0;
3598       //
3599       if (trackp){
3600         for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3601           if (trackp->GetClIndex(ilayer)>0){
3602             chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3603               trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3604           }
3605           else{
3606             chi2p+=9;
3607           }
3608         }
3609       }else{
3610         chi2p = 9*clayer;
3611       }
3612       //
3613       if (trackm){
3614         for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3615           if (trackm->GetClIndex(ilayer)>0){
3616             chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3617               trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3618           }
3619           else{
3620             chi2m+=9;
3621           }
3622         }
3623       }else{
3624         chi2m = 9*clayer;
3625       }
3626       vertex->SetChi2Before(TMath::Min(chi2p,chi2m));
3627       if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10);  // track exist before vertex
3628     }
3629
3630     if (clayer < 5 ){ // calculate chi2 after vertex
3631       Float_t chi2p = 0, chi2m=0;
3632       //
3633       if (trackp&&TMath::Abs(trackp->GetTgl())<1.){
3634         for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3635           if (trackp->GetClIndex(ilayer)>0){
3636             chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3637               trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3638           }
3639           else{
3640             chi2p+=9;
3641           }
3642         }
3643       }else{
3644         chi2p = 0;
3645       }
3646       //
3647       if (trackm&&TMath::Abs(trackm->GetTgl())<1.){
3648         for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3649           if (trackm->GetClIndex(ilayer)>0){
3650             chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3651               trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3652           }
3653           else{
3654             chi2m+=9;
3655           }
3656         }
3657       }else{
3658         chi2m = 0;
3659       }
3660       vertex->SetChi2After(TMath::Max(chi2p,chi2m));
3661       if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20);  // track not found in ITS
3662     }
3663   }
3664   //
3665 }
3666 //------------------------------------------------------------------------
3667 void AliITStrackerMI::FindV02(AliESDEvent *event)
3668 {
3669   //
3670   // V0 finder
3671   //
3672   //  Cuts on DCA -  R dependent
3673   //          max distance DCA between 2 tracks cut
3674   //          maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist);
3675   //
3676   const Float_t kMaxDist0      = 0.1;
3677   const Float_t kMaxDist1      = 0.1;
3678   const Float_t kMaxDist       = 1;
3679   const Float_t kMinPointAngle  = 0.85;
3680   const Float_t kMinPointAngle2 = 0.99;
3681   const Float_t kMinR           = 0.5;
3682   const Float_t kMaxR           = 220;
3683   //const Float_t kCausality0Cut   = 0.19;
3684   //const Float_t kLikelihood01Cut = 0.25;
3685   //const Float_t kPointAngleCut   = 0.9996;
3686   const Float_t kCausality0Cut   = 0.19;
3687   const Float_t kLikelihood01Cut = 0.45;
3688   const Float_t kLikelihood1Cut  = 0.5;
3689   const Float_t kCombinedCut     = 0.55;
3690
3691   //
3692   //
3693   TTreeSRedirector &cstream = *fDebugStreamer;
3694   Int_t ntracks    = event->GetNumberOfTracks();
3695   Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3696   fOriginal.Expand(ntracks);
3697   fTrackHypothesys.Expand(ntracks);
3698   fBestHypothesys.Expand(ntracks);
3699   //
3700   AliHelix * helixes   = new AliHelix[ntracks+2];
3701   TObjArray trackarray(ntracks+2);     //array with tracks - with vertex constrain
3702   TObjArray trackarrayc(ntracks+2);    //array of "best    tracks" - without vertex constrain
3703   TObjArray trackarrayl(ntracks+2);    //array of "longest tracks" - without vertex constrain
3704   Bool_t * forbidden   = new Bool_t [ntracks+2];
3705   Int_t   *itsmap      = new Int_t  [ntracks+2];
3706   Float_t *dist        = new Float_t[ntracks+2];
3707   Float_t *normdist0   = new Float_t[ntracks+2];
3708   Float_t *normdist1   = new Float_t[ntracks+2];
3709   Float_t *normdist    = new Float_t[ntracks+2];
3710   Float_t *norm        = new Float_t[ntracks+2];
3711   Float_t *maxr        = new Float_t[ntracks+2];
3712   Float_t *minr        = new Float_t[ntracks+2];
3713   Float_t *minPointAngle= new Float_t[ntracks+2];
3714   //
3715   AliV0 *pvertex      = new AliV0;
3716   AliITStrackMI * dummy= new AliITStrackMI;
3717   dummy->SetLabel(0);
3718   AliITStrackMI  trackat0;    //temporary track for DCA calculation
3719   //
3720   Float_t primvertex[3]={GetX(),GetY(),GetZ()};
3721   //
3722   // make ITS -  ESD map
3723   //
3724   for (Int_t itrack=0;itrack<ntracks+2;itrack++) {
3725     itsmap[itrack]        = -1;
3726     forbidden[itrack]     = kFALSE;
3727     maxr[itrack]          = kMaxR;
3728     minr[itrack]          = kMinR;
3729     minPointAngle[itrack] = kMinPointAngle;
3730   }
3731   for (Int_t itrack=0;itrack<nitstracks;itrack++){
3732     AliITStrackMI * original =   (AliITStrackMI*)(fOriginal.At(itrack));
3733     Int_t           esdindex =   original->GetESDtrack()->GetID();
3734     itsmap[esdindex]         =   itrack;
3735   }
3736   //
3737   // create ITS tracks from ESD tracks if not done before
3738   //
3739   for (Int_t itrack=0;itrack<ntracks;itrack++){
3740     if (itsmap[itrack]>=0) continue;
3741     AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack)));
3742     //tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY());
3743     //tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ();
3744     tpctrack->GetDZ(GetX(),GetY(),GetZ(),tpctrack->GetDP());   //I.B.
3745     if (tpctrack->GetD(0)<20 && tpctrack->GetD(1)<20){
3746       // tracks which can reach inner part of ITS
3747       // propagate track to outer its volume - with correction for material
3748       CorrectForTPCtoITSDeadZoneMaterial(tpctrack);
3749     }
3750     itsmap[itrack] = nitstracks;
3751     fOriginal.AddAt(tpctrack,nitstracks);
3752     nitstracks++;
3753   }
3754   //
3755   // fill temporary arrays
3756   //
3757   for (Int_t itrack=0;itrack<ntracks;itrack++){
3758     AliESDtrack *  esdtrack = event->GetTrack(itrack);
3759     Int_t          itsindex = itsmap[itrack];
3760     AliITStrackMI *original = (AliITStrackMI*)fOriginal.At(itsmap[itrack]);
3761     if (!original) continue;
3762     AliITStrackMI *bestConst  = 0;
3763     AliITStrackMI *bestLong   = 0;
3764     AliITStrackMI *best       = 0;
3765     //
3766     //
3767     TObjArray * array    = (TObjArray*)  fTrackHypothesys.At(itsindex);
3768     Int_t       hentries = (array==0) ?  0 : array->GetEntriesFast();
3769     // Get best track with vertex constrain
3770     for (Int_t ih=0;ih<hentries;ih++){
3771       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
3772       if (!trackh->GetConstrain()) continue;
3773       if (!bestConst) bestConst = trackh;
3774       if (trackh->GetNumberOfClusters()>5.0){
3775         bestConst  = trackh;                         // full track -  with minimal chi2
3776         break;
3777       }
3778       if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone())  continue;
3779       bestConst = trackh;
3780       break;
3781     }
3782     // Get best long track without vertex constrain and best track without vertex constrain
3783     for (Int_t ih=0;ih<hentries;ih++){
3784       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
3785       if (trackh->GetConstrain()) continue;
3786       if (!best)     best     = trackh;
3787       if (!bestLong) bestLong = trackh;
3788       if (trackh->GetNumberOfClusters()>5.0){
3789         bestLong  = trackh;                         // full track -  with minimal chi2
3790         break;
3791       }
3792       if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone())  continue;
3793       bestLong = trackh;
3794     }
3795     if (!best) {
3796       best     = original;
3797       bestLong = original;
3798     }
3799     //I.B. trackat0 = *bestLong;
3800     new (&trackat0) AliITStrackMI(*bestLong);
3801     Double_t xx,yy,zz,alpha;
3802     bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz);
3803     alpha = TMath::ATan2(yy,xx);
3804     trackat0.Propagate(alpha,0);
3805     // calculate normalized distances to the vertex
3806     //
3807     Float_t ptfac  = (1.+100.*TMath::Abs(trackat0.GetC()));
3808     if ( bestLong->GetNumberOfClusters()>3 ){
3809       dist[itsindex]      = trackat0.GetY();
3810       norm[itsindex]      = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
3811       normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
3812       normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
3813       normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
3814       if (!bestConst){
3815         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<6) normdist[itsindex]*=2.;
3816         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<5) normdist[itsindex]*=2.;
3817         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<4) normdist[itsindex]*=2.;
3818       }else{
3819         if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<6) normdist[itsindex]*=1.5;
3820         if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<5) normdist[itsindex]*=1.5;
3821       }
3822     }
3823     else{
3824       if (bestConst&&bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5){
3825         dist[itsindex] = bestConst->GetD(0);
3826         norm[itsindex] = bestConst->GetDnorm(0);
3827         normdist0[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
3828         normdist1[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
3829         normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
3830       }else{
3831         dist[itsindex]      = trackat0.GetY();
3832         norm[itsindex]      = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
3833         normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
3834         normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
3835         normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
3836         if (TMath::Abs(trackat0.GetTgl())>1.05){
3837           if (normdist[itsindex]<3) forbidden[itsindex]=kTRUE;
3838           if (normdist[itsindex]>3) {
3839             minr[itsindex] = TMath::Max(Float_t(40.),minr[itsindex]);
3840           }
3841         }
3842       }
3843     }
3844     //
3845     //-----------------------------------------------------------
3846     //Forbid primary track candidates -
3847     //
3848     //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5");
3849     //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5");
3850     //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0");
3851     //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0");
3852     //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2");
3853     //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1");
3854     //-----------------------------------------------------------
3855     if (bestConst){
3856       if (bestLong->GetNumberOfClusters()<4       && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5)               forbidden[itsindex]=kTRUE;
3857       if (normdist[itsindex]<3 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>5.5)               forbidden[itsindex]=kTRUE;
3858       if (normdist[itsindex]<2 && bestConst->GetClIndex(0)>0 && bestConst->GetClIndex(1)>0 ) forbidden[itsindex]=kTRUE;
3859       if (normdist[itsindex]<1 && bestConst->GetClIndex(0)>0)                              forbidden[itsindex]=kTRUE;
3860       if (normdist[itsindex]<4 && bestConst->GetNormChi2(0)<2)                             forbidden[itsindex]=kTRUE;
3861       if (normdist[itsindex]<5 && bestConst->GetNormChi2(0)<1)                             forbidden[itsindex]=kTRUE;
3862       if (bestConst->GetNormChi2(0)<2.5) {
3863         minPointAngle[itsindex]= 0.9999;
3864         maxr[itsindex]         = 10;
3865       }
3866     }
3867     //
3868     //forbid daughter kink candidates
3869     //
3870     if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
3871     Bool_t isElectron = kTRUE;
3872     Bool_t isProton   = kTRUE;
3873     Double_t pid[5];
3874     esdtrack->GetESDpid(pid);
3875     for (Int_t i=1;i<5;i++){
3876       if (pid[0]<pid[i]) isElectron= kFALSE;
3877       if (pid[4]<pid[i]) isProton= kFALSE;
3878     }
3879     if (isElectron){
3880       forbidden[itsindex]=kFALSE;
3881       normdist[itsindex]*=-1;
3882     }
3883     if (isProton){
3884       if (normdist[itsindex]>2) forbidden[itsindex]=kFALSE;
3885       normdist[itsindex]*=-1;
3886     }
3887
3888     //
3889     // Causality cuts in TPC volume
3890     //
3891     if (esdtrack->GetTPCdensity(0,10) >0.6)  maxr[itsindex] = TMath::Min(Float_t(110),maxr[itsindex]);
3892     if (esdtrack->GetTPCdensity(10,30)>0.6)  maxr[itsindex] = TMath::Min(Float_t(120),maxr[itsindex]);
3893     if (esdtrack->GetTPCdensity(20,40)>0.6)  maxr[itsindex] = TMath::Min(Float_t(130),maxr[itsindex]);
3894     if (esdtrack->GetTPCdensity(30,50)>0.6)  maxr[itsindex] = TMath::Min(Float_t(140),maxr[itsindex]);
3895     //
3896     if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->GetNumberOfClusters()<3) minr[itsindex]=100;
3897     //
3898     //
3899     if (kFALSE){
3900       cstream<<"Track"<<
3901         "Tr0.="<<best<<
3902         "Tr1.="<<((bestConst)? bestConst:dummy)<<
3903         "Tr2.="<<bestLong<<
3904         "Tr3.="<<&trackat0<<
3905         "Esd.="<<esdtrack<<
3906         "Dist="<<dist[itsindex]<<
3907         "ND0="<<normdist0[itsindex]<<
3908         "ND1="<<normdist1[itsindex]<<
3909         "ND="<<normdist[itsindex]<<
3910         "Pz="<<primvertex[2]<<
3911         "Forbid="<<forbidden[itsindex]<<
3912         "\n";
3913       //
3914     }
3915     trackarray.AddAt(best,itsindex);
3916     trackarrayc.AddAt(bestConst,itsindex);
3917     trackarrayl.AddAt(bestLong,itsindex);
3918     new (&helixes[itsindex]) AliHelix(*best);
3919   }
3920   //
3921   //
3922   //
3923   // first iterration of V0 finder
3924   //
3925   for (Int_t iesd0=0;iesd0<ntracks;iesd0++){
3926     Int_t itrack0 = itsmap[iesd0];
3927     if (forbidden[itrack0]) continue;
3928     AliITStrackMI * btrack0 = (AliITStrackMI*)trackarray.At(itrack0);
3929     if (!btrack0) continue;
3930     if (btrack0->GetSign()>0) continue;
3931     AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
3932     //
3933     for (Int_t iesd1=0;iesd1<ntracks;iesd1++){
3934       Int_t itrack1 = itsmap[iesd1];
3935       if (forbidden[itrack1]) continue;
3936
3937       AliITStrackMI * btrack1 = (AliITStrackMI*)trackarray.At(itrack1);
3938       if (!btrack1) continue;
3939       if (btrack1->GetSign()<0) continue;
3940       Bool_t isGold = kFALSE;
3941       if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){
3942         isGold = kTRUE;
3943       }
3944       AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
3945       AliHelix &h1 = helixes[itrack0];
3946       AliHelix &h2 = helixes[itrack1];
3947       //
3948       // find linear distance
3949       Double_t rmin =0;
3950       //
3951       //
3952       //
3953       Double_t phase[2][2],radius[2];
3954       Int_t  points = h1.GetRPHIintersections(h2, phase, radius);
3955       if    (points==0)  continue;
3956       Double_t delta[2]={1000000,1000000};
3957       rmin = radius[0];
3958       h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]);
3959       if (points==2){
3960         if (radius[1]<rmin) rmin = radius[1];
3961         h1.ParabolicDCA(h2,phase[1][0],phase[1][1],radius[1],delta[1]);
3962       }
3963       rmin = TMath::Sqrt(rmin);
3964       Double_t distance = 0;
3965       Double_t radiusC  = 0;
3966       Int_t    iphase   = 0;
3967       if (points==1 || delta[0]<delta[1]){
3968         distance = TMath::Sqrt(delta[0]);
3969         radiusC  = TMath::Sqrt(radius[0]);
3970       }else{
3971         distance = TMath::Sqrt(delta[1]);
3972         radiusC  = TMath::Sqrt(radius[1]);
3973         iphase=1;
3974       }
3975       if (radiusC<TMath::Max(minr[itrack0],minr[itrack1]))    continue;
3976       if (radiusC>TMath::Min(maxr[itrack0],maxr[itrack1]))     continue;
3977       Float_t maxDist  = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1));
3978       if (distance>maxDist) continue;
3979       Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex);
3980       if (pointAngle<TMath::Max(minPointAngle[itrack0],minPointAngle[itrack1])) continue;
3981       //
3982       //
3983       //      Double_t distance = TestV0(h1,h2,pvertex,rmin);
3984       //
3985       //       if (distance>maxDist)           continue;
3986       //       if (pvertex->GetRr()<kMinR)     continue;
3987       //       if (pvertex->GetRr()>kMaxR)     continue;
3988       AliITStrackMI * track0=btrack0;
3989       AliITStrackMI * track1=btrack1;
3990       //      if (pvertex->GetRr()<3.5){
3991       if (radiusC<3.5){
3992         //use longest tracks inside the pipe
3993         track0 = (AliITStrackMI*)trackarrayl.At(itrack0);
3994         track1 = (AliITStrackMI*)trackarrayl.At(itrack1);
3995       }
3996       //
3997       //
3998       pvertex->SetParamN(*track0);
3999       pvertex->SetParamP(*track1);
4000       pvertex->Update(primvertex);
4001       pvertex->SetClusters(track0->ClIndex(),track1->ClIndex());  // register clusters
4002
4003       if (pvertex->GetRr()<kMinR) continue;
4004       if (pvertex->GetRr()>kMaxR) continue;
4005       if (pvertex->GetV0CosineOfPointingAngle()<kMinPointAngle) continue;
4006 //Bo:      if (pvertex->GetDist2()>maxDist) continue;
4007       if (pvertex->GetDcaV0Daughters()>maxDist) continue;
4008 //Bo:        pvertex->SetLab(0,track0->GetLabel());
4009 //Bo:        pvertex->SetLab(1,track1->GetLabel());
4010       pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4011       pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4012       //
4013       AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
4014       AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
4015
4016       //
4017       //
4018       TObjArray * array0b     = (TObjArray*)fBestHypothesys.At(itrack0);
4019       if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->GetTgl())<1.1) {
4020         fCurrentEsdTrack = itrack0;
4021         FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack0),itrack0, kFALSE);
4022       }
4023       TObjArray * array1b    = (TObjArray*)fBestHypothesys.At(itrack1);
4024       if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->GetTgl())<1.1) {
4025         fCurrentEsdTrack = itrack1;
4026         FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack1),itrack1, kFALSE);
4027       }
4028       //
4029       AliITStrackMI * track0b = (AliITStrackMI*)fOriginal.At(itrack0);
4030       AliITStrackMI * track1b = (AliITStrackMI*)fOriginal.At(itrack1);
4031       AliITStrackMI * track0l = (AliITStrackMI*)fOriginal.At(itrack0);
4032       AliITStrackMI * track1l = (AliITStrackMI*)fOriginal.At(itrack1);
4033
4034       Float_t minchi2before0=16;
4035       Float_t minchi2before1=16;
4036       Float_t minchi2after0 =16;
4037       Float_t minchi2after1 =16;
4038       Double_t xrp[3]; pvertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
4039       Int_t maxLayer = GetNearestLayer(xrp);                   //I.B.
4040
4041       if (array0b) for (Int_t i=0;i<5;i++){
4042         // best track after vertex
4043         AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i);
4044         if (!btrack) continue;
4045         if (btrack->GetNumberOfClusters()>track0l->GetNumberOfClusters()) track0l = btrack;
4046         //      if (btrack->fX<pvertex->GetRr()-2.-0.5/(0.1+pvertex->GetAnglep()[2])) {
4047         if (btrack->GetX()<pvertex->GetRr()-2.) {
4048           if ( (maxLayer>i+2|| (i==0)) && btrack->GetNumberOfClusters()==(6-i)&&i<3){
4049             Float_t sumchi2= 0;
4050             Float_t sumn   = 0;
4051             if (maxLayer<3){   // take prim vertex as additional measurement
4052               if (normdist[itrack0]>0 && htrackc0){
4053                 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack0]*normdist[itrack0],16.);
4054               }else{
4055                 sumchi2 +=  TMath::Min((3.-maxLayer)*(3*normdist[itrack0]*normdist[itrack0]+3.),16.);
4056               }
4057               sumn    +=  3-maxLayer;
4058             }
4059             for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4060               sumn+=1.;
4061               if (!btrack->GetClIndex(ilayer)){
4062                 sumchi2+=25;
4063                 continue;
4064               }else{
4065                 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4066                 for (Int_t itrack=0;itrack<4;itrack++){
4067                   if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack0){
4068                     sumchi2+=18.;  //shared cluster
4069                     break;
4070                   }
4071                 }
4072                 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4073                 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4074               }
4075             }
4076             sumchi2/=sumn;
4077             if (sumchi2<minchi2before0) minchi2before0=sumchi2;
4078           }
4079           continue;   //safety space - Geo manager will give exact layer
4080         }
4081         track0b       = btrack;
4082         minchi2after0 = btrack->GetNormChi2(i);
4083         break;
4084       }
4085       if (array1b) for (Int_t i=0;i<5;i++){
4086         // best track after vertex
4087         AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i);
4088         if (!btrack) continue;
4089         if (btrack->GetNumberOfClusters()>track1l->GetNumberOfClusters()) track1l = btrack;
4090         //      if (btrack->fX<pvertex->GetRr()-2-0.5/(0.1+pvertex->GetAnglep()[2])){
4091         if (btrack->GetX()<pvertex->GetRr()-2){
4092           if ((maxLayer>i+2 || (i==0))&&btrack->GetNumberOfClusters()==(6-i)&&(i<3)){
4093             Float_t sumchi2= 0;
4094             Float_t sumn   = 0;
4095             if (maxLayer<3){   // take prim vertex as additional measurement
4096               if (normdist[itrack1]>0 && htrackc1){
4097                 sumchi2 +=  TMath::Min((3.-maxLayer)*normdist[itrack1]*normdist[itrack1],16.);
4098               }else{
4099                 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack1]*normdist[itrack1]+3.),16.);
4100               }
4101               sumn    +=  3-maxLayer;
4102             }
4103             for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4104               sumn+=1.;
4105               if (!btrack->GetClIndex(ilayer)){
4106                 sumchi2+=30;
4107                 continue;
4108               }else{
4109                 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4110                 for (Int_t itrack=0;itrack<4;itrack++){
4111                   if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack1){
4112                     sumchi2+=18.;  //shared cluster
4113                     break;
4114                   }
4115                 }
4116                 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4117                 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4118               }
4119             }
4120             sumchi2/=sumn;
4121             if (sumchi2<minchi2before1) minchi2before1=sumchi2;
4122           }
4123           continue;   //safety space - Geo manager will give exact layer
4124         }
4125         track1b = btrack;
4126         minchi2after1 = btrack->GetNormChi2(i);
4127         break;
4128       }
4129       //
4130       // position resolution - used for DCA cut
4131       Float_t sigmad = track0b->GetSigmaY2()+track0b->GetSigmaZ2()+track1b->GetSigmaY2()+track1b->GetSigmaZ2()+
4132         (track0b->GetX()-pvertex->GetRr())*(track0b->GetX()-pvertex->GetRr())*(track0b->GetSigmaSnp2()+track0b->GetSigmaTgl2())+
4133         (track1b->GetX()-pvertex->GetRr())*(track1b->GetX()-pvertex->GetRr())*(track1b->GetSigmaSnp2()+track1b->GetSigmaTgl2());
4134       sigmad =TMath::Sqrt(sigmad)+0.04;
4135       if (pvertex->GetRr()>50){
4136         Double_t cov0[15],cov1[15];
4137         track0b->GetESDtrack()->GetInnerExternalCovariance(cov0);
4138         track1b->GetESDtrack()->GetInnerExternalCovariance(cov1);
4139         sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+
4140           (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+
4141           (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]);
4142         sigmad =TMath::Sqrt(sigmad)+0.05;
4143       }
4144       //
4145       AliV0 vertex2;
4146       vertex2.SetParamN(*track0b);
4147       vertex2.SetParamP(*track1b);
4148       vertex2.Update(primvertex);
4149       //Bo:      if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4150       if (vertex2.GetDcaV0Daughters()<=pvertex->GetDcaV0Daughters()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4151         pvertex->SetParamN(*track0b);
4152         pvertex->SetParamP(*track1b);
4153         pvertex->Update(primvertex);
4154         pvertex->SetClusters(track0b->ClIndex(),track1b->ClIndex());  // register clusters
4155         pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4156         pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4157       }
4158       pvertex->SetDistSigma(sigmad);
4159       //Bo:      pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);
4160       pvertex->SetNormDCAPrim(normdist[itrack0],normdist[itrack1]);
4161       //
4162       // define likelihhod and causalities
4163       //
4164       Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;
4165       if (maxLayer<1){
4166         Float_t fnorm0 = normdist[itrack0];
4167         if (fnorm0<0) fnorm0*=-3;
4168         Float_t fnorm1 = normdist[itrack1];
4169         if (fnorm1<0) fnorm1*=-3;
4170         if (pvertex->GetAnglep()[2]>0.1 ||  (pvertex->GetRr()<10.5)&& pvertex->GetAnglep()[2]>0.05 || pvertex->GetRr()<3){
4171           pb0    =  TMath::Exp(-TMath::Min(fnorm0,Float_t(16.))/12.);
4172           pb1    =  TMath::Exp(-TMath::Min(fnorm1,Float_t(16.))/12.);
4173         }
4174         pvertex->SetChi2Before(normdist[itrack0]);
4175         pvertex->SetChi2After(normdist[itrack1]);
4176         pvertex->SetNAfter(0);
4177         pvertex->SetNBefore(0);
4178       }else{
4179         pvertex->SetChi2Before(minchi2before0);
4180         pvertex->SetChi2After(minchi2before1);
4181          if (pvertex->GetAnglep()[2]>0.1 || ( pvertex->GetRr()<10.5 && pvertex->GetAnglep()[2]>0.05) || pvertex->GetRr()<3){
4182            pb0    =  TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.);
4183            pb1    =  TMath::Exp(-TMath::Min(minchi2before1,Float_t(16))/12.);
4184          }
4185          pvertex->SetNAfter(maxLayer);
4186          pvertex->SetNBefore(maxLayer);
4187       }
4188       if (pvertex->GetRr()<90){
4189         pa0  *= TMath::Min(track0->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4190         pa1  *= TMath::Min(track1->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4191       }
4192       if (pvertex->GetRr()<20){
4193         pa0  *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2;
4194         pa1  *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2;
4195       }
4196       //
4197       pvertex->SetCausality(pb0,pb1,pa0,pa1);
4198       //
4199       //  Likelihood calculations  - apply cuts
4200       //
4201       Bool_t v0OK = kTRUE;
4202       Float_t p12 = pvertex->GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4];
4203       p12        += pvertex->GetParamN()->GetParameter()[4]*pvertex->GetParamN()->GetParameter()[4];
4204       p12         = TMath::Sqrt(p12);                             // "mean" momenta
4205       Float_t    sigmap0   = 0.0001+0.001/(0.1+pvertex->GetRr());
4206       Float_t    sigmap    = 0.5*sigmap0*(0.6+0.4*p12);           // "resolution: of point angle - as a function of radius and momenta
4207
4208       Float_t causalityA  = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]);
4209       Float_t causalityB  = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Double_t(0.7))*
4210                                         TMath::Min(pvertex->GetCausalityP()[3],Double_t(0.7)));
4211       //
4212       //Bo:      Float_t likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5);
4213       Float_t lDistNorm = pvertex->GetDcaV0Daughters()/pvertex->GetDistSigma();
4214       Float_t likelihood0 = (TMath::Exp(-lDistNorm)+0.1) *(pvertex->GetDcaV0Daughters()<0.5)*(lDistNorm<5);
4215
4216       Float_t likelihood1 = TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/sigmap)+
4217         0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(4.*sigmap))+
4218         0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(8.*sigmap))+
4219         0.1*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/0.01);
4220       //
4221       if (causalityA<kCausality0Cut)                                          v0OK = kFALSE;
4222       if (TMath::Sqrt(likelihood0*likelihood1)<kLikelihood01Cut)              v0OK = kFALSE;
4223       if (likelihood1<kLikelihood1Cut)                                        v0OK = kFALSE;
4224       if (TMath::Power(likelihood0*likelihood1*causalityB,0.33)<kCombinedCut) v0OK = kFALSE;
4225
4226       //
4227       //
4228       if (kFALSE){
4229         Bool_t gold = TMath::Abs(TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel()))==1;
4230         cstream<<"It0"<<
4231           "Tr0.="<<track0<<                       //best without constrain
4232           "Tr1.="<<track1<<                       //best without constrain
4233           "Tr0B.="<<track0b<<                     //best without constrain  after vertex
4234           "Tr1B.="<<track1b<<                     //best without constrain  after vertex
4235           "Tr0C.="<<htrackc0<<                    //best with constrain     if exist
4236           "Tr1C.="<<htrackc1<<                    //best with constrain     if exist
4237           "Tr0L.="<<track0l<<                     //longest best
4238           "Tr1L.="<<track1l<<                     //longest best
4239           "Esd0.="<<track0->GetESDtrack()<<           // esd track0 params
4240           "Esd1.="<<track1->GetESDtrack()<<           // esd track1 params
4241           "V0.="<<pvertex<<                       //vertex properties
4242           "V0b.="<<&vertex2<<                       //vertex properties at "best" track
4243           "ND0="<<normdist[itrack0]<<             //normalize distance for track0
4244           "ND1="<<normdist[itrack1]<<             //normalize distance for track1
4245           "Gold="<<gold<<                         //
4246           //      "RejectBase="<<rejectBase<<             //rejection in First itteration
4247           "OK="<<v0OK<<
4248           "rmin="<<rmin<<
4249           "sigmad="<<sigmad<<
4250           "\n";
4251       }
4252       //if (rejectBase) continue;
4253       //
4254       pvertex->SetStatus(0);
4255       //      if (rejectBase) {
4256       //        pvertex->SetStatus(-100);
4257       //}
4258       if (pvertex->GetV0CosineOfPointingAngle()>kMinPointAngle2) {
4259         //Bo:     pvertex->SetESDindexes(track0->GetESDtrack()->GetID(),track1->GetESDtrack()->GetID());
4260         pvertex->SetIndex(0,track0->GetESDtrack()->GetID());//Bo: consistency 0 for neg
4261         pvertex->SetIndex(1,track1->GetESDtrack()->GetID());//Bo: consistency 1 for pos
4262         if (v0OK){
4263           //      AliV0vertex vertexjuri(*track0,*track1);
4264           //      vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
4265           //      event->AddV0(&vertexjuri);
4266           pvertex->SetStatus(100);
4267         }
4268         pvertex->SetOnFlyStatus(kTRUE);
4269         pvertex->ChangeMassHypothesis(kK0Short);
4270         event->AddV0(pvertex);
4271       }
4272     }
4273   }
4274   //
4275   //
4276   // delete temporary arrays
4277   //
4278   delete[] forbidden;
4279   delete[] minPointAngle;
4280   delete[] maxr;
4281   delete[] minr;
4282   delete[] norm;
4283   delete[] normdist;
4284   delete[] normdist1;
4285   delete[] normdist0;
4286   delete[] dist;
4287   delete[] itsmap;
4288   delete[] helixes;
4289   delete   pvertex;
4290 }
4291 //------------------------------------------------------------------------
4292 void AliITStrackerMI::RefitV02(AliESDEvent *event)
4293 {
4294   //
4295   //try to refit  V0s in the third path of the reconstruction
4296   //
4297   TTreeSRedirector &cstream = *fDebugStreamer;
4298   //
4299   Int_t  nv0s = event->GetNumberOfV0s();
4300   Float_t primvertex[3]={GetX(),GetY(),GetZ()};
4301   AliV0 v0temp;
4302   for (Int_t iv0 = 0; iv0<nv0s;iv0++){
4303     AliV0 * v0mi = (AliV0*)event->GetV0(iv0);
4304     if (!v0mi) continue;
4305     Int_t     itrack0   = v0mi->GetIndex(0);
4306     Int_t     itrack1   = v0mi->GetIndex(1);
4307     AliESDtrack *esd0   = event->GetTrack(itrack0);
4308     AliESDtrack *esd1   = event->GetTrack(itrack1);
4309     if (!esd0||!esd1) continue;
4310     AliITStrackMI tpc0(*esd0);
4311     AliITStrackMI tpc1(*esd1);
4312     Double_t x,y,z; v0mi->GetXYZ(x,y,z); //I.B.
4313     Double_t alpha =TMath::ATan2(y,x);   //I.B.
4314     if (v0mi->GetRr()>85){
4315       if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4316         v0temp.SetParamN(tpc0);
4317         v0temp.SetParamP(tpc1);
4318         v0temp.Update(primvertex);
4319         if (kFALSE) cstream<<"Refit"<<
4320           "V0.="<<v0mi<<
4321           "V0refit.="<<&v0temp<<
4322           "Tr0.="<<&tpc0<<
4323           "Tr1.="<<&tpc1<<
4324           "\n";
4325         //Bo:   if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4326         if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4327           v0mi->SetParamN(tpc0);
4328           v0mi->SetParamP(tpc1);
4329           v0mi->Update(primvertex);
4330         }
4331       }
4332       continue;
4333     }
4334     if (v0mi->GetRr()>35){
4335       CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4336       CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4337       if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4338         v0temp.SetParamN(tpc0);
4339         v0temp.SetParamP(tpc1);
4340         v0temp.Update(primvertex);
4341         if (kFALSE) cstream<<"Refit"<<
4342           "V0.="<<v0mi<<
4343           "V0refit.="<<&v0temp<<
4344           "Tr0.="<<&tpc0<<
4345           "Tr1.="<<&tpc1<<
4346           "\n";
4347         //Bo:   if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4348         if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4349           v0mi->SetParamN(tpc0);
4350           v0mi->SetParamP(tpc1);
4351           v0mi->Update(primvertex);
4352         }
4353       }
4354       continue;
4355     }
4356     CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4357     CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4358     //    if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4359     if (RefitAt(v0mi->GetRr(),&tpc0, v0mi->GetClusters(0)) && RefitAt(v0mi->GetRr(),&tpc1, v0mi->GetClusters(1))){
4360       v0temp.SetParamN(tpc0);
4361       v0temp.SetParamP(tpc1);
4362       v0temp.Update(primvertex);
4363       if (kFALSE) cstream<<"Refit"<<
4364         "V0.="<<v0mi<<
4365         "V0refit.="<<&v0temp<<
4366         "Tr0.="<<&tpc0<<
4367         "Tr1.="<<&tpc1<<
4368         "\n";
4369       //Bo:      if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4370       if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4371         v0mi->SetParamN(tpc0);
4372         v0mi->SetParamP(tpc1);
4373         v0mi->Update(primvertex);
4374       }
4375     }
4376   }
4377 }
4378 //------------------------------------------------------------------------
4379 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4380   //--------------------------------------------------------------------
4381   // Fill a look-up table with mean material
4382   //--------------------------------------------------------------------
4383
4384   Int_t n=1000;
4385   Double_t mparam[7];
4386   Double_t point1[3],point2[3];
4387   Double_t phi,cosphi,sinphi,z;
4388   // 0-5 layers, 6 pipe, 7-8 shields
4389   Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 7.5,25.0};
4390   Double_t rmax[9]={ 5.5, 7.3,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4391
4392   Int_t ifirst=0,ilast=0;
4393   if(material.Contains("Pipe")) {
4394     ifirst=6; ilast=6;
4395   } else if(material.Contains("Shields")) {
4396     ifirst=7; ilast=8;
4397   } else if(material.Contains("Layers")) {
4398     ifirst=0; ilast=5;
4399   } else {
4400     Error("BuildMaterialLUT","Wrong layer name\n");
4401   }
4402
4403   for(Int_t imat=ifirst; imat<=ilast; imat++) {
4404     Double_t param[5]={0.,0.,0.,0.,0.};
4405     for (Int_t i=0; i<n; i++) {
4406       phi = 2.*TMath::Pi()*gRandom->Rndm();
4407       cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4408       z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4409       point1[0] = rmin[imat]*cosphi;
4410       point1[1] = rmin[imat]*sinphi;
4411       point1[2] = z;
4412       point2[0] = rmax[imat]*cosphi;
4413       point2[1] = rmax[imat]*sinphi;
4414       point2[2] = z;
4415       AliTracker::MeanMaterialBudget(point1,point2,mparam);
4416       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4417     }
4418     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4419     if(imat<=5) {
4420       fxOverX0Layer[imat] = param[1];
4421       fxTimesRhoLayer[imat] = param[0]*param[4];
4422     } else if(imat==6) {
4423       fxOverX0Pipe = param[1];
4424       fxTimesRhoPipe = param[0]*param[4];
4425     } else if(imat==7) {
4426       fxOverX0Shield[0] = param[1];
4427       fxTimesRhoShield[0] = param[0]*param[4];
4428     } else if(imat==8) {
4429       fxOverX0Shield[1] = param[1];
4430       fxTimesRhoShield[1] = param[0]*param[4];
4431     }
4432   }
4433   /*
4434   printf("%s\n",material.Data());
4435   printf("%f  %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4436   printf("%f  %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4437   printf("%f  %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4438   printf("%f  %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4439   printf("%f  %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4440   printf("%f  %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4441   printf("%f  %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4442   printf("%f  %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4443   printf("%f  %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4444   */
4445   return;
4446 }
4447 //------------------------------------------------------------------------
4448 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4449                                               TString direction) {
4450   //-------------------------------------------------------------------
4451   // Propagate beyond beam pipe and correct for material
4452   // (material budget in different ways according to fUseTGeo value)
4453   //-------------------------------------------------------------------
4454
4455   // Define budget mode:
4456   // 0: material from AliITSRecoParam (hard coded)
4457   // 1: material from TGeo (on the fly)
4458   // 2: material from lut
4459   // 3: material from TGeo (same for all hypotheses)
4460   Int_t mode;
4461   switch(fUseTGeo) {
4462   case 0:
4463     mode=0;
4464     break;
4465   case 1:
4466     mode=1;
4467     break;
4468   case 2:
4469     mode=2;
4470     break;
4471   case 3:
4472     if(fTrackingPhase.Contains("Clusters2Tracks"))
4473       { mode=3; } else { mode=1; }
4474     break;
4475   case 4:
4476     if(fTrackingPhase.Contains("Clusters2Tracks"))
4477       { mode=3; } else { mode=2; }
4478     break;
4479   default:
4480     mode=0;
4481     break;
4482   }
4483   if(fTrackingPhase.Contains("Default")) mode=0;
4484
4485   Int_t index=fCurrentEsdTrack;
4486
4487   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4488   Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4489   Double_t xToGo; t->GetLocalXat(rToGo,xToGo);
4490
4491   Double_t xOverX0,x0,lengthTimesMeanDensity;
4492   Bool_t anglecorr=kTRUE;
4493
4494   switch(mode) {
4495   case 0:
4496     xOverX0 = AliITSRecoParam::GetdPipe();
4497     x0 = AliITSRecoParam::GetX0Be();
4498     lengthTimesMeanDensity = xOverX0*x0;
4499     break;
4500   case 1:
4501     if (!t->PropagateToTGeo(xToGo,1)) return 0;
4502     return 1;
4503     break;
4504   case 2:
4505     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4506     xOverX0 = fxOverX0Pipe;
4507     lengthTimesMeanDensity = fxTimesRhoPipe;
4508     break;
4509   case 3:
4510     if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4511     if(fxOverX0PipeTrks[index]<0) {
4512       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4513       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4514                                  (1.-t->GetSnp()*t->GetSnp()));
4515       fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4516       fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4517       return 1;
4518     }
4519     xOverX0 = fxOverX0PipeTrks[index];
4520     lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4521     break;
4522   }
4523
4524   lengthTimesMeanDensity *= dir;
4525
4526   if (!t->AliExternalTrackParam::PropagateTo(xToGo,GetBz())) return 0;
4527   if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4528
4529   return 1;
4530 }
4531 //------------------------------------------------------------------------
4532 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4533                                                 TString shield,
4534                                                 TString direction) {
4535   //-------------------------------------------------------------------
4536   // Propagate beyond SPD or SDD shield and correct for material
4537   // (material budget in different ways according to fUseTGeo value)
4538   //-------------------------------------------------------------------
4539
4540   // Define budget mode:
4541   // 0: material from AliITSRecoParam (hard coded)
4542   // 1: material from TGeo (on the fly)
4543   // 2: material from lut
4544   // 3: material from TGeo (same for all hypotheses)
4545   Int_t mode;
4546   switch(fUseTGeo) {
4547   case 0:
4548     mode=0;
4549     break;
4550   case 1:
4551     mode=1;
4552     break;
4553   case 2:
4554     mode=2;
4555     break;
4556   case 3:
4557     if(fTrackingPhase.Contains("Clusters2Tracks"))
4558       { mode=3; } else { mode=1; }
4559     break;
4560   case 4:
4561     if(fTrackingPhase.Contains("Clusters2Tracks"))
4562       { mode=3; } else { mode=2; }
4563     break;
4564   default:
4565     mode=0;
4566     break;
4567   }
4568   if(fTrackingPhase.Contains("Default")) mode=0;
4569
4570   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4571   Double_t rToGo;
4572   Int_t    shieldindex=0;
4573   if (shield.Contains("SDD")) { // SDDouter
4574     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4575     shieldindex=1;
4576   } else if (shield.Contains("SPD")) {        // SPDouter
4577     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4578     shieldindex=0;
4579   } else {
4580     Error("CorrectForShieldMaterial"," Wrong shield name\n");
4581     return 0;
4582   }
4583   Double_t xToGo; t->GetLocalXat(rToGo,xToGo);
4584
4585   Int_t index=2*fCurrentEsdTrack+shieldindex;
4586
4587   Double_t xOverX0,x0,lengthTimesMeanDensity;
4588   Bool_t anglecorr=kTRUE;
4589
4590   switch(mode) {
4591   case 0:
4592     xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4593     x0 = AliITSRecoParam::GetX0shield(shieldindex);
4594     lengthTimesMeanDensity = xOverX0*x0;
4595     break;
4596   case 1:
4597     if (!t->PropagateToTGeo(xToGo,1)) return 0;
4598     return 1;
4599     break;
4600   case 2:
4601     if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4602     xOverX0 = fxOverX0Shield[shieldindex];
4603     lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4604     break;
4605   case 3:
4606     if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4607     if(fxOverX0ShieldTrks[index]<0) {
4608       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4609       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4610                                  (1.-t->GetSnp()*t->GetSnp()));
4611       fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4612       fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4613       return 1;
4614     }
4615     xOverX0 = fxOverX0ShieldTrks[index];
4616     lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4617     break;
4618   }
4619
4620   lengthTimesMeanDensity *= dir;
4621
4622   if (!t->AliExternalTrackParam::PropagateTo(xToGo,GetBz())) return 0;
4623   if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4624
4625   return 1;
4626 }
4627 //------------------------------------------------------------------------
4628 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4629                                                Int_t layerindex,
4630                                                Double_t oldGlobXYZ[3],
4631                                                TString direction) {
4632   //-------------------------------------------------------------------
4633   // Propagate beyond layer and correct for material
4634   // (material budget in different ways according to fUseTGeo value)
4635   //-------------------------------------------------------------------
4636
4637   // Define budget mode:
4638   // 0: material from AliITSRecoParam (hard coded)
4639   // 1: material from TGeo (on the fly)
4640   // 2: material from lut
4641   // 3: material from TGeo (same for all hypotheses)
4642   Int_t mode;
4643   switch(fUseTGeo) {
4644   case 0:
4645     mode=0;
4646     break;
4647   case 1:
4648     mode=1;
4649     break;
4650   case 2:
4651     mode=2;
4652     break;
4653   case 3:
4654     if(fTrackingPhase.Contains("Clusters2Tracks"))
4655       { mode=3; } else { mode=1; }
4656     break;
4657   case 4:
4658     if(fTrackingPhase.Contains("Clusters2Tracks"))
4659       { mode=3; } else { mode=2; }
4660     break;
4661   default:
4662     mode=0;
4663     break;
4664   }
4665   if(fTrackingPhase.Contains("Default")) mode=0;
4666
4667   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4668
4669   Double_t r=fgLayers[layerindex].GetR();
4670   Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4671
4672   Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4673   Double_t xToGo; t->GetLocalXat(rToGo,xToGo);
4674
4675   Int_t index=6*fCurrentEsdTrack+layerindex;
4676
4677   // Bring the track beyond the material
4678   if (!t->AliExternalTrackParam::PropagateTo(xToGo,GetBz())) return 0;
4679   Double_t globXYZ[3];
4680   t->GetXYZ(globXYZ);
4681
4682   Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4683   Double_t mparam[7];
4684   Bool_t anglecorr=kTRUE;
4685
4686   switch(mode) {
4687   case 0:
4688     xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4689     lengthTimesMeanDensity = xOverX0*x0;
4690     break;
4691   case 1:
4692     AliTracker::MeanMaterialBudget(oldGlobXYZ,globXYZ,mparam);
4693     if(mparam[1]>900000) return 0;
4694     xOverX0=mparam[1];
4695     lengthTimesMeanDensity=mparam[0]*mparam[4];
4696     anglecorr=kFALSE;
4697     break;
4698   case 2:
4699     if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4700     xOverX0 = fxOverX0Layer[layerindex];
4701     lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4702     break;
4703   case 3:
4704     if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4705     if(fxOverX0LayerTrks[index]<0) {
4706       AliTracker::MeanMaterialBudget(oldGlobXYZ,globXYZ,mparam);
4707       if(mparam[1]>900000) return 0;
4708       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4709                                  (1.-t->GetSnp()*t->GetSnp()));
4710       xOverX0=mparam[1]/angle;
4711       lengthTimesMeanDensity=mparam[0]*mparam[4]/angle;
4712       fxOverX0LayerTrks[index] = TMath::Abs(xOverX0);
4713       fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity);
4714     }
4715     xOverX0 = fxOverX0LayerTrks[index];
4716     lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4717     break;
4718   }
4719
4720   lengthTimesMeanDensity *= dir;
4721
4722   if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4723
4724   return 1;
4725 }
4726 //------------------------------------------------------------------------
4727 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4728   //-----------------------------------------------------------------
4729   // Initialize LUT for storing material for each prolonged track
4730   //-----------------------------------------------------------------
4731   fxOverX0PipeTrks = new Float_t[ntracks];
4732   fxTimesRhoPipeTrks = new Float_t[ntracks];
4733   fxOverX0ShieldTrks = new Float_t[ntracks*2];
4734   fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4735   fxOverX0LayerTrks = new Float_t[ntracks*6];
4736   fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4737
4738   for(Int_t i=0; i<ntracks; i++) {
4739     fxOverX0PipeTrks[i] = -1.;
4740     fxTimesRhoPipeTrks[i] = -1.;
4741   }
4742   for(Int_t j=0; j<ntracks*2; j++) {
4743     fxOverX0ShieldTrks[j] = -1.;
4744     fxTimesRhoShieldTrks[j] = -1.;
4745   }
4746   for(Int_t k=0; k<ntracks*6; k++) {
4747     fxOverX0LayerTrks[k] = -1.;
4748     fxTimesRhoLayerTrks[k] = -1.;
4749   }
4750
4751   fNtracks = ntracks;
4752
4753   return;
4754 }
4755 //------------------------------------------------------------------------
4756 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4757   //-----------------------------------------------------------------
4758   // Delete LUT for storing material for each prolonged track
4759   //-----------------------------------------------------------------
4760   if(fxOverX0PipeTrks) {
4761     delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4762   }
4763   if(fxOverX0ShieldTrks) {
4764     delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4765   }
4766
4767   if(fxOverX0LayerTrks) {
4768     delete [] fxOverX0LayerTrks;  fxOverX0LayerTrks = 0;
4769   }
4770   if(fxTimesRhoPipeTrks) {
4771     delete [] fxTimesRhoPipeTrks;  fxTimesRhoPipeTrks = 0;
4772   }
4773   if(fxTimesRhoShieldTrks) {
4774     delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4775   }
4776   if(fxTimesRhoLayerTrks) {
4777     delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4778   }
4779   return;
4780 }
4781 //------------------------------------------------------------------------
4782 Int_t AliITStrackerMI::CheckSkipLayer(AliITStrackMI *track,
4783                                       Int_t ilayer,Int_t idet) const {
4784   //-----------------------------------------------------------------
4785   // This method is used to decide whether to allow a prolongation
4786   // without clusters, because we want to skip the layer.
4787   // In this case the return value is > 0:
4788   // return 1: the user requested to skip a layer
4789   // return 2: track outside z acceptance of SSD/SDD and will cross both SPD
4790   //-----------------------------------------------------------------
4791
4792   if (AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer)) return 1;
4793
4794   if (idet<0 && ilayer>1 && AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4795     // check if track will cross SPD outer layer
4796     Double_t phiAtSPD2,zAtSPD2;
4797     if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4798       if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4799     }
4800   }
4801
4802   return 0;
4803 }
4804 //------------------------------------------------------------------------
4805 Int_t AliITStrackerMI::CheckDeadZone(/*AliITStrackMI *track,*/
4806                                      Int_t ilayer,Int_t idet,
4807                                      Double_t zmin,Double_t zmax/*,Double_t ymin,Double_t ymax*/) const {
4808   //-----------------------------------------------------------------
4809   // This method is used to decide whether to allow a prolongation
4810   // without clusters, because there is a dead zone in the road.
4811   // In this case the return value is > 0:
4812   // return 1: dead zone at z=0,+-7cm in SPD
4813   // return 2: dead area from the OCDB // NOT YET IMPLEMENTED
4814   //-----------------------------------------------------------------
4815
4816   // check dead zones at z=0,+-7cm in the SPD
4817   if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4818     Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4819                           fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4820                           fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4821     Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4822                           fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4823                           fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4824     for (Int_t i=0; i<3; i++)
4825       if (zmin<zmaxdead[i] && zmax>zmindead[i]) return 1;
4826   }
4827
4828   // check dead zones from OCDB
4829   if (!AliITSReconstructor::GetRecoParam()->GetUseDeadZonesFromOCDB()) return 0;
4830
4831   if(idet<0) return 0;
4832
4833   // look in OCDB (only entire dead modules for the moment)
4834   if (ilayer==0 || ilayer==1) { // SPD
4835     AliCDBEntry* cdbEntry = AliCDBManager::Instance()->Get("ITS/Calib/SPDDead");
4836     if (!cdbEntry) {
4837       Error("CheckDeadZone","Cannot get CDB entry for SPD\n");
4838       return 0;
4839     }
4840     TObjArray* spdEntry = (TObjArray*)cdbEntry->GetObject();
4841     if (!spdEntry) {
4842       Error("CheckDeadZone","Cannot get CDB entry for SPD\n");
4843       return 0;
4844     }
4845     if(ilayer==1) idet += AliITSgeomTGeo::GetNLadders(1)*AliITSgeomTGeo::GetNDetectors(1);
4846     //printf("SPD det: %d\n",idet);
4847     AliITSCalibrationSPD *calibSPD = (AliITSCalibrationSPD*)spdEntry->At(idet);
4848     if (calibSPD->IsBad()) return 2;
4849   } else if (ilayer==2 || ilayer==3) { // SDD
4850     AliCDBEntry* cdbEntry = AliCDBManager::Instance()->Get("ITS/Calib/CalibSDD");
4851     if (!cdbEntry) {
4852       Error("CheckDeadZone","Cannot get CDB entry for SDD\n");
4853       return 0;
4854     }
4855     TObjArray* sddEntry = (TObjArray*)cdbEntry->GetObject();
4856     if (!sddEntry) {
4857       Error("CheckDeadZone","Cannot get CDB entry for SDD\n");
4858       return 0;
4859     }
4860     if(ilayer==3) idet += AliITSgeomTGeo::GetNLadders(3)*AliITSgeomTGeo::GetNDetectors(3);
4861     //printf("SDD det: %d\n",idet);
4862     AliITSCalibrationSDD *calibSDD = (AliITSCalibrationSDD*)sddEntry->At(idet);
4863     if (calibSDD->IsDead()) return 2;
4864   } else if (ilayer==4 || ilayer==5) { // SSD
4865   } else {
4866     Error("CheckDeadZone","Wrong layer number\n");
4867     if(ilayer==5) idet += AliITSgeomTGeo::GetNLadders(5)*AliITSgeomTGeo::GetNDetectors(5);
4868     return 0;
4869   }
4870
4871   return 0;
4872 }
4873 //------------------------------------------------------------------------
4874 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4875                                        AliITStrackMI *track,
4876                                        Float_t &xloc,Float_t &zloc) const {
4877   //-----------------------------------------------------------------
4878   // Gives position of track in local module ref. frame
4879   //-----------------------------------------------------------------
4880
4881   xloc=0.;
4882   zloc=0.;
4883
4884   if(idet<0) return kFALSE;
4885
4886   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4887
4888   Int_t lad = Int_t(idet/ndet) + 1;
4889
4890   Int_t det = idet - (lad-1)*ndet + 1;
4891
4892   Double_t xyzGlob[3],xyzLoc[3];
4893
4894   track->GetXYZ(xyzGlob);
4895
4896   AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
4897
4898   xloc = (Float_t)xyzLoc[0];
4899   zloc = (Float_t)xyzLoc[2];
4900
4901   return kTRUE;
4902 }
4903 //------------------------------------------------------------------------