ITS/AliITStrackerMI.cxx

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