ITS/AliITStrackerMI.cxx

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