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