ITS/AliITStrackerMI.cxx

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