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