MuonAcceptor.cxx

Go to the documentation of this file.

// MuonAcceptor.cxx
//
// Implementation of the Muon Acceptor class
//
// Only the accept() method is overridden
//
//
// Authors: S. Hassani, S. Dean
//
//

#include "AtlfastAlgs/MuonAcceptor.h"

#include "CLHEP/Random/JamesRandom.h"
#include "CLHEP/Random/RandGauss.h"
#include "CLHEP/Units/SystemOfUnits.h"

#include <iostream>
#include <fstream>
#include <iostream>
#include <sstream>
#include <cmath>

#include "PathResolver/PathResolver.h"


//-------------------------------------------------------------
namespace Atlfast {

  using std::abs;
  using std::pair;
  using std::make_pair;
  
  //=============================================================
  // muEffdata methods
  //=============================================================
 
  muEffdata::muEffdata(){}
  muEffdata::~muEffdata(){}


  void muEffdata::Print(std::ostream* out){
        *out <<std::setw(12)<<std::setprecision(3)<<GetPt()
         <<std::setw(12)<<std::setprecision(3)<<GetEtaMin()
         <<std::setw(12)<<std::setprecision(3)<<GetEtaMax()
         << std::endl;
  }

  void muEffdata::SetPt (double pt)         { m_pt     = pt ;}
  void muEffdata::SetEtaMin (double etamin) { m_etamin = etamin ;}
  void muEffdata::SetEtaMax (double etamax) { m_etamax = etamax ;}
  void muEffdata::SetPtEtaNum(std::vector<int> ptEtaNum){m_ptEtaNum=ptEtaNum;}

  void muEffdata::SetPtEtaNumeDeno(std::vector<std::pair<int, int> > VecPaire)
  {m_VecPaire = VecPaire;}

  double muEffdata::GetPt()    {return  m_pt    ;}
  double muEffdata::GetEtaMin(){return  m_etamin;}
  double muEffdata::GetEtaMax(){return  m_etamax;}
  int    muEffdata::GetPtEtaNum(std::vector<int> &ptEtaNum){
    
    int SizeOfVectorOfmuEffdata = m_ptEtaNum.size(); 
    ptEtaNum = m_ptEtaNum;
    return SizeOfVectorOfmuEffdata;

  }
  
  void muEffdata::GetPtEtaNumeDeno(double &pt,double &etamin,double &etamax,
                                   std::vector<std::pair<int, int> > &VecPaire){
    
    pt       = m_pt ;
    etamin   = m_etamin ;
    etamax   = m_etamax ;
    VecPaire = m_VecPaire;
    
  }
  
  muEffdata::muEffdata(double &pt,double &etamin,double &etamax,
                       std::vector<std::pair<int, int> > &VecPaire){
    
    pt       = m_pt ;
    etamin   = m_etamin ;
    etamax   = m_etamax ;
    VecPaire = m_VecPaire;

  }
  
  //=============================================================
  // MuonAcceptor methods
  //=============================================================

  MuonAcceptor::MuonAcceptor(int muSmearKey, MsgStream& log) : m_log(log) {
    
    m_log << MSG::DEBUG << "MuonAcceptor::MuonAcceptor" << endreq;
    // Initialize random number generator
    m_randomEngine = new HepJamesRandom(12345);
    
    //Open input file

    std::string filename = muSmearKey == 1 ? 
      "atlfastDatafiles/MuonSpectrometerEff.txt" : 
      "atlfastDatafiles/MuonCombinedEff.txt";
    filename = PathResolver::find_file (filename, "DATAPATH");
    m_log << MSG::DEBUG << "Opening file: " << filename << endreq;
    std::ifstream InputFile;
    InputFile.open(filename.c_str());
    if(!InputFile) { 
      m_log << MSG::ERROR << "Error: file could not be opened" << endreq;
      exit(1);
    }
    // 
    std::vector<muEffdata>VectorOfmuEffdata;
    std::vector<std::pair<int, int> > ptEtaNumeDeno;
    VectorOfmuEffdata.clear();
    ptEtaNumeDeno.clear();
    //
    m_muEffdataObject.clear();
    m_ptValues.clear();
    m_etaminValues.clear();
    m_etamaxValues.clear();


    //Loop on evts
    int Istatus = 1;
    while ( Istatus == 1 ) {
  
      //  Get evt data
      //
      Istatus = muReadEff(InputFile,VectorOfmuEffdata);
      //  Do something
      if (Istatus == 1){

        std::vector<int> ptEtaNum;
        int SizeOfVectorOfmuEffdata = VectorOfmuEffdata.size() ;
        for (int Item=0; Item <SizeOfVectorOfmuEffdata  ; Item++){
          double pt;double etamin;double etamax;
          VectorOfmuEffdata[Item].GetPtEtaNumeDeno(pt,etamin,etamax,ptEtaNumeDeno);
          m_muEffdataObject.push_back(VectorOfmuEffdata[Item]);
          m_ptValues.push_back(pt);
          m_etaminValues.push_back(etamin);
          m_etamaxValues.push_back(etamax);
        }
      }
    }
    //Close input file
    InputFile.close();
    m_log << MSG::DEBUG << "Done with MuonAcceptor constructor" << endreq;
  }

