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