/afs/hep.man.ac.uk/u/markowen/ATLAS/SFrameProof/2011EPS/topUtils/src/muon_Staco_SF_EPS.h

#ifndef muon_Staco_SF_EPS_h
#define muon_Staco_SF_EPS_h

#include <iostream>
#include <stdio.h>
#include <math.h>


double mu_staco_trigger_SF(double eta, double phi, double pt);
double mu_barrel_SF(double eta, double phi);
double mu_endcap_SF(double pt);
double check_Phi_range(double phi);

double mu_staco_trigger_SF(double eta, double phi, double pt){

  const double eta_EC = 1.05;
  double mu_SF = 0.;
  if(fabs(eta) < eta_EC ) mu_SF = mu_barrel_SF(eta,phi);
  else mu_SF = mu_endcap_SF(pt);
  
  return mu_SF;
}


// For trigger scale factors
double mu_endcap_SF(double pt){
  double mu_pt=pt;
  int mu_pt_index=-1;
  const double ptbins[8] = {15.,18.,20.,30.,40.,50.,70.,100.};   // lower edges of pT bins
  const double pt_SFarray[8] = {1.0058817,1.0276103,1.0207309,1.0291394,1.0371242,1.0411885,1.0198005,1.0501196};
  
  for (int i=7; i>=0; i--){    // find phi index
    if ( mu_pt > ptbins[i] ) {
      mu_pt_index = i;
      break;
    }
  }
  
  if(fabs(pt_SFarray[mu_pt_index]) >2.) std::cout << "EC ,pt= " << mu_pt << " sf = " << pt_SFarray[mu_pt_index] << std::endl;
  return  pt_SFarray[mu_pt_index];
  
}
// For trigger scale factors
double mu_barrel_SF(double eta, double phi){
  
  double mu_eta=eta;
  double mu_phi=check_Phi_range(phi);
  if(mu_phi != phi ) std::cout << "Original Phi is " << phi << " changed to " << mu_phi << std::endl;  

  int etaI=-1, phiI=-1;  

  const double etabins[6] = {-1.05,-0.908,-0.476,0.,0.476,0.908}; // lower edges of eta bins
  const double phibins[8] = {-2.945,-2.160,-1.374,-0.589,0.196,0.982,1.767,2.553};// lower edges of phi bins
  const double trig_SFmatrix[6][8] = {{1.1502,1.0660,1.4482,1.2095,1.2142,1.2623,1.2654,1.2445},  // rows are eta, columns are phi
                                      {1.0048,1.0102,1.1890,1.0415,1.0575,1.0710,1.0401,1.0279},
                                      {0.9173,1.0329,1.0243,1.0411,1.0814,1.0160,1.0291,0.9801},
                                      {0.9886,1.1381,1.0670,1.0174,1.0696,0.9710,1.0045,0.9845},
                                      {0.9829,0.9998,0.9519,0.9569,0.9735,1.0454,1.0030,0.9903},
                                      {0.9944,0.6583,1.0112,0.7800,0.8263,1.0898,1.0356,0.8622}};
  
  for (int i=7; i>=0; i--){    // find phi index
    if ( mu_phi > phibins[i] ) {
      phiI = i;
      break;
    }
  }
  for (int i=5; i>=0; i--){    // find eta index
    if ( mu_eta > etabins[i] ) {
      etaI = i;
      break;
    }
  }
  if(phiI==-1){
    phiI=0;
    std::cout << "PHI= " << mu_phi << " eta= " << mu_eta << std::endl;
  }
  

  for(int i=0; i<6 ; ++i){
    for(int j=0; j<8 ; ++j){
      //      std::cout << "trig_SFmatrix[" << i << "][" << j  << "] =" << trig_SFmatrix[i][j] << std::endl;
    }
  }



  if(fabs(trig_SFmatrix[etaI][phiI])>2.)std::cout << "Barrel, phiI=" << phiI << " etatI= " << etaI << "  ,eta= " << mu_eta << " phi= " << mu_phi << "  sf = " << trig_SFmatrix[etaI][phiI] << std::endl;
  
  return trig_SFmatrix[etaI][phiI]; 
}



// Helper function to deal with possible phi ambiguity
double check_Phi_range(double phi)
{

  double newphi = phi;
  const double pi = acos(-1.);
  
  if (newphi > pi) {
    printf("<muon_SF>: WARNING: Muon phi %4.2f > pi! Using (phi-2*pi) \n", phi);
    newphi -= 2*pi;
        }
  if (newphi < -pi) {
    printf("<muon_SF>: WARNING: Muon phi %4.2f < -pi! Using (phi+2*pi) \n", phi);
    newphi += 2*pi;
  }
  
  return newphi;
}


#endif

---