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