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

#ifndef muon_SF_R16_h
#define muon_SF_R16_h

/******************************************************************************
 * muon_SF_R16.h                                                              *
 *                                                                            *
 * Simple functions which return data/MC scale factors                        *
 * given a muon's (or two muons') eta and phi (defined b/w -pi and pi).       *
 * Also functions for corresponding uncertainties.                            *
 *                     .                                                      *
 * double mu_trigger_SF(double eta, double phi)                               * 
 * double mu_trigger_SF_err(double eta, double phi, bool isUp)                * 
 * double mu_trigger_SF_staterr(double eta, double phi)                       * 
 * double mu_trigger_SF_systerr(double eta, double phi, bool isUp)            * 
 *                                                                            * 
 * double mu_recoID_SF()                                                      *
 * double mu_recoID_SF_err()                                                  *
 * double mu_recoID_SF_staterr()                                              *
 * double mu_recoID_SF_systerr()                                              *
 *                                                                            *
 * History                                                                    *
 *         08 Dec 2010 -- created by S. Caughron                              *   
 *****************************************************************************/

#include <iostream>

// forward declaration of functions
double mu_trigger_SF(double eta, double phi);
double mu_trigger_SF_err(double eta, double phi, bool isUp);
double mu_trigger_SF_staterr(double eta, double phi);
double mu_trigger_SF_systerr(double eta, double phi, bool isUp);
double mu_recoID_SF();
double mu_recoID_SF_err();
double mu_recoID_SF_staterr();
double mu_recoID_SF_systerr();
double check_phi_range(double phi);


// For trigger scale factors
double mu_trigger_SF(double eta, double phi)
{

    double mu_eta = eta;
    double mu_phi = check_phi_range(phi);
    int etaI=-1, phiI=-1;   

    const double pi = acos(-1.);

    const double etabins[3] = {-1.05,-0.6,0.6};   // lower edges of eta bins
    const double etabins_P2[4] = {-0.5,-0.4,0.2,0.3};   // defines special B2P2 binning

    const double phibins[4] = {-pi,5*pi/16,pi/2,11*pi/16};   // lower edges of phi bins, first and last bin have same SF

    const double trig_SFmatrix[3][3] = {{1.0264,0.9195,0.9522},   // rows are eta, columns are phi
                                        {1.0089,0.6571,0.9065},
                                        {1.0052,0.8430,1.0458}};

    const double EC_eta = 1.05;
    const double EC_SF = 0.9875;

    if ( fabs(mu_eta) > EC_eta )    // only one SF in endcaps

        return EC_SF;       

    else {                          // if in barrel, look up appropriate SF 
        
        for (int i=3; i>=0; i--){    // find phi index
            if ( mu_phi > phibins[i] ) {
                phiI = i;
                break;
            }
        }
        for (int i=2; i>=0; i--){    // find eta index
            if ( mu_eta > etabins[i] ) {
                etaI = i;
                break;
            }
        }

                if (phiI==1 && etaI == 1) {  // "B1P2" and "B2P2" regions have different binning in rel16
                        if ( !((mu_eta > etabins_P2[0] && mu_eta < etabins_P2[1]) || (mu_eta > etabins_P2[2] && mu_eta < etabins_P2[3])) )
                                etaI = 0;
                }

        return trig_SFmatrix[etaI][phiI % 3];   // if phiI==3, use 0 (same SF)

    } //else

    return 0;
}

// For trigger scale factor uncertainties. 
// These are asymmetric. Passing isUp = kTRUE returns uncertainty in positive direction, kFALSE returns negative direction 
// NB: For both cases, returns *magnitude* of uncertainty, i.e. always a positive number
double mu_trigger_SF_err(double eta, double phi, bool isUp)
{

    double mu_eta = eta;
    double mu_phi = check_phi_range(phi);
    double staterr, systerr, toterr;

    staterr = mu_trigger_SF_staterr(mu_eta,mu_phi);
    systerr = mu_trigger_SF_systerr(mu_eta,mu_phi,isUp);

    toterr = sqrt(pow(staterr,2) + pow(systerr,2));

    return toterr;

}

