ITS/AliITStrackerMI.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id$ */
  16
  17 //-------------------------------------------------------------------------
  18 //               Implementation of the ITS tracker class
  19 //    It reads AliITSRecPoint clusters and creates AliITStrackMI tracks
  20 //                   and fills with them the ESD
  21 //          Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
  22 //          Current support and development:
  23 //                     Andrea Dainese, andrea.dainese@lnl.infn.it
  24 //     dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
  25 //     Params moved to AliITSRecoParam by: Andrea Dainese, INFN
  26 //     Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
  27 //-------------------------------------------------------------------------
  28
  29 #include <TMatrixD.h>
  30 #include <TTree.h>
  31 #include <TDatabasePDG.h>
  32 #include <TString.h>
  33 #include <TRandom.h>
  34 #include <TTreeStream.h>
  35
  36 #include "AliLog.h"
  37 #include "AliGeomManager.h"
  38 #include "AliITSPlaneEff.h"
  39 #include "AliTrackPointArray.h"
  40 #include "AliESDEvent.h"
  41 #include "AliESDtrack.h"
  42 #include "AliV0.h"
  43 #include "AliITSChannelStatus.h"
  44 #include "AliITSDetTypeRec.h"
  45 #include "AliITSRecPoint.h"
  46 #include "AliITSRecPointContainer.h"
  47 #include "AliITSgeomTGeo.h"
  48 #include "AliITSReconstructor.h"
  49 #include "AliITSClusterParam.h"
  50 #include "AliITSsegmentation.h"
  51 #include "AliITSCalibration.h"
  52 #include "AliITSPlaneEffSPD.h"
  53 #include "AliITSPlaneEffSDD.h"
  54 #include "AliITSPlaneEffSSD.h"
  55 #include "AliITSV0Finder.h"
  56 #include "AliITStrackerMI.h"
  57 #include "AliMathBase.h"
  58
  59 ClassImp(AliITStrackerMI)
  60
  61 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
  62
  63 AliITStrackerMI::AliITStrackerMI():AliTracker(),
  64 fI(0),
  65 fBestTrack(),
  66 fTrackToFollow(),
  67 fTrackHypothesys(),
  68 fBestHypothesys(),
  69 fOriginal(),
  70 fCurrentEsdTrack(),
  71 fPass(0),
  72 fAfterV0(kFALSE),
  73 fLastLayerToTrackTo(0),
  74 fCoefficients(0),
  75 fEsd(0),
  76 fTrackingPhase("Default"),
  77 fUseTGeo(3),
  78 fNtracks(0),
  79 fxOverX0Pipe(-1.),
  80 fxTimesRhoPipe(-1.),
  81 fxOverX0PipeTrks(0),
  82 fxTimesRhoPipeTrks(0),
  83 fxOverX0ShieldTrks(0),
  84 fxTimesRhoShieldTrks(0),
  85 fxOverX0LayerTrks(0),
  86 fxTimesRhoLayerTrks(0),
  87 fDebugStreamer(0),
  88 fITSChannelStatus(0),
  89 fkDetTypeRec(0),
  90 fPlaneEff(0) {
  91   //Default constructor
  92   Int_t i;
  93   for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
  94   for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
  95   for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
  96 }
  97 //------------------------------------------------------------------------
  98 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
  99 fI(AliITSgeomTGeo::GetNLayers()),
 100 fBestTrack(),
 101 fTrackToFollow(),
 102 fTrackHypothesys(),
 103 fBestHypothesys(),
 104 fOriginal(),
 105 fCurrentEsdTrack(),
 106 fPass(0),
 107 fAfterV0(kFALSE),
 108 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
 109 fCoefficients(0),
 110 fEsd(0),
 111 fTrackingPhase("Default"),
 112 fUseTGeo(3),
 113 fNtracks(0),
 114 fxOverX0Pipe(-1.),
 115 fxTimesRhoPipe(-1.),
 116 fxOverX0PipeTrks(0),
 117 fxTimesRhoPipeTrks(0),
 118 fxOverX0ShieldTrks(0),
 119 fxTimesRhoShieldTrks(0),
 120 fxOverX0LayerTrks(0),
 121 fxTimesRhoLayerTrks(0),
 122 fDebugStreamer(0),
 123 fITSChannelStatus(0),
 124 fkDetTypeRec(0),
 125 fPlaneEff(0) {
 126   //--------------------------------------------------------------------
 127   //This is the AliITStrackerMI constructor
 128   //--------------------------------------------------------------------
 129   if (geom) {
 130     AliWarning("\"geom\" is actually a dummy argument !");
 131   }
 132
 133   fCoefficients = 0;
 134   fAfterV0     = kFALSE;
 135
 136   for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 137     Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
 138     Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
 139
 140     Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
 141     AliITSgeomTGeo::GetTranslation(i,1,1,xyz);
 142     Double_t poff=TMath::ATan2(y,x);
 143     Double_t zoff=z;
 144
 145     AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
 146     Double_t r=TMath::Sqrt(x*x + y*y);
 147
 148     AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
 149     r += TMath::Sqrt(x*x + y*y);
 150     AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
 151     r += TMath::Sqrt(x*x + y*y);
 152     AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
 153     r += TMath::Sqrt(x*x + y*y);
 154     r*=0.25;
 155
 156     new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
 157
 158     for (Int_t j=1; j<nlad+1; j++) {
 159       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
 160         TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
 161         const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
 162         m.Multiply(tm);
 163         Double_t txyz[3]={0.};
 164         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 165         m.LocalToMaster(txyz,xyz);
 166         r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 167         Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
 168
 169         if (phi<0) phi+=TMath::TwoPi();
 170         else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
 171
 172         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 173         new(&det) AliITSdetector(r,phi);
 174         // compute the real radius (with misalignment)
 175         TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
 176         mmisal.Multiply(tm);
 177         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 178         mmisal.LocalToMaster(txyz,xyz);
 179         Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 180         det.SetRmisal(rmisal);
 181
 182       } // end loop on detectors
 183     } // end loop on ladders
 184     fForceSkippingOfLayer[i] = 0;
 185   } // end loop on layers
 186
 187
 188   fI=AliITSgeomTGeo::GetNLayers();
 189
 190   fPass=0;
 191   fConstraint[0]=1; fConstraint[1]=0;
 192
 193   Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
 194                      AliITSReconstructor::GetRecoParam()->GetYVdef(),
 195                      AliITSReconstructor::GetRecoParam()->GetZVdef()};
 196   Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
 197                      AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
 198                      AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
 199   SetVertex(xyzVtx,ersVtx);
 200
 201   fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
 202   for (Int_t i=0;i<100000;i++){
 203     fBestTrackIndex[i]=0;
 204   }
 205
 206   // store positions of centre of SPD modules (in z)
 207   //  The convetion is that fSPDdetzcentre is ordered from -z to +z
 208   Double_t tr[3];
 209   if (tr[2]<0) { // old geom (up to v5asymmPPR)
 210     AliITSgeomTGeo::GetTranslation(1,1,1,tr);
 211     fSPDdetzcentre[0] = tr[2];
 212     AliITSgeomTGeo::GetTranslation(1,1,2,tr);
 213     fSPDdetzcentre[1] = tr[2];
 214     AliITSgeomTGeo::GetTranslation(1,1,3,tr);
 215     fSPDdetzcentre[2] = tr[2];
 216     AliITSgeomTGeo::GetTranslation(1,1,4,tr);
 217     fSPDdetzcentre[3] = tr[2];
 218   } else { // new geom (from v11Hybrid)
 219     AliITSgeomTGeo::GetTranslation(1,1,4,tr);
 220     fSPDdetzcentre[0] = tr[2];
 221     AliITSgeomTGeo::GetTranslation(1,1,3,tr);
 222     fSPDdetzcentre[1] = tr[2];
 223     AliITSgeomTGeo::GetTranslation(1,1,2,tr);
 224     fSPDdetzcentre[2] = tr[2];
 225     AliITSgeomTGeo::GetTranslation(1,1,1,tr);
 226     fSPDdetzcentre[3] = tr[2];
 227   }
 228
 229   fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
 230   if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
 231     AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
 232     fUseTGeo = 3;
 233   }
 234
 235   for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
 236   for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
 237
 238   fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
 239
 240   // only for plane efficiency evaluation
 241   if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
 242       AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
 243     Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
 244     if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
 245       AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
 246     if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
 247     else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
 248     else fPlaneEff = new AliITSPlaneEffSSD();
 249     if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
 250        if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
 251     if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
 252   }
 253 }
 254 //------------------------------------------------------------------------
 255 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
 256 fI(tracker.fI),
 257 fBestTrack(tracker.fBestTrack),
 258 fTrackToFollow(tracker.fTrackToFollow),
 259 fTrackHypothesys(tracker.fTrackHypothesys),
 260 fBestHypothesys(tracker.fBestHypothesys),
 261 fOriginal(tracker.fOriginal),
 262 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
 263 fPass(tracker.fPass),
 264 fAfterV0(tracker.fAfterV0),
 265 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
 266 fCoefficients(tracker.fCoefficients),
 267 fEsd(tracker.fEsd),
 268 fTrackingPhase(tracker.fTrackingPhase),
 269 fUseTGeo(tracker.fUseTGeo),
 270 fNtracks(tracker.fNtracks),
 271 fxOverX0Pipe(tracker.fxOverX0Pipe),
 272 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
 273 fxOverX0PipeTrks(0),
 274 fxTimesRhoPipeTrks(0),
 275 fxOverX0ShieldTrks(0),
 276 fxTimesRhoShieldTrks(0),
 277 fxOverX0LayerTrks(0),
 278 fxTimesRhoLayerTrks(0),
 279 fDebugStreamer(tracker.fDebugStreamer),
 280 fITSChannelStatus(tracker.fITSChannelStatus),
 281 fkDetTypeRec(tracker.fkDetTypeRec),
 282 fPlaneEff(tracker.fPlaneEff) {
 283   //Copy constructor
 284   Int_t i;
 285   for(i=0;i<4;i++) {
 286     fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
 287   }
 288   for(i=0;i<6;i++) {
 289     fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
 290     fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
 291   }
 292   for(i=0;i<2;i++) {
 293     fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
 294     fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
 295   }
 296 }
 297 //------------------------------------------------------------------------
 298 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
 299   //Assignment operator
 300   this->~AliITStrackerMI();
 301   new(this) AliITStrackerMI(tracker);
 302   return *this;
 303 }
 304 //------------------------------------------------------------------------
 305 AliITStrackerMI::~AliITStrackerMI()
 306 {
 307   //
 308   //destructor
 309   //
 310   if (fCoefficients) delete [] fCoefficients;
 311   DeleteTrksMaterialLUT();
 312   if (fDebugStreamer) {
 313     //fDebugStreamer->Close();
 314     delete fDebugStreamer;
 315   }
 316   if(fITSChannelStatus) delete fITSChannelStatus;
 317   if(fPlaneEff) delete fPlaneEff;
 318 }
 319 //------------------------------------------------------------------------
 320 void AliITStrackerMI::ReadBadFromDetTypeRec() {
 321   //--------------------------------------------------------------------
 322   //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
 323   //i.e. from OCDB
 324   //--------------------------------------------------------------------
 325
 326   if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
 327
 328   Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
 329
 330   if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
 331
 332   // ITS channels map
 333   if(fITSChannelStatus) delete fITSChannelStatus;
 334   fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
 335
 336   // ITS detectors and chips
 337   Int_t i=0,j=0,k=0,ndet=0;
 338   for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 339     Int_t nBadDetsPerLayer=0;
 340     ndet=AliITSgeomTGeo::GetNDetectors(i);
 341     for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
 342       for (k=1; k<ndet+1; k++) {
 343         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 344         det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
 345         if(det.IsBad()) {nBadDetsPerLayer++;}
 346       } // end loop on detectors
 347     } // end loop on ladders
 348     Info("ReadBadFromDetTypeRec",Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
 349   } // end loop on layers
 350
 351   return;
 352 }
 353 //------------------------------------------------------------------------
 354 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
 355   //--------------------------------------------------------------------
 356   //This function loads ITS clusters
 357   //--------------------------------------------------------------------
 358
 359   TClonesArray *clusters = NULL;
 360   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
 361   clusters=rpcont->FetchClusters(0,cTree);
 362   if(!(rpcont->IsSPDActive() || rpcont->IsSDDActive() || rpcont->IsSSDActive())){
 363       AliError("ITS is not in a known running configuration: SPD, SDD and SSD are not active");
 364       return 1;
 365   }
 366   Int_t i=0,j=0,ndet=0;
 367   Int_t detector=0;
 368   for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
 369     ndet=fgLayers[i].GetNdetectors();
 370     Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
 371     for (; j<jmax; j++) {
 372       //      if (!cTree->GetEvent(j)) continue;
 373       clusters = rpcont->UncheckedGetClusters(j);
 374       if(!clusters)continue;
 375       Int_t ncl=clusters->GetEntriesFast();
 376       SignDeltas(clusters,GetZ());
 377
 378       while (ncl--) {
 379         AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
 380         detector=c->GetDetectorIndex();
 381
 382         if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
 383
 384         Int_t retval = fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
 385         if(retval) {
 386           AliWarning(Form("Too many clusters on layer %d!",i));
 387           break;
 388         }
 389       }
 390
 391       // add dead zone "virtual" cluster in SPD, if there is a cluster within
 392       // zwindow cm from the dead zone
 393       //  This method assumes that fSPDdetzcentre is ordered from -z to +z
 394       if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
 395         for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
 396           Int_t lab[4]   = {0,0,0,detector};
 397           Int_t info[3]  = {0,0,i};
 398           Float_t q      = 0.; // this identifies virtual clusters
 399           Float_t hit[5] = {xdead,
 400                             0.,
 401                             AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
 402                             AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
 403                             q};
 404           Bool_t local   = kTRUE;
 405           Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
 406           hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
 407           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 408             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 409           hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
 410           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 411             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 412           hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
 413           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 414             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 415           hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
 416           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 417             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 418           hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
 419           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 420             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 421           hit[1] = fSPDdetzcentre[3]-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         }
 425       } // "virtual" clusters in SPD
 426
 427     }
 428     //
 429     fgLayers[i].ResetRoad(); //road defined by the cluster density
 430     fgLayers[i].SortClusters();
 431   }
 432
 433   // check whether we have to skip some layers
 434   SetForceSkippingOfLayer();
 435
 436   return 0;
 437 }
 438 //------------------------------------------------------------------------
 439 void AliITStrackerMI::UnloadClusters() {
 440   //--------------------------------------------------------------------
 441   //This function unloads ITS clusters
 442   //--------------------------------------------------------------------
 443   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
 444 }
 445 //------------------------------------------------------------------------
 446 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
 447   //--------------------------------------------------------------------
 448   // Publishes all pointers to clusters known to the tracker into the
 449   // passed object array.
 450   // The ownership is not transfered - the caller is not expected to delete
 451   // the clusters.
 452   //--------------------------------------------------------------------
 453
 454   for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
 455     for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
 456       AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
 457       array->AddLast(cl);
 458     }
 459   }
 460
 461   return;
 462 }
 463 //------------------------------------------------------------------------
 464 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
 465   //--------------------------------------------------------------------
 466   // Correction for the material between the TPC and the ITS
 467   //--------------------------------------------------------------------
 468   if (t->GetX() > AliITSRecoParam::Getriw()) {   // inward direction
 469       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
 470       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
 471       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))  return 0;// ITS screen
 472   } else if (t->GetX() < AliITSRecoParam::Getrs()) {  // outward direction
 473       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))        return 0;// ITS screen
 474       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1))       return 0;// TPC central drum
 475       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
 476   } else {
 477     printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
 478     return 0;
 479   }
 480
 481   return 1;
 482 }
 483 //------------------------------------------------------------------------
 484 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
 485   //--------------------------------------------------------------------
 486   // This functions reconstructs ITS tracks
 487   // The clusters must be already loaded !
 488   //--------------------------------------------------------------------
 489
 490   AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
 491
 492   fTrackingPhase="Clusters2Tracks";
 493
 494   TObjArray itsTracks(15000);
 495   fOriginal.Clear();
 496   fEsd = event;         // store pointer to the esd
 497
 498   // temporary (for cosmics)
 499   if(event->GetVertex()) {
 500     TString title = event->GetVertex()->GetTitle();
 501     if(title.Contains("cosmics")) {
 502       Double_t xyz[3]={GetX(),GetY(),GetZ()};
 503       Double_t exyz[3]={0.1,0.1,0.1};
 504       SetVertex(xyz,exyz);
 505     }
 506   }
 507   // temporary
 508   Int_t noesd = 0;
 509   {/* Read ESD tracks */
 510     Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
 511     Int_t nentr=event->GetNumberOfTracks();
 512     noesd=nentr;
 513     //    Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
 514     while (nentr--) {
 515       AliESDtrack *esd=event->GetTrack(nentr);
 516       //  ---- for debugging:
 517       //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); }
 518
 519       if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
 520       if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
 521       if (esd->GetStatus()&AliESDtrack::kITSin) continue;
 522       if (esd->GetKinkIndex(0)>0) continue;   //kink daughter
 523       AliITStrackMI *t=0;
 524       try {
 525         t=new AliITStrackMI(*esd);
 526       } catch (const Char_t *msg) {
 527         //Warning("Clusters2Tracks",msg);
 528         delete t;
 529         continue;
 530       }
 531       t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP());              //I.B.
 532       Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
 533
 534
 535       // look at the ESD mass hypothesys !
 536       if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
 537       // write expected q
 538       t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 539
 540       if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
 541         //track - can be  V0 according to TPC
 542       } else {
 543         if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
 544           delete t;
 545           continue;
 546         }
 547         if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
 548           delete t;
 549           continue;
 550         }
 551         if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
 552           delete t;
 553           continue;
 554         }
 555         if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
 556           delete t;
 557           continue;
 558         }
 559       }
 560       t->SetReconstructed(kFALSE);
 561       itsTracks.AddLast(t);
 562       fOriginal.AddLast(t);
 563     }
 564   } /* End Read ESD tracks */
 565
 566   itsTracks.Sort();
 567   fOriginal.Sort();
 568   Int_t nentr=itsTracks.GetEntriesFast();
 569   fTrackHypothesys.Expand(nentr);
 570   fBestHypothesys.Expand(nentr);
 571   MakeCoefficients(nentr);
 572   if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
 573   Int_t ntrk=0;
 574   // THE TWO TRACKING PASSES
 575   for (fPass=0; fPass<2; fPass++) {
 576      Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
 577      for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
 578        AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
 579        if (t==0) continue;              //this track has been already tracked
 580        //cout<<"========== "<<fPass<<"    "<<fCurrentEsdTrack<<" =========\n";
 581        if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue;  //this track was  already  "succesfully" reconstructed
 582        Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz);              //I.B.
 583        if (fConstraint[fPass]) {
 584          if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
 585              TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
 586        }
 587
 588        Int_t tpcLabel=t->GetLabel(); //save the TPC track label
 589        AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
 590        fI = 6;
 591        ResetTrackToFollow(*t);
 592        ResetBestTrack();
 593
 594        FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
 595
 596
 597        SortTrackHypothesys(fCurrentEsdTrack,20,0);  //MI change
 598        //
 599        AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
 600        if (!besttrack) continue;
 601        besttrack->SetLabel(tpcLabel);
 602        //       besttrack->CookdEdx();
 603        CookdEdx(besttrack);
 604        besttrack->SetFakeRatio(1.);
 605        CookLabel(besttrack,0.); //For comparison only
 606        UpdateESDtrack(besttrack,AliESDtrack::kITSin);
 607
 608        if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue;  //to be tracked also without vertex constrain
 609
 610        t->SetReconstructed(kTRUE);
 611        ntrk++;
 612        AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
 613      }
 614      GetBestHypothesysMIP(itsTracks);
 615   } // end loop on the two tracking passes
 616
 617   if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
 618   if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
 619   fAfterV0 = kTRUE;
 620   //
 621   itsTracks.Delete();
 622   //
 623   Int_t entries = fTrackHypothesys.GetEntriesFast();
 624   for (Int_t ientry=0; ientry<entries; ientry++) {
 625     TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
 626     if (array) array->Delete();
 627     delete fTrackHypothesys.RemoveAt(ientry);
 628   }
 629
 630   fTrackHypothesys.Delete();
 631   fBestHypothesys.Delete();
 632   fOriginal.Clear();
 633   delete [] fCoefficients;
 634   fCoefficients=0;
 635   DeleteTrksMaterialLUT();
 636
 637   AliInfo(Form("Number of prolonged tracks: %d out of %d ESD tracks",ntrk,noesd));
 638
 639   fTrackingPhase="Default";
 640
 641   return 0;
 642 }
 643 //------------------------------------------------------------------------
 644 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
 645   //--------------------------------------------------------------------
 646   // This functions propagates reconstructed ITS tracks back
 647   // The clusters must be loaded !
 648   //--------------------------------------------------------------------
 649   fTrackingPhase="PropagateBack";
 650   Int_t nentr=event->GetNumberOfTracks();
 651   //  Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
 652
 653   Int_t ntrk=0;
 654   for (Int_t i=0; i<nentr; i++) {
 655      AliESDtrack *esd=event->GetTrack(i);
 656
 657      if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
 658      if (esd->GetStatus()&AliESDtrack::kITSout) continue;
 659
 660      AliITStrackMI *t=0;
 661      try {
 662         t=new AliITStrackMI(*esd);
 663      } catch (const Char_t *msg) {
 664        //Warning("PropagateBack",msg);
 665         delete t;
 666         continue;
 667      }
 668      t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 669
 670      ResetTrackToFollow(*t);
 671
 672      /*
 673      // propagate to vertex [SR, GSI 17.02.2003]
 674      // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
 675      if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
 676        if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
 677          fTrackToFollow.StartTimeIntegral();
 678        // from vertex to outside pipe
 679        CorrectForPipeMaterial(&fTrackToFollow,"outward");
 680        }*/
 681      // Start time integral and add distance from current position to vertex
 682      Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
 683      fTrackToFollow.GetXYZ(xyzTrk);
 684      Double_t dst2 = 0.;
 685      for (Int_t icoord=0; icoord<3; icoord++) {
 686        Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
 687        dst2 += di*di;
 688      }
 689      fTrackToFollow.StartTimeIntegral();
 690      fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
 691
 692      fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
 693      if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
 694         if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
 695           delete t;
 696           continue;
 697         }
 698         fTrackToFollow.SetLabel(t->GetLabel());
 699         //fTrackToFollow.CookdEdx();
 700         CookLabel(&fTrackToFollow,0.); //For comparison only
 701         fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
 702         //UseClusters(&fTrackToFollow);
 703         ntrk++;
 704      }
 705      delete t;
 706   }
 707
 708   AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
 709
 710   fTrackingPhase="Default";
 711
 712   return 0;
 713 }
 714 //------------------------------------------------------------------------
 715 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
 716   //--------------------------------------------------------------------
 717   // This functions refits ITS tracks using the
 718   // "inward propagated" TPC tracks
 719   // The clusters must be loaded !
 720   //--------------------------------------------------------------------
 721   fTrackingPhase="RefitInward";
 722
 723   if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
 724
 725   Int_t nentr=event->GetNumberOfTracks();
 726   //  Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
 727
 728   Int_t ntrk=0;
 729   for (Int_t i=0; i<nentr; i++) {
 730     AliESDtrack *esd=event->GetTrack(i);
 731
 732     if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
 733     if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
 734     if (esd->GetStatus()&AliESDtrack::kTPCout)
 735       if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
 736
 737     AliITStrackMI *t=0;
 738     try {
 739         t=new AliITStrackMI(*esd);
 740     } catch (const Char_t *msg) {
 741       //Warning("RefitInward",msg);
 742         delete t;
 743         continue;
 744     }
 745     t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 746     if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
 747        delete t;
 748        continue;
 749     }
 750
 751     ResetTrackToFollow(*t);
 752     fTrackToFollow.ResetClusters();
 753
 754     if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
 755       fTrackToFollow.ResetCovariance(10.);
 756
 757     //Refitting...
 758     Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
 759                AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
 760
 761     AliDebug(2,Form("Refit LABEL %d  %d",t->GetLabel(),t->GetNumberOfClusters()));
 762     if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,kTRUE,pe)) {
 763        AliDebug(2,"  refit OK");
 764        fTrackToFollow.SetLabel(t->GetLabel());
 765        //       fTrackToFollow.CookdEdx();
 766        CookdEdx(&fTrackToFollow);
 767
 768        CookLabel(&fTrackToFollow,0.0); //For comparison only
 769
 770        //The beam pipe
 771        if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
 772          fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
 773          AliESDtrack  *esdTrack =fTrackToFollow.GetESDtrack();
 774          //printf("                                       %d\n",esdTrack->GetITSModuleIndex(0));
 775          //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
 776          Double_t r[3]={0.,0.,0.};
 777          Double_t maxD=3.;
 778          esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
 779          ntrk++;
 780        }
 781     }
 782     delete t;
 783   }
 784
 785   AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
 786
 787   fTrackingPhase="Default";
 788
 789   return 0;
 790 }
 791 //------------------------------------------------------------------------
 792 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
 793   //--------------------------------------------------------------------
 794   //       Return pointer to a given cluster
 795   //--------------------------------------------------------------------
 796   Int_t l=(index & 0xf0000000) >> 28;
 797   Int_t c=(index & 0x0fffffff) >> 00;
 798   return fgLayers[l].GetCluster(c);
 799 }
 800 //------------------------------------------------------------------------
 801 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
 802   //--------------------------------------------------------------------
 803   // Get track space point with index i
 804   //--------------------------------------------------------------------
 805
 806   Int_t l=(index & 0xf0000000) >> 28;
 807   Int_t c=(index & 0x0fffffff) >> 00;
 808   AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 809   Int_t idet = cl->GetDetectorIndex();
 810
 811   Float_t xyz[3];
 812   Float_t cov[6];
 813   cl->GetGlobalXYZ(xyz);
 814   cl->GetGlobalCov(cov);
 815   p.SetXYZ(xyz, cov);
 816   p.SetCharge(cl->GetQ());
 817   p.SetDriftTime(cl->GetDriftTime());
 818   p.SetChargeRatio(cl->GetChargeRatio());
 819   p.SetClusterType(cl->GetClusterType());
 820   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 821   switch (l) {
 822   case 0:
 823     iLayer = AliGeomManager::kSPD1;
 824     break;
 825   case 1:
 826     iLayer = AliGeomManager::kSPD2;
 827     break;
 828   case 2:
 829     iLayer = AliGeomManager::kSDD1;
 830     break;
 831   case 3:
 832     iLayer = AliGeomManager::kSDD2;
 833     break;
 834   case 4:
 835     iLayer = AliGeomManager::kSSD1;
 836     break;
 837   case 5:
 838     iLayer = AliGeomManager::kSSD2;
 839     break;
 840   default:
 841     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 842     break;
 843   };
 844   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 845   p.SetVolumeID((UShort_t)volid);
 846   return kTRUE;
 847 }
 848 //------------------------------------------------------------------------
 849 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
 850                         AliTrackPoint& p, const AliESDtrack *t) {
 851   //--------------------------------------------------------------------
 852   // Get track space point with index i
 853   // (assign error estimated during the tracking)
 854   //--------------------------------------------------------------------
 855
 856   Int_t l=(index & 0xf0000000) >> 28;
 857   Int_t c=(index & 0x0fffffff) >> 00;
 858   const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 859   Int_t idet = cl->GetDetectorIndex();
 860
 861   const AliITSdetector &det=fgLayers[l].GetDetector(idet);
 862
 863   // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
 864   Float_t detxy[2];
 865   detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
 866   detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
 867   Double_t alpha = t->GetAlpha();
 868   Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
 869   Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,GetBz()));
 870   phi += alpha-det.GetPhi();
 871   Float_t tgphi = TMath::Tan(phi);
 872
 873   Float_t tgl = t->GetTgl(); // tgl about const along track
 874   Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
 875
 876   Float_t errtrky,errtrkz,covyz;
 877   Bool_t addMisalErr=kFALSE;
 878   AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
 879
 880   Float_t xyz[3];
 881   Float_t cov[6];
 882   cl->GetGlobalXYZ(xyz);
 883   //  cl->GetGlobalCov(cov);
 884   Float_t pos[3] = {0.,0.,0.};
 885   AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
 886   tmpcl.GetGlobalCov(cov);
 887
 888   p.SetXYZ(xyz, cov);
 889   p.SetCharge(cl->GetQ());
 890   p.SetDriftTime(cl->GetDriftTime());
 891   p.SetChargeRatio(cl->GetChargeRatio());
 892   p.SetClusterType(cl->GetClusterType());
 893
 894   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 895   switch (l) {
 896   case 0:
 897     iLayer = AliGeomManager::kSPD1;
 898     break;
 899   case 1:
 900     iLayer = AliGeomManager::kSPD2;
 901     break;
 902   case 2:
 903     iLayer = AliGeomManager::kSDD1;
 904     break;
 905   case 3:
 906     iLayer = AliGeomManager::kSDD2;
 907     break;
 908   case 4:
 909     iLayer = AliGeomManager::kSSD1;
 910     break;
 911   case 5:
 912     iLayer = AliGeomManager::kSSD2;
 913     break;
 914   default:
 915     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 916     break;
 917   };
 918   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 919
 920   p.SetVolumeID((UShort_t)volid);
 921   return kTRUE;
 922 }
 923 //------------------------------------------------------------------------
 924 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
 925 {
 926   //--------------------------------------------------------------------
 927   // Follow prolongation tree
 928   //--------------------------------------------------------------------
 929   //
 930   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
 931   Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
 932
 933
 934   AliESDtrack * esd = otrack->GetESDtrack();
 935   if (esd->GetV0Index(0)>0) {
 936     // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
 937     //                      mapping of ESD track is different as ITS track in Containers
 938     //                      Need something more stable
 939     //                      Indexes are set back again to the ESD track indexes in UpdateTPCV0
 940     for (Int_t i=0;i<3;i++){
 941       Int_t  index = esd->GetV0Index(i);
 942       if (index==0) break;
 943       AliESDv0 * vertex = fEsd->GetV0(index);
 944       if (vertex->GetStatus()<0) continue;     // rejected V0
 945       //
 946       if (esd->GetSign()>0) {
 947         vertex->SetIndex(0,esdindex);
 948       } else {
 949         vertex->SetIndex(1,esdindex);
 950       }
 951     }
 952   }
 953   TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
 954   if (!bestarray){
 955     bestarray = new TObjArray(5);
 956     fBestHypothesys.AddAt(bestarray,esdindex);
 957   }
 958
 959   //
 960   //setup tree of the prolongations
 961   //
 962   static AliITStrackMI tracks[7][100];
 963   AliITStrackMI *currenttrack;
 964   static AliITStrackMI currenttrack1;
 965   static AliITStrackMI currenttrack2;
 966   static AliITStrackMI backuptrack;
 967   Int_t ntracks[7];
 968   Int_t nindexes[7][100];
 969   Float_t normalizedchi2[100];
 970   for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
 971   otrack->SetNSkipped(0);
 972   new (&(tracks[6][0])) AliITStrackMI(*otrack);
 973   ntracks[6]=1;
 974   for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
 975   Int_t modstatus = 1; // found
 976   Float_t xloc,zloc;
 977   //
 978   //
 979   // follow prolongations
 980   for (Int_t ilayer=5; ilayer>=0; ilayer--) {
 981     AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
 982     fI = ilayer;
 983     //
 984     AliITSlayer &layer=fgLayers[ilayer];
 985     Double_t r = layer.GetR();
 986     ntracks[ilayer]=0;
 987     //
 988     //
 989     Int_t nskipped=0;
 990     Float_t nused =0;
 991     for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
 992       //set current track
 993       if (ntracks[ilayer]>=100) break;
 994       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
 995       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
 996       if (ntracks[ilayer]>15+ilayer){
 997         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
 998         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
 999       }
1000
1001       new(&currenttrack1)  AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
1002
1003       // material between SSD and SDD, SDD and SPD
1004       if (ilayer==3)
1005         if(!CorrectForShieldMaterial(&currenttrack1,"SDD","inward")) continue;
1006       if (ilayer==1)
1007         if(!CorrectForShieldMaterial(&currenttrack1,"SPD","inward")) continue;
1008
1009       // detector number
1010       Double_t phi,z;
1011       if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
1012       Int_t idet=layer.FindDetectorIndex(phi,z);
1013
1014       Double_t trackGlobXYZ1[3];
1015       if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1016
1017       // Get the budget to the primary vertex for the current track being prolonged
1018       Double_t budgetToPrimVertex = GetEffectiveThickness();
1019
1020       // check if we allow a prolongation without point
1021       Int_t skip = CheckSkipLayer(&currenttrack1,ilayer,idet);
1022       if (skip) {
1023         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1024         // propagate to the layer radius
1025         Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1026         if(!vtrack->Propagate(xToGo)) continue;
1027         // apply correction for material of the current layer
1028         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1029         vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1030         vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1031         vtrack->SetClIndex(ilayer,-1);
1032         modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1033         if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1034           vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1035         }
1036         if(constrain) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1037         ntracks[ilayer]++;
1038         continue;
1039       }
1040
1041       // track outside layer acceptance in z
1042       if (idet<0) continue;
1043
1044       //propagate to the intersection with the detector plane
1045       const AliITSdetector &det=layer.GetDetector(idet);
1046       new(&currenttrack2)  AliITStrackMI(currenttrack1);
1047       if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1048       if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1049       currenttrack1.SetDetectorIndex(idet);
1050       currenttrack2.SetDetectorIndex(idet);
1051       if(!LocalModuleCoord(ilayer,idet,&currenttrack1,xloc,zloc)) continue; // local module coords
1052
1053       //***************
1054       // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1055       //
1056       // road in global (rphi,z) [i.e. in tracking ref. system]
1057       Double_t zmin,zmax,ymin,ymax;
1058       if (!ComputeRoad(&currenttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1059
1060       // select clusters in road
1061       layer.SelectClusters(zmin,zmax,ymin,ymax);
1062       //********************
1063
1064       // Define criteria for track-cluster association
1065       Double_t msz = currenttrack1.GetSigmaZ2() +
1066         AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1067         AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1068         AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1069       Double_t msy = currenttrack1.GetSigmaY2() +
1070         AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1071         AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1072         AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1073       if (constrain) {
1074         msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1075         msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1076       }  else {
1077         msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1078         msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1079       }
1080       msz = 1./msz; // 1/RoadZ^2
1081       msy = 1./msy; // 1/RoadY^2
1082
1083       //
1084       //
1085       // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1086       //
1087       const AliITSRecPoint *cl=0;
1088       Int_t clidx=-1;
1089       Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1090       Bool_t deadzoneSPD=kFALSE;
1091       currenttrack = &currenttrack1;
1092
1093       // check if the road contains a dead zone
1094       Bool_t noClusters = kFALSE;
1095       if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1096       if (noClusters) AliDebug(2,"no clusters in road");
1097       Double_t dz=0.5*(zmax-zmin);
1098       Double_t dy=0.5*(ymax-ymin);
1099       Int_t dead = CheckDeadZone(&currenttrack1,ilayer,idet,dz,dy,noClusters);
1100       if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1101       // create a prolongation without clusters (check also if there are no clusters in the road)
1102       if (dead ||
1103           (noClusters &&
1104            AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1105         AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1106         updatetrack->SetClIndex(ilayer,-1);
1107         if (dead==0) {
1108           modstatus = 5; // no cls in road
1109         } else if (dead==1) {
1110           modstatus = 7; // holes in z in SPD
1111         } else if (dead==2 || dead==3 || dead==4) {
1112           modstatus = 2; // dead from OCDB
1113         }
1114         updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1115         // apply correction for material of the current layer
1116         CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1117         if (constrain) { // apply vertex constrain
1118           updatetrack->SetConstrain(constrain);
1119           Bool_t isPrim = kTRUE;
1120           if (ilayer<4) { // check that it's close to the vertex
1121             updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1122             if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1123                 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1124                 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1125                 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1126           }
1127           if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1128         }
1129         updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1130         if (dead) {
1131           if (dead==1) { // dead zone at z=0,+-7cm in SPD
1132             updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1133             deadzoneSPD=kTRUE;
1134           } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1135             updatetrack->SetDeadZoneProbability(ilayer,1.);
1136           } else if (dead==4) { // at least a single dead channel from OCDB
1137             updatetrack->SetDeadZoneProbability(ilayer,0.);
1138           }
1139         }
1140         ntracks[ilayer]++;
1141       }
1142
1143       clidx=-1;
1144       // loop over clusters in the road
1145       while ((cl=layer.GetNextCluster(clidx))!=0) {
1146         if (ntracks[ilayer]>95) break; //space for skipped clusters
1147         Bool_t changedet =kFALSE;
1148         if (TMath::Abs(cl->GetQ())<1.e-13 && deadzoneSPD==kTRUE) continue;
1149         Int_t idetc=cl->GetDetectorIndex();
1150
1151         if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1152           // take into account misalignment (bring track to real detector plane)
1153           Double_t xTrOrig = currenttrack->GetX();
1154           if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1155           // a first cut on track-cluster distance
1156           if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1157                (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1158             {  // cluster not associated to track
1159               AliDebug(2,"not associated");
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) 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         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: check which the innermost layer crossed
2252   // by the track
2253   Int_t innermostlayer=5;
2254   Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2255   for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2256     if(drphi < fgLayers[innermostlayer].GetR()) break;
2257   }
2258   AliDebug(2,Form(" drphi  %f  innermost %d",drphi,innermostlayer));
2259
2260   Int_t modstatus=1; // found
2261   Float_t xloc,zloc;
2262   Int_t from, to, step;
2263   if (xx > track->GetX()) {
2264       from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2265       step=+1;
2266   } else {
2267       from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2268       step=-1;
2269   }
2270   TString dir = (step>0 ? "outward" : "inward");
2271
2272   for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2273      AliITSlayer &layer=fgLayers[ilayer];
2274      Double_t r=layer.GetR();
2275
2276      if (step<0 && xx>r) break;
2277
2278      // material between SSD and SDD, SDD and SPD
2279      Double_t hI=ilayer-0.5*step;
2280      if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2281        if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2282      if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2283        if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2284
2285
2286      Double_t oldGlobXYZ[3];
2287      if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2288
2289      // continue if we are already beyond this layer
2290      Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2291      if(step>0 && oldGlobR > r) continue; // going outward
2292      if(step<0 && oldGlobR < r) continue; // going inward
2293
2294      Double_t phi,z;
2295      if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2296
2297      Int_t idet=layer.FindDetectorIndex(phi,z);
2298
2299      // check if we allow a prolongation without point for large-eta tracks
2300      Int_t skip = CheckSkipLayer(track,ilayer,idet);
2301      if (skip==2) {
2302        modstatus = 4; // out in z
2303        if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2304          track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2305        }
2306        // cross layer
2307        // apply correction for material of the current layer
2308        // add time if going outward
2309        CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2310        continue;
2311      }
2312
2313      if (idet<0) return kFALSE;
2314
2315
2316      const AliITSdetector &det=layer.GetDetector(idet);
2317      // only for ITS-SA tracks refit
2318      if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2319      //
2320      if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2321
2322      track->SetDetectorIndex(idet);
2323      if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2324
2325      Double_t dz,zmin,zmax,dy,ymin,ymax;
2326
2327      const AliITSRecPoint *clAcc=0;
2328      Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2329
2330      Int_t idx=index[ilayer];
2331      if (idx>=0) { // cluster in this layer
2332        modstatus = 6; // no refit
2333        const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2334        if (cl) {
2335          if (idet != cl->GetDetectorIndex()) {
2336            idet=cl->GetDetectorIndex();
2337            const AliITSdetector &detc=layer.GetDetector(idet);
2338            if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2339            track->SetDetectorIndex(idet);
2340            if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2341          }
2342          Int_t cllayer = (idx & 0xf0000000) >> 28;;
2343          Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2344          if (chi2<maxchi2) {
2345            clAcc=cl;
2346            maxchi2=chi2;
2347            modstatus = 1; // found
2348          } else {
2349             return kFALSE; //
2350          }
2351        }
2352      } else { // no cluster in this layer
2353        if (skip==1) {
2354          modstatus = 3; // skipped
2355          // Plane Eff determination:
2356          if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2357            if (IsOKForPlaneEff(track,clusters,ilayer))  // only adequate track for plane eff. evaluation
2358               UseTrackForPlaneEff(track,ilayer);
2359          }
2360        } else {
2361          modstatus = 5; // no cls in road
2362          // check dead
2363          if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2364          dz = 0.5*(zmax-zmin);
2365          dy = 0.5*(ymax-ymin);
2366          Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2367          if (dead==1) modstatus = 7; // holes in z in SPD
2368          if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2369        }
2370      }
2371
2372      if (clAcc) {
2373        if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2374        track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2375      }
2376      track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2377
2378
2379      if (extra) { // search for extra clusters in overlapped modules
2380        AliITStrackV2 tmp(*track);
2381        if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2382        layer.SelectClusters(zmin,zmax,ymin,ymax);
2383
2384        const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2385        Int_t idetExtra=-1;
2386        maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2387        Double_t tolerance=0.1;
2388        while ((clExtra=layer.GetNextCluster(ci))!=0) {
2389          // only clusters in another module! (overlaps)
2390          idetExtra = clExtra->GetDetectorIndex();
2391          if (idet == idetExtra) continue;
2392
2393          const AliITSdetector &detx=layer.GetDetector(idetExtra);
2394
2395          if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2396          if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2397          if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2398          if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2399
2400          Double_t chi2=tmp.GetPredictedChi2(clExtra);
2401          if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2402        }
2403        if (cci>=0) {
2404          track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2405          track->SetExtraModule(ilayer,idetExtra);
2406        }
2407      } // end search for extra clusters in overlapped modules
2408
2409      // Correct for material of the current layer
2410      // cross material
2411      // add time if going outward
2412      if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2413      track->SetCheckInvariant(kTRUE);
2414   } // end loop on layers
2415
2416   if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2417
2418   return kTRUE;
2419 }
2420 //------------------------------------------------------------------------
2421 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2422 {
2423   //
2424   // calculate normalized chi2
2425   //  return NormalizedChi2(track,0);
2426   Float_t chi2 = 0;
2427   Float_t sum=0;
2428   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2429   //  track->fdEdxMismatch=0;
2430   Float_t dedxmismatch =0;
2431   Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2432   if (mode<100){
2433     for (Int_t i = 0;i<6;i++){
2434       if (track->GetClIndex(i)>=0){
2435         Float_t cerry, cerrz;
2436         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2437         else
2438           { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2439         cerry*=cerry;
2440         cerrz*=cerrz;
2441         Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2442         if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2443           Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2444           if (ratio<0.5) {
2445             cchi2+=(0.5-ratio)*10.;
2446             //track->fdEdxMismatch+=(0.5-ratio)*10.;
2447             dedxmismatch+=(0.5-ratio)*10.;
2448           }
2449         }
2450         if (i<2 ||i>3){
2451           AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2452           Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2453           if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2454           if (i<2) chi2+=2*cl->GetDeltaProbability();
2455         }
2456         chi2+=cchi2;
2457         sum++;
2458       }
2459     }
2460     if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2461       track->SetdEdxMismatch(dedxmismatch);
2462     }
2463   }
2464   else{
2465     for (Int_t i = 0;i<4;i++){
2466       if (track->GetClIndex(i)>=0){
2467         Float_t cerry, cerrz;
2468         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2469         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2470         cerry*=cerry;
2471         cerrz*=cerrz;
2472         chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2473         chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2474         sum++;
2475       }
2476     }
2477     for (Int_t i = 4;i<6;i++){
2478       if (track->GetClIndex(i)>=0){
2479         Float_t cerry, cerrz;
2480         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2481         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2482         cerry*=cerry;
2483         cerrz*=cerrz;
2484         Float_t cerryb, cerrzb;
2485         if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2486         else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2487         cerryb*=cerryb;
2488         cerrzb*=cerrzb;
2489         chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2490         chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2491         sum++;
2492       }
2493     }
2494   }
2495   if (track->GetESDtrack()->GetTPCsignal()>85){
2496     Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2497     if (ratio<0.5) {
2498       chi2+=(0.5-ratio)*5.;
2499     }
2500     if (ratio>2){
2501       chi2+=(ratio-2.0)*3;
2502     }
2503   }
2504   //
2505   Double_t match = TMath::Sqrt(track->GetChi22());
2506   if (track->GetConstrain())  match/=track->GetNumberOfClusters();
2507   if (!track->GetConstrain()) {
2508     if (track->GetNumberOfClusters()>2) {
2509       match/=track->GetNumberOfClusters()-2.;
2510     } else {
2511       match=0;
2512     }
2513   }
2514   if (match<0) match=0;
2515
2516   // penalty factor for missing points (NDeadZone>0), but no penalty
2517   // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2518   // or there is a dead from OCDB)
2519   Float_t deadzonefactor = 0.;
2520   if (track->GetNDeadZone()>0.) {
2521     Int_t sumDeadZoneProbability=0;
2522     for(Int_t ilay=0;ilay<6;ilay++) {
2523       if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2524     }
2525     Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2526     if(nDeadZoneWithProbNot1>0) {
2527       Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
2528       AliDebug(2,Form("nDeadZone %f sumDZProbability %f nDZWithProbNot1 %f deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
2529       deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2530       Float_t one = 1.;
2531       deadZoneProbability = TMath::Min(deadZoneProbability,one);
2532       deadzonefactor = 3.*(1.1-deadZoneProbability);
2533     }
2534   }
2535
2536   Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2537     (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2538                                 1./(1.+track->GetNSkipped()));
2539   AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
2540   AliDebug(2,Form("NormChi2 %f  cls %d\n",normchi2,track->GetNumberOfClusters()));
2541   return normchi2;
2542 }
2543 //------------------------------------------------------------------------
2544 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2545 {
2546   //
2547   // return matching chi2 between two tracks
2548   Double_t largeChi2=1000.;
2549
2550   AliITStrackMI track3(*track2);
2551   if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2552   TMatrixD vec(5,1);
2553   vec(0,0)=track1->GetY()   - track3.GetY();
2554   vec(1,0)=track1->GetZ()   - track3.GetZ();
2555   vec(2,0)=track1->GetSnp() - track3.GetSnp();
2556   vec(3,0)=track1->GetTgl() - track3.GetTgl();
2557   vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2558   //
2559   TMatrixD cov(5,5);
2560   cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2561   cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2562   cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2563   cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2564   cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2565
2566   cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2567   cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2568   cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2569   cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2570   //
2571   cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2572   cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2573   cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2574   //
2575   cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2576   cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2577   //
2578   cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2579
2580   cov.Invert();
2581   TMatrixD vec2(cov,TMatrixD::kMult,vec);
2582   TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2583   return chi2(0,0);
2584 }
2585 //------------------------------------------------------------------------
2586 Double_t  AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2587 {
2588   //
2589   //  return probability that given point (characterized by z position and error)
2590   //  is in SPD dead zone
2591   //     This method assumes that fSPDdetzcentre is ordered from -z to +z
2592   //
2593   Double_t probability = 0.;
2594   Double_t nearestz = 0.,distz=0.;
2595   Int_t    nearestzone = -1;
2596   Double_t mindistz = 1000.;
2597
2598   // find closest dead zone
2599   for (Int_t i=0; i<3; i++) {
2600     distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2601     if (distz<mindistz) {
2602       nearestzone=i;
2603       nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2604       mindistz=distz;
2605     }
2606   }
2607
2608   // too far from dead zone
2609   if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2610
2611
2612   Double_t zmin, zmax;
2613   if (nearestzone==0) { // dead zone at z = -7
2614     zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2615     zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2616   } else if (nearestzone==1) { // dead zone at z = 0
2617     zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2618     zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2619   } else if (nearestzone==2) { // dead zone at z = +7
2620     zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2621     zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2622   } else {
2623     zmin = 0.;
2624     zmax = 0.;
2625   }
2626   // probability that the true z is in the range [zmin,zmax] (i.e. inside
2627   // dead zone)
2628   probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2629                       AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2630   AliDebug(2,Form("zpos %f +- %f nearestzone %d  zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
2631   return probability;
2632 }
2633 //------------------------------------------------------------------------
2634 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2635 {
2636   //
2637   // calculate normalized chi2
2638   Float_t chi2[6];
2639   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2640   Float_t ncl = 0;
2641   for (Int_t i = 0;i<6;i++){
2642     if (TMath::Abs(track->GetDy(i))>0){
2643       chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2644       chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2645       ncl++;
2646     }
2647     else{chi2[i]=10000;}
2648   }
2649   Int_t index[6];
2650   TMath::Sort(6,chi2,index,kFALSE);
2651   Float_t max = float(ncl)*fac-1.;
2652   Float_t sumchi=0, sumweight=0;
2653   for (Int_t i=0;i<max+1;i++){
2654     Float_t weight = (i<max)?1.:(max+1.-i);
2655     sumchi+=weight*chi2[index[i]];
2656     sumweight+=weight;
2657   }
2658   Double_t normchi2 = sumchi/sumweight;
2659   return normchi2;
2660 }
2661 //------------------------------------------------------------------------
2662 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2663 {
2664   //
2665   // calculate normalized chi2
2666   //  if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2667   Int_t npoints = 0;
2668   Double_t res =0;
2669   for (Int_t i=0;i<6;i++){
2670     if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2671     Double_t sy1 = forwardtrack->GetSigmaY(i);
2672     Double_t sz1 = forwardtrack->GetSigmaZ(i);
2673     Double_t sy2 = backtrack->GetSigmaY(i);
2674     Double_t sz2 = backtrack->GetSigmaZ(i);
2675     if (i<2){ sy2=1000.;sz2=1000;}
2676     //
2677     Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2678     Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2679     //
2680     Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2681     Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2682     //
2683     res+= nz0*nz0+ny0*ny0;
2684     npoints++;
2685   }
2686   if (npoints>1) return
2687                    TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2688                    //2*forwardtrack->fNUsed+
2689                    res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2690                                   1./(1.+forwardtrack->GetNSkipped()));
2691   return 1000;
2692 }
2693 //------------------------------------------------------------------------
2694 Float_t  *AliITStrackerMI::GetWeight(Int_t index) {
2695   //--------------------------------------------------------------------
2696   //       Return pointer to a given cluster
2697   //--------------------------------------------------------------------
2698   Int_t l=(index & 0xf0000000) >> 28;
2699   Int_t c=(index & 0x0fffffff) >> 00;
2700   return fgLayers[l].GetWeight(c);
2701 }
2702 //------------------------------------------------------------------------
2703 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2704 {
2705   //---------------------------------------------
2706   // register track to the list
2707   //
2708   if (track->GetESDtrack()->GetKinkIndex(0)!=0) return;  //don't register kink tracks
2709   //
2710   //
2711   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2712     Int_t index = track->GetClusterIndex(icluster);
2713     Int_t l=(index & 0xf0000000) >> 28;
2714     Int_t c=(index & 0x0fffffff) >> 00;
2715     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2716     for (Int_t itrack=0;itrack<4;itrack++){
2717       if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2718         fgLayers[l].SetClusterTracks(itrack,c,id);
2719         break;
2720       }
2721     }
2722   }
2723 }
2724 //------------------------------------------------------------------------
2725 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2726 {
2727   //---------------------------------------------
2728   // unregister track from the list
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)==id){
2736         fgLayers[l].SetClusterTracks(itrack,c,-1);
2737       }
2738     }
2739   }
2740 }
2741 //------------------------------------------------------------------------
2742 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2743 {
2744   //-------------------------------------------------------------
2745   //get number of shared clusters
2746   //-------------------------------------------------------------
2747   Float_t shared=0;
2748   for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2749   // mean  number of clusters
2750   Float_t *ny = GetNy(id), *nz = GetNz(id);
2751
2752
2753   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2754     Int_t index = track->GetClusterIndex(icluster);
2755     Int_t l=(index & 0xf0000000) >> 28;
2756     Int_t c=(index & 0x0fffffff) >> 00;
2757     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2758     if (ny[l]<1.e-13){
2759       printf("problem\n");
2760     }
2761     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2762     Float_t weight=1;
2763     //
2764     Float_t deltan = 0;
2765     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2766     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2767       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2768     if (l<2 || l>3){
2769       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2770     }
2771     else{
2772       deltan = (cl->GetNz()-nz[l]);
2773     }
2774     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2775     weight = 2./TMath::Max(3.+deltan,2.);
2776     //
2777     for (Int_t itrack=0;itrack<4;itrack++){
2778       if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2779         list[l]=index;
2780         clist[l] = (AliITSRecPoint*)GetCluster(index);
2781         shared+=weight;
2782         break;
2783       }
2784     }
2785   }
2786   track->SetNUsed(shared);
2787   return shared;
2788 }
2789 //------------------------------------------------------------------------
2790 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2791 {
2792   //
2793   // find first shared track
2794   //
2795   // mean  number of clusters
2796   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2797   //
2798   for (Int_t i=0;i<6;i++) overlist[i]=-1;
2799   Int_t sharedtrack=100000;
2800   Int_t tracks[24],trackindex=0;
2801   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2802   //
2803   for (Int_t icluster=0;icluster<6;icluster++){
2804     if (clusterlist[icluster]<0) continue;
2805     Int_t index = clusterlist[icluster];
2806     Int_t l=(index & 0xf0000000) >> 28;
2807     Int_t c=(index & 0x0fffffff) >> 00;
2808     if (ny[l]<1.e-13){
2809       printf("problem\n");
2810     }
2811     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2812     //if (l>3) continue;
2813     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2814     //
2815     Float_t deltan = 0;
2816     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2817     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2818       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2819     if (l<2 || l>3){
2820       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2821     }
2822     else{
2823       deltan = (cl->GetNz()-nz[l]);
2824     }
2825     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2826     //
2827     for (Int_t itrack=3;itrack>=0;itrack--){
2828       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2829       if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2830        tracks[trackindex]  = fgLayers[l].GetClusterTracks(itrack,c);
2831        trackindex++;
2832       }
2833     }
2834   }
2835   if (trackindex==0) return -1;
2836   if (trackindex==1){
2837     sharedtrack = tracks[0];
2838   }else{
2839     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2840     else{
2841       //
2842       Int_t tracks2[24], cluster[24];
2843       for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2844       Int_t index =0;
2845       //
2846       for (Int_t i=0;i<trackindex;i++){
2847         if (tracks[i]<0) continue;
2848         tracks2[index] = tracks[i];
2849         cluster[index]++;
2850         for (Int_t j=i+1;j<trackindex;j++){
2851           if (tracks[j]<0) continue;
2852           if (tracks[j]==tracks[i]){
2853             cluster[index]++;
2854             tracks[j]=-1;
2855           }
2856         }
2857         index++;
2858       }
2859       Int_t max=0;
2860       for (Int_t i=0;i<index;i++){
2861         if (cluster[index]>max) {
2862           sharedtrack=tracks2[index];
2863           max=cluster[index];
2864         }
2865       }
2866     }
2867   }
2868
2869   if (sharedtrack>=100000) return -1;
2870   //
2871   // make list of overlaps
2872   shared =0;
2873   for (Int_t icluster=0;icluster<6;icluster++){
2874     if (clusterlist[icluster]<0) continue;
2875     Int_t index = clusterlist[icluster];
2876     Int_t l=(index & 0xf0000000) >> 28;
2877     Int_t c=(index & 0x0fffffff) >> 00;
2878     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2879     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2880     if (l==0 || l==1){
2881       if (cl->GetNy()>2) continue;
2882       if (cl->GetNz()>2) continue;
2883     }
2884     if (l==4 || l==5){
2885       if (cl->GetNy()>3) continue;
2886       if (cl->GetNz()>3) continue;
2887     }
2888     //
2889     for (Int_t itrack=3;itrack>=0;itrack--){
2890       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2891       if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2892         overlist[l]=index;
2893         shared++;
2894       }
2895     }
2896   }
2897   return sharedtrack;
2898 }
2899 //------------------------------------------------------------------------
2900 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2901   //
2902   // try to find track hypothesys without conflicts
2903   // with minimal chi2;
2904   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2905   Int_t entries1 = arr1->GetEntriesFast();
2906   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2907   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2908   Int_t entries2 = arr2->GetEntriesFast();
2909   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2910   //
2911   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2912   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2913   if (track10->Pt()>0.5+track20->Pt()) return track10;
2914
2915   for (Int_t itrack=0;itrack<entries1;itrack++){
2916     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2917     UnRegisterClusterTracks(track,trackID1);
2918   }
2919   //
2920   for (Int_t itrack=0;itrack<entries2;itrack++){
2921     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2922     UnRegisterClusterTracks(track,trackID2);
2923   }
2924   Int_t index1=0;
2925   Int_t index2=0;
2926   Float_t maxconflicts=6;
2927   Double_t maxchi2 =1000.;
2928   //
2929   // get weight of hypothesys - using best hypothesy
2930   Double_t w1,w2;
2931
2932   Int_t list1[6],list2[6];
2933   AliITSRecPoint *clist1[6], *clist2[6] ;
2934   RegisterClusterTracks(track10,trackID1);
2935   RegisterClusterTracks(track20,trackID2);
2936   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2937   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2938   UnRegisterClusterTracks(track10,trackID1);
2939   UnRegisterClusterTracks(track20,trackID2);
2940   //
2941   // normalized chi2
2942   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2943   Float_t nerry[6],nerrz[6];
2944   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2945   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2946   for (Int_t i=0;i<6;i++){
2947      if ( (erry1[i]>0) && (erry2[i]>0)) {
2948        nerry[i] = TMath::Min(erry1[i],erry2[i]);
2949        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2950      }else{
2951        nerry[i] = TMath::Max(erry1[i],erry2[i]);
2952        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2953      }
2954      if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2955        chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2956        chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2957        ncl1++;
2958      }
2959      if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2960        chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2961        chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2962        ncl2++;
2963      }
2964   }
2965   chi21/=ncl1;
2966   chi22/=ncl2;
2967   //
2968   //
2969   Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2970   Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2971   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2972   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2973   //
2974   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2975         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2976         +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2977         );
2978   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2979         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2980         +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2981         );
2982
2983   Double_t sumw = w1+w2;
2984   w1/=sumw;
2985   w2/=sumw;
2986   if (w1<w2*0.5) {
2987     w1 = (d2+0.5)/(d1+d2+1.);
2988     w2 = (d1+0.5)/(d1+d2+1.);
2989   }
2990   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2991   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2992   //
2993   // get pair of "best" hypothesys
2994   //
2995   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2996   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2997
2998   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2999     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3000     //if (track1->fFakeRatio>0) continue;
3001     RegisterClusterTracks(track1,trackID1);
3002     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3003       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3004
3005       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3006       //if (track2->fFakeRatio>0) continue;
3007       Float_t nskipped=0;
3008       RegisterClusterTracks(track2,trackID2);
3009       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3010       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3011       UnRegisterClusterTracks(track2,trackID2);
3012       //
3013       if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3014       if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3015       if (nskipped>0.5) continue;
3016       //
3017       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3018       if (conflict1+1<cconflict1) continue;
3019       if (conflict2+1<cconflict2) continue;
3020       Float_t conflict=0;
3021       Float_t sumchi2=0;
3022       Float_t sum=0;
3023       for (Int_t i=0;i<6;i++){
3024         //
3025         Float_t c1 =0.; // conflict coeficients
3026         Float_t c2 =0.;
3027         if (clist1[i]&&clist2[i]){
3028           Float_t deltan = 0;
3029           if (i<2 || i>3){
3030             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3031           }
3032           else{
3033             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3034           }
3035           c1  = 2./TMath::Max(3.+deltan,2.);
3036           c2  = 2./TMath::Max(3.+deltan,2.);
3037         }
3038         else{
3039           if (clist1[i]){
3040             Float_t deltan = 0;
3041             if (i<2 || i>3){
3042               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3043             }
3044             else{
3045               deltan = (clist1[i]->GetNz()-nz1[i]);
3046             }
3047             c1  = 2./TMath::Max(3.+deltan,2.);
3048             c2  = 0;
3049           }
3050
3051           if (clist2[i]){
3052             Float_t deltan = 0;
3053             if (i<2 || i>3){
3054               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3055             }
3056             else{
3057               deltan = (clist2[i]->GetNz()-nz2[i]);
3058             }
3059             c2  = 2./TMath::Max(3.+deltan,2.);
3060             c1  = 0;
3061           }
3062         }
3063         //
3064         chi21=0;chi22=0;
3065         if (TMath::Abs(track1->GetDy(i))>0.) {
3066           chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3067             (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3068           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3069           //  (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3070         }else{
3071           if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3072         }
3073         //
3074         if (TMath::Abs(track2->GetDy(i))>0.) {
3075           chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3076             (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3077           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3078           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3079         }
3080         else{
3081           if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3082         }
3083         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3084         if (chi21>0) sum+=w1;
3085         if (chi22>0) sum+=w2;
3086         conflict+=(c1+c2);
3087       }
3088       Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3089       if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3090       Double_t normchi2 = 2*conflict+sumchi2/sum;
3091       if ( normchi2 <maxchi2 ){
3092         index1 = itrack1;
3093         index2 = itrack2;
3094         maxconflicts = conflict;
3095         maxchi2 = normchi2;
3096       }
3097     }
3098     UnRegisterClusterTracks(track1,trackID1);
3099   }
3100   //
3101   //  if (maxconflicts<4 && maxchi2<th0){
3102   if (maxchi2<th0*2.){
3103     Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3104     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3105     track1->SetChi2MIP(5,maxconflicts);
3106     track1->SetChi2MIP(6,maxchi2);
3107     track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3108     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
3109     track1->SetChi2MIP(8,index1);
3110     fBestTrackIndex[trackID1] =index1;
3111     UpdateESDtrack(track1, AliESDtrack::kITSin);
3112   }
3113   else if (track10->GetChi2MIP(0)<th1){
3114     track10->SetChi2MIP(5,maxconflicts);
3115     track10->SetChi2MIP(6,maxchi2);
3116     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
3117     UpdateESDtrack(track10,AliESDtrack::kITSin);
3118   }
3119
3120   for (Int_t itrack=0;itrack<entries1;itrack++){
3121     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3122     UnRegisterClusterTracks(track,trackID1);
3123   }
3124   //
3125   for (Int_t itrack=0;itrack<entries2;itrack++){
3126     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3127     UnRegisterClusterTracks(track,trackID2);
3128   }
3129
3130   if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3131       &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3132     //  if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3133   //    &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3134     RegisterClusterTracks(track10,trackID1);
3135   }
3136   if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3137       &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3138     //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3139     //  &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3140     RegisterClusterTracks(track20,trackID2);
3141   }
3142   return track10;
3143
3144 }
3145 //------------------------------------------------------------------------
3146 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3147   //--------------------------------------------------------------------
3148   // This function marks clusters assigned to the track
3149   //--------------------------------------------------------------------
3150   AliTracker::UseClusters(t,from);
3151
3152   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3153   //if (c->GetQ()>2) c->Use();
3154   if (c->GetSigmaZ2()>0.1) c->Use();
3155   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3156   //if (c->GetQ()>2) c->Use();
3157   if (c->GetSigmaZ2()>0.1) c->Use();
3158
3159 }
3160 //------------------------------------------------------------------------
3161 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3162 {
3163   //------------------------------------------------------------------
3164   // add track to the list of hypothesys
3165   //------------------------------------------------------------------
3166
3167   if (esdindex>=fTrackHypothesys.GetEntriesFast())
3168     fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3169   //
3170   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3171   if (!array) {
3172     array = new TObjArray(10);
3173     fTrackHypothesys.AddAt(array,esdindex);
3174   }
3175   array->AddLast(track);
3176 }
3177 //------------------------------------------------------------------------
3178 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3179 {
3180   //-------------------------------------------------------------------
3181   // compress array of track hypothesys
3182   // keep only maxsize best hypothesys
3183   //-------------------------------------------------------------------
3184   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3185   if (! (fTrackHypothesys.At(esdindex)) ) return;
3186   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3187   Int_t entries = array->GetEntriesFast();
3188   //
3189   //- find preliminary besttrack as a reference
3190   Float_t minchi2=10000;
3191   Int_t maxn=0;
3192   AliITStrackMI * besttrack=0;
3193   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3194     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3195     if (!track) continue;
3196     Float_t chi2 = NormalizedChi2(track,0);
3197     //
3198     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3199     track->SetLabel(tpcLabel);
3200     CookdEdx(track);
3201     track->SetFakeRatio(1.);
3202     CookLabel(track,0.); //For comparison only
3203     //
3204     //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3205     if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3206       if (track->GetNumberOfClusters()<maxn) continue;
3207       maxn = track->GetNumberOfClusters();
3208       if (chi2<minchi2){
3209         minchi2=chi2;
3210         besttrack=track;
3211       }
3212     }
3213     else{
3214       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3215         delete array->RemoveAt(itrack);
3216       }
3217     }
3218   }
3219   if (!besttrack) return;
3220   //
3221   //
3222   //take errors of best track as a reference
3223   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3224   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3225   for (Int_t j=0;j<6;j++) {
3226     if (besttrack->GetClIndex(j)>=0){
3227       erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3228       errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3229       ny[j]   = besttrack->GetNy(j);
3230       nz[j]   = besttrack->GetNz(j);
3231     }
3232   }
3233   //
3234   // calculate normalized chi2
3235   //
3236   Float_t * chi2        = new Float_t[entries];
3237   Int_t * index         = new Int_t[entries];
3238   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3239   for (Int_t itrack=0;itrack<entries;itrack++){
3240     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3241     if (track){
3242       AliDebug(2,Form("track %d  ncls %d\n",itrack,track->GetNumberOfClusters()));
3243       track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3244       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3245         chi2[itrack] = track->GetChi2MIP(0);
3246       else{
3247         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3248           delete array->RemoveAt(itrack);
3249         }
3250       }
3251     }
3252   }
3253   //
3254   TMath::Sort(entries,chi2,index,kFALSE);
3255   besttrack = (AliITStrackMI*)array->At(index[0]);
3256   if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3257   if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3258     for (Int_t j=0;j<6;j++){
3259       if (besttrack->GetClIndex(j)>=0){
3260         erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3261         errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3262         ny[j]   = besttrack->GetNy(j);
3263         nz[j]   = besttrack->GetNz(j);
3264       }
3265     }
3266   }
3267   //
3268   // calculate one more time with updated normalized errors
3269   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3270   for (Int_t itrack=0;itrack<entries;itrack++){
3271     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3272     if (track){
3273       track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3274       AliDebug(2,Form("track %d  ncls %d\n",itrack,track->GetNumberOfClusters()));
3275       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3276         chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3277       else
3278         {
3279           if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3280             delete array->RemoveAt(itrack);
3281           }
3282         }
3283     }
3284   }
3285   entries = array->GetEntriesFast();
3286   //
3287   //
3288   if (entries>0){
3289     TObjArray * newarray = new TObjArray();
3290     TMath::Sort(entries,chi2,index,kFALSE);
3291     besttrack = (AliITStrackMI*)array->At(index[0]);
3292     if (besttrack){
3293       AliDebug(2,Form("ncls best track %d     %f   %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3294       //
3295       for (Int_t j=0;j<6;j++){
3296         if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3297           erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3298           errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3299           ny[j]   = besttrack->GetNy(j);
3300           nz[j]   = besttrack->GetNz(j);
3301         }
3302       }
3303       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3304       minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3305       Float_t minn = besttrack->GetNumberOfClusters()-3;
3306       Int_t accepted=0;
3307       for (Int_t i=0;i<entries;i++){
3308         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3309         if (!track) continue;
3310         if (accepted>maxcut) break;
3311         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3312         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3313           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3314             delete array->RemoveAt(index[i]);
3315             continue;
3316           }
3317         }
3318         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3319         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3320           if (!shortbest) accepted++;
3321           //
3322           newarray->AddLast(array->RemoveAt(index[i]));
3323           for (Int_t j=0;j<6;j++){
3324             if (nz[j]==0){
3325               erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3326               errz[j] = track->GetSigmaZ(j); errz[j]   = track->GetSigmaZ(j+6);
3327               ny[j]   = track->GetNy(j);
3328               nz[j]   = track->GetNz(j);
3329             }
3330           }
3331         }
3332         else{
3333           delete array->RemoveAt(index[i]);
3334         }
3335       }
3336       array->Delete();
3337       delete fTrackHypothesys.RemoveAt(esdindex);
3338       fTrackHypothesys.AddAt(newarray,esdindex);
3339     }
3340     else{
3341       array->Delete();
3342       delete fTrackHypothesys.RemoveAt(esdindex);
3343     }
3344   }
3345   delete [] chi2;
3346   delete [] index;
3347 }
3348 //------------------------------------------------------------------------
3349 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3350 {
3351   //-------------------------------------------------------------
3352   // try to find best hypothesy
3353   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3354   //-------------------------------------------------------------
3355   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3356   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3357   if (!array) return 0;
3358   Int_t entries = array->GetEntriesFast();
3359   if (!entries) return 0;
3360   Float_t minchi2 = 100000;
3361   AliITStrackMI * besttrack=0;
3362   //
3363   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
3364   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3365   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3366   Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3367   //
3368   for (Int_t i=0;i<entries;i++){
3369     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3370     if (!track) continue;
3371     Float_t sigmarfi,sigmaz;
3372     GetDCASigma(track,sigmarfi,sigmaz);
3373     track->SetDnorm(0,sigmarfi);
3374     track->SetDnorm(1,sigmaz);
3375     //
3376     track->SetChi2MIP(1,1000000);
3377     track->SetChi2MIP(2,1000000);
3378     track->SetChi2MIP(3,1000000);
3379     //
3380     // backtrack
3381     backtrack = new(backtrack) AliITStrackMI(*track);
3382     if (track->GetConstrain()) {
3383       if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3384       if (!backtrack->Improve(0,xyzVtx,ersVtx))         continue;
3385       backtrack->ResetCovariance(10.);
3386     }else{
3387       backtrack->ResetCovariance(10.);
3388     }
3389     backtrack->ResetClusters();
3390
3391     Double_t x = original->GetX();
3392     if (!RefitAt(x,backtrack,track)) continue;
3393     //
3394     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3395     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3396     if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.)  continue;
3397     track->SetChi22(GetMatchingChi2(backtrack,original));
3398
3399     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
3400     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
3401     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
3402
3403
3404     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1))  continue;
3405     //
3406     //forward track - without constraint
3407     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3408     forwardtrack->ResetClusters();
3409     x = track->GetX();
3410     RefitAt(x,forwardtrack,track);
3411     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3412     if  (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0)  continue;
3413     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2))  continue;
3414
3415     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3416     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3417     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
3418     forwardtrack->SetD(0,track->GetD(0));
3419     forwardtrack->SetD(1,track->GetD(1));
3420     {
3421       Int_t list[6];
3422       AliITSRecPoint* clist[6];
3423       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3424       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3425     }
3426
3427     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3428     if  ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3429     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3430       track->SetChi2MIP(3,1000);
3431       continue;
3432     }
3433     Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3434     //
3435     for (Int_t ichi=0;ichi<5;ichi++){
3436       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3437     }
3438     if (chi2 < minchi2){
3439       //besttrack = new AliITStrackMI(*forwardtrack);
3440       besttrack = track;
3441       besttrack->SetLabel(track->GetLabel());
3442       besttrack->SetFakeRatio(track->GetFakeRatio());
3443       minchi2   = chi2;
3444       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3445       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3446       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
3447     }
3448   }
3449   delete backtrack;
3450   delete forwardtrack;
3451   Int_t accepted=0;
3452   for (Int_t i=0;i<entries;i++){
3453     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3454
3455     if (!track) continue;
3456
3457     if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3458         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3459         track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3460       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3461         delete array->RemoveAt(i);
3462         continue;
3463       }
3464     }
3465     else{
3466       accepted++;
3467     }
3468   }
3469   //
3470   array->Compress();
3471   SortTrackHypothesys(esdindex,checkmax,1);
3472
3473   array = (TObjArray*) fTrackHypothesys.At(esdindex);
3474   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3475   besttrack = (AliITStrackMI*)array->At(0);
3476   if (!besttrack)  return 0;
3477   besttrack->SetChi2MIP(8,0);
3478   fBestTrackIndex[esdindex]=0;
3479   entries = array->GetEntriesFast();
3480   AliITStrackMI *longtrack =0;
3481   minchi2 =1000;
3482   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3483   for (Int_t itrack=entries-1;itrack>0;itrack--) {
3484     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3485     if (!track->GetConstrain()) continue;
3486     if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3487     if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3488     if (track->GetChi2MIP(0)>4.) continue;
3489     minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3490     longtrack =track;
3491   }
3492   //if (longtrack) besttrack=longtrack;
3493
3494   Int_t list[6];
3495   AliITSRecPoint * clist[6];
3496   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3497   if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3498       &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3499     RegisterClusterTracks(besttrack,esdindex);
3500   }
3501   //
3502   //
3503   if (shared>0.0){
3504     Int_t nshared;
3505     Int_t overlist[6];
3506     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3507     if (sharedtrack>=0){
3508       //
3509       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3510       if (besttrack){
3511         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3512       }
3513       else return 0;
3514     }
3515   }
3516
3517   if (shared>2.5) return 0;
3518   if (shared>1.0) return besttrack;
3519   //
3520   // Don't sign clusters if not gold track
3521   //
3522   if (!besttrack->IsGoldPrimary()) return besttrack;
3523   if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
3524   //
3525   if (fConstraint[fPass]){
3526     //
3527     // sign clusters
3528     //
3529     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3530     for (Int_t i=0;i<6;i++){
3531       Int_t index = besttrack->GetClIndex(i);
3532       if (index<0) continue;
3533       Int_t ilayer =  (index & 0xf0000000) >> 28;
3534       if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3535       AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3536       if (!c) continue;
3537       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3538       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3539       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
3540       if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3541         if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3542       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
3543
3544       Bool_t cansign = kTRUE;
3545       for (Int_t itrack=0;itrack<entries; itrack++){
3546         AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3547         if (!track) continue;
3548         if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3549         if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3550           cansign = kFALSE;
3551           break;
3552         }
3553       }
3554       if (cansign){
3555         if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3556         if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3557         if (!c->IsUsed()) c->Use();
3558       }
3559     }
3560   }
3561   return besttrack;
3562 }
3563 //------------------------------------------------------------------------
3564 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3565 {
3566   //
3567   // get "best" hypothesys
3568   //
3569
3570   Int_t nentries = itsTracks.GetEntriesFast();
3571   for (Int_t i=0;i<nentries;i++){
3572     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3573     if (!track) continue;
3574     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3575     if (!array) continue;
3576     if (array->GetEntriesFast()<=0) continue;
3577     //
3578     AliITStrackMI* longtrack=0;
3579     Float_t minn=0;
3580     Float_t maxchi2=1000;
3581     for (Int_t j=0;j<array->GetEntriesFast();j++){
3582       AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3583       if (!trackHyp) continue;
3584       if (trackHyp->GetGoldV0()) {
3585         longtrack = trackHyp;   //gold V0 track taken
3586         break;
3587       }
3588       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3589       Float_t chi2 = trackHyp->GetChi2MIP(0);
3590       if (fAfterV0){
3591         if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3592       }
3593       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3594       //
3595       if (chi2 > maxchi2) continue;
3596       minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3597       maxchi2 = chi2;
3598       longtrack=trackHyp;
3599     }
3600     //
3601     //
3602     //
3603     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3604     if (!longtrack) {longtrack = besttrack;}
3605     else besttrack= longtrack;
3606     //
3607     if (besttrack) {
3608       Int_t list[6];
3609       AliITSRecPoint * clist[6];
3610       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3611       //
3612       track->SetNUsed(shared);
3613       track->SetNSkipped(besttrack->GetNSkipped());
3614       track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3615       if (shared>0) {
3616         if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3617         Int_t nshared;
3618         Int_t overlist[6];
3619         //
3620         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3621         //if (sharedtrack==-1) sharedtrack=0;
3622         if (sharedtrack>=0) {
3623           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3624         }
3625       }
3626       if (besttrack&&fAfterV0) {
3627         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3628       }
3629       if (besttrack&&fConstraint[fPass])
3630         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3631       if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3632         if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3633              TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3634       }
3635
3636     }
3637   }
3638 }
3639 //------------------------------------------------------------------------
3640 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3641   //--------------------------------------------------------------------
3642   //This function "cooks" a track label. If label<0, this track is fake.
3643   //--------------------------------------------------------------------
3644   Int_t tpcLabel=-1;
3645
3646   if (track->GetESDtrack()){
3647     tpcLabel = track->GetESDtrack()->GetTPCLabel();
3648     ULong_t trStatus=track->GetESDtrack()->GetStatus();
3649     if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
3650   }
3651    track->SetChi2MIP(9,0);
3652    Int_t nwrong=0;
3653    for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3654      Int_t cindex = track->GetClusterIndex(i);
3655      Int_t l=(cindex & 0xf0000000) >> 28;
3656      AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3657      Int_t isWrong=1;
3658      for (Int_t ind=0;ind<3;ind++){
3659        if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
3660        //AliDebug(2,Form("icl %d  ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3661      }
3662      track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3663      nwrong+=isWrong;
3664    }
3665    Int_t nclusters = track->GetNumberOfClusters();
3666    if (nclusters > 0) //PH Some tracks don't have any cluster
3667      track->SetFakeRatio(double(nwrong)/double(nclusters));
3668    if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
3669      track->SetLabel(-tpcLabel);
3670    } else {
3671      track->SetLabel(tpcLabel);
3672    }
3673    AliDebug(2,Form(" nls %d wrong %d  label %d  tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
3674
3675 }
3676 //------------------------------------------------------------------------
3677 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
3678   //
3679   // Fill the dE/dx in this track
3680   //
3681   track->SetChi2MIP(9,0);
3682   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3683     Int_t cindex = track->GetClusterIndex(i);
3684     Int_t l=(cindex & 0xf0000000) >> 28;
3685     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3686     Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3687     Int_t isWrong=1;
3688     for (Int_t ind=0;ind<3;ind++){
3689       if (cl->GetLabel(ind)==lab) isWrong=0;
3690     }
3691     track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3692   }
3693   Double_t low=0.;
3694   Double_t up=0.51;
3695   track->CookdEdx(low,up);
3696 }
3697 //------------------------------------------------------------------------
3698 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3699   //
3700   // Create some arrays
3701   //
3702   if (fCoefficients) delete []fCoefficients;
3703   fCoefficients = new Float_t[ntracks*48];
3704   for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3705 }
3706 //------------------------------------------------------------------------
3707 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3708 {
3709   //
3710   // Compute predicted chi2
3711   //
3712   Float_t erry,errz,covyz;
3713   Float_t theta = track->GetTgl();
3714   Float_t phi   = track->GetSnp();
3715   phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
3716   AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
3717   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()));
3718   // Take into account the mis-alignment (bring track to cluster plane)
3719   Double_t xTrOrig=track->GetX();
3720   if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
3721   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()));
3722   Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
3723   // Bring the track back to detector plane in ideal geometry
3724   // [mis-alignment will be accounted for in UpdateMI()]
3725   if (!track->Propagate(xTrOrig)) return 1000.;
3726   Float_t ny,nz;
3727   AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3728   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3729   if (delta>1){
3730     chi2+=0.5*TMath::Min(delta/2,2.);
3731     chi2+=2.*cluster->GetDeltaProbability();
3732   }
3733   //
3734   track->SetNy(layer,ny);
3735   track->SetNz(layer,nz);
3736   track->SetSigmaY(layer,erry);
3737   track->SetSigmaZ(layer, errz);
3738   track->SetSigmaYZ(layer,covyz);
3739   //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3740   track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
3741   return chi2;
3742
3743 }
3744 //------------------------------------------------------------------------
3745 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3746 {
3747   //
3748   // Update ITS track
3749   //
3750   Int_t layer = (index & 0xf0000000) >> 28;
3751   track->SetClIndex(layer, index);
3752   if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3753     if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3754       chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3755       track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3756     }
3757   }
3758
3759   if (TMath::Abs(cl->GetQ())<1.e-13) return 0;  // ingore the "virtual" clusters
3760
3761
3762   // Take into account the mis-alignment (bring track to cluster plane)
3763   Double_t xTrOrig=track->GetX();
3764   Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3765   AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3766   AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3767   AliDebug(2,Form(" xtr %f  xcl %f",track->GetX(),cl->GetX()));
3768
3769   if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
3770
3771   AliCluster c(*cl);
3772   c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3773   c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3774   c.SetSigmaYZ(track->GetSigmaYZ(layer));
3775
3776
3777   Int_t updated = track->UpdateMI(&c,chi2,index);
3778   // Bring the track back to detector plane in ideal geometry
3779   if (!track->Propagate(xTrOrig)) return 0;
3780
3781   if(!updated) AliDebug(2,"update failed");
3782   return updated;
3783 }
3784
3785 //------------------------------------------------------------------------
3786 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3787 {
3788   //
3789   //DCA sigmas parameterization
3790   //to be paramterized using external parameters in future
3791   //
3792   //
3793   sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3794   sigmaz   = 0.0110+4.37*TMath::Abs(track->GetC());
3795 }
3796 //------------------------------------------------------------------------
3797 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
3798 {
3799   //
3800   // Clusters from delta electrons?
3801   //
3802   Int_t entries = clusterArray->GetEntriesFast();
3803   if (entries<4) return;
3804   AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
3805   Int_t layer = cluster->GetLayer();
3806   if (layer>1) return;
3807   Int_t index[10000];
3808   Int_t ncandidates=0;
3809   Float_t r = (layer>0)? 7:4;
3810   //
3811   for (Int_t i=0;i<entries;i++){
3812     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
3813     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3814     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3815     index[ncandidates] = i;  //candidate to belong to delta electron track
3816     ncandidates++;
3817     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3818       cl0->SetDeltaProbability(1);
3819     }
3820   }
3821   //
3822   //
3823   //
3824   for (Int_t i=0;i<ncandidates;i++){
3825     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
3826     if (cl0->GetDeltaProbability()>0.8) continue;
3827     //
3828     Int_t ncl = 0;
3829     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3830     sumy=sumz=sumy2=sumyz=sumw=0.0;
3831     for (Int_t j=0;j<ncandidates;j++){
3832       if (i==j) continue;
3833       AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
3834       //
3835       Float_t dz = cl0->GetZ()-cl1->GetZ();
3836       Float_t dy = cl0->GetY()-cl1->GetY();
3837       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3838         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3839         y[ncl] = cl1->GetY();
3840         z[ncl] = cl1->GetZ();
3841         sumy+= y[ncl]*weight;
3842         sumz+= z[ncl]*weight;
3843         sumy2+=y[ncl]*y[ncl]*weight;
3844         sumyz+=y[ncl]*z[ncl]*weight;
3845         sumw+=weight;
3846         ncl++;
3847       }
3848     }
3849     if (ncl<4) continue;
3850     Float_t det = sumw*sumy2  - sumy*sumy;
3851     Float_t delta=1000;
3852     if (TMath::Abs(det)>0.01){
3853       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
3854       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
3855       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3856     }
3857     else{
3858       Float_t z0  = sumyz/sumy;
3859       delta = TMath::Abs(cl0->GetZ()-z0);
3860     }
3861     if ( delta<0.05) {
3862       cl0->SetDeltaProbability(1-20.*delta);
3863     }
3864   }
3865 }
3866 //------------------------------------------------------------------------
3867 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3868 {
3869   //
3870   // Update ESD track
3871   //
3872   track->UpdateESDtrack(flags);
3873   AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3874   if (oldtrack) delete oldtrack;
3875   track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3876   if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3877     printf("Problem\n");
3878   }
3879 }
3880 //------------------------------------------------------------------------
3881 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3882   //
3883   // Get nearest upper layer close to the point xr.
3884   // rough approximation
3885   //
3886   const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3887   Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3888   Int_t res =6;
3889   for (Int_t i=0;i<6;i++){
3890     if (radius<kRadiuses[i]){
3891       res =i;
3892       break;
3893     }
3894   }
3895   return res;
3896 }
3897 //------------------------------------------------------------------------
3898 void AliITStrackerMI::BuildMaterialLUT(TString material) {
3899   //--------------------------------------------------------------------
3900   // Fill a look-up table with mean material
3901   //--------------------------------------------------------------------
3902
3903   Int_t n=1000;
3904   Double_t mparam[7];
3905   Double_t point1[3],point2[3];
3906   Double_t phi,cosphi,sinphi,z;
3907   // 0-5 layers, 6 pipe, 7-8 shields
3908   Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
3909   Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
3910
3911   Int_t ifirst=0,ilast=0;
3912   if(material.Contains("Pipe")) {
3913     ifirst=6; ilast=6;
3914   } else if(material.Contains("Shields")) {
3915     ifirst=7; ilast=8;
3916   } else if(material.Contains("Layers")) {
3917     ifirst=0; ilast=5;
3918   } else {
3919     Error("BuildMaterialLUT","Wrong layer name\n");
3920   }
3921
3922   for(Int_t imat=ifirst; imat<=ilast; imat++) {
3923     Double_t param[5]={0.,0.,0.,0.,0.};
3924     for (Int_t i=0; i<n; i++) {
3925       phi = 2.*TMath::Pi()*gRandom->Rndm();
3926       cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
3927       z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
3928       point1[0] = rmin[imat]*cosphi;
3929       point1[1] = rmin[imat]*sinphi;
3930       point1[2] = z;
3931       point2[0] = rmax[imat]*cosphi;
3932       point2[1] = rmax[imat]*sinphi;
3933       point2[2] = z;
3934       AliTracker::MeanMaterialBudget(point1,point2,mparam);
3935       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
3936     }
3937     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
3938     if(imat<=5) {
3939       fxOverX0Layer[imat] = param[1];
3940       fxTimesRhoLayer[imat] = param[0]*param[4];
3941     } else if(imat==6) {
3942       fxOverX0Pipe = param[1];
3943       fxTimesRhoPipe = param[0]*param[4];
3944     } else if(imat==7) {
3945       fxOverX0Shield[0] = param[1];
3946       fxTimesRhoShield[0] = param[0]*param[4];
3947     } else if(imat==8) {
3948       fxOverX0Shield[1] = param[1];
3949       fxTimesRhoShield[1] = param[0]*param[4];
3950     }
3951   }
3952   /*
3953   printf("%s\n",material.Data());
3954   printf("%f  %f\n",fxOverX0Pipe,fxTimesRhoPipe);
3955   printf("%f  %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
3956   printf("%f  %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
3957   printf("%f  %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
3958   printf("%f  %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
3959   printf("%f  %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
3960   printf("%f  %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
3961   printf("%f  %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
3962   printf("%f  %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
3963   */
3964   return;
3965 }
3966 //------------------------------------------------------------------------
3967 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
3968                                               TString direction) {
3969   //-------------------------------------------------------------------
3970   // Propagate beyond beam pipe and correct for material
3971   // (material budget in different ways according to fUseTGeo value)
3972   // Add time if going outward (PropagateTo or PropagateToTGeo)
3973   //-------------------------------------------------------------------
3974
3975   // Define budget mode:
3976   // 0: material from AliITSRecoParam (hard coded)
3977   // 1: material from TGeo in one step (on the fly)
3978   // 2: material from lut
3979   // 3: material from TGeo in one step (same for all hypotheses)
3980   Int_t mode;
3981   switch(fUseTGeo) {
3982   case 0:
3983     mode=0;
3984     break;
3985   case 1:
3986     mode=1;
3987     break;
3988   case 2:
3989     mode=2;
3990     break;
3991   case 3:
3992     if(fTrackingPhase.Contains("Clusters2Tracks"))
3993       { mode=3; } else { mode=1; }
3994     break;
3995   case 4:
3996     if(fTrackingPhase.Contains("Clusters2Tracks"))
3997       { mode=3; } else { mode=2; }
3998     break;
3999   default:
4000     mode=0;
4001     break;
4002   }
4003   if(fTrackingPhase.Contains("Default")) mode=0;
4004
4005   Int_t index=fCurrentEsdTrack;
4006
4007   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4008   Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4009   Double_t xToGo;
4010   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4011
4012   Double_t xOverX0,x0,lengthTimesMeanDensity;
4013
4014   switch(mode) {
4015   case 0:
4016     xOverX0 = AliITSRecoParam::GetdPipe();
4017     x0 = AliITSRecoParam::GetX0Be();
4018     lengthTimesMeanDensity = xOverX0*x0;
4019     lengthTimesMeanDensity *= dir;
4020     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4021     break;
4022   case 1:
4023     if (!t->PropagateToTGeo(xToGo,1)) return 0;
4024     break;
4025   case 2:
4026     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4027     xOverX0 = fxOverX0Pipe;
4028     lengthTimesMeanDensity = fxTimesRhoPipe;
4029     lengthTimesMeanDensity *= dir;
4030     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4031     break;
4032   case 3:
4033     if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4034     if(fxOverX0PipeTrks[index]<0) {
4035       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4036       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4037                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4038       fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4039       fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4040       return 1;
4041     }
4042     xOverX0 = fxOverX0PipeTrks[index];
4043     lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4044     lengthTimesMeanDensity *= dir;
4045     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4046     break;
4047   }
4048
4049   return 1;
4050 }
4051 //------------------------------------------------------------------------
4052 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4053                                                 TString shield,
4054                                                 TString direction) {
4055   //-------------------------------------------------------------------
4056   // Propagate beyond SPD or SDD shield and correct for material
4057   // (material budget in different ways according to fUseTGeo value)
4058   // Add time if going outward (PropagateTo or PropagateToTGeo)
4059   //-------------------------------------------------------------------
4060
4061   // Define budget mode:
4062   // 0: material from AliITSRecoParam (hard coded)
4063   // 1: material from TGeo in steps of X cm (on the fly)
4064   //    X = AliITSRecoParam::GetStepSizeTGeo()
4065   // 2: material from lut
4066   // 3: material from TGeo in one step (same for all hypotheses)
4067   Int_t mode;
4068   switch(fUseTGeo) {
4069   case 0:
4070     mode=0;
4071     break;
4072   case 1:
4073     mode=1;
4074     break;
4075   case 2:
4076     mode=2;
4077     break;
4078   case 3:
4079     if(fTrackingPhase.Contains("Clusters2Tracks"))
4080       { mode=3; } else { mode=1; }
4081     break;
4082   case 4:
4083     if(fTrackingPhase.Contains("Clusters2Tracks"))
4084       { mode=3; } else { mode=2; }
4085     break;
4086   default:
4087     mode=0;
4088     break;
4089   }
4090   if(fTrackingPhase.Contains("Default")) mode=0;
4091
4092   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4093   Double_t rToGo;
4094   Int_t    shieldindex=0;
4095   if (shield.Contains("SDD")) { // SDDouter
4096     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4097     shieldindex=1;
4098   } else if (shield.Contains("SPD")) {        // SPDouter
4099     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4100     shieldindex=0;
4101   } else {
4102     Error("CorrectForShieldMaterial"," Wrong shield name\n");
4103     return 0;
4104   }
4105
4106   // do nothing if we are already beyond the shield
4107   Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4108   if(dir<0 && rTrack > rToGo) return 1; // going outward
4109   if(dir>0 && rTrack < rToGo) return 1; // going inward
4110
4111
4112   Double_t xToGo;
4113   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4114
4115   Int_t index=2*fCurrentEsdTrack+shieldindex;
4116
4117   Double_t xOverX0,x0,lengthTimesMeanDensity;
4118   Int_t nsteps=1;
4119
4120   switch(mode) {
4121   case 0:
4122     xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4123     x0 = AliITSRecoParam::GetX0shield(shieldindex);
4124     lengthTimesMeanDensity = xOverX0*x0;
4125     lengthTimesMeanDensity *= dir;
4126     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4127     break;
4128   case 1:
4129     nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4130     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4131     break;
4132   case 2:
4133     if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4134     xOverX0 = fxOverX0Shield[shieldindex];
4135     lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4136     lengthTimesMeanDensity *= dir;
4137     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4138     break;
4139   case 3:
4140     if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4141     if(fxOverX0ShieldTrks[index]<0) {
4142       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4143       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4144                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4145       fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4146       fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4147       return 1;
4148     }
4149     xOverX0 = fxOverX0ShieldTrks[index];
4150     lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4151     lengthTimesMeanDensity *= dir;
4152     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4153     break;
4154   }
4155
4156   return 1;
4157 }
4158 //------------------------------------------------------------------------
4159 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4160                                                Int_t layerindex,
4161                                                Double_t oldGlobXYZ[3],
4162                                                TString direction) {
4163   //-------------------------------------------------------------------
4164   // Propagate beyond layer and correct for material
4165   // (material budget in different ways according to fUseTGeo value)
4166   // Add time if going outward (PropagateTo or PropagateToTGeo)
4167   //-------------------------------------------------------------------
4168
4169   // Define budget mode:
4170   // 0: material from AliITSRecoParam (hard coded)
4171   // 1: material from TGeo in stepsof X cm (on the fly)
4172   //    X = AliITSRecoParam::GetStepSizeTGeo()
4173   // 2: material from lut
4174   // 3: material from TGeo in one step (same for all hypotheses)
4175   Int_t mode;
4176   switch(fUseTGeo) {
4177   case 0:
4178     mode=0;
4179     break;
4180   case 1:
4181     mode=1;
4182     break;
4183   case 2:
4184     mode=2;
4185     break;
4186   case 3:
4187     if(fTrackingPhase.Contains("Clusters2Tracks"))
4188       { mode=3; } else { mode=1; }
4189     break;
4190   case 4:
4191     if(fTrackingPhase.Contains("Clusters2Tracks"))
4192       { mode=3; } else { mode=2; }
4193     break;
4194   default:
4195     mode=0;
4196     break;
4197   }
4198   if(fTrackingPhase.Contains("Default")) mode=0;
4199
4200   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4201
4202   Double_t r=fgLayers[layerindex].GetR();
4203   Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4204
4205   Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4206   Double_t xToGo;
4207   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4208
4209   Int_t index=6*fCurrentEsdTrack+layerindex;
4210
4211
4212   Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4213   Int_t nsteps=1;
4214
4215   // back before material (no correction)
4216   Double_t rOld,xOld;
4217   rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4218   if (!t->GetLocalXat(rOld,xOld)) return 0;
4219   if (!t->Propagate(xOld)) return 0;
4220
4221   switch(mode) {
4222   case 0:
4223     xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4224     lengthTimesMeanDensity = xOverX0*x0;
4225     lengthTimesMeanDensity *= dir;
4226     // Bring the track beyond the material
4227     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4228     break;
4229   case 1:
4230     nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4231     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4232     break;
4233   case 2:
4234     if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4235     xOverX0 = fxOverX0Layer[layerindex];
4236     lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4237     lengthTimesMeanDensity *= dir;
4238     // Bring the track beyond the material
4239     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4240     break;
4241   case 3:
4242     if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4243     if(fxOverX0LayerTrks[index]<0) {
4244       nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4245       if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4246       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4247                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4248       fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4249       fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4250       return 1;
4251     }
4252     xOverX0 = fxOverX0LayerTrks[index];
4253     lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4254     lengthTimesMeanDensity *= dir;
4255     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4256     break;
4257   }
4258
4259
4260   return 1;
4261 }
4262 //------------------------------------------------------------------------
4263 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4264   //-----------------------------------------------------------------
4265   // Initialize LUT for storing material for each prolonged track
4266   //-----------------------------------------------------------------
4267   fxOverX0PipeTrks = new Float_t[ntracks];
4268   fxTimesRhoPipeTrks = new Float_t[ntracks];
4269   fxOverX0ShieldTrks = new Float_t[ntracks*2];
4270   fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4271   fxOverX0LayerTrks = new Float_t[ntracks*6];
4272   fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4273
4274   for(Int_t i=0; i<ntracks; i++) {
4275     fxOverX0PipeTrks[i] = -1.;
4276     fxTimesRhoPipeTrks[i] = -1.;
4277   }
4278   for(Int_t j=0; j<ntracks*2; j++) {
4279     fxOverX0ShieldTrks[j] = -1.;
4280     fxTimesRhoShieldTrks[j] = -1.;
4281   }
4282   for(Int_t k=0; k<ntracks*6; k++) {
4283     fxOverX0LayerTrks[k] = -1.;
4284     fxTimesRhoLayerTrks[k] = -1.;
4285   }
4286
4287   fNtracks = ntracks;
4288
4289   return;
4290 }
4291 //------------------------------------------------------------------------
4292 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4293   //-----------------------------------------------------------------
4294   // Delete LUT for storing material for each prolonged track
4295   //-----------------------------------------------------------------
4296   if(fxOverX0PipeTrks) {
4297     delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4298   }
4299   if(fxOverX0ShieldTrks) {
4300     delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4301   }
4302
4303   if(fxOverX0LayerTrks) {
4304     delete [] fxOverX0LayerTrks;  fxOverX0LayerTrks = 0;
4305   }
4306   if(fxTimesRhoPipeTrks) {
4307     delete [] fxTimesRhoPipeTrks;  fxTimesRhoPipeTrks = 0;
4308   }
4309   if(fxTimesRhoShieldTrks) {
4310     delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4311   }
4312   if(fxTimesRhoLayerTrks) {
4313     delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4314   }
4315   return;
4316 }
4317 //------------------------------------------------------------------------
4318 void AliITStrackerMI::SetForceSkippingOfLayer() {
4319   //-----------------------------------------------------------------
4320   // Check if we are forced to skip layers
4321   // either we set to skip them in RecoParam
4322   // or they were off during data-taking
4323   //-----------------------------------------------------------------
4324
4325   const AliEventInfo *eventInfo = GetEventInfo();
4326
4327   for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4328     fForceSkippingOfLayer[l] = 0;
4329     // check reco param
4330     if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4331     // check run info
4332
4333     if(eventInfo &&
4334        AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4335       AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4336       if(l==0 || l==1)  {
4337         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4338       } else if(l==2 || l==3) {
4339         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4340       } else {
4341         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4342       }
4343     }
4344   }
4345   return;
4346 }
4347 //------------------------------------------------------------------------
4348 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4349                                       Int_t ilayer,Int_t idet) const {
4350   //-----------------------------------------------------------------
4351   // This method is used to decide whether to allow a prolongation
4352   // without clusters, because we want to skip the layer.
4353   // In this case the return value is > 0:
4354   // return 1: the user requested to skip a layer
4355   // return 2: track outside z acceptance
4356   //-----------------------------------------------------------------
4357
4358   if (ForceSkippingOfLayer(ilayer)) return 1;
4359
4360   Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4361
4362   if (idet<0 &&  // out in z
4363       ilayer>innerLayCanSkip &&
4364       AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4365     // check if track will cross SPD outer layer
4366     Double_t phiAtSPD2,zAtSPD2;
4367     if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4368       if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4369     }
4370     return 2; // always allow skipping, changed on 05.11.2009
4371   }
4372
4373   return 0;
4374 }
4375 //------------------------------------------------------------------------
4376 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4377                                      Int_t ilayer,Int_t idet,
4378                                      Double_t dz,Double_t dy,
4379                                      Bool_t noClusters) const {
4380   //-----------------------------------------------------------------
4381   // This method is used to decide whether to allow a prolongation
4382   // without clusters, because there is a dead zone in the road.
4383   // In this case the return value is > 0:
4384   // return 1: dead zone at z=0,+-7cm in SPD
4385   //     This method assumes that fSPDdetzcentre is ordered from -z to +z
4386   // return 2: all road is "bad" (dead or noisy) from the OCDB
4387   // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4388   // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4389   //-----------------------------------------------------------------
4390
4391   // check dead zones at z=0,+-7cm in the SPD
4392   if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4393     Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4394                           fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4395                           fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4396     Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4397                           fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4398                           fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4399     for (Int_t i=0; i<3; i++)
4400       if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4401         AliDebug(2,Form("crack SPD %d track z %f   %f   %f  %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4402         if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4403       }
4404   }
4405
4406   // check bad zones from OCDB
4407   if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4408
4409   if (idet<0) return 0;
4410
4411   AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4412
4413   Int_t detType=-1;
4414   Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4415   if (ilayer==0 || ilayer==1) {        // ----------  SPD
4416     detType = 0;
4417   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
4418     detType = 1;
4419     detSizeFactorX *= 2.;
4420   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
4421     detType = 2;
4422   }
4423   AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4424   if (detType==2) segm->SetLayer(ilayer+1);
4425   Float_t detSizeX = detSizeFactorX*segm->Dx();
4426   Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4427
4428   // check if the road overlaps with bad chips
4429   Float_t xloc,zloc;
4430   LocalModuleCoord(ilayer,idet,track,xloc,zloc);
4431   Float_t zlocmin = zloc-dz;
4432   Float_t zlocmax = zloc+dz;
4433   Float_t xlocmin = xloc-dy;
4434   Float_t xlocmax = xloc+dy;
4435   Int_t chipsInRoad[100];
4436
4437   // check if road goes out of detector
4438   Bool_t touchNeighbourDet=kFALSE;
4439   if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4440   if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4441   if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4442   if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4443   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));
4444
4445   // check if this detector is bad
4446   if (det.IsBad()) {
4447     AliDebug(2,Form("lay %d  bad detector %d",ilayer,idet));
4448     if(!touchNeighbourDet) {
4449       return 2; // all detectors in road are bad
4450     } else {
4451       return 3; // at least one is bad
4452     }
4453   }
4454
4455   Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4456   AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4457   if (!nChipsInRoad) return 0;
4458
4459   Bool_t anyBad=kFALSE,anyGood=kFALSE;
4460   for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4461     if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4462     AliDebug(2,Form("  chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4463     if (det.IsChipBad(chipsInRoad[iCh])) {
4464       anyBad=kTRUE;
4465     } else {
4466       anyGood=kTRUE;
4467     }
4468   }
4469
4470   if (!anyGood) {
4471     if(!touchNeighbourDet) {
4472       AliDebug(2,"all bad in road");
4473       return 2;  // all chips in road are bad
4474     } else {
4475       return 3; // at least a bad chip in road
4476     }
4477   }
4478
4479   if (anyBad) {
4480     AliDebug(2,"at least a bad in road");
4481     return 3; // at least a bad chip in road
4482   }
4483
4484
4485   if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4486       || !noClusters) return 0;
4487
4488   // There are no clusters in road: check if there is at least
4489   // a bad SPD pixel or SDD anode or SSD strips on both sides
4490
4491   Int_t idetInITS=idet;
4492   for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4493
4494   if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4495     AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4496     return 4;
4497   }
4498   //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4499
4500   return 0;
4501 }
4502 //------------------------------------------------------------------------
4503 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4504                                        const AliITStrackMI *track,
4505                                        Float_t &xloc,Float_t &zloc) const {
4506   //-----------------------------------------------------------------
4507   // Gives position of track in local module ref. frame
4508   //-----------------------------------------------------------------
4509
4510   xloc=0.;
4511   zloc=0.;
4512
4513   if(idet<0) return kTRUE; // track out of z acceptance of layer
4514
4515   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4516
4517   Int_t lad = Int_t(idet/ndet) + 1;
4518
4519   Int_t det = idet - (lad-1)*ndet + 1;
4520
4521   Double_t xyzGlob[3],xyzLoc[3];
4522
4523   AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4524   // take into account the misalignment: xyz at real detector plane
4525   if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4526
4527   if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4528
4529   xloc = (Float_t)xyzLoc[0];
4530   zloc = (Float_t)xyzLoc[2];
4531
4532   return kTRUE;
4533 }
4534 //------------------------------------------------------------------------
4535 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4536 //
4537 // Method to be optimized further:
4538 // Aim: decide whether a track can be used for PlaneEff evaluation
4539 //      the decision is taken based on the track quality at the layer under study
4540 //      no information on the clusters on this layer has to be used
4541 //      The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4542 //      the cut is done on number of sigmas from the boundaries
4543 //
4544 //  Input: Actual track, layer [0,5] under study
4545 //  Output: none
4546 //  Return: kTRUE if this is a good track
4547 //
4548 // it will apply a pre-selection to obtain good quality tracks.
4549 // Here also  you will have the possibility to put a control on the
4550 // impact point of the track on the basic block, in order to exclude border regions
4551 // this will be done by calling a proper method of the AliITSPlaneEff class.
4552 //
4553 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4554 // return: Bool_t   -> kTRUE if usable track, kFALSE if not usable.
4555 //
4556   Int_t index[AliITSgeomTGeo::kNLayers];
4557   Int_t k;
4558   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4559   //
4560   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4561     index[k]=clusters[k];
4562   }
4563
4564   if(!fPlaneEff)
4565     {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4566   AliITSlayer &layer=fgLayers[ilayer];
4567   Double_t r=layer.GetR();
4568   AliITStrackMI tmp(*track);
4569
4570 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4571   Int_t ncl=0;
4572   for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) {
4573     AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4574                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4575     if (tmp.GetClIndex(lay)>=0) ncl++;
4576   }
4577   Bool_t nextout = kFALSE;
4578   if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4579   else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4580   Bool_t nextin = kFALSE;
4581   if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4582   else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4583   if(ncl<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4584      return kFALSE;
4585   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout)  return kFALSE;
4586   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin)   return kFALSE;
4587   if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4588  //  if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff()  && !tmp.GetConstrain()) return kFALSE;
4589
4590 // detector number
4591   Double_t phi,z;
4592   if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4593   Int_t idet=layer.FindDetectorIndex(phi,z);
4594   if(idet<0) { AliInfo(Form("cannot find detector"));
4595     return kFALSE;}
4596
4597   // here check if it has good Chi Square.
4598
4599   //propagate to the intersection with the detector plane
4600   const AliITSdetector &det=layer.GetDetector(idet);
4601   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4602
4603   Float_t locx; //
4604   Float_t locz; //
4605   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4606   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4607   if(key>fPlaneEff->Nblock()) return kFALSE;
4608   Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4609   if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4610   //***************
4611   // DEFINITION OF SEARCH ROAD FOR accepting a track
4612   //
4613   //For the time being they are hard-wired, later on from AliITSRecoParam
4614   // Double_t nsigx=AliITSRecoParam::GetNSigXFarFromBoundary();
4615   // Double_t nsigz=AliITSRecoParam::GetNSigZFarFromBoundary();
4616   Double_t nsigz=4;
4617   Double_t nsigx=4;
4618   Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
4619   Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
4620                                                 // done for RecPoints
4621
4622   // exclude tracks at boundary between detectors
4623   //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4624   Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4625   AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4626   AliDebug(2,Form("Local:    track impact x=%f, z=%f",locx,locz));
4627   AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4628
4629   if ( (locx-dx < blockXmn+boundaryWidth) ||
4630        (locx+dx > blockXmx-boundaryWidth) ||
4631        (locz-dz < blockZmn+boundaryWidth) ||
4632        (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4633   return kTRUE;
4634 }
4635 //------------------------------------------------------------------------
4636 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4637 //
4638 // This Method has to be optimized! For the time-being it uses the same criteria
4639 // as those used in the search of extra clusters for overlapping modules.
4640 //
4641 // Method Purpose: estabilish whether a track has produced a recpoint or not
4642 //                 in the layer under study (For Plane efficiency)
4643 //
4644 // inputs: AliITStrackMI* track  (pointer to a usable track)
4645 // outputs: none
4646 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4647 //               traversed by this very track. In details:
4648 //               - if a cluster can be associated to the track then call
4649 //                  AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4650 //               - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4651 //
4652   if(!fPlaneEff)
4653     {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4654   AliITSlayer &layer=fgLayers[ilayer];
4655   Double_t r=layer.GetR();
4656   AliITStrackMI tmp(*track);
4657
4658 // detector number
4659   Double_t phi,z;
4660   if (!tmp.GetPhiZat(r,phi,z)) return;
4661   Int_t idet=layer.FindDetectorIndex(phi,z);
4662
4663   if(idet<0) { AliInfo(Form("cannot find detector"));
4664     return;}
4665
4666
4667 //propagate to the intersection with the detector plane
4668   const AliITSdetector &det=layer.GetDetector(idet);
4669   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
4670
4671
4672 //***************
4673 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
4674 //
4675   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
4676                     TMath::Sqrt(tmp.GetSigmaZ2() +
4677                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4678                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4679                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
4680   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
4681                     TMath::Sqrt(tmp.GetSigmaY2() +
4682                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4683                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4684                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
4685
4686 // road in global (rphi,z) [i.e. in tracking ref. system]
4687   Double_t zmin = tmp.GetZ() - dz;
4688   Double_t zmax = tmp.GetZ() + dz;
4689   Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
4690   Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
4691
4692 // select clusters in road
4693   layer.SelectClusters(zmin,zmax,ymin,ymax);
4694
4695 // Define criteria for track-cluster association
4696   Double_t msz = tmp.GetSigmaZ2() +
4697   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4698   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4699   AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
4700   Double_t msy = tmp.GetSigmaY2() +
4701   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4702   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4703   AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
4704   if (tmp.GetConstrain()) {
4705     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
4706     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
4707   }  else {
4708     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
4709     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
4710   }
4711   msz = 1./msz; // 1/RoadZ^2
4712   msy = 1./msy; // 1/RoadY^2
4713 //
4714
4715   const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
4716   Int_t idetc=-1;
4717   Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
4718   //Double_t  tolerance=0.2;
4719   /*while ((cl=layer.GetNextCluster(clidx))!=0) {
4720     idetc = cl->GetDetectorIndex();
4721     if(idet!=idetc) continue;
4722     //Int_t ilay = cl->GetLayer();
4723
4724     if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
4725     if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
4726
4727     Double_t chi2=tmp.GetPredictedChi2(cl);
4728     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
4729   }*/
4730   Float_t locx; //
4731   Float_t locz; //
4732   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
4733 //
4734   AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
4735   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4736   if(key>fPlaneEff->Nblock()) return;
4737   Bool_t found=kFALSE;
4738   //if (ci>=0) {
4739   Double_t chi2;
4740   while ((cl=layer.GetNextCluster(clidx))!=0) {
4741     idetc = cl->GetDetectorIndex();
4742     if(idet!=idetc) continue;
4743     // here real control to see whether the cluster can be associated to the track.
4744     // cluster not associated to track
4745     if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
4746          (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy   > 1. ) continue;
4747     // calculate track-clusters chi2
4748     chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
4749                                                // in particular, the error associated to the cluster
4750     //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
4751     // chi2 cut
4752     if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
4753     found=kTRUE;
4754     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
4755    // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
4756    // track->SetExtraModule(ilayer,idetExtra);
4757   }
4758   if(!fPlaneEff->UpDatePlaneEff(found,key))
4759        AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
4760   if(fPlaneEff->GetCreateHistos()&&  AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
4761     Float_t tr[4]={99999.,99999.,9999.,9999.};    // initialize to high values
4762     Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
4763     Int_t cltype[2]={-999,-999};
4764
4765     tr[0]=locx;
4766     tr[1]=locz;
4767     tr[2]=TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
4768     tr[3]=TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
4769
4770     if (found){
4771       clu[0]=layer.GetCluster(ci)->GetDetLocalX();
4772       clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
4773       cltype[0]=layer.GetCluster(ci)->GetNy();
4774       cltype[1]=layer.GetCluster(ci)->GetNz();
4775
4776      // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
4777      //  X->50/sqrt(12)=14 micron   Z->450/sqrt(12)= 120 micron)
4778      // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
4779      // It is computed properly by calling the method
4780      // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4781      // T
4782      //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
4783       //if (tmp.PropagateTo(x,0.,0.)) {
4784         chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
4785         AliCluster c(*layer.GetCluster(ci));
4786         c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
4787         c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
4788         //if (layer.GetCluster(ci)->GetGlobalCov(cov))  // by using this, instead, you got nominal cluster errors
4789         clu[2]=TMath::Sqrt(c.GetSigmaY2());
4790         clu[3]=TMath::Sqrt(c.GetSigmaZ2());
4791       //}
4792     }
4793     fPlaneEff->FillHistos(key,found,tr,clu,cltype);
4794   }
4795 return;
4796 }