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