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