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