ITS/AliITStrackerMI.cxx

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