]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ITS/AliITStrackerMI.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / ITS / AliITStrackerMI.cxx
index b916c2b28c854a8e5a32e85198f20dc572188919..d7df282aa40886f7144165587f3707f1f827661b 100644 (file)
 #include <TRandom.h>
 #include <TTreeStream.h>
 #include <TVector3.h>
+#include <TBits.h>
 
 #include "AliLog.h"
 #include "AliGeomManager.h"
 #include "AliITSPlaneEff.h"
 #include "AliTrackPointArray.h"
 #include "AliESDEvent.h"
+#include "AliESDV0Params.h"
 #include "AliESDtrack.h"
 #include "AliV0.h"
 #include "AliITSChannelStatus.h"
@@ -56,6 +58,7 @@
 #include "AliITSV0Finder.h"
 #include "AliITStrackerMI.h"
 #include "AliMathBase.h"
+#include "AliPID.h"
 
 
 ClassImp(AliITStrackerMI)
@@ -93,6 +96,7 @@ fDebugStreamer(0),
 fITSChannelStatus(0),
 fkDetTypeRec(0),
 fPlaneEff(0),
+fSPDChipIntPlaneEff(0),
 fITSPid(0)
 
  {
@@ -143,6 +147,7 @@ fDebugStreamer(0),
 fITSChannelStatus(0),
 fkDetTypeRec(0),
 fPlaneEff(0),
+fSPDChipIntPlaneEff(0),
 fITSPid(0) {
   //--------------------------------------------------------------------
   //This is the AliITStrackerMI constructor
@@ -265,9 +270,13 @@ fITSPid(0) {
   if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
       AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
     Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
-    if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
+    if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)!=1)
       AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
-    if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
+    if (iplane<2) {
+      fPlaneEff = new AliITSPlaneEffSPD();
+      fSPDChipIntPlaneEff = new Bool_t[AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip];
+      for (UInt_t i=0; i<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip; i++) fSPDChipIntPlaneEff[i]=kFALSE;
+    }
     else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
     else fPlaneEff = new AliITSPlaneEffSSD();
     if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
@@ -353,6 +362,7 @@ AliITStrackerMI::~AliITStrackerMI()
   if(fITSChannelStatus) delete fITSChannelStatus;
   if(fPlaneEff) delete fPlaneEff;
   if(fITSPid) delete fITSPid;
+  if (fSPDChipIntPlaneEff) delete [] fSPDChipIntPlaneEff;
 
 }
 //------------------------------------------------------------------------
@@ -439,8 +449,8 @@ Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
          Float_t q      = 0.; // this identifies virtual clusters
          Float_t hit[6] = {xdead,
                            0.,
-                           AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
-                           AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
+                           static_cast<Float_t>(AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2()),
+                           static_cast<Float_t>(AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2()),
                            q,
                            0.};
          Bool_t local   = kTRUE;
@@ -554,7 +564,6 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
   // temporary
   Int_t noesd = 0;
   {/* Read ESD tracks */
-    Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
     Int_t nentr=event->GetNumberOfTracks();
     noesd=nentr;
     //    Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
@@ -571,9 +580,6 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
       t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP());              //I.B.
       Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
 
-
-      // look at the ESD mass hypothesys !
-      if (t->GetMass()<0.9*pimass) t->SetMass(pimass); 
       // write expected q
       t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 
@@ -698,54 +704,60 @@ Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
   fTrackingPhase="PropagateBack";
   Int_t nentr=event->GetNumberOfTracks();
   //  Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
