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