ITS/AliITStrackerMI.cxx

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