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