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