ITS/AliITStrackerMI.cxx

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