-
+  double bz0 = GetBz();
+  const double kWatchStep=10.; // for larger steps watch arc vs segment difference
+  //
   Int_t ntrk=0;
   for (Int_t i=0; i<nentr; i++) {
      AliESDtrack *esd=event->GetTrack(i);
 
-     if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
-     if (esd->GetStatus()&AliESDtrack::kITSout) continue;
-
-     AliITStrackMI *t = new AliITStrackMI(*esd);
-
-     t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
-
-     ResetTrackToFollow(*t);
-
-     /*
-     // propagate to vertex [SR, GSI 17.02.2003]
-     // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
-     if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
-       if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
-        fTrackToFollow.StartTimeIntegral();
-       // from vertex to outside pipe
-       CorrectForPipeMaterial(&fTrackToFollow,"outward");
-       }*/
      // Start time integral and add distance from current position to vertex 
+     if (esd->GetStatus()&AliESDtrack::kITSout) continue;
+     AliITStrackMI t(*esd);
      Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
-     fTrackToFollow.GetXYZ(xyzTrk);
+     t.GetXYZ(xyzTrk); 
      Double_t dst2 = 0.;
-     for (Int_t icoord=0; icoord<3; icoord++) {  
-       Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
-       dst2 += di*di; 
+     {
+       double dxs = xyzTrk[0] - xyzVtx[0];
+       double dys = xyzTrk[1] - xyzVtx[1];
+       double dzs = xyzTrk[2] - xyzVtx[2];
+       // RS: for large segment steps d use approximation of cicrular arc s by
+       // s = 2R*asin(d/2R) = d/p asin(p) \approx d/p (p + 1/6 p^3) = d (1+1/6 p^2)
+       // where R is the track radius, p = d/2R = 2C*d (C is the abs curvature)
+       // Hence s^2/d^2 = (1+1/6 p^2)^2
+       dst2 = dxs*dxs + dys*dys;
+       if (dst2 > kWatchStep*kWatchStep) { // correct circular part for arc/segment factor
+        double crv = TMath::Abs(esd->GetC(bz0));
+        double fcarc = 1.+crv*crv*dst2/6.;
+        dst2 *= fcarc*fcarc;
+       }
+       dst2 += dzs*dzs;
      }
-     fTrackToFollow.StartTimeIntegral();
-     fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
+     t.StartTimeIntegral();
+     t.AddTimeStep(TMath::Sqrt(dst2));
+     //
+     // transfer the time integral to ESD track
+     esd->SetStatus(AliESDtrack::kTIME);
+     Double_t times[AliPID::kSPECIESC];
+     t.GetIntegratedTimes(times,AliPID::kSPECIESC); 
+     esd->SetIntegratedTimes(times);
+     esd->SetIntegratedLength(t.GetIntegratedLength());
+     //
+     if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
 
+     t.SetExpQ(TMath::Max(0.8*t.GetESDtrack()->GetTPCsignal(),30.));
+     ResetTrackToFollow(t);
+     //
      fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
-     if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
-        if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
-          delete t;
-          continue;
-        }
-        fTrackToFollow.SetLabel(t->GetLabel());
-        //fTrackToFollow.CookdEdx();
-        CookLabel(&fTrackToFollow,0.); //For comparison only
-        fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
-        //UseClusters(&fTrackToFollow);
-        ntrk++;
+     if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,&t)) {
+       if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) continue;
+       //       fTrackToFollow.SetLabel(t.GetLabel());              // why do we neet this
+       //fTrackToFollow.CookdEdx();
+       CookLabel(&fTrackToFollow,0.); //For comparison only // why do we need this?
+       fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
+       //UseClusters(&fTrackToFollow);
+       ntrk++;
      }
-     delete t;
   }
 
   AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
@@ -771,6 +783,13 @@ Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
   Int_t nentr=event->GetNumberOfTracks();
   //  Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
 
