ITS/AliITStrackerMI.cxx

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