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      if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2307
2308      track->SetDetectorIndex(idet);
2309      if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2310
2311      Double_t dz,zmin,zmax,dy,ymin,ymax;
2312
2313      const AliITSRecPoint *clAcc=0;
2314      Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2315
2316      Int_t idx=index[ilayer];
2317      if (idx>=0) { // cluster in this layer
2318        modstatus = 6; // no refit
2319        const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2320        if (cl) {
2321          if (idet != cl->GetDetectorIndex()) {
2322            idet=cl->GetDetectorIndex();
2323            const AliITSdetector &detc=layer.GetDetector(idet);
2324            if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2325            track->SetDetectorIndex(idet);
2326            if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2327          }
2328          Int_t cllayer = (idx & 0xf0000000) >> 28;;
2329          Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2330          if (chi2<maxchi2) {
2331            clAcc=cl;
2332            maxchi2=chi2;
2333            modstatus = 1; // found
2334          } else {
2335             return kFALSE; //
2336          }
2337        }
2338      } else { // no cluster in this layer
2339        if (skip==1) {
2340          modstatus = 3; // skipped
2341          // Plane Eff determination:
2342          if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2343            if (IsOKForPlaneEff(track,clusters,ilayer))  // only adequate track for plane eff. evaluation
2344               UseTrackForPlaneEff(track,ilayer);
2345          }
2346        } else {
2347          modstatus = 5; // no cls in road
2348          // check dead
2349          if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2350          dz = 0.5*(zmax-zmin);
2351          dy = 0.5*(ymax-ymin);
2352          Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2353          if (dead==1) modstatus = 7; // holes in z in SPD
2354          if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2355        }
2356      }
2357
2358      if (clAcc) {
2359        if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2360        track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2361      }
2362      track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2363
2364
2365      if (extra) { // search for extra clusters in overlapped modules
2366        AliITStrackV2 tmp(*track);
2367        if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2368        layer.SelectClusters(zmin,zmax,ymin,ymax);
2369
2370        const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2371        Int_t idetExtra=-1;
2372        maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2373        Double_t tolerance=0.1;
2374        while ((clExtra=layer.GetNextCluster(ci))!=0) {
2375          // only clusters in another module! (overlaps)
2376          idetExtra = clExtra->GetDetectorIndex();
2377          if (idet == idetExtra) continue;
2378
2379          const AliITSdetector &detx=layer.GetDetector(idetExtra);
2380
2381          if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2382          if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2383          if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2384          if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2385
2386          Double_t chi2=tmp.GetPredictedChi2(clExtra);
2387          if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2388        }
2389        if (cci>=0) {
2390          track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2391          track->SetExtraModule(ilayer,idetExtra);
2392        }
2393      } // end search for extra clusters in overlapped modules
2394
2395      // Correct for material of the current layer
2396      // cross material
2397      // add time if going outward
2398      if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2399
2400   } // end loop on layers
2401
2402   if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2403
2404   return kTRUE;
2405 }
2406 //------------------------------------------------------------------------
2407 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2408 {
2409   //
2410   // calculate normalized chi2
2411   //  return NormalizedChi2(track,0);
2412   Float_t chi2 = 0;
2413   Float_t sum=0;
2414   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2415   //  track->fdEdxMismatch=0;
2416   Float_t dedxmismatch =0;
2417   Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2418   if (mode<100){
2419     for (Int_t i = 0;i<6;i++){
2420       if (track->GetClIndex(i)>=0){
2421         Float_t cerry, cerrz;
2422         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2423         else
2424           { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2425         cerry*=cerry;
2426         cerrz*=cerrz;
2427         Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2428         if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2429           Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2430           if (ratio<0.5) {
2431             cchi2+=(0.5-ratio)*10.;
2432             //track->fdEdxMismatch+=(0.5-ratio)*10.;
2433             dedxmismatch+=(0.5-ratio)*10.;
2434           }
2435         }
2436         if (i<2 ||i>3){
2437           AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2438           Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2439           if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2440           if (i<2) chi2+=2*cl->GetDeltaProbability();
2441         }
2442         chi2+=cchi2;
2443         sum++;
2444       }
2445     }
2446     if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2447       track->SetdEdxMismatch(dedxmismatch);
2448     }
2449   }
2450   else{
2451     for (Int_t i = 0;i<4;i++){
2452       if (track->GetClIndex(i)>=0){
2453         Float_t cerry, cerrz;
2454         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2455         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2456         cerry*=cerry;
2457         cerrz*=cerrz;
2458         chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2459         chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2460         sum++;
2461       }
2462     }
2463     for (Int_t i = 4;i<6;i++){
2464       if (track->GetClIndex(i)>=0){
2465         Float_t cerry, cerrz;
2466         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2467         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2468         cerry*=cerry;
2469         cerrz*=cerrz;
2470         Float_t cerryb, cerrzb;
2471         if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2472         else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2473         cerryb*=cerryb;
2474         cerrzb*=cerrzb;
2475         chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2476         chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2477         sum++;
2478       }
2479     }
2480   }
2481   if (track->GetESDtrack()->GetTPCsignal()>85){
2482     Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2483     if (ratio<0.5) {
2484       chi2+=(0.5-ratio)*5.;
2485     }
2486     if (ratio>2){
2487       chi2+=(ratio-2.0)*3;
2488     }
2489   }
2490   //
2491   Double_t match = TMath::Sqrt(track->GetChi22());
2492   if (track->GetConstrain())  match/=track->GetNumberOfClusters();
2493   if (!track->GetConstrain()) {
2494     if (track->GetNumberOfClusters()>2) {
2495       match/=track->GetNumberOfClusters()-2.;
2496     } else {
2497       match=0;
2498     }
2499   }
2500   if (match<0) match=0;
2501
2502   // penalty factor for missing points (NDeadZone>0), but no penalty
2503   // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2504   // or there is a dead from OCDB)
2505   Float_t deadzonefactor = 0.;
2506   if (track->GetNDeadZone()>0.) {
2507     Float_t sumDeadZoneProbability=0;
2508     for(Int_t ilay=0;ilay<6;ilay++) sumDeadZoneProbability+=track->GetDeadZoneProbability(ilay);
2509     Float_t nDeadZoneWithProbNot1=(Float_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2510     if(nDeadZoneWithProbNot1>0.) {
2511       Float_t deadZoneProbability = sumDeadZoneProbability - (Float_t)((Int_t)sumDeadZoneProbability);
2512       deadZoneProbability /= nDeadZoneWithProbNot1;
2513       deadzonefactor = 3.*(1.1-deadZoneProbability);
2514     }
2515   }
2516
2517   Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2518     (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2519                                 1./(1.+track->GetNSkipped()));
2520
2521  return normchi2;
2522 }
2523 //------------------------------------------------------------------------
2524 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2525 {
2526   //
2527   // return matching chi2 between two tracks
2528   Double_t largeChi2=1000.;
2529
2530   AliITStrackMI track3(*track2);
2531   if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2532   TMatrixD vec(5,1);
2533   vec(0,0)=track1->GetY()   - track3.GetY();
2534   vec(1,0)=track1->GetZ()   - track3.GetZ();
2535   vec(2,0)=track1->GetSnp() - track3.GetSnp();
2536   vec(3,0)=track1->GetTgl() - track3.GetTgl();
2537   vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2538   //
2539   TMatrixD cov(5,5);
2540   cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2541   cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2542   cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2543   cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2544   cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2545
2546   cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2547   cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2548   cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2549   cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2550   //
2551   cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2552   cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2553   cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2554   //
2555   cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2556   cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2557   //
2558   cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2559
2560   cov.Invert();
2561   TMatrixD vec2(cov,TMatrixD::kMult,vec);
2562   TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2563   return chi2(0,0);
2564 }
2565 //------------------------------------------------------------------------
2566 Double_t  AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2567 {
2568   //
2569   //  return probability that given point (characterized by z position and error)
2570   //  is in SPD dead zone
2571   //
2572   Double_t probability = 0.;
2573   Double_t absz = TMath::Abs(zpos);
2574   Double_t nearestz = (absz<2.) ? 0.5*(fSPDdetzcentre[1]+fSPDdetzcentre[2]) :
2575                                   0.5*(fSPDdetzcentre[2]+fSPDdetzcentre[3]);
2576   if (TMath::Abs(absz-nearestz)>0.25+3.*zerr) return probability;
2577   Double_t zmin, zmax;
2578   if (zpos<-6.) { // dead zone at z = -7
2579     zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2580     zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2581   } else if (zpos>6.) { // dead zone at z = +7
2582     zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2583     zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2584   } else if (absz<2.) { // dead zone at z = 0
2585     zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2586     zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2587   } else {
2588     zmin = 0.;
2589     zmax = 0.;
2590   }
2591   // probability that the true z is in the range [zmin,zmax] (i.e. inside
2592   // dead zone)
2593   probability = 0.5*( TMath::Erf((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2594                       TMath::Erf((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2595   return probability;
2596 }
2597 //------------------------------------------------------------------------
2598 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2599 {
2600   //
2601   // calculate normalized chi2
2602   Float_t chi2[6];
2603   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2604   Float_t ncl = 0;
2605   for (Int_t i = 0;i<6;i++){
2606     if (TMath::Abs(track->GetDy(i))>0){
2607       chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2608       chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2609       ncl++;
2610     }
2611     else{chi2[i]=10000;}
2612   }
2613   Int_t index[6];
2614   TMath::Sort(6,chi2,index,kFALSE);
2615   Float_t max = float(ncl)*fac-1.;
2616   Float_t sumchi=0, sumweight=0;
2617   for (Int_t i=0;i<max+1;i++){
2618     Float_t weight = (i<max)?1.:(max+1.-i);
2619     sumchi+=weight*chi2[index[i]];
2620     sumweight+=weight;
2621   }
2622   Double_t normchi2 = sumchi/sumweight;
2623   return normchi2;
2624 }
2625 //------------------------------------------------------------------------
2626 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2627 {
2628   //
2629   // calculate normalized chi2
2630   //  if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2631   Int_t npoints = 0;
2632   Double_t res =0;
2633   for (Int_t i=0;i<6;i++){
2634     if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2635     Double_t sy1 = forwardtrack->GetSigmaY(i);
2636     Double_t sz1 = forwardtrack->GetSigmaZ(i);
2637     Double_t sy2 = backtrack->GetSigmaY(i);
2638     Double_t sz2 = backtrack->GetSigmaZ(i);
2639     if (i<2){ sy2=1000.;sz2=1000;}
2640     //
2641     Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2642     Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2643     //
2644     Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2645     Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2646     //
2647     res+= nz0*nz0+ny0*ny0;
2648     npoints++;
2649   }
2650   if (npoints>1) return
2651                    TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2652                    //2*forwardtrack->fNUsed+
2653                    res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2654                                   1./(1.+forwardtrack->GetNSkipped()));
2655   return 1000;
2656 }
2657 //------------------------------------------------------------------------
2658 Float_t  *AliITStrackerMI::GetWeight(Int_t index) {
2659   //--------------------------------------------------------------------
2660   //       Return pointer to a given cluster
2661   //--------------------------------------------------------------------
2662   Int_t l=(index & 0xf0000000) >> 28;
2663   Int_t c=(index & 0x0fffffff) >> 00;
2664   return fgLayers[l].GetWeight(c);
2665 }
2666 //------------------------------------------------------------------------
2667 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2668 {
2669   //---------------------------------------------
2670   // register track to the list
2671   //
2672   if (track->GetESDtrack()->GetKinkIndex(0)!=0) return;  //don't register kink tracks
2673   //
2674   //
2675   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2676     Int_t index = track->GetClusterIndex(icluster);
2677     Int_t l=(index & 0xf0000000) >> 28;
2678     Int_t c=(index & 0x0fffffff) >> 00;
2679     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2680     for (Int_t itrack=0;itrack<4;itrack++){
2681       if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2682         fgLayers[l].SetClusterTracks(itrack,c,id);
2683         break;
2684       }
2685     }
2686   }
2687 }
2688 //------------------------------------------------------------------------
2689 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2690 {
2691   //---------------------------------------------
2692   // unregister track from the list
2693   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2694     Int_t index = track->GetClusterIndex(icluster);
2695     Int_t l=(index & 0xf0000000) >> 28;
2696     Int_t c=(index & 0x0fffffff) >> 00;
2697     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2698     for (Int_t itrack=0;itrack<4;itrack++){
2699       if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2700         fgLayers[l].SetClusterTracks(itrack,c,-1);
2701       }
2702     }
2703   }
2704 }
2705 //------------------------------------------------------------------------
2706 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2707 {
2708   //-------------------------------------------------------------
2709   //get number of shared clusters
2710   //-------------------------------------------------------------
2711   Float_t shared=0;
2712   for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2713   // mean  number of clusters
2714   Float_t *ny = GetNy(id), *nz = GetNz(id);
2715
2716
2717   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2718     Int_t index = track->GetClusterIndex(icluster);
2719     Int_t l=(index & 0xf0000000) >> 28;
2720     Int_t c=(index & 0x0fffffff) >> 00;
2721     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2722     if (ny[l]==0){
2723       printf("problem\n");
2724     }
2725     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2726     Float_t weight=1;
2727     //
2728     Float_t deltan = 0;
2729     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2730     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2731       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2732     if (l<2 || l>3){
2733       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2734     }
2735     else{
2736       deltan = (cl->GetNz()-nz[l]);
2737     }
2738     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2739     weight = 2./TMath::Max(3.+deltan,2.);
2740     //
2741     for (Int_t itrack=0;itrack<4;itrack++){
2742       if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2743         list[l]=index;
2744         clist[l] = (AliITSRecPoint*)GetCluster(index);
2745         shared+=weight;
2746         break;
2747       }
2748     }
2749   }
2750   track->SetNUsed(shared);
2751   return shared;
2752 }
2753 //------------------------------------------------------------------------
2754 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2755 {
2756   //
2757   // find first shared track
2758   //
2759   // mean  number of clusters
2760   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2761   //
2762   for (Int_t i=0;i<6;i++) overlist[i]=-1;
2763   Int_t sharedtrack=100000;
2764   Int_t tracks[24],trackindex=0;
2765   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2766   //
2767   for (Int_t icluster=0;icluster<6;icluster++){
2768     if (clusterlist[icluster]<0) continue;
2769     Int_t index = clusterlist[icluster];
2770     Int_t l=(index & 0xf0000000) >> 28;
2771     Int_t c=(index & 0x0fffffff) >> 00;
2772     if (ny[l]==0){
2773       printf("problem\n");
2774     }
2775     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2776     //if (l>3) continue;
2777     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2778     //
2779     Float_t deltan = 0;
2780     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2781     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2782       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2783     if (l<2 || l>3){
2784       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2785     }
2786     else{
2787       deltan = (cl->GetNz()-nz[l]);
2788     }
2789     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2790     //
2791     for (Int_t itrack=3;itrack>=0;itrack--){
2792       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2793       if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2794        tracks[trackindex]  = fgLayers[l].GetClusterTracks(itrack,c);
2795        trackindex++;
2796       }
2797     }
2798   }
2799   if (trackindex==0) return -1;
2800   if (trackindex==1){
2801     sharedtrack = tracks[0];
2802   }else{
2803     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2804     else{
2805       //
2806       Int_t tracks2[24], cluster[24];
2807       for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2808       Int_t index =0;
2809       //
2810       for (Int_t i=0;i<trackindex;i++){
2811         if (tracks[i]<0) continue;
2812         tracks2[index] = tracks[i];
2813         cluster[index]++;
2814         for (Int_t j=i+1;j<trackindex;j++){
2815           if (tracks[j]<0) continue;
2816           if (tracks[j]==tracks[i]){
2817             cluster[index]++;
2818             tracks[j]=-1;
2819           }
2820         }
2821         index++;
2822       }
2823       Int_t max=0;
2824       for (Int_t i=0;i<index;i++){
2825         if (cluster[index]>max) {
2826           sharedtrack=tracks2[index];
2827           max=cluster[index];
2828         }
2829       }
2830     }
2831   }
2832
2833   if (sharedtrack>=100000) return -1;
2834   //
2835   // make list of overlaps
2836   shared =0;
2837   for (Int_t icluster=0;icluster<6;icluster++){
2838     if (clusterlist[icluster]<0) continue;
2839     Int_t index = clusterlist[icluster];
2840     Int_t l=(index & 0xf0000000) >> 28;
2841     Int_t c=(index & 0x0fffffff) >> 00;
2842     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2843     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2844     if (l==0 || l==1){
2845       if (cl->GetNy()>2) continue;
2846       if (cl->GetNz()>2) continue;
2847     }
2848     if (l==4 || l==5){
2849       if (cl->GetNy()>3) continue;
2850       if (cl->GetNz()>3) continue;
2851     }
2852     //
2853     for (Int_t itrack=3;itrack>=0;itrack--){
2854       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2855       if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2856         overlist[l]=index;
2857         shared++;
2858       }
2859     }
2860   }
2861   return sharedtrack;
2862 }
2863 //------------------------------------------------------------------------
2864 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2865   //
2866   // try to find track hypothesys without conflicts
2867   // with minimal chi2;
2868   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2869   Int_t entries1 = arr1->GetEntriesFast();
2870   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2871   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2872   Int_t entries2 = arr2->GetEntriesFast();
2873   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2874   //
2875   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2876   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2877   if (track10->Pt()>0.5+track20->Pt()) return track10;
2878
2879   for (Int_t itrack=0;itrack<entries1;itrack++){
2880     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2881     UnRegisterClusterTracks(track,trackID1);
2882   }
2883   //
2884   for (Int_t itrack=0;itrack<entries2;itrack++){
2885     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2886     UnRegisterClusterTracks(track,trackID2);
2887   }
2888   Int_t index1=0;
2889   Int_t index2=0;
2890   Float_t maxconflicts=6;
2891   Double_t maxchi2 =1000.;
2892   //
2893   // get weight of hypothesys - using best hypothesy
2894   Double_t w1,w2;
2895
2896   Int_t list1[6],list2[6];
2897   AliITSRecPoint *clist1[6], *clist2[6] ;
2898   RegisterClusterTracks(track10,trackID1);
2899   RegisterClusterTracks(track20,trackID2);
2900   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2901   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2902   UnRegisterClusterTracks(track10,trackID1);
2903   UnRegisterClusterTracks(track20,trackID2);
2904   //
2905   // normalized chi2
2906   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2907   Float_t nerry[6],nerrz[6];
2908   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2909   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2910   for (Int_t i=0;i<6;i++){
2911      if ( (erry1[i]>0) && (erry2[i]>0)) {
2912        nerry[i] = TMath::Min(erry1[i],erry2[i]);
2913        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2914      }else{
2915        nerry[i] = TMath::Max(erry1[i],erry2[i]);
2916        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2917      }
2918      if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2919        chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2920        chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2921        ncl1++;
2922      }
2923      if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2924        chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2925        chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2926        ncl2++;
2927      }
2928   }
2929   chi21/=ncl1;
2930   chi22/=ncl2;
2931   //
2932   //
2933   Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2934   Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2935   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2936   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2937   //
2938   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2939         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2940         +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2941         );
2942   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2943         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2944         +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2945         );
2946
2947   Double_t sumw = w1+w2;
2948   w1/=sumw;
2949   w2/=sumw;
2950   if (w1<w2*0.5) {
2951     w1 = (d2+0.5)/(d1+d2+1.);
2952     w2 = (d1+0.5)/(d1+d2+1.);
2953   }
2954   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2955   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2956   //
2957   // get pair of "best" hypothesys
2958   //
2959   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2960   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2961
2962   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2963     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2964     //if (track1->fFakeRatio>0) continue;
2965     RegisterClusterTracks(track1,trackID1);
2966     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2967       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2968
2969       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2970       //if (track2->fFakeRatio>0) continue;
2971       Float_t nskipped=0;
2972       RegisterClusterTracks(track2,trackID2);
2973       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2974       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2975       UnRegisterClusterTracks(track2,trackID2);
2976       //
2977       if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2978       if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2979       if (nskipped>0.5) continue;
2980       //
2981       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2982       if (conflict1+1<cconflict1) continue;
2983       if (conflict2+1<cconflict2) continue;
2984       Float_t conflict=0;
2985       Float_t sumchi2=0;
2986       Float_t sum=0;
2987       for (Int_t i=0;i<6;i++){
2988         //
2989         Float_t c1 =0.; // conflict coeficients
2990         Float_t c2 =0.;
2991         if (clist1[i]&&clist2[i]){
2992           Float_t deltan = 0;
2993           if (i<2 || i>3){
2994             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2995           }
2996           else{
2997             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2998           }
2999           c1  = 2./TMath::Max(3.+deltan,2.);
3000           c2  = 2./TMath::Max(3.+deltan,2.);
3001         }
3002         else{
3003           if (clist1[i]){
3004             Float_t deltan = 0;
3005             if (i<2 || i>3){
3006               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3007             }
3008             else{
3009               deltan = (clist1[i]->GetNz()-nz1[i]);
3010             }
3011             c1  = 2./TMath::Max(3.+deltan,2.);
3012             c2  = 0;
3013           }
3014
3015           if (clist2[i]){
3016             Float_t deltan = 0;
3017             if (i<2 || i>3){
3018               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3019             }
3020             else{
3021               deltan = (clist2[i]->GetNz()-nz2[i]);
3022             }
3023             c2  = 2./TMath::Max(3.+deltan,2.);
3024             c1  = 0;
3025           }
3026         }
3027         //
3028         chi21=0;chi22=0;
3029         if (TMath::Abs(track1->GetDy(i))>0.) {
3030           chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3031             (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3032           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3033           //  (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3034         }else{
3035           if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3036         }
3037         //
3038         if (TMath::Abs(track2->GetDy(i))>0.) {
3039           chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3040             (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3041           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3042           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3043         }
3044         else{
3045           if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3046         }
3047         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3048         if (chi21>0) sum+=w1;
3049         if (chi22>0) sum+=w2;
3050         conflict+=(c1+c2);
3051       }
3052       Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3053       if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3054       Double_t normchi2 = 2*conflict+sumchi2/sum;
3055       if ( normchi2 <maxchi2 ){
3056         index1 = itrack1;
3057         index2 = itrack2;
3058         maxconflicts = conflict;
3059         maxchi2 = normchi2;
3060       }
3061     }
3062     UnRegisterClusterTracks(track1,trackID1);
3063   }
3064   //
3065   //  if (maxconflicts<4 && maxchi2<th0){
3066   if (maxchi2<th0*2.){
3067     Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3068     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3069     track1->SetChi2MIP(5,maxconflicts);
3070     track1->SetChi2MIP(6,maxchi2);
3071     track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3072     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
3073     track1->SetChi2MIP(8,index1);
3074     fBestTrackIndex[trackID1] =index1;
3075     UpdateESDtrack(track1, AliESDtrack::kITSin);
3076   }
3077   else if (track10->GetChi2MIP(0)<th1){
3078     track10->SetChi2MIP(5,maxconflicts);
3079     track10->SetChi2MIP(6,maxchi2);
3080     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
3081     UpdateESDtrack(track10,AliESDtrack::kITSin);
3082   }
3083
3084   for (Int_t itrack=0;itrack<entries1;itrack++){
3085     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3086     UnRegisterClusterTracks(track,trackID1);
3087   }
3088   //
3089   for (Int_t itrack=0;itrack<entries2;itrack++){
3090     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3091     UnRegisterClusterTracks(track,trackID2);
3092   }
3093
3094   if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3095       &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3096     //  if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3097   //    &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3098     RegisterClusterTracks(track10,trackID1);
3099   }
3100   if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3101       &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3102     //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3103     //  &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3104     RegisterClusterTracks(track20,trackID2);
3105   }
3106   return track10;
3107
3108 }
3109 //------------------------------------------------------------------------
3110 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3111   //--------------------------------------------------------------------
3112   // This function marks clusters assigned to the track
3113   //--------------------------------------------------------------------
3114   AliTracker::UseClusters(t,from);
3115
3116   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3117   //if (c->GetQ()>2) c->Use();
3118   if (c->GetSigmaZ2()>0.1) c->Use();
3119   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3120   //if (c->GetQ()>2) c->Use();
3121   if (c->GetSigmaZ2()>0.1) c->Use();
3122
3123 }
3124 //------------------------------------------------------------------------
3125 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3126 {
3127   //------------------------------------------------------------------
3128   // add track to the list of hypothesys
3129   //------------------------------------------------------------------
3130
3131   if (esdindex>=fTrackHypothesys.GetEntriesFast())
3132     fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3133   //
3134   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3135   if (!array) {
3136     array = new TObjArray(10);
3137     fTrackHypothesys.AddAt(array,esdindex);
3138   }
3139   array->AddLast(track);
3140 }
3141 //------------------------------------------------------------------------
3142 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3143 {
3144   //-------------------------------------------------------------------
3145   // compress array of track hypothesys
3146   // keep only maxsize best hypothesys
3147   //-------------------------------------------------------------------
3148   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3149   if (! (fTrackHypothesys.At(esdindex)) ) return;
3150   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3151   Int_t entries = array->GetEntriesFast();
3152   //
3153   //- find preliminary besttrack as a reference
3154   Float_t minchi2=10000;
3155   Int_t maxn=0;
3156   AliITStrackMI * besttrack=0;
3157   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3158     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3159     if (!track) continue;
3160     Float_t chi2 = NormalizedChi2(track,0);
3161     //
3162     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3163     track->SetLabel(tpcLabel);
3164     CookdEdx(track);
3165     track->SetFakeRatio(1.);
3166     CookLabel(track,0.); //For comparison only
3167     //
3168     //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3169     if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3170       if (track->GetNumberOfClusters()<maxn) continue;
3171       maxn = track->GetNumberOfClusters();
3172       if (chi2<minchi2){
3173         minchi2=chi2;
3174         besttrack=track;
3175       }
3176     }
3177     else{
3178       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3179         delete array->RemoveAt(itrack);
3180       }
3181     }
3182   }
3183   if (!besttrack) return;
3184   //
3185   //
3186   //take errors of best track as a reference
3187   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3188   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3189   for (Int_t j=0;j<6;j++) {
3190     if (besttrack->GetClIndex(j)>=0){
3191       erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3192       errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3193       ny[j]   = besttrack->GetNy(j);
3194       nz[j]   = besttrack->GetNz(j);
3195     }
3196   }
3197   //
3198   // calculate normalized chi2
3199   //
3200   Float_t * chi2        = new Float_t[entries];
3201   Int_t * index         = new Int_t[entries];
3202   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3203   for (Int_t itrack=0;itrack<entries;itrack++){
3204     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3205     if (track){
3206       track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3207       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3208         chi2[itrack] = track->GetChi2MIP(0);
3209       else{
3210         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3211           delete array->RemoveAt(itrack);
3212         }
3213       }
3214     }
3215   }
3216   //
3217   TMath::Sort(entries,chi2,index,kFALSE);
3218   besttrack = (AliITStrackMI*)array->At(index[0]);
3219   if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3220     for (Int_t j=0;j<6;j++){
3221       if (besttrack->GetClIndex(j)>=0){
3222         erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3223         errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3224         ny[j]   = besttrack->GetNy(j);
3225         nz[j]   = besttrack->GetNz(j);
3226       }
3227     }
3228   }
3229   //
3230   // calculate one more time with updated normalized errors
3231   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3232   for (Int_t itrack=0;itrack<entries;itrack++){
3233     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3234     if (track){
3235       track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3236       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3237         chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3238       else
3239         {
3240           if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3241             delete array->RemoveAt(itrack);
3242           }
3243         }
3244     }
3245   }
3246   entries = array->GetEntriesFast();
3247   //
3248   //
3249   if (entries>0){
3250     TObjArray * newarray = new TObjArray();
3251     TMath::Sort(entries,chi2,index,kFALSE);
3252     besttrack = (AliITStrackMI*)array->At(index[0]);
3253     if (besttrack){
3254       //
3255       for (Int_t j=0;j<6;j++){
3256         if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3257           erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3258           errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3259           ny[j]   = besttrack->GetNy(j);
3260           nz[j]   = besttrack->GetNz(j);
3261         }
3262       }
3263       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3264       minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3265       Float_t minn = besttrack->GetNumberOfClusters()-3;
3266       Int_t accepted=0;
3267       for (Int_t i=0;i<entries;i++){
3268         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3269         if (!track) continue;
3270         if (accepted>maxcut) break;
3271         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3272         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3273           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3274             delete array->RemoveAt(index[i]);
3275             continue;
3276           }
3277         }
3278         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3279         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3280           if (!shortbest) accepted++;
3281           //
3282           newarray->AddLast(array->RemoveAt(index[i]));
3283           for (Int_t j=0;j<6;j++){
3284             if (nz[j]==0){
3285               erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3286               errz[j] = track->GetSigmaZ(j); errz[j]   = track->GetSigmaZ(j+6);
3287               ny[j]   = track->GetNy(j);
3288               nz[j]   = track->GetNz(j);
3289             }
3290           }
3291         }
3292         else{
3293           delete array->RemoveAt(index[i]);
3294         }
3295       }
3296       array->Delete();
3297       delete fTrackHypothesys.RemoveAt(esdindex);
3298       fTrackHypothesys.AddAt(newarray,esdindex);
3299     }
3300     else{
3301       array->Delete();
3302       delete fTrackHypothesys.RemoveAt(esdindex);
3303     }
3304   }
3305   delete [] chi2;
3306   delete [] index;
3307 }
3308 //------------------------------------------------------------------------
3309 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3310 {
3311   //-------------------------------------------------------------
3312   // try to find best hypothesy
3313   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3314   //-------------------------------------------------------------
3315   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3316   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3317   if (!array) return 0;
3318   Int_t entries = array->GetEntriesFast();
3319   if (!entries) return 0;
3320   Float_t minchi2 = 100000;
3321   AliITStrackMI * besttrack=0;
3322   //
3323   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
3324   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3325   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3326   Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3327   //
3328   for (Int_t i=0;i<entries;i++){
3329     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3330     if (!track) continue;
3331     Float_t sigmarfi,sigmaz;
3332     GetDCASigma(track,sigmarfi,sigmaz);
3333     track->SetDnorm(0,sigmarfi);
3334     track->SetDnorm(1,sigmaz);
3335     //
3336     track->SetChi2MIP(1,1000000);
3337     track->SetChi2MIP(2,1000000);
3338     track->SetChi2MIP(3,1000000);
3339     //
3340     // backtrack
3341     backtrack = new(backtrack) AliITStrackMI(*track);
3342     if (track->GetConstrain()) {
3343       if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3344       if (!backtrack->Improve(0,xyzVtx,ersVtx))         continue;
3345       backtrack->ResetCovariance(10.);
3346     }else{
3347       backtrack->ResetCovariance(10.);
3348     }
3349     backtrack->ResetClusters();
3350
3351     Double_t x = original->GetX();
3352     if (!RefitAt(x,backtrack,track)) continue;
3353     //
3354     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3355     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3356     if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.)  continue;
3357     track->SetChi22(GetMatchingChi2(backtrack,original));
3358
3359     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
3360     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
3361     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
3362
3363
3364     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1))  continue;
3365     //
3366     //forward track - without constraint
3367     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3368     forwardtrack->ResetClusters();
3369     x = track->GetX();
3370     RefitAt(x,forwardtrack,track);
3371     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3372     if  (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0)  continue;
3373     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2))  continue;
3374
3375     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3376     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3377     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
3378     forwardtrack->SetD(0,track->GetD(0));
3379     forwardtrack->SetD(1,track->GetD(1));
3380     {
3381       Int_t list[6];
3382       AliITSRecPoint* clist[6];
3383       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3384       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3385     }
3386
3387     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3388     if  ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3389     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3390       track->SetChi2MIP(3,1000);
3391       continue;
3392     }
3393     Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3394     //
3395     for (Int_t ichi=0;ichi<5;ichi++){
3396       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3397     }
3398     if (chi2 < minchi2){
3399       //besttrack = new AliITStrackMI(*forwardtrack);
3400       besttrack = track;
3401       besttrack->SetLabel(track->GetLabel());
3402       besttrack->SetFakeRatio(track->GetFakeRatio());
3403       minchi2   = chi2;
3404       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3405       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3406       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
3407     }
3408   }
3409   delete backtrack;
3410   delete forwardtrack;
3411   Int_t accepted=0;
3412   for (Int_t i=0;i<entries;i++){
3413     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3414
3415     if (!track) continue;
3416
3417     if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3418         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3419         track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3420       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3421         delete array->RemoveAt(i);
3422         continue;
3423       }
3424     }
3425     else{
3426       accepted++;
3427     }
3428   }
3429   //
3430   array->Compress();
3431   SortTrackHypothesys(esdindex,checkmax,1);
3432
3433   array = (TObjArray*) fTrackHypothesys.At(esdindex);
3434   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3435   besttrack = (AliITStrackMI*)array->At(0);
3436   if (!besttrack)  return 0;
3437   besttrack->SetChi2MIP(8,0);
3438   fBestTrackIndex[esdindex]=0;
3439   entries = array->GetEntriesFast();
3440   AliITStrackMI *longtrack =0;
3441   minchi2 =1000;
3442   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3443   for (Int_t itrack=entries-1;itrack>0;itrack--) {
3444     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3445     if (!track->GetConstrain()) continue;
3446     if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3447     if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3448     if (track->GetChi2MIP(0)>4.) continue;
3449     minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3450     longtrack =track;
3451   }
3452   //if (longtrack) besttrack=longtrack;
3453
3454   Int_t list[6];
3455   AliITSRecPoint * clist[6];
3456   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3457   if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3458       &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3459     RegisterClusterTracks(besttrack,esdindex);
3460   }
3461   //
3462   //
3463   if (shared>0.0){
3464     Int_t nshared;
3465     Int_t overlist[6];
3466     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3467     if (sharedtrack>=0){
3468       //
3469       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3470       if (besttrack){
3471         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3472       }
3473       else return 0;
3474     }
3475   }
3476
3477   if (shared>2.5) return 0;
3478   if (shared>1.0) return besttrack;
3479   //
3480   // Don't sign clusters if not gold track
3481   //
3482   if (!besttrack->IsGoldPrimary()) return besttrack;
3483   if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
3484   //
3485   if (fConstraint[fPass]){
3486     //
3487     // sign clusters
3488     //
3489     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3490     for (Int_t i=0;i<6;i++){
3491       Int_t index = besttrack->GetClIndex(i);
3492       if (index<0) continue;
3493       Int_t ilayer =  (index & 0xf0000000) >> 28;
3494       if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3495       AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3496       if (!c) continue;
3497       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3498       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3499       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
3500       if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3501         if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3502       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
3503
3504       Bool_t cansign = kTRUE;
3505       for (Int_t itrack=0;itrack<entries; itrack++){
3506         AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3507         if (!track) continue;
3508         if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3509         if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3510           cansign = kFALSE;
3511           break;
3512         }
3513       }
3514       if (cansign){
3515         if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3516         if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3517         if (!c->IsUsed()) c->Use();
3518       }
3519     }
3520   }
3521   return besttrack;
3522 }
3523 //------------------------------------------------------------------------
3524 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3525 {
3526   //
3527   // get "best" hypothesys
3528   //
3529
3530   Int_t nentries = itsTracks.GetEntriesFast();
3531   for (Int_t i=0;i<nentries;i++){
3532     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3533     if (!track) continue;
3534     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3535     if (!array) continue;
3536     if (array->GetEntriesFast()<=0) continue;
3537     //
3538     AliITStrackMI* longtrack=0;
3539     Float_t minn=0;
3540     Float_t maxchi2=1000;
3541     for (Int_t j=0;j<array->GetEntriesFast();j++){
3542       AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3543       if (!trackHyp) continue;
3544       if (trackHyp->GetGoldV0()) {
3545         longtrack = trackHyp;   //gold V0 track taken
3546         break;
3547       }
3548       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3549       Float_t chi2 = trackHyp->GetChi2MIP(0);
3550       if (fAfterV0){
3551         if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3552       }
3553       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3554       //
3555       if (chi2 > maxchi2) continue;
3556       minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3557       maxchi2 = chi2;
3558       longtrack=trackHyp;
3559     }
3560     //
3561     //
3562     //
3563     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3564     if (!longtrack) {longtrack = besttrack;}
3565     else besttrack= longtrack;
3566     //
3567     if (besttrack) {
3568       Int_t list[6];
3569       AliITSRecPoint * clist[6];
3570       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3571       //
3572       track->SetNUsed(shared);
3573       track->SetNSkipped(besttrack->GetNSkipped());
3574       track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3575       if (shared>0) {
3576         if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3577         Int_t nshared;
3578         Int_t overlist[6];
3579         //
3580         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3581         //if (sharedtrack==-1) sharedtrack=0;
3582         if (sharedtrack>=0) {
3583           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3584         }
3585       }
3586       if (besttrack&&fAfterV0) {
3587         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3588       }
3589       if (besttrack&&fConstraint[fPass])
3590         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3591       if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3592         if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3593              TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3594       }
3595
3596     }
3597   }
3598 }
3599 //------------------------------------------------------------------------
3600 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3601   //--------------------------------------------------------------------
3602   //This function "cooks" a track label. If label<0, this track is fake.
3603   //--------------------------------------------------------------------
3604   Int_t tpcLabel=-1;
3605
3606   if ( track->GetESDtrack())   tpcLabel =  TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3607
3608    track->SetChi2MIP(9,0);
3609    Int_t nwrong=0;
3610    for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3611      Int_t cindex = track->GetClusterIndex(i);
3612      Int_t l=(cindex & 0xf0000000) >> 28;
3613      AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3614      Int_t isWrong=1;
3615      for (Int_t ind=0;ind<3;ind++){
3616        if (tpcLabel>0)
3617          if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
3618        AliDebug(2,Form("icl %d  ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3619      }
3620      track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3621      nwrong+=isWrong;
3622    }
3623    Int_t nclusters = track->GetNumberOfClusters();
3624    if (nclusters > 0) //PH Some tracks don't have any cluster
3625      track->SetFakeRatio(double(nwrong)/double(nclusters));
3626    if (tpcLabel>0){
3627      if (track->GetFakeRatio()>wrong) track->SetLabel(-tpcLabel);
3628      else
3629        track->SetLabel(tpcLabel);
3630    }
3631    AliDebug(2,Form(" nls %d wrong %d  label %d  tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
3632
3633 }
3634 //------------------------------------------------------------------------
3635 void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
3636 {
3637   //
3638   track->SetChi2MIP(9,0);
3639   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3640     Int_t cindex = track->GetClusterIndex(i);
3641     Int_t l=(cindex & 0xf0000000) >> 28;
3642     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3643     Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3644     Int_t isWrong=1;
3645     for (Int_t ind=0;ind<3;ind++){
3646       if (cl->GetLabel(ind)==lab) isWrong=0;
3647     }
3648     track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3649   }
3650   Double_t low=0.;
3651   Double_t up=0.51;
3652   track->CookdEdx(low,up);
3653 }
3654 //------------------------------------------------------------------------
3655 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3656   //
3657   // Create some arrays
3658   //
3659   if (fCoefficients) delete []fCoefficients;
3660   fCoefficients = new Float_t[ntracks*48];
3661   for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3662 }
3663 //------------------------------------------------------------------------
3664 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3665 {
3666   //
3667   // Compute predicted chi2
3668   //
3669   Float_t erry,errz,covyz;
3670   Float_t theta = track->GetTgl();
3671   Float_t phi   = track->GetSnp();
3672   phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
3673   AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
3674   AliDebug(2,Form(" chi2: tr-cl   %f  %f   tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
3675   // Take into account the mis-alignment (bring track to cluster plane)
3676   Double_t xTrOrig=track->GetX();
3677   if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
3678   AliDebug(2,Form(" chi2: tr-cl   %f  %f   tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
3679   Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
3680   // Bring the track back to detector plane in ideal geometry
3681   // [mis-alignment will be accounted for in UpdateMI()]
3682   if (!track->Propagate(xTrOrig)) return 1000.;
3683   Float_t ny,nz;
3684   AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3685   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3686   if (delta>1){
3687     chi2+=0.5*TMath::Min(delta/2,2.);
3688     chi2+=2.*cluster->GetDeltaProbability();
3689   }
3690   //
3691   track->SetNy(layer,ny);
3692   track->SetNz(layer,nz);
3693   track->SetSigmaY(layer,erry);
3694   track->SetSigmaZ(layer, errz);
3695   track->SetSigmaYZ(layer,covyz);
3696   //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3697   track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
3698   return chi2;
3699
3700 }
3701 //------------------------------------------------------------------------
3702 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3703 {
3704   //
3705   // Update ITS track
3706   //
3707   Int_t layer = (index & 0xf0000000) >> 28;
3708   track->SetClIndex(layer, index);
3709   if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3710     if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3711       chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3712       track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3713     }
3714   }
3715
3716   if (cl->GetQ()<=0) return 0;  // ingore the "virtual" clusters
3717
3718
3719   // Take into account the mis-alignment (bring track to cluster plane)
3720   Double_t xTrOrig=track->GetX();
3721   Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3722   AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3723   AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3724   AliDebug(2,Form(" xtr %f  xcl %f",track->GetX(),cl->GetX()));
3725
3726   if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
3727
3728   AliCluster c(*cl);
3729   c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3730   c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3731   c.SetSigmaYZ(track->GetSigmaYZ(layer));
3732
3733
3734   Int_t updated = track->UpdateMI(&c,chi2,index);
3735   // Bring the track back to detector plane in ideal geometry
3736   if (!track->Propagate(xTrOrig)) return 0;
3737
3738   if(!updated) AliDebug(2,"update failed");
3739   return updated;
3740 }
3741
3742 //------------------------------------------------------------------------
3743 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3744 {
3745   //
3746   //DCA sigmas parameterization
3747   //to be paramterized using external parameters in future
3748   //
3749   //
3750   sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3751   sigmaz   = 0.0110+4.37*TMath::Abs(track->GetC());
3752 }
3753 //------------------------------------------------------------------------
3754 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
3755 {
3756   //
3757   // Clusters from delta electrons?
3758   //
3759   Int_t entries = clusterArray->GetEntriesFast();
3760   if (entries<4) return;
3761   AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
3762   Int_t layer = cluster->GetLayer();
3763   if (layer>1) return;
3764   Int_t index[10000];
3765   Int_t ncandidates=0;
3766   Float_t r = (layer>0)? 7:4;
3767   //
3768   for (Int_t i=0;i<entries;i++){
3769     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
3770     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3771     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3772     index[ncandidates] = i;  //candidate to belong to delta electron track
3773     ncandidates++;
3774     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3775       cl0->SetDeltaProbability(1);
3776     }
3777   }
3778   //
3779   //
3780   //
3781   for (Int_t i=0;i<ncandidates;i++){
3782     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
3783     if (cl0->GetDeltaProbability()>0.8) continue;
3784     //
3785     Int_t ncl = 0;
3786     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3787     sumy=sumz=sumy2=sumyz=sumw=0.0;
3788     for (Int_t j=0;j<ncandidates;j++){
3789       if (i==j) continue;
3790       AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
3791       //
3792       Float_t dz = cl0->GetZ()-cl1->GetZ();
3793       Float_t dy = cl0->GetY()-cl1->GetY();
3794       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3795         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3796         y[ncl] = cl1->GetY();
3797         z[ncl] = cl1->GetZ();
3798         sumy+= y[ncl]*weight;
3799         sumz+= z[ncl]*weight;
3800         sumy2+=y[ncl]*y[ncl]*weight;
3801         sumyz+=y[ncl]*z[ncl]*weight;
3802         sumw+=weight;
3803         ncl++;
3804       }
3805     }
3806     if (ncl<4) continue;
3807     Float_t det = sumw*sumy2  - sumy*sumy;
3808     Float_t delta=1000;
3809     if (TMath::Abs(det)>0.01){
3810       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
3811       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
3812       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3813     }
3814     else{
3815       Float_t z0  = sumyz/sumy;
3816       delta = TMath::Abs(cl0->GetZ()-z0);
3817     }
3818     if ( delta<0.05) {
3819       cl0->SetDeltaProbability(1-20.*delta);
3820     }
3821   }
3822 }
3823 //------------------------------------------------------------------------
3824 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3825 {
3826   //
3827   // Update ESD track
3828   //
3829   track->UpdateESDtrack(flags);
3830   AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3831   if (oldtrack) delete oldtrack;
3832   track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3833   if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3834     printf("Problem\n");
3835   }
3836 }
3837 //------------------------------------------------------------------------
3838 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3839   //
3840   // Get nearest upper layer close to the point xr.
3841   // rough approximation
3842   //
3843   const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3844   Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3845   Int_t res =6;
3846   for (Int_t i=0;i<6;i++){
3847     if (radius<kRadiuses[i]){
3848       res =i;
3849       break;
3850     }
3851   }
3852   return res;
3853 }
3854 //------------------------------------------------------------------------
3855 void AliITStrackerMI::BuildMaterialLUT(TString material) {
3856   //--------------------------------------------------------------------
3857   // Fill a look-up table with mean material
3858   //--------------------------------------------------------------------
3859
3860   Int_t n=1000;
3861   Double_t mparam[7];
3862   Double_t point1[3],point2[3];
3863   Double_t phi,cosphi,sinphi,z;
3864   // 0-5 layers, 6 pipe, 7-8 shields
3865   Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
3866   Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
3867
3868   Int_t ifirst=0,ilast=0;
3869   if(material.Contains("Pipe")) {
3870     ifirst=6; ilast=6;
3871   } else if(material.Contains("Shields")) {
3872     ifirst=7; ilast=8;
3873   } else if(material.Contains("Layers")) {
3874     ifirst=0; ilast=5;
3875   } else {
3876     Error("BuildMaterialLUT","Wrong layer name\n");
3877   }
3878
3879   for(Int_t imat=ifirst; imat<=ilast; imat++) {
3880     Double_t param[5]={0.,0.,0.,0.,0.};
3881     for (Int_t i=0; i<n; i++) {
3882       phi = 2.*TMath::Pi()*gRandom->Rndm();
3883       cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
3884       z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
3885       point1[0] = rmin[imat]*cosphi;
3886       point1[1] = rmin[imat]*sinphi;
3887       point1[2] = z;
3888       point2[0] = rmax[imat]*cosphi;
3889       point2[1] = rmax[imat]*sinphi;
3890       point2[2] = z;
3891       AliTracker::MeanMaterialBudget(point1,point2,mparam);
3892       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
3893     }
3894     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
3895     if(imat<=5) {
3896       fxOverX0Layer[imat] = param[1];
3897       fxTimesRhoLayer[imat] = param[0]*param[4];
3898     } else if(imat==6) {
3899       fxOverX0Pipe = param[1];
3900       fxTimesRhoPipe = param[0]*param[4];
3901     } else if(imat==7) {
3902       fxOverX0Shield[0] = param[1];
3903       fxTimesRhoShield[0] = param[0]*param[4];
3904     } else if(imat==8) {
3905       fxOverX0Shield[1] = param[1];
3906       fxTimesRhoShield[1] = param[0]*param[4];
3907     }
3908   }
3909   /*
3910   printf("%s\n",material.Data());
3911   printf("%f  %f\n",fxOverX0Pipe,fxTimesRhoPipe);
3912   printf("%f  %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
3913   printf("%f  %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
3914   printf("%f  %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
3915   printf("%f  %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
3916   printf("%f  %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
3917   printf("%f  %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
3918   printf("%f  %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
3919   printf("%f  %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
3920   */
3921   return;
3922 }
3923 //------------------------------------------------------------------------
3924 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
3925                                               TString direction) {
3926   //-------------------------------------------------------------------
3927   // Propagate beyond beam pipe and correct for material
3928   // (material budget in different ways according to fUseTGeo value)
3929   // Add time if going outward (PropagateTo or PropagateToTGeo)
3930   //-------------------------------------------------------------------
3931
3932   // Define budget mode:
3933   // 0: material from AliITSRecoParam (hard coded)
3934   // 1: material from TGeo in one step (on the fly)
3935   // 2: material from lut
3936   // 3: material from TGeo in one step (same for all hypotheses)
3937   Int_t mode;
3938   switch(fUseTGeo) {
3939   case 0:
3940     mode=0;
3941     break;
3942   case 1:
3943     mode=1;
3944     break;
3945   case 2:
3946     mode=2;
3947     break;
3948   case 3:
3949     if(fTrackingPhase.Contains("Clusters2Tracks"))
3950       { mode=3; } else { mode=1; }
3951     break;
3952   case 4:
3953     if(fTrackingPhase.Contains("Clusters2Tracks"))
3954       { mode=3; } else { mode=2; }
3955     break;
3956   default:
3957     mode=0;
3958     break;
3959   }
3960   if(fTrackingPhase.Contains("Default")) mode=0;
3961
3962   Int_t index=fCurrentEsdTrack;
3963
3964   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
3965   Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
3966   Double_t xToGo;
3967   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
3968
3969   Double_t xOverX0,x0,lengthTimesMeanDensity;
3970
3971   switch(mode) {
3972   case 0:
3973     xOverX0 = AliITSRecoParam::GetdPipe();
3974     x0 = AliITSRecoParam::GetX0Be();
3975     lengthTimesMeanDensity = xOverX0*x0;
3976     lengthTimesMeanDensity *= dir;
3977     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
3978     break;
3979   case 1:
3980     if (!t->PropagateToTGeo(xToGo,1)) return 0;
3981     break;
3982   case 2:
3983     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
3984     xOverX0 = fxOverX0Pipe;
3985     lengthTimesMeanDensity = fxTimesRhoPipe;
3986     lengthTimesMeanDensity *= dir;
3987     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
3988     break;
3989   case 3:
3990     if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
3991     if(fxOverX0PipeTrks[index]<0) {
3992       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
3993       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
3994                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
3995       fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
3996       fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
3997       return 1;
3998     }
3999     xOverX0 = fxOverX0PipeTrks[index];
4000     lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4001     lengthTimesMeanDensity *= dir;
4002     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4003     break;
4004   }
4005
4006   return 1;
4007 }
4008 //------------------------------------------------------------------------
4009 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4010                                                 TString shield,
4011                                                 TString direction) {
4012   //-------------------------------------------------------------------
4013   // Propagate beyond SPD or SDD shield and correct for material
4014   // (material budget in different ways according to fUseTGeo value)
4015   // Add time if going outward (PropagateTo or PropagateToTGeo)
4016   //-------------------------------------------------------------------
4017
4018   // Define budget mode:
4019   // 0: material from AliITSRecoParam (hard coded)
4020   // 1: material from TGeo in steps of X cm (on the fly)
4021   //    X = AliITSRecoParam::GetStepSizeTGeo()
4022   // 2: material from lut
4023   // 3: material from TGeo in one step (same for all hypotheses)
4024   Int_t mode;
4025   switch(fUseTGeo) {
4026   case 0:
4027     mode=0;
4028     break;
4029   case 1:
4030     mode=1;
4031     break;
4032   case 2:
4033     mode=2;
4034     break;
4035   case 3:
4036     if(fTrackingPhase.Contains("Clusters2Tracks"))
4037       { mode=3; } else { mode=1; }
4038     break;
4039   case 4:
4040     if(fTrackingPhase.Contains("Clusters2Tracks"))
4041       { mode=3; } else { mode=2; }
4042     break;
4043   default:
4044     mode=0;
4045     break;
4046   }
4047   if(fTrackingPhase.Contains("Default")) mode=0;
4048
4049   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4050   Double_t rToGo;
4051   Int_t    shieldindex=0;
4052   if (shield.Contains("SDD")) { // SDDouter
4053     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4054     shieldindex=1;
4055   } else if (shield.Contains("SPD")) {        // SPDouter
4056     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4057     shieldindex=0;
4058   } else {
4059     Error("CorrectForShieldMaterial"," Wrong shield name\n");
4060     return 0;
4061   }
4062
4063   // do nothing if we are already beyond the shield
4064   Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4065   if(dir<0 && rTrack > rToGo) return 1; // going outward
4066   if(dir>0 && rTrack < rToGo) return 1; // going inward
4067
4068
4069   Double_t xToGo;
4070   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4071
4072   Int_t index=2*fCurrentEsdTrack+shieldindex;
4073
4074   Double_t xOverX0,x0,lengthTimesMeanDensity;
4075   Int_t nsteps=1;
4076
4077   switch(mode) {
4078   case 0:
4079     xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4080     x0 = AliITSRecoParam::GetX0shield(shieldindex);
4081     lengthTimesMeanDensity = xOverX0*x0;
4082     lengthTimesMeanDensity *= dir;
4083     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4084     break;
4085   case 1:
4086     nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4087     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4088     break;
4089   case 2:
4090     if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4091     xOverX0 = fxOverX0Shield[shieldindex];
4092     lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4093     lengthTimesMeanDensity *= dir;
4094     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4095     break;
4096   case 3:
4097     if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4098     if(fxOverX0ShieldTrks[index]<0) {
4099       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4100       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4101                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4102       fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4103       fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4104       return 1;
4105     }
4106     xOverX0 = fxOverX0ShieldTrks[index];
4107     lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4108     lengthTimesMeanDensity *= dir;
4109     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4110     break;
4111   }
4112
4113   return 1;
4114 }
4115 //------------------------------------------------------------------------
4116 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4117                                                Int_t layerindex,
4118                                                Double_t oldGlobXYZ[3],
4119                                                TString direction) {
4120   //-------------------------------------------------------------------
4121   // Propagate beyond layer and correct for material
4122   // (material budget in different ways according to fUseTGeo value)
4123   // Add time if going outward (PropagateTo or PropagateToTGeo)
4124   //-------------------------------------------------------------------
4125
4126   // Define budget mode:
4127   // 0: material from AliITSRecoParam (hard coded)
4128   // 1: material from TGeo in stepsof X cm (on the fly)
4129   //    X = AliITSRecoParam::GetStepSizeTGeo()
4130   // 2: material from lut
4131   // 3: material from TGeo in one step (same for all hypotheses)
4132   Int_t mode;
4133   switch(fUseTGeo) {
4134   case 0:
4135     mode=0;
4136     break;
4137   case 1:
4138     mode=1;
4139     break;
4140   case 2:
4141     mode=2;
4142     break;
4143   case 3:
4144     if(fTrackingPhase.Contains("Clusters2Tracks"))
4145       { mode=3; } else { mode=1; }
4146     break;
4147   case 4:
4148     if(fTrackingPhase.Contains("Clusters2Tracks"))
4149       { mode=3; } else { mode=2; }
4150     break;
4151   default:
4152     mode=0;
4153     break;
4154   }
4155   if(fTrackingPhase.Contains("Default")) mode=0;
4156
4157   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4158
4159   Double_t r=fgLayers[layerindex].GetR();
4160   Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4161
4162   Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4163   Double_t xToGo;
4164   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4165
4166   Int_t index=6*fCurrentEsdTrack+layerindex;
4167
4168
4169   Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4170   Int_t nsteps=1;
4171
4172   // back before material (no correction)
4173   Double_t rOld,xOld;
4174   rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4175   if (!t->GetLocalXat(rOld,xOld)) return 0;
4176   if (!t->Propagate(xOld)) return 0;
4177
4178   switch(mode) {
4179   case 0:
4180     xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4181     lengthTimesMeanDensity = xOverX0*x0;
4182     lengthTimesMeanDensity *= dir;
4183     // Bring the track beyond the material
4184     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4185     break;
4186   case 1:
4187     nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4188     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4189     break;
4190   case 2:
4191     if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4192     xOverX0 = fxOverX0Layer[layerindex];
4193     lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4194     lengthTimesMeanDensity *= dir;
4195     // Bring the track beyond the material
4196     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4197     break;
4198   case 3:
4199     if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4200     if(fxOverX0LayerTrks[index]<0) {
4201       nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4202       if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4203       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4204                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4205       fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4206       fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4207       return 1;
4208     }
4209     xOverX0 = fxOverX0LayerTrks[index];
4210     lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4211     lengthTimesMeanDensity *= dir;
4212     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4213     break;
4214   }
4215
4216
4217   return 1;
4218 }
4219 //------------------------------------------------------------------------
4220 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4221   //-----------------------------------------------------------------
4222   // Initialize LUT for storing material for each prolonged track
4223   //-----------------------------------------------------------------
4224   fxOverX0PipeTrks = new Float_t[ntracks];
4225   fxTimesRhoPipeTrks = new Float_t[ntracks];
4226   fxOverX0ShieldTrks = new Float_t[ntracks*2];
4227   fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4228   fxOverX0LayerTrks = new Float_t[ntracks*6];
4229   fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4230
4231   for(Int_t i=0; i<ntracks; i++) {
4232     fxOverX0PipeTrks[i] = -1.;
4233     fxTimesRhoPipeTrks[i] = -1.;
4234   }
4235   for(Int_t j=0; j<ntracks*2; j++) {
4236     fxOverX0ShieldTrks[j] = -1.;
4237     fxTimesRhoShieldTrks[j] = -1.;
4238   }
4239   for(Int_t k=0; k<ntracks*6; k++) {
4240     fxOverX0LayerTrks[k] = -1.;
4241     fxTimesRhoLayerTrks[k] = -1.;
4242   }
4243
4244   fNtracks = ntracks;
4245
4246   return;
4247 }
4248 //------------------------------------------------------------------------
4249 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4250   //-----------------------------------------------------------------
4251   // Delete LUT for storing material for each prolonged track
4252   //-----------------------------------------------------------------
4253   if(fxOverX0PipeTrks) {
4254     delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4255   }
4256   if(fxOverX0ShieldTrks) {
4257     delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4258   }
4259
4260   if(fxOverX0LayerTrks) {
4261     delete [] fxOverX0LayerTrks;  fxOverX0LayerTrks = 0;
4262   }
4263   if(fxTimesRhoPipeTrks) {
4264     delete [] fxTimesRhoPipeTrks;  fxTimesRhoPipeTrks = 0;
4265   }
4266   if(fxTimesRhoShieldTrks) {
4267     delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4268   }
4269   if(fxTimesRhoLayerTrks) {
4270     delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4271   }
4272   return;
4273 }
4274 //------------------------------------------------------------------------
4275 void AliITStrackerMI::SetForceSkippingOfLayer() {
4276   //-----------------------------------------------------------------
4277   // Check if we are forced to skip layers
4278   // either we set to skip them in RecoParam
4279   // or they were off during data-taking
4280   //-----------------------------------------------------------------
4281
4282   const AliEventInfo *eventInfo = GetEventInfo();
4283
4284   for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4285     fForceSkippingOfLayer[l] = 0;
4286     // check reco param
4287     if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4288     // check run info
4289
4290     if(eventInfo) {
4291       AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4292       if(l==0 || l==1)  {
4293         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4294       } else if(l==2 || l==3) {
4295         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4296       } else {
4297         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4298       }
4299     }
4300   }
4301   return;
4302 }
4303 //------------------------------------------------------------------------
4304 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4305                                       Int_t ilayer,Int_t idet) const {
4306   //-----------------------------------------------------------------
4307   // This method is used to decide whether to allow a prolongation
4308   // without clusters, because we want to skip the layer.
4309   // In this case the return value is > 0:
4310   // return 1: the user requested to skip a layer
4311   // return 2: track outside z acceptance
4312   //-----------------------------------------------------------------
4313
4314   if (ForceSkippingOfLayer(ilayer)) return 1;
4315
4316   Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4317
4318   if (idet<0 &&  // out in z
4319       ilayer>innerLayCanSkip &&
4320       AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4321     // check if track will cross SPD outer layer
4322     Double_t phiAtSPD2,zAtSPD2;
4323     if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4324       if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4325     }
4326     return 2; // always allow skipping, changed on 05.11.2009
4327   }
4328
4329   return 0;
4330 }
4331 //------------------------------------------------------------------------
4332 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4333                                      Int_t ilayer,Int_t idet,
4334                                      Double_t dz,Double_t dy,
4335                                      Bool_t noClusters) const {
4336   //-----------------------------------------------------------------
4337   // This method is used to decide whether to allow a prolongation
4338   // without clusters, because there is a dead zone in the road.
4339   // In this case the return value is > 0:
4340   // return 1: dead zone at z=0,+-7cm in SPD
4341   // return 2: all road is "bad" (dead or noisy) from the OCDB
4342   // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4343   // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4344   //-----------------------------------------------------------------
4345
4346   // check dead zones at z=0,+-7cm in the SPD
4347   if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4348     Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4349                           fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4350                           fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4351     Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4352                           fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4353                           fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4354     for (Int_t i=0; i<3; i++)
4355       if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4356         AliDebug(2,Form("crack SPD %d",ilayer));
4357         return 1;
4358       }
4359   }
4360
4361   // check bad zones from OCDB
4362   if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4363
4364   if (idet<0) return 0;
4365
4366   AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4367
4368   Int_t detType=-1;
4369   Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4370   if (ilayer==0 || ilayer==1) {        // ----------  SPD
4371     detType = 0;
4372   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
4373     detType = 1;
4374     detSizeFactorX *= 2.;
4375   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
4376     detType = 2;
4377   }
4378   AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4379   if (detType==2) segm->SetLayer(ilayer+1);
4380   Float_t detSizeX = detSizeFactorX*segm->Dx();
4381   Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4382
4383   // check if the road overlaps with bad chips
4384   Float_t xloc,zloc;
4385   LocalModuleCoord(ilayer,idet,track,xloc,zloc);
4386   Float_t zlocmin = zloc-dz;
4387   Float_t zlocmax = zloc+dz;
4388   Float_t xlocmin = xloc-dy;
4389   Float_t xlocmax = xloc+dy;
4390   Int_t chipsInRoad[100];
4391
4392   // check if road goes out of detector
4393   Bool_t touchNeighbourDet=kFALSE;
4394   if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4395   if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4396   if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4397   if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4398   AliDebug(2,Form("layer %d det %d zmim zmax %f %f xmin xmax %f %f   %f %f",ilayer,idet,zlocmin,zlocmax,xlocmin,xlocmax,detSizeZ,detSizeX));
4399
4400   // check if this detector is bad
4401   if (det.IsBad()) {
4402     AliDebug(2,Form("lay %d  bad detector %d",ilayer,idet));
4403     if(!touchNeighbourDet) {
4404       return 2; // all detectors in road are bad
4405     } else {
4406       return 3; // at least one is bad
4407     }
4408   }
4409
4410   Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4411   AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4412   if (!nChipsInRoad) return 0;
4413
4414   Bool_t anyBad=kFALSE,anyGood=kFALSE;
4415   for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4416     if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4417     AliDebug(2,Form("  chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4418     if (det.IsChipBad(chipsInRoad[iCh])) {
4419       anyBad=kTRUE;
4420     } else {
4421       anyGood=kTRUE;
4422     }
4423   }
4424
4425   if (!anyGood) {
4426     if(!touchNeighbourDet) {
4427       AliDebug(2,"all bad in road");
4428       return 2;  // all chips in road are bad
4429     } else {
4430       return 3; // at least a bad chip in road
4431     }
4432   }
4433
4434   if (anyBad) {
4435     AliDebug(2,"at least a bad in road");
4436     return 3; // at least a bad chip in road
4437   }
4438
4439
4440   if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4441       || !noClusters) return 0;
4442
4443   // There are no clusters in road: check if there is at least
4444   // a bad SPD pixel or SDD anode or SSD strips on both sides
4445
4446   Int_t idetInITS=idet;
4447   for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4448
4449   if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4450     AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4451     return 4;
4452   }
4453   //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4454
4455   return 0;
4456 }
4457 //------------------------------------------------------------------------
4458 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4459                                        const AliITStrackMI *track,
4460                                        Float_t &xloc,Float_t &zloc) const {
4461   //-----------------------------------------------------------------
4462   // Gives position of track in local module ref. frame
4463   //-----------------------------------------------------------------
4464
4465   xloc=0.;
4466   zloc=0.;
4467
4468   if(idet<0) return kTRUE; // track out of z acceptance of layer
4469
4470   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4471
4472   Int_t lad = Int_t(idet/ndet) + 1;
4473
4474   Int_t det = idet - (lad-1)*ndet + 1;
4475
4476   Double_t xyzGlob[3],xyzLoc[3];
4477
4478   AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4479   // take into account the misalignment: xyz at real detector plane
4480   if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4481
4482   if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4483
4484   xloc = (Float_t)xyzLoc[0];
4485   zloc = (Float_t)xyzLoc[2];
4486
4487   return kTRUE;
4488 }
4489 //------------------------------------------------------------------------
4490 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4491 //
4492 // Method to be optimized further:
4493 // Aim: decide whether a track can be used for PlaneEff evaluation
4494 //      the decision is taken based on the track quality at the layer under study
4495 //      no information on the clusters on this layer has to be used
4496 //      The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4497 //      the cut is done on number of sigmas from the boundaries
4498 //
4499 //  Input: Actual track, layer [0,5] under study
4500 //  Output: none
4501 //  Return: kTRUE if this is a good track
4502 //
4503 // it will apply a pre-selection to obtain good quality tracks.
4504 // Here also  you will have the possibility to put a control on the
4505 // impact point of the track on the basic block, in order to exclude border regions
4506 // this will be done by calling a proper method of the AliITSPlaneEff class.
4507 //
4508 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4509 // return: Bool_t   -> kTRUE if usable track, kFALSE if not usable.
4510 //
4511   Int_t index[AliITSgeomTGeo::kNLayers];
4512   Int_t k;
4513   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4514   //
4515   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4516     index[k]=clusters[k];
4517   }
4518
4519   if(!fPlaneEff)
4520     {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4521   AliITSlayer &layer=fgLayers[ilayer];
4522   Double_t r=layer.GetR();
4523   AliITStrackMI tmp(*track);
4524
4525 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4526   Int_t ncl=0;
4527   for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) {
4528     AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4529                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4530     if (tmp.GetClIndex(lay)>=0) ncl++;
4531   }
4532   Bool_t nextout = kFALSE;
4533   if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4534   else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4535   Bool_t nextin = kFALSE;
4536   if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4537   else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4538   if(ncl<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4539      return kFALSE;
4540   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout)  return kFALSE;
4541   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin)   return kFALSE;
4542   if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4543  //  if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff()  && !tmp.GetConstrain()) return kFALSE;
4544
4545 // detector number
4546   Double_t phi,z;
4547   if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4548   Int_t idet=layer.FindDetectorIndex(phi,z);
4549   if(idet<0) { AliInfo(Form("cannot find detector"));
4550     return kFALSE;}
4551
4552   // here check if it has good Chi Square.
4553
4554   //propagate to the intersection with the detector plane
4555   const AliITSdetector &det=layer.GetDetector(idet);
4556   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4557
4558   Float_t locx; //
4559   Float_t locz; //
4560   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4561   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4562   if(key>fPlaneEff->Nblock()) return kFALSE;
4563   Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4564   if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4565   //***************
4566   // DEFINITION OF SEARCH ROAD FOR accepting a track
4567   //
4568   //For the time being they are hard-wired, later on from AliITSRecoParam
4569   // Double_t nsigx=AliITSRecoParam::GetNSigXFarFromBoundary();
4570   // Double_t nsigz=AliITSRecoParam::GetNSigZFarFromBoundary();
4571   Double_t nsigz=4;
4572   Double_t nsigx=4;
4573   Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
4574   Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
4575                                                 // done for RecPoints
4576
4577   // exclude tracks at boundary between detectors
4578   //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4579   Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4580   AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4581   AliDebug(2,Form("Local:    track impact x=%f, z=%f",locx,locz));
4582   AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4583
4584   if ( (locx-dx < blockXmn+boundaryWidth) ||
4585        (locx+dx > blockXmx-boundaryWidth) ||
4586        (locz-dz < blockZmn+boundaryWidth) ||
4587        (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4588   return kTRUE;
4589 }
4590 //------------------------------------------------------------------------
4591 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4592 //
4593 // This Method has to be optimized! For the time-being it uses the same criteria
4594 // as those used in the search of extra clusters for overlapping modules.
4595 //
4596 // Method Purpose: estabilish whether a track has produced a recpoint or not
4597 //                 in the layer under study (For Plane efficiency)
4598 //
4599 // inputs: AliITStrackMI* track  (pointer to a usable track)
4600 // outputs: none
4601 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4602 //               traversed by this very track. In details:
4603 //               - if a cluster can be associated to the track then call
4604 //                  AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4605 //               - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4606 //
4607   if(!fPlaneEff)
4608     {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4609   AliITSlayer &layer=fgLayers[ilayer];
4610   Double_t r=layer.GetR();
4611   AliITStrackMI tmp(*track);
4612
4613 // detector number
4614   Double_t phi,z;
4615   if (!tmp.GetPhiZat(r,phi,z)) return;
4616   Int_t idet=layer.FindDetectorIndex(phi,z);
4617
4618   if(idet<0) { AliInfo(Form("cannot find detector"));
4619     return;}
4620
4621
4622 //propagate to the intersection with the detector plane
4623   const AliITSdetector &det=layer.GetDetector(idet);
4624   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
4625
4626
4627 //***************
4628 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
4629 //
4630   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
4631                     TMath::Sqrt(tmp.GetSigmaZ2() +
4632                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4633                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4634                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
4635   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
4636                     TMath::Sqrt(tmp.GetSigmaY2() +
4637                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4638                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4639                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
4640
4641 // road in global (rphi,z) [i.e. in tracking ref. system]
4642   Double_t zmin = tmp.GetZ() - dz;
4643   Double_t zmax = tmp.GetZ() + dz;
4644   Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
4645   Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
4646
4647 // select clusters in road
4648   layer.SelectClusters(zmin,zmax,ymin,ymax);
4649
4650 // Define criteria for track-cluster association
4651   Double_t msz = tmp.GetSigmaZ2() +
4652   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4653   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4654   AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
4655   Double_t msy = tmp.GetSigmaY2() +
4656   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4657   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4658   AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
4659   if (tmp.GetConstrain()) {
4660     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
4661     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
4662   }  else {
4663     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
4664     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
4665   }
4666   msz = 1./msz; // 1/RoadZ^2
4667   msy = 1./msy; // 1/RoadY^2
4668 //
4669
4670   const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
4671   Int_t idetc=-1;
4672   Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
4673   //Double_t  tolerance=0.2;
4674   /*while ((cl=layer.GetNextCluster(clidx))!=0) {
4675     idetc = cl->GetDetectorIndex();
4676     if(idet!=idetc) continue;
4677     //Int_t ilay = cl->GetLayer();
4678
4679     if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
4680     if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
4681
4682     Double_t chi2=tmp.GetPredictedChi2(cl);
4683     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
4684   }*/
4685   Float_t locx; //
4686   Float_t locz; //
4687   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
4688 //
4689   AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
4690   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4691   if(key>fPlaneEff->Nblock()) return;
4692   Bool_t found=kFALSE;
4693   //if (ci>=0) {
4694   Double_t chi2;
4695   while ((cl=layer.GetNextCluster(clidx))!=0) {
4696     idetc = cl->GetDetectorIndex();
4697     if(idet!=idetc) continue;
4698     // here real control to see whether the cluster can be associated to the track.
4699     // cluster not associated to track
4700     if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
4701          (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy   > 1. ) continue;
4702     // calculate track-clusters chi2
4703     chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
4704                                                // in particular, the error associated to the cluster
4705     //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
4706     // chi2 cut
4707     if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
4708     found=kTRUE;
4709     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
4710    // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
4711    // track->SetExtraModule(ilayer,idetExtra);
4712   }
4713   if(!fPlaneEff->UpDatePlaneEff(found,key))
4714        AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
4715   if(fPlaneEff->GetCreateHistos()&&  AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
4716     Float_t tr[4]={99999.,99999.,9999.,9999.};    // initialize to high values
4717     Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
4718     Int_t cltype[2]={-999,-999};
4719
4720     tr[0]=locx;
4721     tr[1]=locz;
4722     tr[2]=TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
4723     tr[3]=TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
4724
4725     if (found){
4726       clu[0]=layer.GetCluster(ci)->GetDetLocalX();
4727       clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
4728       cltype[0]=layer.GetCluster(ci)->GetNy();
4729       cltype[1]=layer.GetCluster(ci)->GetNz();
4730
4731      // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
4732      //  X->50/sqrt(12)=14 micron   Z->450/sqrt(12)= 120 micron)
4733      // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
4734      // It is computed properly by calling the method
4735      // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4736      // T
4737      //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
4738       //if (tmp.PropagateTo(x,0.,0.)) {
4739         chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
4740         AliCluster c(*layer.GetCluster(ci));
4741         c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
4742         c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
4743         //if (layer.GetCluster(ci)->GetGlobalCov(cov))  // by using this, instead, you got nominal cluster errors
4744         clu[2]=TMath::Sqrt(c.GetSigmaY2());
4745         clu[3]=TMath::Sqrt(c.GetSigmaZ2());
4746       //}
4747     }
4748     fPlaneEff->FillHistos(key,found,tr,clu,cltype);
4749   }
4750 return;
4751 }
4752