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