  bool MuonAcceptor::accept( const ReconstructedParticle &particle, MsgStream &log ){
    
    log << MSG::DEBUG << "MuonAcceptor::accept" << endreq;

    // The efficiencies are a function of pT in GeV
    double pTinGeV = particle.pT() / GeV;
 
    double muEff = muEfficiency(pTinGeV,
                                particle.eta(),
                                particle.phi()).first;
    log << MSG::DEBUG << "pT,eta,phi,efficiency = " << pTinGeV
        << "," << particle.eta() << "," 
        << particle.phi() << "," << muEff << endreq;
    double random = m_randomEngine->flat();
    log << MSG::DEBUG << "random = " << random << endreq;
    return ( random < muEff) ? true : false;
  }
  
  //-------------------------------------
  int MuonAcceptor::muReadEff(std::ifstream& InputFile,
                              std::vector<muEffdata>& VectorOfmuEffdata){
    
    VectorOfmuEffdata.clear();
    //
    std::vector<int> ptEtaNume;
    std::vector<int> ptEtaDeno;
 
    std::vector<int> ptEtaNume1;
    std::vector<int> ptEtaDeno1;

    std::vector<int> ptEtaNume2;
    std::vector<int> ptEtaDeno2; 

    std::vector<int> ptEtaNume3;
    std::vector<int> ptEtaDeno3;

    std::vector<int> ptEtaNume4;
    std::vector<int> ptEtaDeno4; 
    //
    std::vector<std::pair<int, int> > VecPaire; 
    std::pair<int, int> mapaire;
    //
    std::vector<std::pair<int, int> > VecPaire1; 
    std::pair<int, int> mapaire1;
    //
    std::vector<std::pair<int, int> > VecPaire2; 
    std::pair<int, int> mapaire2;
    //
    std::vector<std::pair<int, int> > VecPaire3; 
    std::pair<int, int> mapaire3;
    //
    std::vector<std::pair<int, int> > VecPaire4; 
    std::pair<int, int> mapaire4;
    //
    ptEtaNume.clear();
    ptEtaDeno.clear();
    ptEtaNume1.clear();
    ptEtaDeno1.clear();
    ptEtaNume2.clear();
    ptEtaDeno2.clear();
    ptEtaNume3.clear();
    ptEtaDeno3.clear();
    ptEtaNume4.clear();
    ptEtaDeno4.clear();
    //
    muEffdata pmuEffdata;
    //  
    if(!InputFile)      return 0;
    if(InputFile.eof()) return 0;
  
    double pt;
    double etamin,etamax;
    int nBin;
    int nnum,ndeno;  
    //
    if(InputFile.eof()) return 0;
    InputFile >> pt;
    if(InputFile.eof()) return 0;
    InputFile >>etamin>>etamax>>nBin;
    if(InputFile.eof()) return 0;
    //
    pmuEffdata.SetPt(pt);
    pmuEffdata.SetEtaMin(etamin );
    pmuEffdata.SetEtaMax(etamax);
    //
    for (int ibin=0; ibin < nBin ; ibin++){
      if(InputFile.eof()) return 0;
      InputFile >> nnum ;
      ptEtaNume.push_back(nnum); 
    }
    if(InputFile.eof()) return 0;
    for (int ibin=0; ibin < nBin ; ibin++){
      if(InputFile.eof()) return 0;
      InputFile >> ndeno ;
      ptEtaDeno.push_back(ndeno); 
    }
    pmuEffdata.SetPtEtaNum(ptEtaNume);
    pmuEffdata.SetPtEtaNum(ptEtaDeno);
    //
    int SizeOfptEtaNume =ptEtaNume.size();
    for(int i=0;i<SizeOfptEtaNume;i++){
      mapaire.first=ptEtaNume[i];
      mapaire.second = ptEtaDeno[i];
      VecPaire.push_back(mapaire);
    }
    pmuEffdata.SetPtEtaNumeDeno(VecPaire);
   
    //End first eta part
    if(InputFile.eof()) return 0; 
    // 
    VectorOfmuEffdata.push_back(pmuEffdata);
    // Re-Initialization
    int nBin1 =0;
    double etamin1=0;
    double etamax1=0;
    int ndeno1=0;
    int nnum1=0;
    //
    InputFile >>etamin1>>etamax1>>nBin1;
    if(InputFile.eof()) return 0; 
    //
    pmuEffdata.SetEtaMin(etamin1 );
    pmuEffdata.SetEtaMax(etamax1);
    //
    for (int Item=0; Item < nBin1 ; Item++){
      if(InputFile.eof()) return 0;
      InputFile >> nnum1 ;
      ptEtaNume1.push_back(nnum1); 
    }
    if(InputFile.eof()) return 0;
    for (int Item=0; Item < nBin1 ; Item++){
      if(InputFile.eof()) return 0;
      InputFile >> ndeno1 ;
      ptEtaDeno1.push_back(ndeno1);
    }
    pmuEffdata.SetPtEtaNum(ptEtaNume1);
    pmuEffdata.SetPtEtaNum(ptEtaDeno1);
    //  
    //End second eta part
    if(InputFile.eof()) return 0;
    //
    int SizeOfptEtaNume1 =ptEtaNume1.size();
    for(int i=0;i<SizeOfptEtaNume1;i++){
      mapaire1.first=ptEtaNume1[i];
      mapaire1.second = ptEtaDeno1[i];
      VecPaire1.push_back(mapaire1);
    }
    pmuEffdata.SetPtEtaNumeDeno(VecPaire1);
    VectorOfmuEffdata.push_back(pmuEffdata);
    //
    // Re-Initialization
    int nBin2 =0;
    double etamin2=0;
    double etamax2=0;
    int ndeno2=0;
    int nnum2=0;
    //
    InputFile >>etamin2>>etamax2>>nBin2;
    if(InputFile.eof()) return 0; 
    //
    pmuEffdata.SetEtaMin(etamin2 );
    pmuEffdata.SetEtaMax(etamax2);
    //
    for (int Item=0; Item < nBin2 ; Item++){
      if(InputFile.eof()) return 0;
      InputFile >> nnum2 ;
      ptEtaNume2.push_back(nnum2); 
    }
    if(InputFile.eof()) return 0;
    for (int Item=0; Item < nBin2 ; Item++){
      if(InputFile.eof()) return 0;
      InputFile >> ndeno2 ;
      ptEtaDeno2.push_back(ndeno2);
    }
    pmuEffdata.SetPtEtaNum(ptEtaNume2);
    pmuEffdata.SetPtEtaNum(ptEtaDeno2);
    //  
    //End second eta part
    if(InputFile.eof()) return 0;
    //
    int SizeOfptEtaNume2 =ptEtaNume2.size();
    for(int i=0;i<SizeOfptEtaNume2;i++){
      mapaire2.first=ptEtaNume2[i];
      mapaire2.second = ptEtaDeno2[i];
      VecPaire2.push_back(mapaire2);
    }
    pmuEffdata.SetPtEtaNumeDeno(VecPaire2);
    VectorOfmuEffdata.push_back(pmuEffdata);
    //
    //
    // Re-Initialization
    int nBin3 =0;
    double etamin3=0;
    double etamax3=0;
    int ndeno3=0;
    int nnum3=0;
    //
    InputFile >>etamin3>>etamax3>>nBin3;
    if(InputFile.eof()) return 0; 
    //
    pmuEffdata.SetEtaMin(etamin3 );
    pmuEffdata.SetEtaMax(etamax3);
    //
    for (int Item=0; Item < nBin3 ; Item++){
      if(InputFile.eof()) return 0;
      InputFile >> nnum3 ;
      ptEtaNume3.push_back(nnum3); 
    }
    if(InputFile.eof()) return 0;
    for (int Item=0; Item < nBin3 ; Item++){
      if(InputFile.eof()) return 0;
      InputFile >> ndeno3 ;
      ptEtaDeno3.push_back(ndeno3);
    }
    pmuEffdata.SetPtEtaNum(ptEtaNume3);
    pmuEffdata.SetPtEtaNum(ptEtaDeno3);
    //  
    //End second eta part
    if(InputFile.eof()) return 0;
    //
    int SizeOfptEtaNume3 =ptEtaNume3.size();
    for(int i=0;i<SizeOfptEtaNume3;i++){
      mapaire3.first=ptEtaNume3[i];
      mapaire3.second = ptEtaDeno3[i];
      VecPaire3.push_back(mapaire3);
    }
    pmuEffdata.SetPtEtaNumeDeno(VecPaire3);
    VectorOfmuEffdata.push_back(pmuEffdata);
    //
    //
    //
    // Re-Initialization
    int nBin4 =0;
    double etamin4=0;
    double etamax4=0;
    int ndeno4=0;
    int nnum4=0;
    //
    InputFile >>etamin4>>etamax4>>nBin4;
    if(InputFile.eof()) return 0; 
    //
    pmuEffdata.SetEtaMin(etamin4 );
    pmuEffdata.SetEtaMax(etamax4);
    //
    for (int Item=0; Item < nBin4 ; Item++){
      if(InputFile.eof()) return 0;
      InputFile >> nnum4 ;
      ptEtaNume4.push_back(nnum4); 
    }
    if(InputFile.eof()) return 0;
    for (int Item=0; Item < nBin4 ; Item++){
      if(InputFile.eof()) return 0;
      InputFile >> ndeno4 ;
      ptEtaDeno4.push_back(ndeno4);
    }
    pmuEffdata.SetPtEtaNum(ptEtaNume4);
    pmuEffdata.SetPtEtaNum(ptEtaDeno4);
    //  
    //End second eta part
    if(InputFile.eof()) return 0;
    //
    int SizeOfptEtaNume4 =ptEtaNume4.size();
    for(int i=0;i<SizeOfptEtaNume4;i++){
      mapaire4.first=ptEtaNume4[i];
      mapaire4.second = ptEtaDeno4[i];
      VecPaire4.push_back(mapaire4);
    }
    pmuEffdata.SetPtEtaNumeDeno(VecPaire4);
    VectorOfmuEffdata.push_back(pmuEffdata);
    //
    //
    return 1;
  }
  
  
  //-------------------------------------
  std::pair<double,double>
  MuonAcceptor::muEfficiency(double ptIn,double etaIn, double phi)
  {

    double efficiency    =0;
    double errEfficiency =0; 
   
    double eta = fabs(etaIn);
    double pt  = fabs(ptIn);

    // Interpolate between pt points
    double pt1 = 0.;
    double pt2 = 0.;

    // Get the two pt points
    std::pair<double,double>ptLowHigh;
    ptLowHigh= muGetPtLowHigh(pt);
    pt1 = ptLowHigh.first;
    pt2 = ptLowHigh.second;

    if(pt >= pt1 && pt < pt2){    
      std::pair<double,double> eff1_errff1;
      std::pair<double,double> eff2_errff2;

      eff1_errff1 = muEff(pt1,eta,phi);
      eff2_errff2 = muEff(pt2,eta,phi); 

      double eff1  =  eff1_errff1.first ;
      double eeff1 =  eff1_errff1.second;

      double eff2  =  eff2_errff2.first;
      double eeff2 =  eff2_errff2.second;

      efficiency    =              eff2 *log(pt/pt1)/log(pt2/pt1);     
      efficiency    = efficiency + eff1 *log(pt2/pt)/log(pt2/pt1);     
      errEfficiency = pow((eeff2*log(pt/pt1)/log(pt2/pt1)),2); 
      errEfficiency = errEfficiency +pow((eeff1*log(pt2/pt)/log(pt2/pt1)),2); 
      errEfficiency = sqrt(errEfficiency); 
    }

    if(efficiency < 0){
      efficiency   =0.;
      errEfficiency=0.;
    }
    if(efficiency > 1){
      efficiency   =1.;
      errEfficiency=0.;
    }

    return make_pair(efficiency,errEfficiency);
  }
  //-------------------------------------
  std::pair<double,double>
  MuonAcceptor::muEff(double ptIn,double etaIn, double phi)
  {
    double eta = fabs(etaIn);
    double pt  = fabs(ptIn);
 
    double xEtaMinOut = 0.;
    double xEtaMaxOut = 2.5;

    int nNumOut  = 0;            
    int nDenoOut = 0;  

    int nNumTotOut = 0;
    int nDenoTotOut= 0;
    //
    int iEtaBin = muGetEtaBin(eta);
    //
    std::vector<std::pair<int, int> > ptEtaNumeDeno;
    ptEtaNumeDeno.clear();
    ptEtaNumeDeno= muEffGetPaire(pt,xEtaMinOut,xEtaMaxOut);

    nNumOut=ptEtaNumeDeno[iEtaBin].first;
    nDenoOut=ptEtaNumeDeno[iEtaBin].second;

    int nBinOut = ptEtaNumeDeno.size();

    for(unsigned int i=0 ; i<ptEtaNumeDeno.size() ; i++){
      nNumTotOut = nNumTotOut+ptEtaNumeDeno[i].first;
      nDenoTotOut = nDenoTotOut+ptEtaNumeDeno[i].second;
    }

    double eff1    = double(nNumOut)/double(nDenoOut);                     
    double eeff1   = sqrt(eff1*(1.-eff1)/double(nDenoOut));                
    double delta   = (xEtaMaxOut-xEtaMinOut)/double(nBinOut) ;                 
    double xcbin1  = xEtaMinOut+delta/2.+delta*double(iEtaBin) ; 
    //
    muPhiEfficiency(pt,xcbin1,phi,eff1,eeff1);
    //
    return make_pair(eff1,eeff1);
  }
  //-------------------------------------
  int MuonAcceptor::muGetEtaBin(double etaIn){

    int iEtaBin = -1;
    double eta = fabs(etaIn);

    double xMaxOut =2.5;
    double xMinOut = 0.;
    double pt=50.;
    //
    std::vector<std::pair<int, int> > ptEtaNumeDeno;
    ptEtaNumeDeno.clear();
    ptEtaNumeDeno= muEffGetPaire(pt,xMinOut,xMaxOut);
    //
    int  NbinOut   = ptEtaNumeDeno.size();

    double delta=double(xMaxOut-xMinOut)/double(NbinOut); 
    iEtaBin = int(floor(eta/delta));                         
    if( iEtaBin >= NbinOut) iEtaBin= NbinOut-1;  

    return iEtaBin;
  }
  //-------------------------------------
  void MuonAcceptor::muPhiEfficiency(double pt,double etaIn,double phi,double &eff,double &eeff){

    double pi = 3.1415927;
    double phig = phi;
    if (phi < 0.) {
      phig = pi*2 + phi;
    }

    double effin  = eff;
    //double eeffin = eeff;

    double eta = fabs(etaIn);

    //---------------------------
    //Eta region 0. <|eta|<0.15
    //---------------------------

    double etamin1 = 0.;
    double etamax1 = 0.15;
    if (eta >= etamin1 && eta < etamax1) {
      double pion4  = pi/4.;
      double pion8  = pi/8.;
      double phig2  = fmod(phig,pion4);
      if (phig2 >= pion8) phig2 = pion4-phig2;
      //
      std::pair<double,double> eff_erreff;
      eff_erreff = getPhiEff(pt,effin,phig2,etamin1,etamax1);
      eff  =  eff_erreff.first ;
      eeff =  eff_erreff.second;
    }

    //------------------------------------------------------
    //Feets: Eta region 0.3 <|eta|<0.4 && 0.8 <|eta|<0.9
    //------------------------------------------------------

    double etamin2 = 0.3;
    double etamax2 = 0.4;
    double etamin22 = 0.8;
    double etamax22 = 0.9;
    if ((eta >= etamin2 && eta < etamax2) ||
        (eta >= etamin22 && eta < etamax22) ) {
      std::pair<double,double> eff_erreff;
      eff_erreff = getPhiEff(pt,effin,phig,etamin2,etamax2);
      eff  =  eff_erreff.first ;
      eeff =  eff_erreff.second;
    }
    //---------------------------
    // Eta region 1.05 <|eta|<1.15
    //---------------------------
    double etamin3 =1.05;
    double etamax3 =1.15;
    if (eta >= etamin3 && eta < etamax3) {

      std::pair<double,double> eff_erreff;
      eff_erreff = getPhiEff(pt,effin,phig,etamin3,etamax3);
      eff  =  eff_erreff.first ;
      eeff =  eff_erreff.second;
    }
    //
    //---------------------------
    // Eta region 1.20 <|eta|<1.25
    //---------------------------
    double etamin4 =1.20;
    double etamax4 =1.25;

    if (eta >= etamin4 && eta < etamax4) {

      std::pair<double,double> eff_erreff;
      eff_erreff = getPhiEff(pt,effin,phig,etamin4,etamax4);
      eff  =  eff_erreff.first ;
      eeff =  eff_erreff.second;
    }

  }

