ITS/AliITStrackerMI.cxx

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