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