  //------------------------------------- 
  int MuonAcceptor::muGetPhiBin(double phi){

    double pi = 3.1415927;
    int iPhiBin = -1;

    double xMaxOut =360;
    double xMinOut = 0.;
    int  NbinOut   = 20;

    double phig = phi;
    if (phi < 0.) {
      phig = pi*2 + phi;
    }
    double phideg = phig*180./ pi;

    double delta=double(xMaxOut-xMinOut)/double(NbinOut); 
    iPhiBin = int(floor(phideg/delta));                         
    if( iPhiBin >= NbinOut) iPhiBin= NbinOut-1;

    return iPhiBin;
  }
  //------------------------------------- 
  int MuonAcceptor::muGetPhiBin1(double phi){

    double pi = 3.1415927;
    int iPhiBin = -1;

    double xMaxOut = 22.5;
    double xMinOut = 0.;
    int  NbinOut   = 20;

    double phig = phi;
    if (phi < 0.) {
      phig = pi*2 + phi;
    }
    double phideg = phig*180./ pi;

    double delta=double(xMaxOut-xMinOut)/double(NbinOut); 
    iPhiBin = int(floor(phideg/delta)); 
    if( iPhiBin >= NbinOut) iPhiBin= NbinOut-1;               
 
    return iPhiBin;
  }