+  // only for PlaneEff and in case of SPD (for FO studies)
+  if( AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
+      AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0 && 
+      AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()<2) {
+      for (UInt_t i=0; i<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip; i++) fSPDChipIntPlaneEff[i]=kFALSE;     
+  }
+
   Int_t ntrk=0;
   for (Int_t i=0; i<nentr; i++) {
     AliESDtrack *esd=event->GetTrack(i);
@@ -813,7 +832,7 @@ Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
        //       fTrackToFollow.CookdEdx();
        CookdEdx(&fTrackToFollow);
 
-       CookLabel(&fTrackToFollow,0.0); //For comparison only
+       CookLabel(&fTrackToFollow,0.0); //For comparison only // RS why do we need this?
 
        //The beam pipe
        if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
@@ -1428,7 +1447,7 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin
   // update TPC V0 information
   //
   if (otrack->GetESDtrack()->GetV0Index(0)>0){    
-    Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
+    Float_t fprimvertex[3]={static_cast<Float_t>(GetX()),static_cast<Float_t>(GetY()),static_cast<Float_t>(GetZ())};
     for (Int_t i=0;i<3;i++){
       Int_t  index = otrack->GetESDtrack()->GetV0Index(i); 
       if (index==0) break;
@@ -1529,7 +1548,7 @@ fNMaxSigmaCl(3)
   fYB[0]=0;
   fYB[1]=0;
 
-  for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer5; j++) {
+  for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
     for (Int_t j1=0; j1<6; j1++) {
       fClusters5[j1][j]=0;
       fClusterIndex5[j1][j]=-1;
@@ -1541,7 +1560,7 @@ fNMaxSigmaCl(3)
     }
   }
 
-  for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer10; j++) {
+  for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
     for (Int_t j1=0; j1<11; j1++) {
       fClusters10[j1][j]=0;
       fClusterIndex10[j1][j]=-1;
@@ -1553,7 +1572,7 @@ fNMaxSigmaCl(3)
     }
   }
 
-  for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer20; j++) {
+  for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
     for (Int_t j1=0; j1<21; j1++) {
       fClusters20[j1][j]=0;
       fClusterIndex20[j1][j]=-1;
@@ -1564,7 +1583,7 @@ fNMaxSigmaCl(3)
       fBy20[j1][1]=0;
     }
   }
-  for(Int_t i=0;i<AliITSRecoParam::fgkMaxClusterPerLayer;i++){
+  for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
     fClusters[i]=NULL;
     fClusterIndex[i]=0;
   }
@@ -1619,7 +1638,7 @@ fNMaxSigmaCl(3) {
   fYB[0]=0;
   fYB[1]=0;
 
 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer5; j++) {
for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
     for (Int_t j1=0; j1<6; j1++) {
       fClusters5[j1][j]=0;
       fClusterIndex5[j1][j]=-1;
@@ -1631,7 +1650,7 @@ fNMaxSigmaCl(3) {
     }
   }
 
-  for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer10; j++) {
+  for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
     for (Int_t j1=0; j1<11; j1++) {
       fClusters10[j1][j]=0;
       fClusterIndex10[j1][j]=-1;
@@ -1643,7 +1662,7 @@ fNMaxSigmaCl(3) {
     }
   }
 
-  for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer20; j++) {
+  for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
     for (Int_t j1=0; j1<21; j1++) {
       fClusters20[j1][j]=0;
       fClusterIndex20[j1][j]=-1;
@@ -1654,7 +1673,7 @@ fNMaxSigmaCl(3) {
       fBy20[j1][1]=0;
     }
   }
-  for(Int_t i=0;i<AliITSRecoParam::fgkMaxClusterPerLayer;i++){
+  for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
     fClusters[i]=NULL;
     fClusterIndex[i]=0;
   }
@@ -3071,7 +3090,7 @@ Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID
     else{
       //
       Int_t tracks2[24], cluster[24];
-      for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
+      for (Int_t i=0;i<24;i++){ tracks2[i]=-1; cluster[i]=0;}
       Int_t index =0;
       //
       for (Int_t i=0;i<trackindex;i++){
@@ -3397,13 +3416,11 @@ void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
   // add track to the list of hypothesys
   //------------------------------------------------------------------
 
-  if (esdindex>=fTrackHypothesys.GetEntriesFast()) 
-    fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
   //
   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
   if (!array) {
     array = new TObjArray(10);
-    fTrackHypothesys.AddAt(array,esdindex);
+    fTrackHypothesys.AddAtAndExpand(array,esdindex);
   }
   array->AddLast(track);
 }
@@ -3898,7 +3915,7 @@ void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
 } 
 
 //------------------------------------------------------------------------
-void AliITStrackerMI::FlagFakes(TObjArray &itsTracks)
+void AliITStrackerMI::FlagFakes(const TObjArray &itsTracks)
 {
   //
   // RS: flag those tracks which are suxpected to have fake clusters
@@ -3968,6 +3985,64 @@ void AliITStrackerMI::FlagFakes(TObjArray &itsTracks)
   for (int i=6;i--;) delete refArr[i];
 }
 
+
+
+//------------------------------------------------------------------------
+void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
+  //--------------------------------------------------------------------
+  //This function "cooks" a track label. If label<0, this track is fake.
+  //--------------------------------------------------------------------
+  const int kMaxLbPerCl = 3;
+  int lbID[36],lbStat[36];
+  Int_t nLab=0, nCl = track->GetNumberOfClusters();
+  //
+  //  track->SetLabel(-1);
+  //  track->SetFakeRatio(0);
+  //
+  for (Int_t i=0;i<nCl;i++) { // fill all labels
+    Int_t cindex = track->GetClusterIndex(i);
+    //    Int_t l=(cindex & 0xf0000000) >> 28;
+    AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
+    //
+    for (int imc=0;imc<kMaxLbPerCl;imc++) { // labels within single cluster
+      int trLb = cl->GetLabel(imc);
+      if (trLb<0) break;
+      // search this mc track in already accounted ones
+      int iLab;
+      for (iLab=0;iLab<nLab;iLab++) if (lbID[iLab]==trLb) break;
+      if (iLab<nLab) lbStat[iLab]++;
+      else {
+       lbID[nLab] = trLb;
+       lbStat[nLab++] = 1;
+      }
+    } // loop over given cluster's labels   
+  } // loop over clusters
+  //
+  if (nLab<1) return; // no labels at all
+  //
+  Int_t tpcLabel=-1; 
+  if (track->GetESDtrack() && track->GetESDtrack()->IsOn(AliESDtrack::kTPCin)){
+    tpcLabel = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
+  }
+  //
+  // find majority label
+  if (nCl && nLab) {
+    int maxLab=0,tpcLabID=-1;
+    for (int ilb=nLab;ilb--;) {
+      int st = lbStat[ilb];
+      if (lbStat[maxLab]<st) maxLab = ilb;
+      if (lbID[ilb] == tpcLabel) tpcLabID = ilb;
+    }
+    // if there is an equal choice, prefer ITS label consistent with TPC label
+    if (tpcLabel>0 && (tpcLabID!=maxLab) && lbStat[maxLab]==lbStat[tpcLabID]) maxLab=tpcLabID;
+                                                                              
+    track->SetFakeRatio(1.-float(lbStat[maxLab])/nCl);
+    track->SetLabel( lbStat[maxLab]>=nCl-wrong ? lbID[maxLab] : -lbID[maxLab]);
+  }
+  //
+}
+
+/*
 //------------------------------------------------------------------------
 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
   //--------------------------------------------------------------------
@@ -4002,9 +4077,10 @@ void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
    } else {
      track->SetLabel(tpcLabel);
    }
-   AliDebug(2,Form(" nls %d wrong %d  label %d  tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
-   
+   AliDebug(2,Form(" nls %d wrong %d  label %d  tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel)); 
 }
+*/
+
 //------------------------------------------------------------------------
 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
   //
