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