  //-------------------------------------
  std::pair<double,double>
  MuonAcceptor::getPhiEff(double pt,double effin,double phig,double etamin,double etamax)
  {
    double eff = 1.;
    double eeff =0.;

    double phivalue[21];
    double efficiency[21];
    double errefficiency[21];

    for (int j = 0; j < 21; j++) {
      efficiency[j]    = 0.;
      errefficiency[j] = 0.;
      phivalue[j]      = 0.;
    }
    //
    double pi = 3.1415927;
    double degrad  = pi/180.;
    //
    for (int j = 0; j < 21; j++) {
      phivalue[j] = j*18.*degrad;
      if(etamin ==0. && etamax ==0.15)phivalue[j] = j*1.125*degrad;
    }
    //
    // Mean Efficiency in phi

    std::pair<double,double> eff_errff;
    double phiLow  = phivalue[0];
    double phiHigh = phivalue[20];
    eff_errff         = muPhiEff(pt,etamin,etamax,phiLow,phiHigh);
    efficiency[20]    = eff_errff.first ;
    errefficiency[20] = eff_errff.second;

    // The other bins
    double delta2 = 8.9*degrad;
    if(etamin ==0. && etamax ==0.15)delta2 = 0.55*degrad;
    for (int j = 0; j < 20; j++) {
      double phicenter = (phivalue[j+1]+phivalue[j])/2.;
      double phimin    =  phicenter - delta2;
      double phimax    =  phicenter + delta2;
      eff_errff = muPhiEff(pt,etamin,etamax,phimin,phimax);
      efficiency[j] =  eff_errff.first ;
      errefficiency[j] = eff_errff.second;
    }
    //
    double averaefficiency = effin/efficiency[20];
    double factornorm = averaefficiency;
    if(factornorm > 1. || factornorm < 0.9)factornorm =1;
    //
    if (phig >= phivalue[0] && phig <= phivalue[1]) {
      eff  = efficiency[0] * factornorm;
      eeff = errefficiency[0];
    }
    for (int j = 1; j < 20; j++) {
      if (phig > phivalue[j] && phig <= phivalue[j+1]) {
        eff = efficiency[j] * factornorm;
        eeff= errefficiency[j];
      }
    }
    if (phig > phivalue[20] && phig <= phivalue[21]) {
      eff = efficiency[19] * factornorm;
      eeff= errefficiency[19];
    }
    //
    if (eff < 0.) {
      eff = 0.;
      eeff = 0.;
    }
    if (eff > 1.) {
      eff = 1.;
      eeff =0.;
    }
    //
    return make_pair(eff,eeff);
  }
                            