@@ -4251,7 +4327,7 @@ void AliITStrackerMI::BuildMaterialLUT(TString material) {
   } else {
     Error("BuildMaterialLUT","Wrong layer name\n");
   }
-
+  const double kAngEps = 1e-4; // tiny slope to avoid tracks strictly normal to Z axis
   for(Int_t imat=ifirst; imat<=ilast; imat++) {
     Double_t param[5]={0.,0.,0.,0.,0.};
     for (Int_t i=0; i<n; i++) {
@@ -4263,8 +4339,9 @@ void AliITStrackerMI::BuildMaterialLUT(TString material) {
       point1[2] = z;
       point2[0] = rmax[imat]*cosphi;
       point2[1] = rmax[imat]*sinphi;
-      point2[2] = z;
+      point2[2] = z+(rmax[imat]-rmin[imat])*kAngEps;
       AliTracker::MeanMaterialBudget(point1,point2,mparam);
+      if (mparam[1]>999) {n--; continue;}  // skip anomalous values in failed propagation
       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
     }
     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
@@ -4785,6 +4862,7 @@ Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
     }
   }
 
+  if(zlocmin>zlocmax)return 0;
   Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
   AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
   if (!nChipsInRoad) return 0;
@@ -4866,7 +4944,7 @@ Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
 }
 //------------------------------------------------------------------------
 //------------------------------------------------------------------------
-Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
+Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer){
 //
 // Method to be optimized further: 
 // Aim: decide whether a track can be used for PlaneEff evaluation
@@ -4902,19 +4980,19 @@ Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t
   AliITStrackMI tmp(*track);
 
 // require a minimal number of cluster in other layers and eventually clusters in closest layers 
-  Int_t ncl_out=0; Int_t ncl_in=0;
+  Int_t nclout=0; Int_t nclin=0;
   for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
  AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
-   // if (tmp.GetClIndex(lay)>=0) ncl_out++;
-if(index[lay]>=0)ncl_out++;
+   // if (tmp.GetClIndex(lay)>=0) nclout++;
+if(index[lay]>=0)nclout++;
   }
   for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
     AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
-   if (index[lay]>=0) ncl_in++; 
+   if (index[lay]>=0) nclin++; 
   }
-  Int_t ncl=ncl_out+ncl_in;
+  Int_t ncl=nclout+nclin;
   Bool_t nextout = kFALSE;
   if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
   else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
@@ -4922,7 +5000,7 @@ if(index[lay]>=0)ncl_out++;
   if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
   else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
   // maximum number of missing clusters allowed in outermost layers
-  if(ncl_out<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff()) 
+  if(nclout<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff()) 
      return kFALSE; 
   // maximum number of missing clusters allowed (both in innermost and in outermost layers)
   if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff()) 
@@ -4957,9 +5035,20 @@ if(index[lay]>=0)ncl_out++;
   //
   Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
   Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
-  Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
-  Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
-                                                    // done for RecPoints
+  Double_t distx=AliITSReconstructor::GetRecoParam()->GetDistXFromBoundaryPlaneEff();
+  Double_t distz=AliITSReconstructor::GetRecoParam()->GetDistZFromBoundaryPlaneEff();
+  Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()) + distx;
+  // those are precisions in the tracking reference system
+  Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()) + distz;
+  // Use it also for the module reference system, as it is done for RecPoints
+
+  if(AliITSReconstructor::GetRecoParam()->GetSwitchOnMaxDistNSigFrmBndPlaneEff()){
+   if(nsigx*TMath::Sqrt(tmp.GetSigmaY2())<=distx) dx -= nsigx*TMath::Sqrt(tmp.GetSigmaY2());
+   else dx -= distx;
+
+   if(nsigz*TMath::Sqrt(tmp.GetSigmaZ2())<=distz) dz -= nsigz*TMath::Sqrt(tmp.GetSigmaZ2());
+   else dz -= distz;
+  }
 
   // exclude tracks at boundary between detectors
   //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
@@ -4971,6 +5060,28 @@ if(index[lay]>=0)ncl_out++;
        (locx+dx > blockXmx-boundaryWidth) ||
        (locz-dz < blockZmn+boundaryWidth) ||
        (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
+
+  if(ilayer==0){
+       const AliESDEvent *myesd = ((AliESDtrack*)tmp.GetESDtrack())->GetESDEvent();
+       //The beam pipe
+       if (CorrectForPipeMaterial(&tmp,"inward")) {
+          const AliESDVertex* vtx = myesd->GetVertex();
+          if(!vtx) return kFALSE;
+          Double_t ddz[2],cov[3];
+          Double_t maxD=3.;
+          if(!tmp.PropagateToDCA(vtx,AliTracker::GetBz(),maxD,ddz,cov)) return kFALSE;
+          if(TMath::Abs(ddz[0])>=AliITSReconstructor::GetRecoParam()->GetDCACutPlaneEff()) return kFALSE;
+
+          Double_t covar[6]; vtx->GetCovMatrix(covar);
+          Double_t p[2]={tmp.GetParameter()[0]-ddz[0],tmp.GetParameter()[1]-ddz[1]};
+          Double_t c[3]={covar[2],0.,covar[5]};
+          Double_t chi2= ((AliESDtrack*)tmp.GetESDtrack())->GetPredictedChi2(p,c);
+          if (chi2>AliITSReconstructor::GetRecoParam()->GetVertexChi2CutPlaneEff()) return kFALSE;  // Use this to cut on chi^2
+
+       } else return kFALSE;
+    }
+
+
   return kTRUE;
 }
 //------------------------------------------------------------------------
@@ -5049,9 +5160,32 @@ void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilay
     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
   }
