ITS/AliITStrackerMI.cxx

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