  //------------------------------------- 
  std::pair<double,double>
  MuonAcceptor::muPhiEff(double pt,double etamin,double etamax,double phimin,double phimax){

    int NnumTot=0 ;                                  
    int NdenoTot=0;  

    double eff =0;
    double eeff=0;  

    int ibinmin;
    int ibinmax;

    if(etamin==0. && etamax==0.15){
      ibinmin = muGetPhiBin1(phimin) ;
      ibinmax = muGetPhiBin1(phimax) ; 
    }
    else{
      ibinmin = muGetPhiBin(phimin) ;
      ibinmax = muGetPhiBin(phimax) ; 
    }
    std::vector<std::pair<int, int> > ptEtaNumeDeno;
    ptEtaNumeDeno.clear();
    ptEtaNumeDeno= muEffGetPaire(pt,etamin,etamax);

    for (int ibin=ibinmin;ibin<=ibinmax;ibin++){
      NnumTot =NnumTot+ptEtaNumeDeno[ibin].first;
      NdenoTot=NdenoTot+ptEtaNumeDeno[ibin].second;
    }
 
    eff =double(NnumTot)/double(NdenoTot); 
    eeff=sqrt(eff*(1.-eff)/double(NdenoTot)); 
    //
    return make_pair(eff,eeff);
  }