// For trigger scale factor statistical uncertainties. (symmetric)
double mu_trigger_SF_staterr(double eta, double phi)
{

    double mu_eta = eta;
    double mu_phi = check_phi_range(phi);
    int etaI=-1, phiI=-1;   
    //bool isEC = kFALSE;

    const double pi = acos(-1.);

    const double etabins[3] = {-1.05,-0.6,0.6};   // lower edges of eta bins
    const double etabins_P2[4] = {-0.5,-0.4,0.2,0.3};   // defines special B2P2 binning

    const double phibins[4] = {-pi,5*pi/16,pi/2,11*pi/16};   // lower edges of phi bins, first and last bin have same SF

    const double errmatrix[3][3] = {{0.0104,0.0167,0.0300},   // rows are eta, columns are phi
                                        {0.0056,0.0501,0.0188},
                                        {0.0101,0.0530,0.0288}};

    const double EC_eta = 1.05;
    const double EC_err= 0.004;

    if ( fabs(mu_eta) > EC_eta ) {   // only one SF in endcaps

                return EC_err;

    } 
    else {                          // if in barrel, look up appropriate SF 
        
        for (int i=2; i>=0; i--){    // find eta index
            if ( mu_eta > etabins[i] ) {
                etaI = i;
                break;
            }
        }

        for (int i=3; i>=0; i--){    // find phi index
            if ( mu_phi > phibins[i] ) {
                phiI = i;
                break;
            }
        }
                
                if (phiI==1 && etaI == 1) {  // "B1P2" and "B2P2" regions have different binning in rel16
                        if ( !((mu_eta > etabins_P2[0] && mu_eta < etabins_P2[1]) || (mu_eta > etabins_P2[2] && mu_eta < etabins_P2[3])) )
                                etaI = 0;
                }

        return errmatrix[etaI][phiI % 3];   // if phiI==3, use 0 (same SF)

    } //else

    return 0;
}

// For trigger scale factor systematic uncertainties. 
// These are asymmetric. Passing isUp = kTRUE returns uncertainty in positive direction, kFALSE returns negative direction 
// NB: For both cases, returns *magnitude* of uncertainty, i.e. always a positive number
double mu_trigger_SF_systerr(double eta, double phi, bool isUp)
{

    double mu_eta = eta;
    double mu_phi = check_phi_range(phi);
    int etaI=-1, phiI=-1;   
    //bool isEC = kFALSE;

    const double pi = acos(-1.);

    const double etabins[3] = {-1.05,-0.6,0.6};   // lower edges of eta bins
    const double etabins_P2[4] = {-0.5,-0.4,0.2,0.3};   // defines special B2P2 binning

    const double phibins[4] = {-pi,5*pi/16,pi/2,11*pi/16};   // lower edges of phi bins, first and last bin have same SF

    const double pos_errmatrix[3][3] = {{0.003,0.007,0.002},   // rows are eta, columns are phi
                                        {0.001,0.010,0.000},
                                        {0.002,0.000,0.011}};

    const double neg_errmatrix[3][3] = {{0.002,0.000,0.003},   // rows are eta, columns are phi
                                        {0.002,0.000,0.004},
                                        {0.003,0.013,0.009}};

    const double EC_eta = 1.05;
    const double EC_err_up = 0.001;
    const double EC_err_down = 0.001;

    if ( fabs(mu_eta) > EC_eta ) {   // only one SF in endcaps

        if (isUp)
                return EC_err_up;
        else
                return EC_err_down;

    } 
    else {                          // if in barrel, look up appropriate SF 
        
        for (int i=2; i>=0; i--){    // find eta index
            if ( mu_eta > etabins[i] ) {
                etaI = i;
                break;
            }
        }

        for (int i=3; i>=0; i--){    // find phi index
            if ( mu_phi > phibins[i] ) {
                phiI = i;
                break;
            }
        }
                
                if (phiI==1 && etaI == 1) {  // "B1P2" and "B2P2" regions have different binning in rel16
                        if ( !((mu_eta > etabins_P2[0] && mu_eta < etabins_P2[1]) || (mu_eta > etabins_P2[2] && mu_eta < etabins_P2[3])) )
                                etaI = 0;
                }

        if (isUp)
                return pos_errmatrix[etaI][phiI % 3];   // if phiI==3, use 0 (same SF)
        else
                return neg_errmatrix[etaI][phiI % 3];   

    } //else

    return 0;
}


// For reco+ID scale factors
double mu_recoID_SF()
{

        const double recoID_SF = 0.999;  // in rel16, just a number
    return recoID_SF;

}

// For reco+ID scale factor total uncertainties. 
double mu_recoID_SF_err()
{

    double staterr, systerr, toterr;

    staterr = mu_recoID_SF_staterr();
    systerr = mu_recoID_SF_systerr();

    toterr = sqrt(pow(staterr,2) + pow(systerr,2));

    return toterr;

}

// For reco+ID scale factor statistical uncertainties (symmetric)
double mu_recoID_SF_staterr()
{

        const double SFerr = 0.002;  //in rel16, just a number
        return SFerr;    

}

// For reco+ID scale factor systematic uncertainties 
double mu_recoID_SF_systerr()
{

        const double SFerr = 0.003;
        return SFerr;

}

// 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

---