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