  //-------------------------------------
  std::vector<std::pair<int, int> >
  MuonAcceptor::muEffGetPaire(double pt,double etamin,double etamax){
    
    std::vector<std::pair<int, int> > vecPaireNumeDeno;
    vecPaireNumeDeno.clear();
    //
    std::vector<std::pair<int, int> > ptEtaNumeDeno;
    ptEtaNumeDeno.clear();
    // 
    int iposition = muGetPaireNumeDenoPosi(pt,etamin,etamax);
    
    if(iposition != -1){
      int i=iposition;
      m_muEffdataObject[i].GetPtEtaNumeDeno(pt,etamin,etamax,ptEtaNumeDeno);
      vecPaireNumeDeno=ptEtaNumeDeno;
    }
    //
    return vecPaireNumeDeno;
  }
  //-------------------------------------  //-------------------------------------
  int MuonAcceptor::muGetPaireNumeDenoPosi(double pt,double etamin,double etamax){
    
    int iposition =-1;
    const double epsilon=1.e-6;
    for(unsigned int i= 0;i<m_ptValues.size();i++){
      if( fabs(pt-m_ptValues[i]) < epsilon &&
          fabs(etamin - m_etaminValues[i]) < epsilon &&
          fabs( etamax - m_etamaxValues[i]) < epsilon 
          )     
        iposition=i;
    }
    return iposition;
  }
  //-------------------------------------
  std::pair<double,double> 
  MuonAcceptor::muGetPtLowHigh(double pt){

    std::pair<double,double> ptLowHigh;

    int sizeOfptvalues   = m_ptValues.size();
    double ptmin= m_ptValues.at(0);
    double ptmax= m_ptValues.at(sizeOfptvalues-1);

    for (int i=0; i <sizeOfptvalues-1 ; i++){     
      ptmin= m_ptValues[i];
      ptmax= m_ptValues[i+1];
      if(pt >=ptmin && ptmin<ptmax){
        double ptlow  =m_ptValues[i] ;
        double pthigh =m_ptValues[i+1] ;
        ptLowHigh=make_pair(ptlow,pthigh);     
      } 
    }  
    if(pt>=m_ptValues.at(sizeOfptvalues-1)){
      double ptlow  = m_ptValues.at(sizeOfptvalues-2);
      double pthigh = m_ptValues.at(sizeOfptvalues-1);
      ptLowHigh=make_pair(ptlow,pthigh);
    }
    if(pt <=m_ptValues.at(0)){
      double ptlow  = m_ptValues.at(0);
      double pthigh = m_ptValues.at(0) ;
      ptLowHigh=make_pair(ptlow,pthigh);
    }
    return ptLowHigh ;
  }
  
} // end of namespace

---