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