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<trackindex;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   if (esdindex>=fTrackHypothesys.GetEntriesFast())
3404     fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3405   //
3406   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3407   if (!array) {
3408     array = new TObjArray(10);
3409     fTrackHypothesys.AddAt(array,esdindex);
3410   }
3411   array->AddLast(track);
3412 }
3413 //------------------------------------------------------------------------
3414 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3415 {
3416   //-------------------------------------------------------------------
3417   // compress array of track hypothesys
3418   // keep only maxsize best hypothesys
3419   //-------------------------------------------------------------------
3420   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3421   if (! (fTrackHypothesys.At(esdindex)) ) return;
3422   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3423   Int_t entries = array->GetEntriesFast();
3424   //
3425   //- find preliminary besttrack as a reference
3426   Float_t minchi2=10000;
3427   Int_t maxn=0;
3428   AliITStrackMI * besttrack=0;
3429   //
3430   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3431     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3432     if (!track) continue;
3433     Float_t chi2 = NormalizedChi2(track,0);
3434     //
3435     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3436     track->SetLabel(tpcLabel);
3437     CookdEdx(track);
3438     track->SetFakeRatio(1.);
3439     CookLabel(track,0.); //For comparison only
3440     //
3441     //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3442     if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3443       if (track->GetNumberOfClusters()<maxn) continue;
3444       maxn = track->GetNumberOfClusters();
3445       //      if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2 *= track->GetChi2MIP(3); // RS
3446       if (chi2<minchi2){
3447         minchi2=chi2;
3448         besttrack=track;
3449       }
3450     }
3451     else{
3452       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3453         delete array->RemoveAt(itrack);
3454       }
3455     }
3456   }
3457   if (!besttrack) return;
3458   //
3459   //
3460   //take errors of best track as a reference
3461   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3462   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3463   for (Int_t j=0;j<6;j++) {
3464     if (besttrack->GetClIndex(j)>=0){
3465       erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3466       errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3467       ny[j]   = besttrack->GetNy(j);
3468       nz[j]   = besttrack->GetNz(j);
3469     }
3470   }
3471   //
3472   // calculate normalized chi2
3473   //
3474   Float_t * chi2        = new Float_t[entries];
3475   Int_t * index         = new Int_t[entries];
3476   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3477   for (Int_t itrack=0;itrack<entries;itrack++){
3478     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3479     if (track){
3480       AliDebug(2,Form("track %d  ncls %d\n",itrack,track->GetNumberOfClusters()));
3481       double chi2t = GetNormalizedChi2(track, mode);
3482       track->SetChi2MIP(0,chi2t);
3483       if (chi2t<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3484         if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3485         chi2[itrack] = chi2t;
3486       }
3487       else{
3488         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3489           delete array->RemoveAt(itrack);
3490         }
3491       }
3492     }
3493   }
3494   //
3495   TMath::Sort(entries,chi2,index,kFALSE);
3496   besttrack = (AliITStrackMI*)array->At(index[0]);
3497   if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3498   if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3499     for (Int_t j=0;j<6;j++){
3500       if (besttrack->GetClIndex(j)>=0){
3501         erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3502         errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3503         ny[j]   = besttrack->GetNy(j);
3504         nz[j]   = besttrack->GetNz(j);
3505       }
3506     }
3507   }
3508   //
3509   // calculate one more time with updated normalized errors
3510   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3511   for (Int_t itrack=0;itrack<entries;itrack++){
3512     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3513     if (track){
3514       double chi2t = GetNormalizedChi2(track, mode);
3515       track->SetChi2MIP(0,chi2t);
3516       AliDebug(2,Form("track %d  ncls %d\n",itrack,track->GetNumberOfClusters()));
3517       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3518         if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3519         chi2[itrack] = chi2t;  //-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3520       }
3521       else {
3522         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3523           delete array->RemoveAt(itrack);
3524         }
3525       }
3526     }
3527   }
3528   entries = array->GetEntriesFast();
3529   //
3530   //
3531   if (entries>0){
3532     TObjArray * newarray = new TObjArray();
3533     TMath::Sort(entries,chi2,index,kFALSE);
3534     besttrack = (AliITStrackMI*)array->At(index[0]);
3535     if (besttrack){
3536       AliDebug(2,Form("ncls best track %d     %f   %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3537       //
3538       for (Int_t j=0;j<6;j++){
3539         if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3540           erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3541           errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3542           ny[j]   = besttrack->GetNy(j);
3543           nz[j]   = besttrack->GetNz(j);
3544         }
3545       }
3546       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3547       minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3548       Float_t minn = besttrack->GetNumberOfClusters()-3;
3549       Int_t accepted=0;
3550       for (Int_t i=0;i<entries;i++){
3551         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3552         if (!track) continue;
3553         if (accepted>maxcut) break;
3554         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3555         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3556           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3557             delete array->RemoveAt(index[i]);
3558             continue;
3559           }
3560         }
3561         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3562         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3563           if (!shortbest) accepted++;
3564           //
3565           newarray->AddLast(array->RemoveAt(index[i]));
3566           for (Int_t j=0;j<6;j++){
3567             if (nz[j]==0){
3568               erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3569               errz[j] = track->GetSigmaZ(j); errz[j]   = track->GetSigmaZ(j+6);
3570               ny[j]   = track->GetNy(j);
3571               nz[j]   = track->GetNz(j);
3572             }
3573           }
3574         }
3575         else{
3576           delete array->RemoveAt(index[i]);
3577         }
3578       }
3579       array->Delete();
3580       delete fTrackHypothesys.RemoveAt(esdindex);
3581       fTrackHypothesys.AddAt(newarray,esdindex);
3582     }
3583     else{
3584       array->Delete();
3585       delete fTrackHypothesys.RemoveAt(esdindex);
3586     }
3587   }
3588   delete [] chi2;
3589   delete [] index;
3590 }
3591 //------------------------------------------------------------------------
3592 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3593 {
3594   //-------------------------------------------------------------
3595   // try to find best hypothesy
3596   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3597   // RS: optionally changing this to product of Chi2MIP(0)*Chi2MIP(3) == (chi2*chi2_interpolated)
3598   //-------------------------------------------------------------
3599   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3600   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3601   if (!array) return 0;
3602   Int_t entries = array->GetEntriesFast();
3603   if (!entries) return 0;
3604   Float_t minchi2 = 100000;
3605   AliITStrackMI * besttrack=0;
3606   //
3607   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
3608   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3609   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3610   Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3611   //
3612   for (Int_t i=0;i<entries;i++){
3613     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3614     if (!track) continue;
3615     Float_t sigmarfi,sigmaz;
3616     GetDCASigma(track,sigmarfi,sigmaz);
3617     track->SetDnorm(0,sigmarfi);
3618     track->SetDnorm(1,sigmaz);
3619     //
3620     track->SetChi2MIP(1,1000000);
3621     track->SetChi2MIP(2,1000000);
3622     track->SetChi2MIP(3,1000000);
3623     //
3624     // backtrack
3625     backtrack = new(backtrack) AliITStrackMI(*track);
3626     if (track->GetConstrain()) {
3627       if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3628       if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3629         if (fUseImproveKalman) {if (!backtrack->ImproveKalman(xyzVtx,ersVtx,0,0,0)) continue;}
3630         else                   {if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;}
3631       }
3632       backtrack->ResetCovariance(10.);
3633     }else{
3634       backtrack->ResetCovariance(10.);
3635     }
3636     backtrack->ResetClusters();
3637
3638     Double_t x = original->GetX();
3639     if (!RefitAt(x,backtrack,track)) continue;
3640     //
3641     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3642     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3643     if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.)  continue;
3644     track->SetChi22(GetMatchingChi2(backtrack,original));
3645
3646     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
3647     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
3648     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
3649
3650
3651     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1))  continue;
3652     //
3653     //forward track - without constraint
3654     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3655     forwardtrack->ResetClusters();
3656     x = track->GetX();
3657     if (!RefitAt(x,forwardtrack,track) && fSelectBestMIP03) continue;  // w/o fwd track MIP03 is meaningless
3658     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3659     if  (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0)  continue;
3660     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2))  continue;
3661
3662     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3663     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3664     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
3665     forwardtrack->SetD(0,track->GetD(0));
3666     forwardtrack->SetD(1,track->GetD(1));
3667     {
3668       Int_t list[6];
3669       AliITSRecPoint* clist[6];
3670       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3671       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3672     }
3673
3674     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3675     if  ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3676     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3677       track->SetChi2MIP(3,1000);
3678       continue;
3679     }
3680     Double_t chi2 = track->GetChi2MIP(0); // +track->GetNUsed();    //RS
3681     if (fSelectBestMIP03) chi2 *= track->GetChi2MIP(3);
3682     else chi2 += track->GetNUsed();
3683     //
3684     for (Int_t ichi=0;ichi<5;ichi++){
3685       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3686     }
3687     if (chi2 < minchi2){
3688       //besttrack = new AliITStrackMI(*forwardtrack);
3689       besttrack = track;
3690       besttrack->SetLabel(track->GetLabel());
3691       besttrack->SetFakeRatio(track->GetFakeRatio());
3692       minchi2   = chi2;
3693       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3694       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3695       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
3696     }
3697   }
3698   delete backtrack;
3699   delete forwardtrack;
3700
3701   if (!besttrack)  return 0;
3702
3703   Int_t accepted=0;
3704   for (Int_t i=0;i<entries;i++){
3705     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3706
3707     if (!track) continue;
3708
3709     if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3710         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)
3711         // RS: don't apply this cut when fSelectBestMIP03 is on
3712         || (!fSelectBestMIP03 && (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.))
3713         ){
3714       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3715         delete array->RemoveAt(i);
3716         continue;
3717       }
3718     }
3719     else{
3720       accepted++;
3721     }
3722   }
3723   //
3724   array->Compress();
3725   SortTrackHypothesys(esdindex,checkmax,1);
3726
3727   array = (TObjArray*) fTrackHypothesys.At(esdindex);
3728   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3729   besttrack = (AliITStrackMI*)array->At(0);
3730   if (!besttrack)  return 0;
3731   besttrack->SetChi2MIP(8,0);
3732   fBestTrackIndex[esdindex]=0;
3733   entries = array->GetEntriesFast();
3734   AliITStrackMI *longtrack =0;
3735   minchi2 =1000;
3736   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3737   for (Int_t itrack=entries-1;itrack>0;itrack--) {
3738     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3739     if (!track->GetConstrain()) continue;
3740     if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3741     if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3742     if (track->GetChi2MIP(0)>4.) continue;
3743     minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3744     longtrack =track;
3745   }
3746   //if (longtrack) besttrack=longtrack;
3747   //
3748   // RS do shared cluster analysis here only if the new sharing analysis is not requested
3749   //RRR if (fFlagFakes) return besttrack;
3750
3751   Int_t list[6];
3752   AliITSRecPoint * clist[6];
3753   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3754   if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3755       &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3756     RegisterClusterTracks(besttrack,esdindex);
3757   }
3758   //
3759   //
3760   if (shared>0.0){
3761     Int_t nshared;
3762     Int_t overlist[6];
3763     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3764     if (sharedtrack>=0){
3765       //
3766       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5,original);
3767       if (besttrack){
3768         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3769       }
3770       else return 0;
3771     }
3772   }
3773
3774   if (shared>2.5) return 0;
3775   if (shared>1.0) return besttrack;
3776   //
3777   // Don't sign clusters if not gold track
3778   //
3779   if (!besttrack->IsGoldPrimary()) return besttrack;
3780   if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
3781   //
3782   if (fConstraint[fPass]){
3783     //
3784     // sign clusters
3785     //
3786     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3787     for (Int_t i=0;i<6;i++){
3788       Int_t index = besttrack->GetClIndex(i);
3789       if (index<0) continue;
3790       Int_t ilayer =  (index & 0xf0000000) >> 28;
3791       if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3792       AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3793       if (!c) continue;
3794       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3795       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3796       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
3797       if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3798         if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3799       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
3800
3801       Bool_t cansign = kTRUE;
3802       for (Int_t itrack=0;itrack<entries; itrack++){
3803         AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3804         if (!track) continue;
3805         if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3806         if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3807           cansign = kFALSE;
3808           break;
3809         }
3810       }
3811       if (cansign){
3812         if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3813         if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3814         if (!c->IsUsed()) c->Use();
3815       }
3816     }
3817   }
3818   return besttrack;
3819 }
3820 //------------------------------------------------------------------------
3821 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3822 {
3823   //
3824   // get "best" hypothesys
3825   //
3826
3827   Int_t nentries = itsTracks.GetEntriesFast();
3828   for (Int_t i=0;i<nentries;i++){
3829     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3830     if (!track) continue;
3831     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3832     if (!array) continue;
3833     if (array->GetEntriesFast()<=0) continue;
3834     //
3835     AliITStrackMI* longtrack=0;
3836     Float_t minn=0;
3837     Float_t maxchi2=1000;
3838     for (Int_t j=0;j<array->GetEntriesFast();j++){
3839       AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3840       if (!trackHyp) continue;
3841       if (trackHyp->GetGoldV0()) {
3842         longtrack = trackHyp;   //gold V0 track taken
3843         break;
3844       }
3845       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3846       Float_t chi2 = trackHyp->GetChi2MIP(0);
3847       if (fSelectBestMIP03) chi2 *= trackHyp->GetChi2MIP(3);
3848       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = chi2; //trackHyp->GetChi2MIP(0);
3849       //
3850       if (fAfterV0){ // ??? RS
3851         if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3852       }
3853       if (chi2 > maxchi2) continue;
3854       minn = trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3855       if (fSelectBestMIP03) minn++; // allow next to longest to win
3856       maxchi2 = chi2;
3857       longtrack=trackHyp;
3858     }
3859     //
3860     //
3861     //
3862     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3863     if (!longtrack) {longtrack = besttrack;}
3864     else besttrack= longtrack;
3865     //
3866     if (besttrack) {
3867       Int_t list[6];
3868       AliITSRecPoint * clist[6];
3869       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3870       //
3871       track->SetNUsed(shared);
3872       track->SetNSkipped(besttrack->GetNSkipped());
3873       track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3874       if (shared>0) {
3875         if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3876         Int_t nshared;
3877         Int_t overlist[6];
3878         //
3879         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3880         //if (sharedtrack==-1) sharedtrack=0;
3881         if (sharedtrack>=0) {
3882           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5,track);
3883         }
3884       }
3885       if (besttrack&&fAfterV0) {
3886         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3887         track->SetWinner(besttrack);
3888       }
3889       if (besttrack) {
3890         if (fConstraint[fPass]) {
3891           UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3892           track->SetWinner(besttrack);
3893         }
3894         if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3895           if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3896                TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3897         }
3898       }
3899     }
3900   }
3901 }
3902
3903 //------------------------------------------------------------------------
3904 void AliITStrackerMI::FlagFakes(const TObjArray &itsTracks)
3905 {
3906   //
3907   // RS: flag those tracks which are suxpected to have fake clusters
3908   //
3909   const double kThreshPt = 0.5;
3910   AliRefArray *refArr[6];
3911   //
3912   for (int i=0;i<6;i++) {
3913     int ncl = fgLayers[i].GetNumberOfClusters();
3914     refArr[i] = new AliRefArray(ncl,TMath::Min(ncl,1000));
3915   }
3916   Int_t nentries = itsTracks.GetEntriesFast();
3917   //
3918   // fill cluster->track associations
3919   for (Int_t itr=0;itr<nentries;itr++){
3920     AliITStrackMI* track = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3921     if (!track) continue;
3922     AliITStrackMI* trackITS = track->GetWinner();
3923     if (!trackITS) continue;
3924     for (int il=trackITS->GetNumberOfClusters();il--;) {
3925       int idx = trackITS->GetClusterIndex(il);
3926       Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3927       //      if (c>fgLayers[l].GetNumberOfClusters()) continue;
3928       refArr[l]->AddReference(c, itr);
3929     }
3930   }
3931   //
3932   const UInt_t kMaxRef = 100;
3933   UInt_t crefs[kMaxRef];
3934   Int_t ncrefs=0;
3935   // process tracks with shared clusters
3936   for (int itr=0;itr<nentries;itr++){
3937     AliITStrackMI* track0 = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3938     AliITStrackMI* trackH0 = track0->GetWinner();
3939     if (!trackH0) continue;
3940     AliESDtrack* esd0 = track0->GetESDtrack();
3941     //
3942     for (int il=0;il<trackH0->GetNumberOfClusters();il++) {
3943       int idx = trackH0->GetClusterIndex(il);
3944       Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3945       ncrefs = refArr[l]->GetReferences(c,crefs,kMaxRef);
3946       if (ncrefs<2) continue; // there will be always self-reference, for sharing needs at least 2
3947       esd0->SetITSSharedFlag(l);
3948       for (int ir=ncrefs;ir--;) {
3949         if (int(crefs[ir]) <= itr) continue; // ==:selfreference, <: the same pair will be checked with >
3950         AliITStrackMI* track1 = (AliITStrackMI*)itsTracks.UncheckedAt(crefs[ir]);
3951         AliITStrackMI* trackH1 = track1->GetWinner();
3952         AliESDtrack* esd1 = track1->GetESDtrack();
3953         esd1->SetITSSharedFlag(l);
3954         //
3955         double pt0 = trackH0->Pt(), pt1 = trackH1->Pt(), res = 0.;
3956         if      (pt0>kThreshPt && pt0-pt1>0.2+0.2*(pt0-kThreshPt) ) res = -100;
3957         else if (pt1>kThreshPt && pt1-pt0>0.2+0.2*(pt1-kThreshPt) ) res = 100;
3958
3959         // select the one with smallest chi2's product
3960         res += trackH0->GetChi2MIP(0)*trackH0->GetChi2MIP(3);
3961         res -= trackH1->GetChi2MIP(0)*trackH1->GetChi2MIP(3);
3962         //
3963         if (res<0) esd1->SetITSFakeFlag();  // esd0 is winner
3964         else       esd0->SetITSFakeFlag();  // esd1 is winner
3965       }
3966       //
3967     }
3968     //
3969   }
3970   //
3971   for (int i=6;i--;) delete refArr[i];
3972 }
3973
3974 //------------------------------------------------------------------------
3975 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3976   //--------------------------------------------------------------------
3977   //This function "cooks" a track label. If label<0, this track is fake.
3978   //--------------------------------------------------------------------
3979   Int_t tpcLabel=-1;
3980
3981   if (track->GetESDtrack()){
3982     tpcLabel = track->GetESDtrack()->GetTPCLabel();
3983     ULong_t trStatus=track->GetESDtrack()->GetStatus();
3984     if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
3985   }
3986    track->SetChi2MIP(9,0);
3987    Int_t nwrong=0;
3988    for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3989      Int_t cindex = track->GetClusterIndex(i);
3990      Int_t l=(cindex & 0xf0000000) >> 28;
3991      AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3992      Int_t isWrong=1;
3993      for (Int_t ind=0;ind<3;ind++){
3994        if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
3995        //AliDebug(2,Form("icl %d  ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3996      }
3997      track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3998      nwrong+=isWrong;
3999    }
4000    Int_t nclusters = track->GetNumberOfClusters();
4001    if (nclusters > 0) //PH Some tracks don't have any cluster
4002      track->SetFakeRatio(double(nwrong)/double(nclusters));
4003    if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
4004      track->SetLabel(-tpcLabel);
4005    } else {
4006      track->SetLabel(tpcLabel);
4007    }
4008    AliDebug(2,Form(" nls %d wrong %d  label %d  tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
4009
4010 }
4011 //------------------------------------------------------------------------
4012 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
4013   //
4014   // Fill the dE/dx in this track
4015   //
4016   track->SetChi2MIP(9,0);
4017   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
4018     Int_t cindex = track->GetClusterIndex(i);
4019     Int_t l=(cindex & 0xf0000000) >> 28;
4020     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4021     Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
4022     Int_t isWrong=1;
4023     for (Int_t ind=0;ind<3;ind++){
4024       if (cl->GetLabel(ind)==lab) isWrong=0;
4025     }
4026     track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4027   }
4028   Double_t low=0.;
4029   Double_t up=0.51;
4030   track->CookdEdx(low,up);
4031 }
4032 //------------------------------------------------------------------------
4033 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
4034   //
4035   // Create some arrays
4036   //
4037   if (fCoefficients) delete []fCoefficients;
4038   fCoefficients = new Float_t[ntracks*48];
4039   for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
4040 }
4041 //------------------------------------------------------------------------
4042 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4043 {
4044   //
4045   // Compute predicted chi2
4046   //
4047   // Take into account the mis-alignment (bring track to cluster plane)
4048   Double_t xTrOrig=track->GetX();
4049   if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
4050   Float_t erry,errz,covyz;
4051   Float_t theta = track->GetTgl();
4052   Float_t phi   = track->GetSnp();
4053   phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
4054   AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
4055   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()));
4056   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()));
4057   Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
4058   // Bring the track back to detector plane in ideal geometry
4059   // [mis-alignment will be accounted for in UpdateMI()]
4060   if (!track->Propagate(xTrOrig)) return 1000.;
4061   Float_t ny,nz;
4062   AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
4063   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
4064   if (delta>1){
4065     chi2+=0.5*TMath::Min(delta/2,2.);
4066     chi2+=2.*cluster->GetDeltaProbability();
4067   }
4068   //
4069   track->SetNy(layer,ny);
4070   track->SetNz(layer,nz);
4071   track->SetSigmaY(layer,erry);
4072   track->SetSigmaZ(layer, errz);
4073   track->SetSigmaYZ(layer,covyz);
4074   //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
4075   track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
4076   return chi2;
4077
4078 }
4079 //------------------------------------------------------------------------
4080 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
4081 {
4082   //
4083   // Update ITS track
4084   //
4085   Int_t layer = (index & 0xf0000000) >> 28;
4086   track->SetClIndex(layer, index);
4087   if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
4088     if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
4089       chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
4090       track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
4091     }
4092   }
4093
4094   if (TMath::Abs(cl->GetQ())<1.e-13) return 0;  // ingore the "virtual" clusters
4095
4096
4097   // Take into account the mis-alignment (bring track to cluster plane)
4098   Double_t xTrOrig=track->GetX();
4099   Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
4100   AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
4101   AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
4102   AliDebug(2,Form(" xtr %f  xcl %f",track->GetX(),cl->GetX()));
4103
4104   if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
4105
4106   AliCluster c(*cl);
4107   c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
4108   c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
4109   c.SetSigmaYZ(track->GetSigmaYZ(layer));
4110
4111
4112   Int_t updated = track->UpdateMI(&c,chi2,index);
4113   // Bring the track back to detector plane in ideal geometry
4114   if (!track->Propagate(xTrOrig)) return 0;
4115
4116   if(!updated) AliDebug(2,"update failed");
4117   return updated;
4118 }
4119
4120 //------------------------------------------------------------------------
4121 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
4122 {
4123   //
4124   //DCA sigmas parameterization
4125   //to be paramterized using external parameters in future
4126   //
4127   //
4128   Double_t curv=track->GetC();
4129   sigmarfi = 0.0040+1.4 *TMath::Abs(curv)+332.*curv*curv;
4130   sigmaz   = 0.0110+4.37*TMath::Abs(curv);
4131 }
4132 //------------------------------------------------------------------------
4133 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
4134 {
4135   //
4136   // Clusters from delta electrons?
4137   //
4138   Int_t entries = clusterArray->GetEntriesFast();
4139   if (entries<4) return;
4140   AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
4141   Int_t layer = cluster->GetLayer();
4142   if (layer>1) return;
4143   Int_t index[10000];
4144   Int_t ncandidates=0;
4145   Float_t r = (layer>0)? 7:4;
4146   //
4147   for (Int_t i=0;i<entries;i++){
4148     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
4149     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
4150     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
4151     index[ncandidates] = i;  //candidate to belong to delta electron track
4152     ncandidates++;
4153     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
4154       cl0->SetDeltaProbability(1);
4155     }
4156   }
4157   //
4158   //
4159   //
4160   for (Int_t i=0;i<ncandidates;i++){
4161     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
4162     if (cl0->GetDeltaProbability()>0.8) continue;
4163     //
4164     Int_t ncl = 0;
4165     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
4166     sumy=sumz=sumy2=sumyz=sumw=0.0;
4167     for (Int_t j=0;j<ncandidates;j++){
4168       if (i==j) continue;
4169       AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
4170       //
4171       Float_t dz = cl0->GetZ()-cl1->GetZ();
4172       Float_t dy = cl0->GetY()-cl1->GetY();
4173       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
4174         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
4175         y[ncl] = cl1->GetY();
4176         z[ncl] = cl1->GetZ();
4177         sumy+= y[ncl]*weight;
4178         sumz+= z[ncl]*weight;
4179         sumy2+=y[ncl]*y[ncl]*weight;
4180         sumyz+=y[ncl]*z[ncl]*weight;
4181         sumw+=weight;
4182         ncl++;
4183       }
4184     }
4185     if (ncl<4) continue;
4186     Float_t det = sumw*sumy2  - sumy*sumy;
4187     Float_t delta=1000;
4188     if (TMath::Abs(det)>0.01){
4189       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
4190       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
4191       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
4192     }
4193     else{
4194       Float_t z0  = sumyz/sumy;
4195       delta = TMath::Abs(cl0->GetZ()-z0);
4196     }
4197     if ( delta<0.05) {
4198       cl0->SetDeltaProbability(1-20.*delta);
4199     }
4200   }
4201 }
4202 //------------------------------------------------------------------------
4203 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
4204 {
4205   //
4206   // Update ESD track
4207   //
4208   track->UpdateESDtrack(flags);
4209   AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
4210   if (oldtrack) delete oldtrack;
4211   track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
4212   // if (TMath::Abs(track->GetDnorm(1))<0.000000001){
4213   //   printf("Problem\n");
4214   // }
4215 }
4216 //------------------------------------------------------------------------
4217 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
4218   //
4219   // Get nearest upper layer close to the point xr.
4220   // rough approximation
4221   //
4222   const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
4223   Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
4224   Int_t res =6;
4225   for (Int_t i=0;i<6;i++){
4226     if (radius<kRadiuses[i]){
4227       res =i;
4228       break;
4229     }
4230   }
4231   return res;
4232 }
4233 //------------------------------------------------------------------------
4234 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4235   //--------------------------------------------------------------------
4236   // Fill a look-up table with mean material
4237   //--------------------------------------------------------------------
4238
4239   Int_t n=1000;
4240   Double_t mparam[7];
4241   Double_t point1[3],point2[3];
4242   Double_t phi,cosphi,sinphi,z;
4243   // 0-5 layers, 6 pipe, 7-8 shields
4244   Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4245   Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4246
4247   Int_t ifirst=0,ilast=0;
4248   if(material.Contains("Pipe")) {
4249     ifirst=6; ilast=6;
4250   } else if(material.Contains("Shields")) {
4251     ifirst=7; ilast=8;
4252   } else if(material.Contains("Layers")) {
4253     ifirst=0; ilast=5;
4254   } else {
4255     Error("BuildMaterialLUT","Wrong layer name\n");
4256   }
4257
4258   for(Int_t imat=ifirst; imat<=ilast; imat++) {
4259     Double_t param[5]={0.,0.,0.,0.,0.};
4260     for (Int_t i=0; i<n; i++) {
4261       phi = 2.*TMath::Pi()*gRandom->Rndm();
4262       cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4263       z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4264       point1[0] = rmin[imat]*cosphi;
4265       point1[1] = rmin[imat]*sinphi;
4266       point1[2] = z;
4267       point2[0] = rmax[imat]*cosphi;
4268       point2[1] = rmax[imat]*sinphi;
4269       point2[2] = z;
4270       AliTracker::MeanMaterialBudget(point1,point2,mparam);
4271       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4272     }
4273     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4274     if(imat<=5) {
4275       fxOverX0Layer[imat] = param[1];
4276       fxTimesRhoLayer[imat] = param[0]*param[4];
4277     } else if(imat==6) {
4278       fxOverX0Pipe = param[1];
4279       fxTimesRhoPipe = param[0]*param[4];
4280     } else if(imat==7) {
4281       fxOverX0Shield[0] = param[1];
4282       fxTimesRhoShield[0] = param[0]*param[4];
4283     } else if(imat==8) {
4284       fxOverX0Shield[1] = param[1];
4285       fxTimesRhoShield[1] = param[0]*param[4];
4286     }
4287   }
4288   /*
4289   printf("%s\n",material.Data());
4290   printf("%f  %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4291   printf("%f  %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4292   printf("%f  %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4293   printf("%f  %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4294   printf("%f  %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4295   printf("%f  %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4296   printf("%f  %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4297   printf("%f  %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4298   printf("%f  %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4299   */
4300   return;
4301 }
4302 //------------------------------------------------------------------------
4303 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4304                                               TString direction) {
4305   //-------------------------------------------------------------------
4306   // Propagate beyond beam pipe and correct for material
4307   // (material budget in different ways according to fUseTGeo value)
4308   // Add time if going outward (PropagateTo or PropagateToTGeo)
4309   //-------------------------------------------------------------------
4310
4311   // Define budget mode:
4312   // 0: material from AliITSRecoParam (hard coded)
4313   // 1: material from TGeo in one step (on the fly)
4314   // 2: material from lut
4315   // 3: material from TGeo in one step (same for all hypotheses)
4316   Int_t mode;
4317   switch(fUseTGeo) {
4318   case 0:
4319     mode=0;
4320     break;
4321   case 1:
4322     mode=1;
4323     break;
4324   case 2:
4325     mode=2;
4326     break;
4327   case 3:
4328     if(fTrackingPhase.Contains("Clusters2Tracks"))
4329       { mode=3; } else { mode=1; }
4330     break;
4331   case 4:
4332     if(fTrackingPhase.Contains("Clusters2Tracks"))
4333       { mode=3; } else { mode=2; }
4334     break;
4335   default:
4336     mode=0;
4337     break;
4338   }
4339   if(fTrackingPhase.Contains("Default")) mode=0;
4340
4341   Int_t index=fCurrentEsdTrack;
4342
4343   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4344   Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4345   Double_t xToGo;
4346   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4347
4348   Double_t xOverX0,x0,lengthTimesMeanDensity;
4349
4350   switch(mode) {
4351   case 0:
4352     xOverX0 = AliITSRecoParam::GetdPipe();
4353     x0 = AliITSRecoParam::GetX0Be();
4354     lengthTimesMeanDensity = xOverX0*x0;
4355     lengthTimesMeanDensity *= dir;
4356     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4357     break;
4358   case 1:
4359     if (!t->PropagateToTGeo(xToGo,1)) return 0;
4360     break;
4361   case 2:
4362     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4363     xOverX0 = fxOverX0Pipe;
4364     lengthTimesMeanDensity = fxTimesRhoPipe;
4365     lengthTimesMeanDensity *= dir;
4366     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4367     break;
4368   case 3:
4369     if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4370     if(fxOverX0PipeTrks[index]<0) {
4371       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4372       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4373                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4374       fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4375       fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4376       return 1;
4377     }
4378     xOverX0 = fxOverX0PipeTrks[index];
4379     lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4380     lengthTimesMeanDensity *= dir;
4381     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4382     break;
4383   }
4384
4385   return 1;
4386 }
4387 //------------------------------------------------------------------------
4388 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4389                                                 TString shield,
4390                                                 TString direction) {
4391   //-------------------------------------------------------------------
4392   // Propagate beyond SPD or SDD shield and correct for material
4393   // (material budget in different ways according to fUseTGeo value)
4394   // Add time if going outward (PropagateTo or PropagateToTGeo)
4395   //-------------------------------------------------------------------
4396
4397   // Define budget mode:
4398   // 0: material from AliITSRecoParam (hard coded)
4399   // 1: material from TGeo in steps of X cm (on the fly)
4400   //    X = AliITSRecoParam::GetStepSizeTGeo()
4401   // 2: material from lut
4402   // 3: material from TGeo in one step (same for all hypotheses)
4403   Int_t mode;
4404   switch(fUseTGeo) {
4405   case 0:
4406     mode=0;
4407     break;
4408   case 1:
4409     mode=1;
4410     break;
4411   case 2:
4412     mode=2;
4413     break;
4414   case 3:
4415     if(fTrackingPhase.Contains("Clusters2Tracks"))
4416       { mode=3; } else { mode=1; }
4417     break;
4418   case 4:
4419     if(fTrackingPhase.Contains("Clusters2Tracks"))
4420       { mode=3; } else { mode=2; }
4421     break;
4422   default:
4423     mode=0;
4424     break;
4425   }
4426   if(fTrackingPhase.Contains("Default")) mode=0;
4427
4428   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4429   Double_t rToGo;
4430   Int_t    shieldindex=0;
4431   if (shield.Contains("SDD")) { // SDDouter
4432     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4433     shieldindex=1;
4434   } else if (shield.Contains("SPD")) {        // SPDouter
4435     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4436     shieldindex=0;
4437   } else {
4438     Error("CorrectForShieldMaterial"," Wrong shield name\n");
4439     return 0;
4440   }
4441
4442   // do nothing if we are already beyond the shield
4443   Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4444   if(dir<0 && rTrack > rToGo) return 1; // going outward
4445   if(dir>0 && rTrack < rToGo) return 1; // going inward
4446
4447
4448   Double_t xToGo;
4449   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4450
4451   Int_t index=2*fCurrentEsdTrack+shieldindex;
4452
4453   Double_t xOverX0,x0,lengthTimesMeanDensity;
4454   Int_t nsteps=1;
4455
4456   switch(mode) {
4457   case 0:
4458     xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4459     x0 = AliITSRecoParam::GetX0shield(shieldindex);
4460     lengthTimesMeanDensity = xOverX0*x0;
4461     lengthTimesMeanDensity *= dir;
4462     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4463     break;
4464   case 1:
4465     nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4466     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4467     break;
4468   case 2:
4469     if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4470     xOverX0 = fxOverX0Shield[shieldindex];
4471     lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4472     lengthTimesMeanDensity *= dir;
4473     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4474     break;
4475   case 3:
4476     if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4477     if(fxOverX0ShieldTrks[index]<0) {
4478       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4479       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4480                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4481       fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4482       fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4483       return 1;
4484     }
4485     xOverX0 = fxOverX0ShieldTrks[index];
4486     lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4487     lengthTimesMeanDensity *= dir;
4488     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4489     break;
4490   }
4491
4492   return 1;
4493 }
4494 //------------------------------------------------------------------------
4495 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4496                                                Int_t layerindex,
4497                                                Double_t oldGlobXYZ[3],
4498                                                TString direction) {
4499   //-------------------------------------------------------------------
4500   // Propagate beyond layer and correct for material
4501   // (material budget in different ways according to fUseTGeo value)
4502   // Add time if going outward (PropagateTo or PropagateToTGeo)
4503   //-------------------------------------------------------------------
4504
4505   // Define budget mode:
4506   // 0: material from AliITSRecoParam (hard coded)
4507   // 1: material from TGeo in stepsof X cm (on the fly)
4508   //    X = AliITSRecoParam::GetStepSizeTGeo()
4509   // 2: material from lut
4510   // 3: material from TGeo in one step (same for all hypotheses)
4511   Int_t mode;
4512   switch(fUseTGeo) {
4513   case 0:
4514     mode=0;
4515     break;
4516   case 1:
4517     mode=1;
4518     break;
4519   case 2:
4520     mode=2;
4521     break;
4522   case 3:
4523     if(fTrackingPhase.Contains("Clusters2Tracks"))
4524       { mode=3; } else { mode=1; }
4525     break;
4526   case 4:
4527     if(fTrackingPhase.Contains("Clusters2Tracks"))
4528       { mode=3; } else { mode=2; }
4529     break;
4530   default:
4531     mode=0;
4532     break;
4533   }
4534   if(fTrackingPhase.Contains("Default")) mode=0;
4535
4536   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4537
4538   Double_t r=fgLayers[layerindex].GetR();
4539   Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4540
4541   Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4542   Double_t xToGo;
4543   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4544
4545   Int_t index=6*fCurrentEsdTrack+layerindex;
4546
4547
4548   Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4549   Int_t nsteps=1;
4550
4551   // back before material (no correction)
4552   Double_t rOld,xOld;
4553   rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4554   if (!t->GetLocalXat(rOld,xOld)) return 0;
4555   if (!t->Propagate(xOld)) return 0;
4556
4557   switch(mode) {
4558   case 0:
4559     xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4560     lengthTimesMeanDensity = xOverX0*x0;
4561     lengthTimesMeanDensity *= dir;
4562     // Bring the track beyond the material
4563     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4564     break;
4565   case 1:
4566     nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4567     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4568     break;
4569   case 2:
4570     if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4571     xOverX0 = fxOverX0Layer[layerindex];
4572     lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4573     lengthTimesMeanDensity *= dir;
4574     // Bring the track beyond the material
4575     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4576     break;
4577   case 3:
4578     if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4579     if(fxOverX0LayerTrks[index]<0) {
4580       nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4581       if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4582       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4583                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4584       fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4585       fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4586       return 1;
4587     }
4588     xOverX0 = fxOverX0LayerTrks[index];
4589     lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4590     lengthTimesMeanDensity *= dir;
4591     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4592     break;
4593   }
4594
4595
4596   return 1;
4597 }
4598 //------------------------------------------------------------------------
4599 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4600   //-----------------------------------------------------------------
4601   // Initialize LUT for storing material for each prolonged track
4602   //-----------------------------------------------------------------
4603   fxOverX0PipeTrks = new Float_t[ntracks];
4604   fxTimesRhoPipeTrks = new Float_t[ntracks];
4605   fxOverX0ShieldTrks = new Float_t[ntracks*2];
4606   fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4607   fxOverX0LayerTrks = new Float_t[ntracks*6];
4608   fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4609
4610   for(Int_t i=0; i<ntracks; i++) {
4611     fxOverX0PipeTrks[i] = -1.;
4612     fxTimesRhoPipeTrks[i] = -1.;
4613   }
4614   for(Int_t j=0; j<ntracks*2; j++) {
4615     fxOverX0ShieldTrks[j] = -1.;
4616     fxTimesRhoShieldTrks[j] = -1.;
4617   }
4618   for(Int_t k=0; k<ntracks*6; k++) {
4619     fxOverX0LayerTrks[k] = -1.;
4620     fxTimesRhoLayerTrks[k] = -1.;
4621   }
4622
4623   fNtracks = ntracks;
4624
4625   return;
4626 }
4627 //------------------------------------------------------------------------
4628 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4629   //-----------------------------------------------------------------
4630   // Delete LUT for storing material for each prolonged track
4631   //-----------------------------------------------------------------
4632   if(fxOverX0PipeTrks) {
4633     delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4634   }
4635   if(fxOverX0ShieldTrks) {
4636     delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4637   }
4638
4639   if(fxOverX0LayerTrks) {
4640     delete [] fxOverX0LayerTrks;  fxOverX0LayerTrks = 0;
4641   }
4642   if(fxTimesRhoPipeTrks) {
4643     delete [] fxTimesRhoPipeTrks;  fxTimesRhoPipeTrks = 0;
4644   }
4645   if(fxTimesRhoShieldTrks) {
4646     delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4647   }
4648   if(fxTimesRhoLayerTrks) {
4649     delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4650   }
4651   return;
4652 }
4653 //------------------------------------------------------------------------
4654 void AliITStrackerMI::SetForceSkippingOfLayer() {
4655   //-----------------------------------------------------------------
4656   // Check if we are forced to skip layers
4657   // either we set to skip them in RecoParam
4658   // or they were off during data-taking
4659   //-----------------------------------------------------------------
4660
4661   const AliEventInfo *eventInfo = GetEventInfo();
4662
4663   for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4664     fForceSkippingOfLayer[l] = 0;
4665     // check reco param
4666     if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4667     // check run info
4668
4669     if(eventInfo &&
4670        AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4671       AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4672       if(l==0 || l==1)  {
4673         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4674       } else if(l==2 || l==3) {
4675         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4676       } else {
4677         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4678       }
4679     }
4680   }
4681   return;
4682 }
4683 //------------------------------------------------------------------------
4684 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4685                                       Int_t ilayer,Int_t idet) const {
4686   //-----------------------------------------------------------------
4687   // This method is used to decide whether to allow a prolongation
4688   // without clusters, because we want to skip the layer.
4689   // In this case the return value is > 0:
4690   // return 1: the user requested to skip a layer
4691   // return 2: track outside z acceptance
4692   //-----------------------------------------------------------------
4693
4694   if (ForceSkippingOfLayer(ilayer)) return 1;
4695
4696   Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4697
4698   if (idet<0 &&  // out in z
4699       ilayer>innerLayCanSkip &&
4700       AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4701     // check if track will cross SPD outer layer
4702     Double_t phiAtSPD2,zAtSPD2;
4703     if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4704       if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4705     }
4706     return 2; // always allow skipping, changed on 05.11.2009
4707   }
4708
4709   return 0;
4710 }
4711 //------------------------------------------------------------------------
4712 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4713                                      Int_t ilayer,Int_t idet,
4714                                      Double_t dz,Double_t dy,
4715                                      Bool_t noClusters) const {
4716   //-----------------------------------------------------------------
4717   // This method is used to decide whether to allow a prolongation
4718   // without clusters, because there is a dead zone in the road.
4719   // In this case the return value is > 0:
4720   // return 1: dead zone at z=0,+-7cm in SPD
4721   //     This method assumes that fSPDdetzcentre is ordered from -z to +z
4722   // return 2: all road is "bad" (dead or noisy) from the OCDB
4723   // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4724   // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4725   //-----------------------------------------------------------------
4726
4727   // check dead zones at z=0,+-7cm in the SPD
4728   if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4729     Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4730                           fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4731                           fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4732     Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4733                           fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4734                           fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4735     for (Int_t i=0; i<3; i++)
4736       if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4737         AliDebug(2,Form("crack SPD %d track z %f   %f   %f  %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4738         if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4739       }
4740   }
4741
4742   // check bad zones from OCDB
4743   if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4744
4745   if (idet<0) return 0;
4746
4747   AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4748
4749   Int_t detType=-1;
4750   Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4751   if (ilayer==0 || ilayer==1) {        // ----------  SPD
4752     detType = 0;
4753   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
4754     detType = 1;
4755     detSizeFactorX *= 2.;
4756   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
4757     detType = 2;
4758   }
4759   AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4760   if (detType==2) segm->SetLayer(ilayer+1);
4761   Float_t detSizeX = detSizeFactorX*segm->Dx();
4762   Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4763
4764   // check if the road overlaps with bad chips
4765   Float_t xloc,zloc;
4766   if(!(LocalModuleCoord(ilayer,idet,track,xloc,zloc)))return 0;
4767   Float_t zlocmin = zloc-dz;
4768   Float_t zlocmax = zloc+dz;
4769   Float_t xlocmin = xloc-dy;
4770   Float_t xlocmax = xloc+dy;
4771   Int_t chipsInRoad[100];
4772
4773   // check if road goes out of detector
4774   Bool_t touchNeighbourDet=kFALSE;
4775   if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4776   if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4777   if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4778   if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4779   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));
4780
4781   // check if this detector is bad
4782   if (det.IsBad()) {
4783     AliDebug(2,Form("lay %d  bad detector %d",ilayer,idet));
4784     if(!touchNeighbourDet) {
4785       return 2; // all detectors in road are bad
4786     } else {
4787       return 3; // at least one is bad
4788     }
4789   }
4790
4791   if(zlocmin>zlocmax)return 0;
4792   Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4793   AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4794   if (!nChipsInRoad) return 0;
4795
4796   Bool_t anyBad=kFALSE,anyGood=kFALSE;
4797   for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4798     if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4799     AliDebug(2,Form("  chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4800     if (det.IsChipBad(chipsInRoad[iCh])) {
4801       anyBad=kTRUE;
4802     } else {
4803       anyGood=kTRUE;
4804     }
4805   }
4806
4807   if (!anyGood) {
4808     if(!touchNeighbourDet) {
4809       AliDebug(2,"all bad in road");
4810       return 2;  // all chips in road are bad
4811     } else {
4812       return 3; // at least a bad chip in road
4813     }
4814   }
4815
4816   if (anyBad) {
4817     AliDebug(2,"at least a bad in road");
4818     return 3; // at least a bad chip in road
4819   }
4820
4821
4822   if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4823       || !noClusters) return 0;
4824
4825   // There are no clusters in road: check if there is at least
4826   // a bad SPD pixel or SDD anode or SSD strips on both sides
4827
4828   Int_t idetInITS=idet;
4829   for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4830
4831   if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4832     AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4833     return 4;
4834   }
4835   //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4836
4837   return 0;
4838 }
4839 //------------------------------------------------------------------------
4840 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4841                                        const AliITStrackMI *track,
4842                                        Float_t &xloc,Float_t &zloc) const {
4843   //-----------------------------------------------------------------
4844   // Gives position of track in local module ref. frame
4845   //-----------------------------------------------------------------
4846
4847   xloc=0.;
4848   zloc=0.;
4849
4850   if(idet<0) return kTRUE; // track out of z acceptance of layer
4851
4852   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4853
4854   Int_t lad = Int_t(idet/ndet) + 1;
4855
4856   Int_t det = idet - (lad-1)*ndet + 1;
4857
4858   Double_t xyzGlob[3],xyzLoc[3];
4859
4860   AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4861   // take into account the misalignment: xyz at real detector plane
4862   if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4863
4864   if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4865
4866   xloc = (Float_t)xyzLoc[0];
4867   zloc = (Float_t)xyzLoc[2];
4868
4869   return kTRUE;
4870 }
4871 //------------------------------------------------------------------------
4872 //------------------------------------------------------------------------
4873 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4874 //
4875 // Method to be optimized further:
4876 // Aim: decide whether a track can be used for PlaneEff evaluation
4877 //      the decision is taken based on the track quality at the layer under study
4878 //      no information on the clusters on this layer has to be used
4879 //      The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4880 //      the cut is done on number of sigmas from the boundaries
4881 //
4882 //  Input: Actual track, layer [0,5] under study
4883 //  Output: none
4884 //  Return: kTRUE if this is a good track
4885 //
4886 // it will apply a pre-selection to obtain good quality tracks.
4887 // Here also  you will have the possibility to put a control on the
4888 // impact point of the track on the basic block, in order to exclude border regions
4889 // this will be done by calling a proper method of the AliITSPlaneEff class.
4890 //
4891 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4892 // return: Bool_t   -> kTRUE if usable track, kFALSE if not usable.
4893 //
4894   Int_t index[AliITSgeomTGeo::kNLayers];
4895   Int_t k;
4896   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4897   //
4898   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4899     index[k]=clusters[k];
4900   }
4901
4902   if(!fPlaneEff)
4903     {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4904   AliITSlayer &layer=fgLayers[ilayer];
4905   Double_t r=layer.GetR();
4906   AliITStrackMI tmp(*track);
4907
4908 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4909   Int_t nclout=0; Int_t nclin=0;
4910   for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
4911  AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4912                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4913    // if (tmp.GetClIndex(lay)>=0) nclout++;
4914 if(index[lay]>=0)nclout++;
4915   }
4916   for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
4917     AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4918                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4919    if (index[lay]>=0) nclin++;
4920   }
4921   Int_t ncl=nclout+nclin;
4922   Bool_t nextout = kFALSE;
4923   if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4924   else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4925   Bool_t nextin = kFALSE;
4926   if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4927   else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4928   // maximum number of missing clusters allowed in outermost layers
4929   if(nclout<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff())
4930      return kFALSE;
4931   // maximum number of missing clusters allowed (both in innermost and in outermost layers)
4932   if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4933      return kFALSE;
4934   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout)  return kFALSE;
4935   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin)   return kFALSE;
4936   if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4937  //  if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff()  && !tmp.GetConstrain()) return kFALSE;
4938
4939 // detector number
4940   Double_t phi,z;
4941   if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4942   Int_t idet=layer.FindDetectorIndex(phi,z);
4943   if(idet<0) { AliInfo(Form("cannot find detector"));
4944     return kFALSE;}
4945
4946   // here check if it has good Chi Square.
4947
4948   //propagate to the intersection with the detector plane
4949   const AliITSdetector &det=layer.GetDetector(idet);
4950   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4951
4952   Float_t locx; //
4953   Float_t locz; //
4954   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4955   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4956   if(key>fPlaneEff->Nblock()) return kFALSE;
4957   Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4958   if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4959   //***************
4960   // DEFINITION OF SEARCH ROAD FOR accepting a track
4961   //
4962   Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
4963   Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
4964   Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
4965   Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
4966                                                     // done for RecPoints
4967
4968   // exclude tracks at boundary between detectors
4969   //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4970   Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4971   AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4972   AliDebug(2,Form("Local:    track impact x=%f, z=%f",locx,locz));
4973   AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4974   if ( (locx-dx < blockXmn+boundaryWidth) ||
4975        (locx+dx > blockXmx-boundaryWidth) ||
4976        (locz-dz < blockZmn+boundaryWidth) ||
4977        (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4978   return kTRUE;
4979 }
4980 //------------------------------------------------------------------------
4981 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4982 //
4983 // This Method has to be optimized! For the time-being it uses the same criteria
4984 // as those used in the search of extra clusters for overlapping modules.
4985 //
4986 // Method Purpose: estabilish whether a track has produced a recpoint or not
4987 //                 in the layer under study (For Plane efficiency)
4988 //
4989 // inputs: AliITStrackMI* track  (pointer to a usable track)
4990 // outputs: none
4991 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4992 //               traversed by this very track. In details:
4993 //               - if a cluster can be associated to the track then call
4994 //                  AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4995 //               - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4996 //
4997   if(!fPlaneEff)
4998     {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4999   AliITSlayer &layer=fgLayers[ilayer];
5000   Double_t r=layer.GetR();
5001   AliITStrackMI tmp(*track);
5002
5003 // detector number
5004   Double_t phi,z;
5005   if (!tmp.GetPhiZat(r,phi,z)) return;
5006   Int_t idet=layer.FindDetectorIndex(phi,z);
5007
5008   if(idet<0) { AliInfo(Form("cannot find detector"));
5009     return;}
5010
5011
5012 //propagate to the intersection with the detector plane
5013   const AliITSdetector &det=layer.GetDetector(idet);
5014   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5015
5016
5017 //***************
5018 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5019 //
5020   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5021                     TMath::Sqrt(tmp.GetSigmaZ2() +
5022                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5023                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5024                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5025   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5026                     TMath::Sqrt(tmp.GetSigmaY2() +
5027                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5028                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5029                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5030
5031 // road in global (rphi,z) [i.e. in tracking ref. system]
5032   Double_t zmin = tmp.GetZ() - dz;
5033   Double_t zmax = tmp.GetZ() + dz;
5034   Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5035   Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5036
5037 // select clusters in road
5038   layer.SelectClusters(zmin,zmax,ymin,ymax);
5039
5040 // Define criteria for track-cluster association
5041   Double_t msz = tmp.GetSigmaZ2() +
5042   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5043   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5044   AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5045   Double_t msy = tmp.GetSigmaY2() +
5046   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5047   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5048   AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5049   if (tmp.GetConstrain()) {
5050     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5051     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5052   }  else {
5053     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5054     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5055   }
5056   msz = 1./msz; // 1/RoadZ^2
5057   msy = 1./msy; // 1/RoadY^2
5058 //
5059
5060   const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5061   Int_t idetc=-1;
5062   Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5063   //Double_t  tolerance=0.2;
5064   /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5065     idetc = cl->GetDetectorIndex();
5066     if(idet!=idetc) continue;
5067     //Int_t ilay = cl->GetLayer();
5068
5069     if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5070     if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5071
5072     Double_t chi2=tmp.GetPredictedChi2(cl);
5073     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5074   }*/
5075   Float_t locx; //
5076   Float_t locz; //
5077   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5078 //
5079   AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5080   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5081   if(key>fPlaneEff->Nblock()) return;
5082   Bool_t found=kFALSE;
5083   //if (ci>=0) {
5084   Double_t chi2;
5085   while ((cl=layer.GetNextCluster(clidx))!=0) {
5086     idetc = cl->GetDetectorIndex();
5087     if(idet!=idetc) continue;
5088     // here real control to see whether the cluster can be associated to the track.
5089     // cluster not associated to track
5090     if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5091          (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy   > 1. ) continue;
5092     // calculate track-clusters chi2
5093     chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5094                                                // in particular, the error associated to the cluster
5095     //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5096     // chi2 cut
5097     if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5098     found=kTRUE;
5099     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5100    // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5101    // track->SetExtraModule(ilayer,idetExtra);
5102   }
5103   if(!fPlaneEff->UpDatePlaneEff(found,key))
5104        AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5105
5106 // this for FO efficiency studies (only for SPD) //
5107    UInt_t keyFO=999999;
5108    Bool_t foundFO=kFALSE;
5109    if(ilayer<2){ //ONLY SPD layers for FastOr studies
5110     TBits mapFO = fkDetTypeRec->GetFastOrFiredMap();
5111     Int_t phase = (fEsd->GetBunchCrossNumber())%4;
5112     if(!fSPDChipIntPlaneEff[key]){
5113       AliITSPlaneEffSPD spd;
5114       keyFO = spd.SwitchChipKeyNumbering(key);
5115       if(mapFO.TestBitNumber(keyFO))foundFO=kTRUE;
5116        keyFO = key + (AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip)*(phase+1);
5117        if(keyFO<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip) {
5118          AliWarning(Form("UseTrackForPlaneEff: too small keyF0 (= %d), setting it to 999999",keyFO));
5119          keyFO=999999;
5120        }
5121        if(!fPlaneEff->UpDatePlaneEff(foundFO,keyFO))
5122           AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for FastOR for key=%d",keyFO));
5123      }
5124   }
5125
5126
5127
5128   if(fPlaneEff->GetCreateHistos()&&  AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5129     Float_t tr[4]={99999.,99999.,9999.,9999.};    // initialize to high values
5130     Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5131     Int_t cltype[2]={-999,-999};
5132                                                           // and the module
5133
5134     Float_t angleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
5135
5136     tr[0]=locx;
5137     tr[1]=locz;
5138     tr[2]=TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
5139     tr[3]=TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
5140
5141     if (found){
5142       clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5143       clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5144       cltype[0]=layer.GetCluster(ci)->GetNy();
5145       cltype[1]=layer.GetCluster(ci)->GetNz();
5146
5147      // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5148      //  X->50/sqrt(12)=14 micron   Z->450/sqrt(12)= 120 micron)
5149      // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5150      // It is computed properly by calling the method
5151      // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5152      // T
5153      //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5154       //if (tmp.PropagateTo(x,0.,0.)) {
5155         chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5156         AliCluster c(*layer.GetCluster(ci));
5157         c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5158         c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5159         //if (layer.GetCluster(ci)->GetGlobalCov(cov))  // by using this, instead, you got nominal cluster errors
5160         clu[2]=TMath::Sqrt(c.GetSigmaY2());
5161         clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5162       //}
5163     }
5164   // Compute the angles between the track and the module
5165       // compute the angle "in phi direction", i.e. the angle in the transverse plane
5166       // between the normal to the module and the projection (in the transverse plane) of the
5167       // track trajectory
5168     // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
5169     Float_t tgl = tmp.GetTgl();
5170     Float_t phitr   = tmp.GetSnp();
5171     phitr = TMath::ASin(phitr);
5172     Int_t volId = AliGeomManager::LayerToVolUIDSafe(ilayer+1 ,idet );
5173
5174     Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
5175     Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
5176    Double_t alpha =0.;
5177     alpha = tmp.GetAlpha();
5178     Double_t phiglob = alpha+phitr;
5179     Double_t p[3];
5180     p[0] = TMath::Cos(phiglob);
5181     p[1] = TMath::Sin(phiglob);
5182     p[2] = tgl;
5183     TVector3 pvec(p[0],p[1],p[2]);
5184     TVector3 normvec(rot[1],rot[4],rot[7]);
5185     Double_t angle = pvec.Angle(normvec);
5186
5187     if(angle>0.5*TMath::Pi()) angle = (TMath::Pi()-angle);
5188     angle *= 180./TMath::Pi();
5189
5190     //Trasverse Plane
5191     TVector3 pt(p[0],p[1],0);
5192     TVector3 normt(rot[1],rot[4],0);
5193     Double_t anglet = pt.Angle(normt);
5194
5195     Double_t phiPt = TMath::ATan2(p[1],p[0]);
5196     if(phiPt<0)phiPt+=2.*TMath::Pi();
5197     Double_t phiNorm = TMath::ATan2(rot[4],rot[1]);
5198     if(phiNorm<0) phiNorm+=2.*TMath::Pi();
5199     if(anglet>0.5*TMath::Pi()) anglet = (TMath::Pi()-anglet);
5200     if(phiNorm>phiPt) anglet*=-1.;// pt-->normt  clockwise: anglet>0
5201     if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
5202     anglet *= 180./TMath::Pi();
5203
5204      angleModTrack[2]=(Float_t) angle;
5205      angleModTrack[0]=(Float_t) anglet;
5206      // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
5207     angleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
5208     angleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
5209     angleModTrack[1]*=180./TMath::Pi(); // in degree
5210
5211     fPlaneEff->FillHistos(key,found,tr,clu,cltype,angleModTrack);
5212
5213     // For FO efficiency studies of SPD
5214     if(ilayer<2 && !fSPDChipIntPlaneEff[key]) fPlaneEff->FillHistos(keyFO,foundFO,tr,clu,cltype,angleModTrack);
5215   }
5216   if(ilayer<2) fSPDChipIntPlaneEff[key]=kTRUE;
5217 return;
5218 }
5219
5220 Int_t AliITStrackerMI::FindClusterOfTrack(int label, int lr, int* store) const //RS
5221 {
5222   // find the MC cluster for the label
5223   return fgLayers[lr].FindClusterForLabel(label,store);
5224 }
5225
5226 /*
5227 Int_t AliITStrackerMI::GetPattern(const AliITStrackMI* track, char* patt)
5228 {
5229   // creates pattarn of hits marking fake/corrects by f/c. Used for debugging (RS)
5230   strncpy(patt,"......",6);
5231   int tpcLabel = 0;
5232   if (track->GetESDtrack()) tpcLabel = track->GetESDtrack()->GetTPCLabel();
5233   //
5234   int nwrong = 0;
5235   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
5236     Int_t cindex = track->GetClusterIndex(i);
5237     Int_t l=(cindex & 0xf0000000) >> 28;
5238     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
5239     Int_t isWrong=1;
5240     for (Int_t ind=0;ind<3;ind++) if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
5241     patt[l] = isWrong ? 'f':'c';
5242     nwrong+=isWrong;
5243   }
5244   return nwrong;
5245 }
5246 */
5247 //------------------------------------------------------------------------
5248 Int_t AliITStrackerMI::AliITSlayer::FindClusterForLabel(Int_t label, Int_t *store) const
5249 { //RS
5250   //--------------------------------------------------------------------
5251
5252   int nfound = 0;
5253   for (int ic=0;ic<fN;ic++) {
5254     const AliITSRecPoint *cl = GetCluster(ic);
5255     for (int i=0;i<3;i++) if (cl->GetLabel(i)==label) {
5256         if (nfound<50) {
5257           if (store) store[nfound] = ic;
5258           nfound++;
5259         }
5260         break;
5261       }
5262   }
5263   return nfound;
5264 }
5265