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