+
+  if(AliITSReconstructor::GetRecoParam()->GetSwitchOffStdSearchClusPlaneEff()){
+     Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXSearchClusterPlaneEff();
+     Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZSearchClusterPlaneEff();
+     Double_t distx=AliITSReconstructor::GetRecoParam()->GetDistXSearchClusterPlaneEff();
+     Double_t distz=AliITSReconstructor::GetRecoParam()->GetDistZSearchClusterPlaneEff();
+     msy = nsigx*TMath::Sqrt(tmp.GetSigmaY2()) + distx;
+     msz = nsigz*TMath::Sqrt(tmp.GetSigmaZ2()) + distz;
+
+  if(AliITSReconstructor::GetRecoParam()->GetSwitchOnMaxDistNSigSrhClusPlaneEff()){
+     if(nsigx*TMath::Sqrt(tmp.GetSigmaY2())<=distx) msy -= nsigx*TMath::Sqrt(tmp.GetSigmaY2());
+     else msy -= distx;
+
+     if(nsigz*TMath::Sqrt(tmp.GetSigmaZ2())<=distz) msz -= nsigz*TMath::Sqrt(tmp.GetSigmaZ2());
+     else msz -= distz;
+  }
+
+  msy *= msy;
+  msz *= msz;
+
+  }
+
+  if(msz==0 || msy==0){AliWarning("UseTrackForPlaneEff: null search frame"); return;}
+    
   msz = 1./msz; // 1/RoadZ^2
   msy = 1./msy; // 1/RoadY^2
-//
 
   const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
   Int_t idetc=-1;
@@ -5098,13 +5232,36 @@ void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilay
   }
   if(!fPlaneEff->UpDatePlaneEff(found,key))
        AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
+
+// this for FO efficiency studies (only for SPD) // 
+   UInt_t keyFO=999999;
+   Bool_t foundFO=kFALSE;
+   if(ilayer<2){ //ONLY SPD layers for FastOr studies
+    TBits mapFO = fkDetTypeRec->GetFastOrFiredMap();
+    Int_t phase = (fEsd->GetBunchCrossNumber())%4;
+    if(!fSPDChipIntPlaneEff[key]){
+      AliITSPlaneEffSPD spd; 
+      keyFO = spd.SwitchChipKeyNumbering(key);
+      if(mapFO.TestBitNumber(keyFO))foundFO=kTRUE;
+       keyFO = key + (AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip)*(phase+1);
+       if(keyFO<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip) {
+         AliWarning(Form("UseTrackForPlaneEff: too small keyF0 (= %d), setting it to 999999",keyFO));
+         keyFO=999999;
+       }
+       if(!fPlaneEff->UpDatePlaneEff(foundFO,keyFO))
+          AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for FastOR for key=%d",keyFO));
+     }
+  }
+  
+
+
   if(fPlaneEff->GetCreateHistos()&&  AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
     Float_t tr[4]={99999.,99999.,9999.,9999.};    // initialize to high values 
     Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails) 
     Int_t cltype[2]={-999,-999};
                                                           // and the module
 
-Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
+    Float_t angleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
 
     tr[0]=locx;
     tr[1]=locz;
@@ -5174,15 +5331,19 @@ Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute
     if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
     anglet *= 180./TMath::Pi();
 
-     AngleModTrack[2]=(Float_t) angle;
-     AngleModTrack[0]=(Float_t) anglet;
+     angleModTrack[2]=(Float_t) angle;
+     angleModTrack[0]=(Float_t) anglet;
      // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
-    AngleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
-    AngleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
-    AngleModTrack[1]*=180./TMath::Pi(); // in degree
+    angleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
+    angleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
+    angleModTrack[1]*=180./TMath::Pi(); // in degree
+
+    fPlaneEff->FillHistos(key,found,tr,clu,cltype,angleModTrack);
 
-      fPlaneEff->FillHistos(key,found,tr,clu,cltype,AngleModTrack);
+    // For FO efficiency studies of SPD 
+    if(ilayer<2 && !fSPDChipIntPlaneEff[key]) fPlaneEff->FillHistos(keyFO,foundFO,tr,clu,cltype,angleModTrack);
   }
+  if(ilayer<2) fSPDChipIntPlaneEff[key]=kTRUE;
 return;
 }