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