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         track->SetSharedWeight(l,weight);
2900         shared+=weight;
2901         break;
2902       }
2903     }
2904   }
2905   track->SetNUsed(shared);
2906   return shared;
2907 }
2908 //------------------------------------------------------------------------
2909 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2910 {
2911   //
2912   // find first shared track
2913   //
2914   // mean  number of clusters
2915   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2916   //
2917   for (Int_t i=0;i<6;i++) overlist[i]=-1;
2918   Int_t sharedtrack=100000;
2919   Int_t tracks[24],trackindex=0;
2920   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2921   //
2922   for (Int_t icluster=0;icluster<6;icluster++){
2923     if (clusterlist[icluster]<0) continue;
2924     Int_t index = clusterlist[icluster];
2925     Int_t l=(index & 0xf0000000) >> 28;
2926     Int_t c=(index & 0x0fffffff) >> 00;
2927     // if (ny[l]<1.e-13){
2928     //   printf("problem\n");
2929     // }
2930     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2931     //if (l>3) continue;
2932     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2933     //
2934     Float_t deltan = 0;
2935     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2936     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2937       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2938     if (l<2 || l>3){
2939       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2940     }
2941     else{
2942       deltan = (cl->GetNz()-nz[l]);
2943     }
2944     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2945     //
2946     for (Int_t itrack=3;itrack>=0;itrack--){
2947       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2948       if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2949        tracks[trackindex]  = fgLayers[l].GetClusterTracks(itrack,c);
2950        trackindex++;
2951       }
2952     }
2953   }
2954   if (trackindex==0) return -1;
2955   if (trackindex==1){
2956     sharedtrack = tracks[0];
2957   }else{
2958     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2959     else{
2960       //
2961       Int_t tracks2[24], cluster[24];
2962       for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2963       Int_t index =0;
2964       //
2965       for (Int_t i=0;i<trackindex;i++){
2966         if (tracks[i]<0) continue;
2967         tracks2[index] = tracks[i];
2968         cluster[index]++;
2969         for (Int_t j=i+1;j<trackindex;j++){
2970           if (tracks[j]<0) continue;
2971           if (tracks[j]==tracks[i]){
2972             cluster[index]++;
2973             tracks[j]=-1;
2974           }
2975         }
2976         index++;
2977       }
2978       Int_t max=0;
2979       for (Int_t i=0;i<index;i++){
2980         if (cluster[index]>max) {
2981           sharedtrack=tracks2[index];
2982           max=cluster[index];
2983         }
2984       }
2985     }
2986   }
2987
2988   if (sharedtrack>=100000) return -1;
2989   //
2990   // make list of overlaps
2991   shared =0;
2992   for (Int_t icluster=0;icluster<6;icluster++){
2993     if (clusterlist[icluster]<0) continue;
2994     Int_t index = clusterlist[icluster];
2995     Int_t l=(index & 0xf0000000) >> 28;
2996     Int_t c=(index & 0x0fffffff) >> 00;
2997     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2998     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2999     if (l==0 || l==1){
3000       if (cl->GetNy()>2) continue;
3001       if (cl->GetNz()>2) continue;
3002     }
3003     if (l==4 || l==5){
3004       if (cl->GetNy()>3) continue;
3005       if (cl->GetNz()>3) continue;
3006     }
3007     //
3008     for (Int_t itrack=3;itrack>=0;itrack--){
3009       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3010       if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
3011         overlist[l]=index;
3012         shared++;
3013       }
3014     }
3015   }
3016   return sharedtrack;
3017 }
3018 //------------------------------------------------------------------------
3019 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
3020   //
3021   // try to find track hypothesys without conflicts
3022   // with minimal chi2;
3023   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
3024   Int_t entries1 = arr1->GetEntriesFast();
3025   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
3026   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
3027   Int_t entries2 = arr2->GetEntriesFast();
3028   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
3029   //
3030   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
3031   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
3032   if (track10->Pt()>0.5+track20->Pt()) return track10;
3033
3034   for (Int_t itrack=0;itrack<entries1;itrack++){
3035     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3036     UnRegisterClusterTracks(track,trackID1);
3037   }
3038   //
3039   for (Int_t itrack=0;itrack<entries2;itrack++){
3040     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3041     UnRegisterClusterTracks(track,trackID2);
3042   }
3043   Int_t index1=0;
3044   Int_t index2=0;
3045   Float_t maxconflicts=6;
3046   Double_t maxchi2 =1000.;
3047   //
3048   // get weight of hypothesys - using best hypothesy
3049   Double_t w1,w2;
3050
3051   Int_t list1[6],list2[6];
3052   AliITSRecPoint *clist1[6], *clist2[6] ;
3053   RegisterClusterTracks(track10,trackID1);
3054   RegisterClusterTracks(track20,trackID2);
3055   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
3056   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
3057   UnRegisterClusterTracks(track10,trackID1);
3058   UnRegisterClusterTracks(track20,trackID2);
3059   //
3060   // normalized chi2
3061   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
3062   Float_t nerry[6],nerrz[6];
3063   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
3064   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
3065   for (Int_t i=0;i<6;i++){
3066      if ( (erry1[i]>0) && (erry2[i]>0)) {
3067        nerry[i] = TMath::Min(erry1[i],erry2[i]);
3068        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
3069      }else{
3070        nerry[i] = TMath::Max(erry1[i],erry2[i]);
3071        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
3072      }
3073      if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
3074        chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
3075        chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
3076        ncl1++;
3077      }
3078      if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
3079        chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
3080        chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
3081        ncl2++;
3082      }
3083   }
3084   chi21/=ncl1;
3085   chi22/=ncl2;
3086   //
3087   //
3088   Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
3089   Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
3090   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
3091   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
3092   //
3093   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
3094         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
3095         +1.*track10->Pt()/(track10->Pt()+track20->Pt())
3096         );
3097   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
3098         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
3099         +1.*track20->Pt()/(track10->Pt()+track20->Pt())
3100         );
3101
3102   Double_t sumw = w1+w2;
3103   w1/=sumw;
3104   w2/=sumw;
3105   if (w1<w2*0.5) {
3106     w1 = (d2+0.5)/(d1+d2+1.);
3107     w2 = (d1+0.5)/(d1+d2+1.);
3108   }
3109   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
3110   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
3111   //
3112   // get pair of "best" hypothesys
3113   //
3114   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
3115   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
3116
3117   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3118     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3119     //if (track1->fFakeRatio>0) continue;
3120     RegisterClusterTracks(track1,trackID1);
3121     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3122       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3123
3124       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3125       //if (track2->fFakeRatio>0) continue;
3126       Float_t nskipped=0;
3127       RegisterClusterTracks(track2,trackID2);
3128       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3129       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3130       UnRegisterClusterTracks(track2,trackID2);
3131       //
3132       if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3133       if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3134       if (nskipped>0.5) continue;
3135       //
3136       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3137       if (conflict1+1<cconflict1) continue;
3138       if (conflict2+1<cconflict2) continue;
3139       Float_t conflict=0;
3140       Float_t sumchi2=0;
3141       Float_t sum=0;
3142       for (Int_t i=0;i<6;i++){
3143         //
3144         Float_t c1 =0.; // conflict coeficients
3145         Float_t c2 =0.;
3146         if (clist1[i]&&clist2[i]){
3147           Float_t deltan = 0;
3148           if (i<2 || i>3){
3149             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3150           }
3151           else{
3152             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3153           }
3154           c1  = 2./TMath::Max(3.+deltan,2.);
3155           c2  = 2./TMath::Max(3.+deltan,2.);
3156         }
3157         else{
3158           if (clist1[i]){
3159             Float_t deltan = 0;
3160             if (i<2 || i>3){
3161               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3162             }
3163             else{
3164               deltan = (clist1[i]->GetNz()-nz1[i]);
3165             }
3166             c1  = 2./TMath::Max(3.+deltan,2.);
3167             c2  = 0;
3168           }
3169
3170           if (clist2[i]){
3171             Float_t deltan = 0;
3172             if (i<2 || i>3){
3173               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3174             }
3175             else{
3176               deltan = (clist2[i]->GetNz()-nz2[i]);
3177             }
3178             c2  = 2./TMath::Max(3.+deltan,2.);
3179             c1  = 0;
3180           }
3181         }
3182         //
3183         chi21=0;chi22=0;
3184         if (TMath::Abs(track1->GetDy(i))>0.) {
3185           chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3186             (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3187           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3188           //  (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3189         }else{
3190           if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3191         }
3192         //
3193         if (TMath::Abs(track2->GetDy(i))>0.) {
3194           chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3195             (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3196           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3197           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3198         }
3199         else{
3200           if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3201         }
3202         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3203         if (chi21>0) sum+=w1;
3204         if (chi22>0) sum+=w2;
3205         conflict+=(c1+c2);
3206       }
3207       Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3208       if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3209       Double_t normchi2 = 2*conflict+sumchi2/sum;
3210       if ( normchi2 <maxchi2 ){
3211         index1 = itrack1;
3212         index2 = itrack2;
3213         maxconflicts = conflict;
3214         maxchi2 = normchi2;
3215       }
3216     }
3217     UnRegisterClusterTracks(track1,trackID1);
3218   }
3219   //
3220   //  if (maxconflicts<4 && maxchi2<th0){
3221   if (maxchi2<th0*2.){
3222     Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3223     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3224     track1->SetChi2MIP(5,maxconflicts);
3225     track1->SetChi2MIP(6,maxchi2);
3226     track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3227     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
3228     track1->SetChi2MIP(8,index1);
3229     fBestTrackIndex[trackID1] =index1;
3230     UpdateESDtrack(track1, AliESDtrack::kITSin);
3231   }
3232   else if (track10->GetChi2MIP(0)<th1){
3233     track10->SetChi2MIP(5,maxconflicts);
3234     track10->SetChi2MIP(6,maxchi2);
3235     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
3236     UpdateESDtrack(track10,AliESDtrack::kITSin);
3237   }
3238
3239   for (Int_t itrack=0;itrack<entries1;itrack++){
3240     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3241     UnRegisterClusterTracks(track,trackID1);
3242   }
3243   //
3244   for (Int_t itrack=0;itrack<entries2;itrack++){
3245     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3246     UnRegisterClusterTracks(track,trackID2);
3247   }
3248
3249   if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3250       &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3251     //  if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3252   //    &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3253     RegisterClusterTracks(track10,trackID1);
3254   }
3255   if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3256       &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3257     //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3258     //  &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3259     RegisterClusterTracks(track20,trackID2);
3260   }
3261   return track10;
3262
3263 }
3264 //------------------------------------------------------------------------
3265 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3266   //--------------------------------------------------------------------
3267   // This function marks clusters assigned to the track
3268   //--------------------------------------------------------------------
3269   AliTracker::UseClusters(t,from);
3270
3271   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3272   //if (c->GetQ()>2) c->Use();
3273   if (c->GetSigmaZ2()>0.1) c->Use();
3274   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3275   //if (c->GetQ()>2) c->Use();
3276   if (c->GetSigmaZ2()>0.1) c->Use();
3277
3278 }
3279 //------------------------------------------------------------------------
3280 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3281 {
3282   //------------------------------------------------------------------
3283   // add track to the list of hypothesys
3284   //------------------------------------------------------------------
3285
3286   if (esdindex>=fTrackHypothesys.GetEntriesFast())
3287     fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3288   //
3289   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3290   if (!array) {
3291     array = new TObjArray(10);
3292     fTrackHypothesys.AddAt(array,esdindex);
3293   }
3294   array->AddLast(track);
3295 }
3296 //------------------------------------------------------------------------
3297 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3298 {
3299   //-------------------------------------------------------------------
3300   // compress array of track hypothesys
3301   // keep only maxsize best hypothesys
3302   //-------------------------------------------------------------------
3303   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3304   if (! (fTrackHypothesys.At(esdindex)) ) return;
3305   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3306   Int_t entries = array->GetEntriesFast();
3307   //
3308   //- find preliminary besttrack as a reference
3309   Float_t minchi2=10000;
3310   Int_t maxn=0;
3311   AliITStrackMI * besttrack=0;
3312   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3313     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3314     if (!track) continue;
3315     Float_t chi2 = NormalizedChi2(track,0);
3316     //
3317     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3318     track->SetLabel(tpcLabel);
3319     CookdEdx(track);
3320     track->SetFakeRatio(1.);
3321     CookLabel(track,0.); //For comparison only
3322     //
3323     //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3324     if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3325       if (track->GetNumberOfClusters()<maxn) continue;
3326       maxn = track->GetNumberOfClusters();
3327       if (chi2<minchi2){
3328         minchi2=chi2;
3329         besttrack=track;
3330       }
3331     }
3332     else{
3333       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3334         delete array->RemoveAt(itrack);
3335       }
3336     }
3337   }
3338   if (!besttrack) return;
3339   //
3340   //
3341   //take errors of best track as a reference
3342   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3343   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3344   for (Int_t j=0;j<6;j++) {
3345     if (besttrack->GetClIndex(j)>=0){
3346       erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3347       errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3348       ny[j]   = besttrack->GetNy(j);
3349       nz[j]   = besttrack->GetNz(j);
3350     }
3351   }
3352   //
3353   // calculate normalized chi2
3354   //
3355   Float_t * chi2        = new Float_t[entries];
3356   Int_t * index         = new Int_t[entries];
3357   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3358   for (Int_t itrack=0;itrack<entries;itrack++){
3359     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3360     if (track){
3361       AliDebug(2,Form("track %d  ncls %d\n",itrack,track->GetNumberOfClusters()));
3362       track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3363       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3364         chi2[itrack] = track->GetChi2MIP(0);
3365       else{
3366         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3367           delete array->RemoveAt(itrack);
3368         }
3369       }
3370     }
3371   }
3372   //
3373   TMath::Sort(entries,chi2,index,kFALSE);
3374   besttrack = (AliITStrackMI*)array->At(index[0]);
3375   if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3376   if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3377     for (Int_t j=0;j<6;j++){
3378       if (besttrack->GetClIndex(j)>=0){
3379         erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3380         errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3381         ny[j]   = besttrack->GetNy(j);
3382         nz[j]   = besttrack->GetNz(j);
3383       }
3384     }
3385   }
3386   //
3387   // calculate one more time with updated normalized errors
3388   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3389   for (Int_t itrack=0;itrack<entries;itrack++){
3390     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3391     if (track){
3392       track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3393       AliDebug(2,Form("track %d  ncls %d\n",itrack,track->GetNumberOfClusters()));
3394       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3395         chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3396       else
3397         {
3398           if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3399             delete array->RemoveAt(itrack);
3400           }
3401         }
3402     }
3403   }
3404   entries = array->GetEntriesFast();
3405   //
3406   //
3407   if (entries>0){
3408     TObjArray * newarray = new TObjArray();
3409     TMath::Sort(entries,chi2,index,kFALSE);
3410     besttrack = (AliITStrackMI*)array->At(index[0]);
3411     if (besttrack){
3412       AliDebug(2,Form("ncls best track %d     %f   %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3413       //
3414       for (Int_t j=0;j<6;j++){
3415         if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3416           erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3417           errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3418           ny[j]   = besttrack->GetNy(j);
3419           nz[j]   = besttrack->GetNz(j);
3420         }
3421       }
3422       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3423       minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3424       Float_t minn = besttrack->GetNumberOfClusters()-3;
3425       Int_t accepted=0;
3426       for (Int_t i=0;i<entries;i++){
3427         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3428         if (!track) continue;
3429         if (accepted>maxcut) break;
3430         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3431         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3432           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3433             delete array->RemoveAt(index[i]);
3434             continue;
3435           }
3436         }
3437         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3438         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3439           if (!shortbest) accepted++;
3440           //
3441           newarray->AddLast(array->RemoveAt(index[i]));
3442           for (Int_t j=0;j<6;j++){
3443             if (nz[j]==0){
3444               erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3445               errz[j] = track->GetSigmaZ(j); errz[j]   = track->GetSigmaZ(j+6);
3446               ny[j]   = track->GetNy(j);
3447               nz[j]   = track->GetNz(j);
3448             }
3449           }
3450         }
3451         else{
3452           delete array->RemoveAt(index[i]);
3453         }
3454       }
3455       array->Delete();
3456       delete fTrackHypothesys.RemoveAt(esdindex);
3457       fTrackHypothesys.AddAt(newarray,esdindex);
3458     }
3459     else{
3460       array->Delete();
3461       delete fTrackHypothesys.RemoveAt(esdindex);
3462     }
3463   }
3464   delete [] chi2;
3465   delete [] index;
3466 }
3467 //------------------------------------------------------------------------
3468 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3469 {
3470   //-------------------------------------------------------------
3471   // try to find best hypothesy
3472   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3473   //-------------------------------------------------------------
3474   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3475   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3476   if (!array) return 0;
3477   Int_t entries = array->GetEntriesFast();
3478   if (!entries) return 0;
3479   Float_t minchi2 = 100000;
3480   AliITStrackMI * besttrack=0;
3481   //
3482   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
3483   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3484   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3485   Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3486   //
3487   for (Int_t i=0;i<entries;i++){
3488     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3489     if (!track) continue;
3490     Float_t sigmarfi,sigmaz;
3491     GetDCASigma(track,sigmarfi,sigmaz);
3492     track->SetDnorm(0,sigmarfi);
3493     track->SetDnorm(1,sigmaz);
3494     //
3495     track->SetChi2MIP(1,1000000);
3496     track->SetChi2MIP(2,1000000);
3497     track->SetChi2MIP(3,1000000);
3498     //
3499     // backtrack
3500     backtrack = new(backtrack) AliITStrackMI(*track);
3501     if (track->GetConstrain()) {
3502       if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3503       if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3504         if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;
3505       }
3506       backtrack->ResetCovariance(10.);
3507     }else{
3508       backtrack->ResetCovariance(10.);
3509     }
3510     backtrack->ResetClusters();
3511
3512     Double_t x = original->GetX();
3513     if (!RefitAt(x,backtrack,track)) continue;
3514     //
3515     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3516     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3517     if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.)  continue;
3518     track->SetChi22(GetMatchingChi2(backtrack,original));
3519
3520     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
3521     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
3522     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
3523
3524
3525     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1))  continue;
3526     //
3527     //forward track - without constraint
3528     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3529     forwardtrack->ResetClusters();
3530     x = track->GetX();
3531     RefitAt(x,forwardtrack,track);
3532     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3533     if  (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0)  continue;
3534     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2))  continue;
3535
3536     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3537     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3538     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
3539     forwardtrack->SetD(0,track->GetD(0));
3540     forwardtrack->SetD(1,track->GetD(1));
3541     {
3542       Int_t list[6];
3543       AliITSRecPoint* clist[6];
3544       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3545       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3546     }
3547
3548     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3549     if  ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3550     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3551       track->SetChi2MIP(3,1000);
3552       continue;
3553     }
3554     Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3555     //
3556     for (Int_t ichi=0;ichi<5;ichi++){
3557       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3558     }
3559     if (chi2 < minchi2){
3560       //besttrack = new AliITStrackMI(*forwardtrack);
3561       besttrack = track;
3562       besttrack->SetLabel(track->GetLabel());
3563       besttrack->SetFakeRatio(track->GetFakeRatio());
3564       minchi2   = chi2;
3565       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3566       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3567       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
3568     }
3569   }
3570   delete backtrack;
3571   delete forwardtrack;
3572
3573   if (!besttrack)  return 0;
3574
3575   Int_t accepted=0;
3576   for (Int_t i=0;i<entries;i++){
3577     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3578
3579     if (!track) continue;
3580
3581     if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3582         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3583         track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3584       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3585         delete array->RemoveAt(i);
3586         continue;
3587       }
3588     }
3589     else{
3590       accepted++;
3591     }
3592   }
3593   //
3594   array->Compress();
3595   SortTrackHypothesys(esdindex,checkmax,1);
3596
3597   array = (TObjArray*) fTrackHypothesys.At(esdindex);
3598   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3599   besttrack = (AliITStrackMI*)array->At(0);
3600   if (!besttrack)  return 0;
3601   besttrack->SetChi2MIP(8,0);
3602   fBestTrackIndex[esdindex]=0;
3603   entries = array->GetEntriesFast();
3604   AliITStrackMI *longtrack =0;
3605   minchi2 =1000;
3606   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3607   for (Int_t itrack=entries-1;itrack>0;itrack--) {
3608     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3609     if (!track->GetConstrain()) continue;
3610     if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3611     if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3612     if (track->GetChi2MIP(0)>4.) continue;
3613     minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3614     longtrack =track;
3615   }
3616   //if (longtrack) besttrack=longtrack;
3617
3618   Int_t list[6];
3619   AliITSRecPoint * clist[6];
3620   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3621   if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3622       &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3623     RegisterClusterTracks(besttrack,esdindex);
3624   }
3625   //
3626   //
3627   if (shared>0.0){
3628     Int_t nshared;
3629     Int_t overlist[6];
3630     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3631     if (sharedtrack>=0){
3632       //
3633       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3634       if (besttrack){
3635         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3636       }
3637       else return 0;
3638     }
3639   }
3640
3641   if (shared>2.5) return 0;
3642   if (shared>1.0) return besttrack;
3643   //
3644   // Don't sign clusters if not gold track
3645   //
3646   if (!besttrack->IsGoldPrimary()) return besttrack;
3647   if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
3648   //
3649   if (fConstraint[fPass]){
3650     //
3651     // sign clusters
3652     //
3653     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3654     for (Int_t i=0;i<6;i++){
3655       Int_t index = besttrack->GetClIndex(i);
3656       if (index<0) continue;
3657       Int_t ilayer =  (index & 0xf0000000) >> 28;
3658       if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3659       AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3660       if (!c) continue;
3661       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3662       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3663       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
3664       if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3665         if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3666       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
3667
3668       Bool_t cansign = kTRUE;
3669       for (Int_t itrack=0;itrack<entries; itrack++){
3670         AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3671         if (!track) continue;
3672         if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3673         if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3674           cansign = kFALSE;
3675           break;
3676         }
3677       }
3678       if (cansign){
3679         if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3680         if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3681         if (!c->IsUsed()) c->Use();
3682       }
3683     }
3684   }
3685   return besttrack;
3686 }
3687 //------------------------------------------------------------------------
3688 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3689 {
3690   //
3691   // get "best" hypothesys
3692   //
3693
3694   Int_t nentries = itsTracks.GetEntriesFast();
3695   for (Int_t i=0;i<nentries;i++){
3696     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3697     if (!track) continue;
3698     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3699     if (!array) continue;
3700     if (array->GetEntriesFast()<=0) continue;
3701     //
3702     AliITStrackMI* longtrack=0;
3703     Float_t minn=0;
3704     Float_t maxchi2=1000;
3705     for (Int_t j=0;j<array->GetEntriesFast();j++){
3706       AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3707       if (!trackHyp) continue;
3708       if (trackHyp->GetGoldV0()) {
3709         longtrack = trackHyp;   //gold V0 track taken
3710         break;
3711       }
3712       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3713       Float_t chi2 = trackHyp->GetChi2MIP(0);
3714       if (fAfterV0){
3715         if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3716       }
3717       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3718       //
3719       if (chi2 > maxchi2) continue;
3720       minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3721       maxchi2 = chi2;
3722       longtrack=trackHyp;
3723     }
3724     //
3725     //
3726     //
3727     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3728     if (!longtrack) {longtrack = besttrack;}
3729     else besttrack= longtrack;
3730     //
3731     if (besttrack) {
3732       Int_t list[6];
3733       AliITSRecPoint * clist[6];
3734       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3735       //
3736       track->SetNUsed(shared);
3737       track->SetNSkipped(besttrack->GetNSkipped());
3738       track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3739       if (shared>0) {
3740         if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3741         Int_t nshared;
3742         Int_t overlist[6];
3743         //
3744         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3745         //if (sharedtrack==-1) sharedtrack=0;
3746         if (sharedtrack>=0) {
3747           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3748         }
3749       }
3750       if (besttrack&&fAfterV0) {
3751         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3752       }
3753       if (besttrack) {
3754         if (fConstraint[fPass]) UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3755         if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3756           if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3757                TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3758         }
3759       }
3760     }
3761   }
3762 }
3763 //------------------------------------------------------------------------
3764 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3765   //--------------------------------------------------------------------
3766   //This function "cooks" a track label. If label<0, this track is fake.
3767   //--------------------------------------------------------------------
3768   Int_t tpcLabel=-1;
3769
3770   if (track->GetESDtrack()){
3771     tpcLabel = track->GetESDtrack()->GetTPCLabel();
3772     ULong_t trStatus=track->GetESDtrack()->GetStatus();
3773     if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
3774   }
3775    track->SetChi2MIP(9,0);
3776    Int_t nwrong=0;
3777    for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3778      Int_t cindex = track->GetClusterIndex(i);
3779      Int_t l=(cindex & 0xf0000000) >> 28;
3780      AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3781      Int_t isWrong=1;
3782      for (Int_t ind=0;ind<3;ind++){
3783        if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
3784        //AliDebug(2,Form("icl %d  ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3785      }
3786      track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3787      nwrong+=isWrong;
3788    }
3789    Int_t nclusters = track->GetNumberOfClusters();
3790    if (nclusters > 0) //PH Some tracks don't have any cluster
3791      track->SetFakeRatio(double(nwrong)/double(nclusters));
3792    if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
3793      track->SetLabel(-tpcLabel);
3794    } else {
3795      track->SetLabel(tpcLabel);
3796    }
3797    AliDebug(2,Form(" nls %d wrong %d  label %d  tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
3798
3799 }
3800 //------------------------------------------------------------------------
3801 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
3802   //
3803   // Fill the dE/dx in this track
3804   //
3805   track->SetChi2MIP(9,0);
3806   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3807     Int_t cindex = track->GetClusterIndex(i);
3808     Int_t l=(cindex & 0xf0000000) >> 28;
3809     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3810     Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3811     Int_t isWrong=1;
3812     for (Int_t ind=0;ind<3;ind++){
3813       if (cl->GetLabel(ind)==lab) isWrong=0;
3814     }
3815     track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3816   }
3817   Double_t low=0.;
3818   Double_t up=0.51;
3819   track->CookdEdx(low,up);
3820 }
3821 //------------------------------------------------------------------------
3822 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3823   //
3824   // Create some arrays
3825   //
3826   if (fCoefficients) delete []fCoefficients;
3827   fCoefficients = new Float_t[ntracks*48];
3828   for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3829 }
3830 //------------------------------------------------------------------------
3831 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3832 {
3833   //
3834   // Compute predicted chi2
3835   //
3836   // Take into account the mis-alignment (bring track to cluster plane)
3837   Double_t xTrOrig=track->GetX();
3838   if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
3839   Float_t erry,errz,covyz;
3840   Float_t theta = track->GetTgl();
3841   Float_t phi   = track->GetSnp();
3842   phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
3843   AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
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   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()));
3846   Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
3847   // Bring the track back to detector plane in ideal geometry
3848   // [mis-alignment will be accounted for in UpdateMI()]
3849   if (!track->Propagate(xTrOrig)) return 1000.;
3850   Float_t ny,nz;
3851   AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3852   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3853   if (delta>1){
3854     chi2+=0.5*TMath::Min(delta/2,2.);
3855     chi2+=2.*cluster->GetDeltaProbability();
3856   }
3857   //
3858   track->SetNy(layer,ny);
3859   track->SetNz(layer,nz);
3860   track->SetSigmaY(layer,erry);
3861   track->SetSigmaZ(layer, errz);
3862   track->SetSigmaYZ(layer,covyz);
3863   //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3864   track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
3865   return chi2;
3866
3867 }
3868 //------------------------------------------------------------------------
3869 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3870 {
3871   //
3872   // Update ITS track
3873   //
3874   Int_t layer = (index & 0xf0000000) >> 28;
3875   track->SetClIndex(layer, index);
3876   if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3877     if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3878       chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3879       track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3880     }
3881   }
3882
3883   if (TMath::Abs(cl->GetQ())<1.e-13) return 0;  // ingore the "virtual" clusters
3884
3885
3886   // Take into account the mis-alignment (bring track to cluster plane)
3887   Double_t xTrOrig=track->GetX();
3888   Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3889   AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3890   AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3891   AliDebug(2,Form(" xtr %f  xcl %f",track->GetX(),cl->GetX()));
3892
3893   if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
3894
3895   AliCluster c(*cl);
3896   c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3897   c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3898   c.SetSigmaYZ(track->GetSigmaYZ(layer));
3899
3900
3901   Int_t updated = track->UpdateMI(&c,chi2,index);
3902   // Bring the track back to detector plane in ideal geometry
3903   if (!track->Propagate(xTrOrig)) return 0;
3904
3905   if(!updated) AliDebug(2,"update failed");
3906   return updated;
3907 }
3908
3909 //------------------------------------------------------------------------
3910 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3911 {
3912   //
3913   //DCA sigmas parameterization
3914   //to be paramterized using external parameters in future
3915   //
3916   //
3917   Double_t curv=track->GetC();
3918   sigmarfi = 0.0040+1.4 *TMath::Abs(curv)+332.*curv*curv;
3919   sigmaz   = 0.0110+4.37*TMath::Abs(curv);
3920 }
3921 //------------------------------------------------------------------------
3922 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
3923 {
3924   //
3925   // Clusters from delta electrons?
3926   //
3927   Int_t entries = clusterArray->GetEntriesFast();
3928   if (entries<4) return;
3929   AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
3930   Int_t layer = cluster->GetLayer();
3931   if (layer>1) return;
3932   Int_t index[10000];
3933   Int_t ncandidates=0;
3934   Float_t r = (layer>0)? 7:4;
3935   //
3936   for (Int_t i=0;i<entries;i++){
3937     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
3938     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3939     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3940     index[ncandidates] = i;  //candidate to belong to delta electron track
3941     ncandidates++;
3942     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3943       cl0->SetDeltaProbability(1);
3944     }
3945   }
3946   //
3947   //
3948   //
3949   for (Int_t i=0;i<ncandidates;i++){
3950     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
3951     if (cl0->GetDeltaProbability()>0.8) continue;
3952     //
3953     Int_t ncl = 0;
3954     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3955     sumy=sumz=sumy2=sumyz=sumw=0.0;
3956     for (Int_t j=0;j<ncandidates;j++){
3957       if (i==j) continue;
3958       AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
3959       //
3960       Float_t dz = cl0->GetZ()-cl1->GetZ();
3961       Float_t dy = cl0->GetY()-cl1->GetY();
3962       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3963         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3964         y[ncl] = cl1->GetY();
3965         z[ncl] = cl1->GetZ();
3966         sumy+= y[ncl]*weight;
3967         sumz+= z[ncl]*weight;
3968         sumy2+=y[ncl]*y[ncl]*weight;
3969         sumyz+=y[ncl]*z[ncl]*weight;
3970         sumw+=weight;
3971         ncl++;
3972       }
3973     }
3974     if (ncl<4) continue;
3975     Float_t det = sumw*sumy2  - sumy*sumy;
3976     Float_t delta=1000;
3977     if (TMath::Abs(det)>0.01){
3978       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
3979       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
3980       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3981     }
3982     else{
3983       Float_t z0  = sumyz/sumy;
3984       delta = TMath::Abs(cl0->GetZ()-z0);
3985     }
3986     if ( delta<0.05) {
3987       cl0->SetDeltaProbability(1-20.*delta);
3988     }
3989   }
3990 }
3991 //------------------------------------------------------------------------
3992 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3993 {
3994   //
3995   // Update ESD track
3996   //
3997   track->UpdateESDtrack(flags);
3998   AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3999   if (oldtrack) delete oldtrack;
4000   track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
4001   // if (TMath::Abs(track->GetDnorm(1))<0.000000001){
4002   //   printf("Problem\n");
4003   // }
4004 }
4005 //------------------------------------------------------------------------
4006 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
4007   //
4008   // Get nearest upper layer close to the point xr.
4009   // rough approximation
4010   //
4011   const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
4012   Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
4013   Int_t res =6;
4014   for (Int_t i=0;i<6;i++){
4015     if (radius<kRadiuses[i]){
4016       res =i;
4017       break;
4018     }
4019   }
4020   return res;
4021 }
4022 //------------------------------------------------------------------------
4023 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4024   //--------------------------------------------------------------------
4025   // Fill a look-up table with mean material
4026   //--------------------------------------------------------------------
4027
4028   Int_t n=1000;
4029   Double_t mparam[7];
4030   Double_t point1[3],point2[3];
4031   Double_t phi,cosphi,sinphi,z;
4032   // 0-5 layers, 6 pipe, 7-8 shields
4033   Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4034   Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4035
4036   Int_t ifirst=0,ilast=0;
4037   if(material.Contains("Pipe")) {
4038     ifirst=6; ilast=6;
4039   } else if(material.Contains("Shields")) {
4040     ifirst=7; ilast=8;
4041   } else if(material.Contains("Layers")) {
4042     ifirst=0; ilast=5;
4043   } else {
4044     Error("BuildMaterialLUT","Wrong layer name\n");
4045   }
4046
4047   for(Int_t imat=ifirst; imat<=ilast; imat++) {
4048     Double_t param[5]={0.,0.,0.,0.,0.};
4049     for (Int_t i=0; i<n; i++) {
4050       phi = 2.*TMath::Pi()*gRandom->Rndm();
4051       cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4052       z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4053       point1[0] = rmin[imat]*cosphi;
4054       point1[1] = rmin[imat]*sinphi;
4055       point1[2] = z;
4056       point2[0] = rmax[imat]*cosphi;
4057       point2[1] = rmax[imat]*sinphi;
4058       point2[2] = z;
4059       AliTracker::MeanMaterialBudget(point1,point2,mparam);
4060       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4061     }
4062     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4063     if(imat<=5) {
4064       fxOverX0Layer[imat] = param[1];
4065       fxTimesRhoLayer[imat] = param[0]*param[4];
4066     } else if(imat==6) {
4067       fxOverX0Pipe = param[1];
4068       fxTimesRhoPipe = param[0]*param[4];
4069     } else if(imat==7) {
4070       fxOverX0Shield[0] = param[1];
4071       fxTimesRhoShield[0] = param[0]*param[4];
4072     } else if(imat==8) {
4073       fxOverX0Shield[1] = param[1];
4074       fxTimesRhoShield[1] = param[0]*param[4];
4075     }
4076   }
4077   /*
4078   printf("%s\n",material.Data());
4079   printf("%f  %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4080   printf("%f  %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4081   printf("%f  %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4082   printf("%f  %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4083   printf("%f  %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4084   printf("%f  %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4085   printf("%f  %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4086   printf("%f  %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4087   printf("%f  %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4088   */
4089   return;
4090 }
4091 //------------------------------------------------------------------------
4092 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4093                                               TString direction) {
4094   //-------------------------------------------------------------------
4095   // Propagate beyond beam pipe and correct for material
4096   // (material budget in different ways according to fUseTGeo value)
4097   // Add time if going outward (PropagateTo or PropagateToTGeo)
4098   //-------------------------------------------------------------------
4099
4100   // Define budget mode:
4101   // 0: material from AliITSRecoParam (hard coded)
4102   // 1: material from TGeo in one step (on the fly)
4103   // 2: material from lut
4104   // 3: material from TGeo in one step (same for all hypotheses)
4105   Int_t mode;
4106   switch(fUseTGeo) {
4107   case 0:
4108     mode=0;
4109     break;
4110   case 1:
4111     mode=1;
4112     break;
4113   case 2:
4114     mode=2;
4115     break;
4116   case 3:
4117     if(fTrackingPhase.Contains("Clusters2Tracks"))
4118       { mode=3; } else { mode=1; }
4119     break;
4120   case 4:
4121     if(fTrackingPhase.Contains("Clusters2Tracks"))
4122       { mode=3; } else { mode=2; }
4123     break;
4124   default:
4125     mode=0;
4126     break;
4127   }
4128   if(fTrackingPhase.Contains("Default")) mode=0;
4129
4130   Int_t index=fCurrentEsdTrack;
4131
4132   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4133   Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4134   Double_t xToGo;
4135   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4136
4137   Double_t xOverX0,x0,lengthTimesMeanDensity;
4138
4139   switch(mode) {
4140   case 0:
4141     xOverX0 = AliITSRecoParam::GetdPipe();
4142     x0 = AliITSRecoParam::GetX0Be();
4143     lengthTimesMeanDensity = xOverX0*x0;
4144     lengthTimesMeanDensity *= dir;
4145     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4146     break;
4147   case 1:
4148     if (!t->PropagateToTGeo(xToGo,1)) return 0;
4149     break;
4150   case 2:
4151     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4152     xOverX0 = fxOverX0Pipe;
4153     lengthTimesMeanDensity = fxTimesRhoPipe;
4154     lengthTimesMeanDensity *= dir;
4155     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4156     break;
4157   case 3:
4158     if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4159     if(fxOverX0PipeTrks[index]<0) {
4160       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4161       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4162                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4163       fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4164       fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4165       return 1;
4166     }
4167     xOverX0 = fxOverX0PipeTrks[index];
4168     lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4169     lengthTimesMeanDensity *= dir;
4170     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4171     break;
4172   }
4173
4174   return 1;
4175 }
4176 //------------------------------------------------------------------------
4177 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4178                                                 TString shield,
4179                                                 TString direction) {
4180   //-------------------------------------------------------------------
4181   // Propagate beyond SPD or SDD shield and correct for material
4182   // (material budget in different ways according to fUseTGeo value)
4183   // Add time if going outward (PropagateTo or PropagateToTGeo)
4184   //-------------------------------------------------------------------
4185
4186   // Define budget mode:
4187   // 0: material from AliITSRecoParam (hard coded)
4188   // 1: material from TGeo in steps of X cm (on the fly)
4189   //    X = AliITSRecoParam::GetStepSizeTGeo()
4190   // 2: material from lut
4191   // 3: material from TGeo in one step (same for all hypotheses)
4192   Int_t mode;
4193   switch(fUseTGeo) {
4194   case 0:
4195     mode=0;
4196     break;
4197   case 1:
4198     mode=1;
4199     break;
4200   case 2:
4201     mode=2;
4202     break;
4203   case 3:
4204     if(fTrackingPhase.Contains("Clusters2Tracks"))
4205       { mode=3; } else { mode=1; }
4206     break;
4207   case 4:
4208     if(fTrackingPhase.Contains("Clusters2Tracks"))
4209       { mode=3; } else { mode=2; }
4210     break;
4211   default:
4212     mode=0;
4213     break;
4214   }
4215   if(fTrackingPhase.Contains("Default")) mode=0;
4216
4217   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4218   Double_t rToGo;
4219   Int_t    shieldindex=0;
4220   if (shield.Contains("SDD")) { // SDDouter
4221     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4222     shieldindex=1;
4223   } else if (shield.Contains("SPD")) {        // SPDouter
4224     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4225     shieldindex=0;
4226   } else {
4227     Error("CorrectForShieldMaterial"," Wrong shield name\n");
4228     return 0;
4229   }
4230
4231   // do nothing if we are already beyond the shield
4232   Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4233   if(dir<0 && rTrack > rToGo) return 1; // going outward
4234   if(dir>0 && rTrack < rToGo) return 1; // going inward
4235
4236
4237   Double_t xToGo;
4238   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4239
4240   Int_t index=2*fCurrentEsdTrack+shieldindex;
4241
4242   Double_t xOverX0,x0,lengthTimesMeanDensity;
4243   Int_t nsteps=1;
4244
4245   switch(mode) {
4246   case 0:
4247     xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4248     x0 = AliITSRecoParam::GetX0shield(shieldindex);
4249     lengthTimesMeanDensity = xOverX0*x0;
4250     lengthTimesMeanDensity *= dir;
4251     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4252     break;
4253   case 1:
4254     nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4255     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4256     break;
4257   case 2:
4258     if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4259     xOverX0 = fxOverX0Shield[shieldindex];
4260     lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4261     lengthTimesMeanDensity *= dir;
4262     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4263     break;
4264   case 3:
4265     if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4266     if(fxOverX0ShieldTrks[index]<0) {
4267       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4268       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4269                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4270       fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4271       fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4272       return 1;
4273     }
4274     xOverX0 = fxOverX0ShieldTrks[index];
4275     lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4276     lengthTimesMeanDensity *= dir;
4277     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4278     break;
4279   }
4280
4281   return 1;
4282 }
4283 //------------------------------------------------------------------------
4284 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4285                                                Int_t layerindex,
4286                                                Double_t oldGlobXYZ[3],
4287                                                TString direction) {
4288   //-------------------------------------------------------------------
4289   // Propagate beyond layer and correct for material
4290   // (material budget in different ways according to fUseTGeo value)
4291   // Add time if going outward (PropagateTo or PropagateToTGeo)
4292   //-------------------------------------------------------------------
4293
4294   // Define budget mode:
4295   // 0: material from AliITSRecoParam (hard coded)
4296   // 1: material from TGeo in stepsof X cm (on the fly)
4297   //    X = AliITSRecoParam::GetStepSizeTGeo()
4298   // 2: material from lut
4299   // 3: material from TGeo in one step (same for all hypotheses)
4300   Int_t mode;
4301   switch(fUseTGeo) {
4302   case 0:
4303     mode=0;
4304     break;
4305   case 1:
4306     mode=1;
4307     break;
4308   case 2:
4309     mode=2;
4310     break;
4311   case 3:
4312     if(fTrackingPhase.Contains("Clusters2Tracks"))
4313       { mode=3; } else { mode=1; }
4314     break;
4315   case 4:
4316     if(fTrackingPhase.Contains("Clusters2Tracks"))
4317       { mode=3; } else { mode=2; }
4318     break;
4319   default:
4320     mode=0;
4321     break;
4322   }
4323   if(fTrackingPhase.Contains("Default")) mode=0;
4324
4325   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4326
4327   Double_t r=fgLayers[layerindex].GetR();
4328   Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4329
4330   Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4331   Double_t xToGo;
4332   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4333
4334   Int_t index=6*fCurrentEsdTrack+layerindex;
4335
4336
4337   Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4338   Int_t nsteps=1;
4339
4340   // back before material (no correction)
4341   Double_t rOld,xOld;
4342   rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4343   if (!t->GetLocalXat(rOld,xOld)) return 0;
4344   if (!t->Propagate(xOld)) return 0;
4345
4346   switch(mode) {
4347   case 0:
4348     xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4349     lengthTimesMeanDensity = xOverX0*x0;
4350     lengthTimesMeanDensity *= dir;
4351     // Bring the track beyond the material
4352     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4353     break;
4354   case 1:
4355     nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4356     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4357     break;
4358   case 2:
4359     if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4360     xOverX0 = fxOverX0Layer[layerindex];
4361     lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4362     lengthTimesMeanDensity *= dir;
4363     // Bring the track beyond the material
4364     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4365     break;
4366   case 3:
4367     if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4368     if(fxOverX0LayerTrks[index]<0) {
4369       nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4370       if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4371       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4372                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4373       fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4374       fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4375       return 1;
4376     }
4377     xOverX0 = fxOverX0LayerTrks[index];
4378     lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4379     lengthTimesMeanDensity *= dir;
4380     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4381     break;
4382   }
4383
4384
4385   return 1;
4386 }
4387 //------------------------------------------------------------------------
4388 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4389   //-----------------------------------------------------------------
4390   // Initialize LUT for storing material for each prolonged track
4391   //-----------------------------------------------------------------
4392   fxOverX0PipeTrks = new Float_t[ntracks];
4393   fxTimesRhoPipeTrks = new Float_t[ntracks];
4394   fxOverX0ShieldTrks = new Float_t[ntracks*2];
4395   fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4396   fxOverX0LayerTrks = new Float_t[ntracks*6];
4397   fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4398
4399   for(Int_t i=0; i<ntracks; i++) {
4400     fxOverX0PipeTrks[i] = -1.;
4401     fxTimesRhoPipeTrks[i] = -1.;
4402   }
4403   for(Int_t j=0; j<ntracks*2; j++) {
4404     fxOverX0ShieldTrks[j] = -1.;
4405     fxTimesRhoShieldTrks[j] = -1.;
4406   }
4407   for(Int_t k=0; k<ntracks*6; k++) {
4408     fxOverX0LayerTrks[k] = -1.;
4409     fxTimesRhoLayerTrks[k] = -1.;
4410   }
4411
4412   fNtracks = ntracks;
4413
4414   return;
4415 }
4416 //------------------------------------------------------------------------
4417 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4418   //-----------------------------------------------------------------
4419   // Delete LUT for storing material for each prolonged track
4420   //-----------------------------------------------------------------
4421   if(fxOverX0PipeTrks) {
4422     delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4423   }
4424   if(fxOverX0ShieldTrks) {
4425     delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4426   }
4427
4428   if(fxOverX0LayerTrks) {
4429     delete [] fxOverX0LayerTrks;  fxOverX0LayerTrks = 0;
4430   }
4431   if(fxTimesRhoPipeTrks) {
4432     delete [] fxTimesRhoPipeTrks;  fxTimesRhoPipeTrks = 0;
4433   }
4434   if(fxTimesRhoShieldTrks) {
4435     delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4436   }
4437   if(fxTimesRhoLayerTrks) {
4438     delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4439   }
4440   return;
4441 }
4442 //------------------------------------------------------------------------
4443 void AliITStrackerMI::SetForceSkippingOfLayer() {
4444   //-----------------------------------------------------------------
4445   // Check if we are forced to skip layers
4446   // either we set to skip them in RecoParam
4447   // or they were off during data-taking
4448   //-----------------------------------------------------------------
4449
4450   const AliEventInfo *eventInfo = GetEventInfo();
4451
4452   for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4453     fForceSkippingOfLayer[l] = 0;
4454     // check reco param
4455     if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4456     // check run info
4457
4458     if(eventInfo &&
4459        AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4460       AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4461       if(l==0 || l==1)  {
4462         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4463       } else if(l==2 || l==3) {
4464         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4465       } else {
4466         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4467       }
4468     }
4469   }
4470   return;
4471 }
4472 //------------------------------------------------------------------------
4473 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4474                                       Int_t ilayer,Int_t idet) const {
4475   //-----------------------------------------------------------------
4476   // This method is used to decide whether to allow a prolongation
4477   // without clusters, because we want to skip the layer.
4478   // In this case the return value is > 0:
4479   // return 1: the user requested to skip a layer
4480   // return 2: track outside z acceptance
4481   //-----------------------------------------------------------------
4482
4483   if (ForceSkippingOfLayer(ilayer)) return 1;
4484
4485   Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4486
4487   if (idet<0 &&  // out in z
4488       ilayer>innerLayCanSkip &&
4489       AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4490     // check if track will cross SPD outer layer
4491     Double_t phiAtSPD2,zAtSPD2;
4492     if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4493       if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4494     }
4495     return 2; // always allow skipping, changed on 05.11.2009
4496   }
4497
4498   return 0;
4499 }
4500 //------------------------------------------------------------------------
4501 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4502                                      Int_t ilayer,Int_t idet,
4503                                      Double_t dz,Double_t dy,
4504                                      Bool_t noClusters) const {
4505   //-----------------------------------------------------------------
4506   // This method is used to decide whether to allow a prolongation
4507   // without clusters, because there is a dead zone in the road.
4508   // In this case the return value is > 0:
4509   // return 1: dead zone at z=0,+-7cm in SPD
4510   //     This method assumes that fSPDdetzcentre is ordered from -z to +z
4511   // return 2: all road is "bad" (dead or noisy) from the OCDB
4512   // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4513   // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4514   //-----------------------------------------------------------------
4515
4516   // check dead zones at z=0,+-7cm in the SPD
4517   if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4518     Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4519                           fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4520                           fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4521     Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4522                           fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4523                           fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4524     for (Int_t i=0; i<3; i++)
4525       if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4526         AliDebug(2,Form("crack SPD %d track z %f   %f   %f  %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4527         if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4528       }
4529   }
4530
4531   // check bad zones from OCDB
4532   if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4533
4534   if (idet<0) return 0;
4535
4536   AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4537
4538   Int_t detType=-1;
4539   Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4540   if (ilayer==0 || ilayer==1) {        // ----------  SPD
4541     detType = 0;
4542   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
4543     detType = 1;
4544     detSizeFactorX *= 2.;
4545   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
4546     detType = 2;
4547   }
4548   AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4549   if (detType==2) segm->SetLayer(ilayer+1);
4550   Float_t detSizeX = detSizeFactorX*segm->Dx();
4551   Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4552
4553   // check if the road overlaps with bad chips
4554   Float_t xloc,zloc;
4555   if(!(LocalModuleCoord(ilayer,idet,track,xloc,zloc)))return 0;
4556   Float_t zlocmin = zloc-dz;
4557   Float_t zlocmax = zloc+dz;
4558   Float_t xlocmin = xloc-dy;
4559   Float_t xlocmax = xloc+dy;
4560   Int_t chipsInRoad[100];
4561
4562   // check if road goes out of detector
4563   Bool_t touchNeighbourDet=kFALSE;
4564   if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4565   if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4566   if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4567   if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4568   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));
4569
4570   // check if this detector is bad
4571   if (det.IsBad()) {
4572     AliDebug(2,Form("lay %d  bad detector %d",ilayer,idet));
4573     if(!touchNeighbourDet) {
4574       return 2; // all detectors in road are bad
4575     } else {
4576       return 3; // at least one is bad
4577     }
4578   }
4579
4580   Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4581   AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4582   if (!nChipsInRoad) return 0;
4583
4584   Bool_t anyBad=kFALSE,anyGood=kFALSE;
4585   for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4586     if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4587     AliDebug(2,Form("  chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4588     if (det.IsChipBad(chipsInRoad[iCh])) {
4589       anyBad=kTRUE;
4590     } else {
4591       anyGood=kTRUE;
4592     }
4593   }
4594
4595   if (!anyGood) {
4596     if(!touchNeighbourDet) {
4597       AliDebug(2,"all bad in road");
4598       return 2;  // all chips in road are bad
4599     } else {
4600       return 3; // at least a bad chip in road
4601     }
4602   }
4603
4604   if (anyBad) {
4605     AliDebug(2,"at least a bad in road");
4606     return 3; // at least a bad chip in road
4607   }
4608
4609
4610   if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4611       || !noClusters) return 0;
4612
4613   // There are no clusters in road: check if there is at least
4614   // a bad SPD pixel or SDD anode or SSD strips on both sides
4615
4616   Int_t idetInITS=idet;
4617   for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4618
4619   if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4620     AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4621     return 4;
4622   }
4623   //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4624
4625   return 0;
4626 }
4627 //------------------------------------------------------------------------
4628 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4629                                        const AliITStrackMI *track,
4630                                        Float_t &xloc,Float_t &zloc) const {
4631   //-----------------------------------------------------------------
4632   // Gives position of track in local module ref. frame
4633   //-----------------------------------------------------------------
4634
4635   xloc=0.;
4636   zloc=0.;
4637
4638   if(idet<0) return kTRUE; // track out of z acceptance of layer
4639
4640   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4641
4642   Int_t lad = Int_t(idet/ndet) + 1;
4643
4644   Int_t det = idet - (lad-1)*ndet + 1;
4645
4646   Double_t xyzGlob[3],xyzLoc[3];
4647
4648   AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4649   // take into account the misalignment: xyz at real detector plane
4650   if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4651
4652   if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4653
4654   xloc = (Float_t)xyzLoc[0];
4655   zloc = (Float_t)xyzLoc[2];
4656
4657   return kTRUE;
4658 }
4659 //------------------------------------------------------------------------
4660 //------------------------------------------------------------------------
4661 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4662 //
4663 // Method to be optimized further:
4664 // Aim: decide whether a track can be used for PlaneEff evaluation
4665 //      the decision is taken based on the track quality at the layer under study
4666 //      no information on the clusters on this layer has to be used
4667 //      The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4668 //      the cut is done on number of sigmas from the boundaries
4669 //
4670 //  Input: Actual track, layer [0,5] under study
4671 //  Output: none
4672 //  Return: kTRUE if this is a good track
4673 //
4674 // it will apply a pre-selection to obtain good quality tracks.
4675 // Here also  you will have the possibility to put a control on the
4676 // impact point of the track on the basic block, in order to exclude border regions
4677 // this will be done by calling a proper method of the AliITSPlaneEff class.
4678 //
4679 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4680 // return: Bool_t   -> kTRUE if usable track, kFALSE if not usable.
4681 //
4682   Int_t index[AliITSgeomTGeo::kNLayers];
4683   Int_t k;
4684   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4685   //
4686   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4687     index[k]=clusters[k];
4688   }
4689
4690   if(!fPlaneEff)
4691     {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4692   AliITSlayer &layer=fgLayers[ilayer];
4693   Double_t r=layer.GetR();
4694   AliITStrackMI tmp(*track);
4695
4696 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4697   Int_t ncl_out=0; Int_t ncl_in=0;
4698   for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
4699  AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4700                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4701    // if (tmp.GetClIndex(lay)>=0) ncl_out++;
4702 if(index[lay]>=0)ncl_out++;
4703   }
4704   for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
4705     AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4706                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4707    if (index[lay]>=0) ncl_in++;
4708   }
4709   Int_t ncl=ncl_out+ncl_in;
4710   Bool_t nextout = kFALSE;
4711   if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4712   else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4713   Bool_t nextin = kFALSE;
4714   if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4715   else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4716   // maximum number of missing clusters allowed in outermost layers
4717   if(ncl_out<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff())
4718      return kFALSE;
4719   // maximum number of missing clusters allowed (both in innermost and in outermost layers)
4720   if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4721      return kFALSE;
4722   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout)  return kFALSE;
4723   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin)   return kFALSE;
4724   if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4725  //  if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff()  && !tmp.GetConstrain()) return kFALSE;
4726
4727 // detector number
4728   Double_t phi,z;
4729   if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4730   Int_t idet=layer.FindDetectorIndex(phi,z);
4731   if(idet<0) { AliInfo(Form("cannot find detector"));
4732     return kFALSE;}
4733
4734   // here check if it has good Chi Square.
4735
4736   //propagate to the intersection with the detector plane
4737   const AliITSdetector &det=layer.GetDetector(idet);
4738   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4739
4740   Float_t locx; //
4741   Float_t locz; //
4742   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4743   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4744   if(key>fPlaneEff->Nblock()) return kFALSE;
4745   Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4746   if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4747   //***************
4748   // DEFINITION OF SEARCH ROAD FOR accepting a track
4749   //
4750   Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
4751   Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
4752   Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
4753   Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
4754                                                     // done for RecPoints
4755
4756   // exclude tracks at boundary between detectors
4757   //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4758   Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4759   AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4760   AliDebug(2,Form("Local:    track impact x=%f, z=%f",locx,locz));
4761   AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4762   if ( (locx-dx < blockXmn+boundaryWidth) ||
4763        (locx+dx > blockXmx-boundaryWidth) ||
4764        (locz-dz < blockZmn+boundaryWidth) ||
4765        (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4766   return kTRUE;
4767 }
4768 //------------------------------------------------------------------------
4769 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4770 //
4771 // This Method has to be optimized! For the time-being it uses the same criteria
4772 // as those used in the search of extra clusters for overlapping modules.
4773 //
4774 // Method Purpose: estabilish whether a track has produced a recpoint or not
4775 //                 in the layer under study (For Plane efficiency)
4776 //
4777 // inputs: AliITStrackMI* track  (pointer to a usable track)
4778 // outputs: none
4779 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4780 //               traversed by this very track. In details:
4781 //               - if a cluster can be associated to the track then call
4782 //                  AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4783 //               - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4784 //
4785   if(!fPlaneEff)
4786     {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4787   AliITSlayer &layer=fgLayers[ilayer];
4788   Double_t r=layer.GetR();
4789   AliITStrackMI tmp(*track);
4790
4791 // detector number
4792   Double_t phi,z;
4793   if (!tmp.GetPhiZat(r,phi,z)) return;
4794   Int_t idet=layer.FindDetectorIndex(phi,z);
4795
4796   if(idet<0) { AliInfo(Form("cannot find detector"));
4797     return;}
4798
4799
4800 //propagate to the intersection with the detector plane
4801   const AliITSdetector &det=layer.GetDetector(idet);
4802   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
4803
4804
4805 //***************
4806 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
4807 //
4808   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
4809                     TMath::Sqrt(tmp.GetSigmaZ2() +
4810                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4811                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4812                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
4813   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
4814                     TMath::Sqrt(tmp.GetSigmaY2() +
4815                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4816                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4817                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
4818
4819 // road in global (rphi,z) [i.e. in tracking ref. system]
4820   Double_t zmin = tmp.GetZ() - dz;
4821   Double_t zmax = tmp.GetZ() + dz;
4822   Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
4823   Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
4824
4825 // select clusters in road
4826   layer.SelectClusters(zmin,zmax,ymin,ymax);
4827
4828 // Define criteria for track-cluster association
4829   Double_t msz = tmp.GetSigmaZ2() +
4830   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4831   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4832   AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
4833   Double_t msy = tmp.GetSigmaY2() +
4834   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4835   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4836   AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
4837   if (tmp.GetConstrain()) {
4838     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
4839     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
4840   }  else {
4841     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
4842     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
4843   }
4844   msz = 1./msz; // 1/RoadZ^2
4845   msy = 1./msy; // 1/RoadY^2
4846 //
4847
4848   const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
4849   Int_t idetc=-1;
4850   Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
4851   //Double_t  tolerance=0.2;
4852   /*while ((cl=layer.GetNextCluster(clidx))!=0) {
4853     idetc = cl->GetDetectorIndex();
4854     if(idet!=idetc) continue;
4855     //Int_t ilay = cl->GetLayer();
4856
4857     if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
4858     if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
4859
4860     Double_t chi2=tmp.GetPredictedChi2(cl);
4861     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
4862   }*/
4863   Float_t locx; //
4864   Float_t locz; //
4865   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
4866 //
4867   AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
4868   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4869   if(key>fPlaneEff->Nblock()) return;
4870   Bool_t found=kFALSE;
4871   //if (ci>=0) {
4872   Double_t chi2;
4873   while ((cl=layer.GetNextCluster(clidx))!=0) {
4874     idetc = cl->GetDetectorIndex();
4875     if(idet!=idetc) continue;
4876     // here real control to see whether the cluster can be associated to the track.
4877     // cluster not associated to track
4878     if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
4879          (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy   > 1. ) continue;
4880     // calculate track-clusters chi2
4881     chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
4882                                                // in particular, the error associated to the cluster
4883     //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
4884     // chi2 cut
4885     if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
4886     found=kTRUE;
4887     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
4888    // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
4889    // track->SetExtraModule(ilayer,idetExtra);
4890   }
4891   if(!fPlaneEff->UpDatePlaneEff(found,key))
4892        AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
4893   if(fPlaneEff->GetCreateHistos()&&  AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
4894     Float_t tr[4]={99999.,99999.,9999.,9999.};    // initialize to high values
4895     Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
4896     Int_t cltype[2]={-999,-999};
4897                                                           // and the module
4898
4899 Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
4900
4901     tr[0]=locx;
4902     tr[1]=locz;
4903     tr[2]=TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
4904     tr[3]=TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
4905
4906     if (found){
4907       clu[0]=layer.GetCluster(ci)->GetDetLocalX();
4908       clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
4909       cltype[0]=layer.GetCluster(ci)->GetNy();
4910       cltype[1]=layer.GetCluster(ci)->GetNz();
4911
4912      // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
4913      //  X->50/sqrt(12)=14 micron   Z->450/sqrt(12)= 120 micron)
4914      // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
4915      // It is computed properly by calling the method
4916      // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4917      // T
4918      //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
4919       //if (tmp.PropagateTo(x,0.,0.)) {
4920         chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
4921         AliCluster c(*layer.GetCluster(ci));
4922         c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
4923         c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
4924         //if (layer.GetCluster(ci)->GetGlobalCov(cov))  // by using this, instead, you got nominal cluster errors
4925         clu[2]=TMath::Sqrt(c.GetSigmaY2());
4926         clu[3]=TMath::Sqrt(c.GetSigmaZ2());
4927       //}
4928     }
4929   // Compute the angles between the track and the module
4930       // compute the angle "in phi direction", i.e. the angle in the transverse plane
4931       // between the normal to the module and the projection (in the transverse plane) of the
4932       // track trajectory
4933     // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
4934     Float_t tgl = tmp.GetTgl();
4935     Float_t phitr   = tmp.GetSnp();
4936     phitr = TMath::ASin(phitr);
4937     Int_t volId = AliGeomManager::LayerToVolUIDSafe(ilayer+1 ,idet );
4938
4939     Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
4940     Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
4941    Double_t alpha =0.;
4942     alpha = tmp.GetAlpha();
4943     Double_t phiglob = alpha+phitr;
4944     Double_t p[3];
4945     p[0] = TMath::Cos(phiglob);
4946     p[1] = TMath::Sin(phiglob);
4947     p[2] = tgl;
4948     TVector3 pvec(p[0],p[1],p[2]);
4949     TVector3 normvec(rot[1],rot[4],rot[7]);
4950     Double_t angle = pvec.Angle(normvec);
4951
4952     if(angle>0.5*TMath::Pi()) angle = (TMath::Pi()-angle);
4953     angle *= 180./TMath::Pi();
4954
4955     //Trasverse Plane
4956     TVector3 pt(p[0],p[1],0);
4957     TVector3 normt(rot[1],rot[4],0);
4958     Double_t anglet = pt.Angle(normt);
4959
4960     Double_t phiPt = TMath::ATan2(p[1],p[0]);
4961     if(phiPt<0)phiPt+=2.*TMath::Pi();
4962     Double_t phiNorm = TMath::ATan2(rot[4],rot[1]);
4963     if(phiNorm<0) phiNorm+=2.*TMath::Pi();
4964     if(anglet>0.5*TMath::Pi()) anglet = (TMath::Pi()-anglet);
4965     if(phiNorm>phiPt) anglet*=-1.;// pt-->normt  clockwise: anglet>0
4966     if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
4967     anglet *= 180./TMath::Pi();
4968
4969      AngleModTrack[2]=(Float_t) angle;
4970      AngleModTrack[0]=(Float_t) anglet;
4971      // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
4972     AngleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
4973     AngleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
4974     AngleModTrack[1]*=180./TMath::Pi(); // in degree
4975
4976       fPlaneEff->FillHistos(key,found,tr,clu,cltype,AngleModTrack);
4977   }
4978 return;
4979 }