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