//////////////////////////////////////////////////////////////////////////
// LhoodCalc.cpp
//
// Created: Oct 6, 2003  Joe Kozminski
//
//////////////////////////////////////////////////////////////////////////

#include "emreco/LhoodCalc.hpp"
using namespace std;

namespace emid {

LhoodCalc::LhoodCalc(){

  emf_min=0.9;
  emf_max=1.0;
  emf_max_not42=1.1;
  hmx_min=0.0;
  hmx_max=50.0;
  dca_min=-0.1;
  dca_max=0.1;
  spchi_min=-25.0;
  spchi_max=0.0;
  etpt_min=0.0;
  etpt_max=4.5;
  ndr05_min=1;
  ndr05_max=9;
  tpt_min=0;
  tpt_max=5;
  cps_min=0;
  cps_max=60.0;
  
  eta_type=0;

}
LhoodCalc::~LhoodCalc(){}

double LhoodCalc::calc_lh(double deta, double emf, double hmx, double etpt, double dca, double spchi,double ndr05, double tpt, double cps, int cps_flag, int t42_flag, int hmx_flag){
  LHOOD = -1;
  if(fabs(deta)<1.1)eta_type=1; //CC
  else {
    eta_type=2;
    // JS: per request by NP Physics group, also calculate likelihood for
    //     1.1<|deta|<1.5 and |deta|>2.5
  }


  if(eta_type==1 || eta_type==2){
    if(t42_flag==1){
      emfprobS=emfSigLookup(eta_type,emf);
      emfprobB=emfBkgLookup(eta_type,emf);
      if(hmx_flag==1){   //hmx_flag == 1 use hmx!!!
	hmxprobS=hmxSigLookup(eta_type,hmx);
	hmxprobB=hmxBkgLookup(eta_type,hmx);
      }
      else{
	hmxprobS=1.0;
	hmxprobB=1.0;
      }	
    }
    if(t42_flag==0){
      emfprobS=emfSigLookup_not42(eta_type,emf);
      emfprobB=emfBkgLookup_not42(eta_type,emf);
      if(hmx_flag==1){
	hmxprobS=hmxSigLookup_not42(eta_type,hmx);
	hmxprobB=hmxBkgLookup_not42(eta_type,hmx);
      }
      else{
	hmxprobS=1.0;
	hmxprobB=1.0;
      }
    }
    dcaprobS=dcaSigLookup(eta_type,dca);
    dcaprobB=dcaBkgLookup(eta_type,dca);
    etptprobS=etptSigLookup(eta_type,etpt);
    etptprobB=etptBkgLookup(eta_type,etpt);
    spchiprobS=spchiSigLookup(eta_type,spchi);
    spchiprobB=spchiBkgLookup(eta_type,spchi);
    ndr05probS=ndr05SigLookup(eta_type,ndr05);
    ndr05probB=ndr05BkgLookup(eta_type,ndr05);
    tptprobS=tptSigLookup(eta_type,tpt);
    tptprobB=tptBkgLookup(eta_type,tpt);
    if(eta_type==1 && cps_flag ==1){
      cpsprobS=tptSigLookup(eta_type,cps); // should be cpsSigLookup
      cpsprobB=tptBkgLookup(eta_type,cps); // should be cpsBkgLookup
    }
    else{
      cpsprobS=1.0;
      cpsprobB=1.0;
    }        

    probS = emfprobS*hmxprobS*dcaprobS*etptprobS*spchiprobS*ndr05probS*tptprobS*cpsprobS;
    probB = emfprobB*hmxprobB*dcaprobB*etptprobB*spchiprobB*ndr05probB*tptprobB*cpsprobB;

    if((probS+probB)>0)LHOOD=probS/(probS+probB);
    return LHOOD;
  }
  
  return LHOOD;

  
}

double LhoodCalc::emfSigLookup(int region, double emf){
  
  //const int nemf=22;
  const int nemf=12;
  
  //EMF Vector CC
  float emf_cc[nemf] = {0, 0.00027467, 0.000683096, 0.00135263, 0.00249991, 0.00424062, 0.00653142, 0.00915567, 0.0218821, 0.0727876, 0.880592, 0};
  
  //EMF Vector EC
  float emf_ec[nemf] = {0, 0, 7.83786e-05, 0.000293326, 0.000852665, 0.00199034, 0.00379543, 0.00606009, 0.0147174, 0.0582923, 0.91392, 0};
  
  int bin = 0;
  if( emf == emf_max) {
    bin = nemf-2;
  } else if( emf > emf_max) {
    bin = nemf-1;
  } else if (emf < emf_min) {
    bin = 0;
  } else {
    bin=(int)(1+((emf-emf_min)/(emf_max-emf_min))*(double)(nemf-2));
  }
  if(region==1) return emf_cc[bin];
  if(region==2) return emf_ec[bin];
  return 0;
}

double LhoodCalc::emfBkgLookup(int region, double emf){
  
  //const int nemf=22;
  const int nemf=12;
  
  //EMF Vector CC
  float emf_cc[nemf] = {0, 0.00298351, 0.00988157, 0.0183231, 0.029408, 0.0434928, 0.0601263, 0.0782082, 0.0981876, 0.135096, 0.524293, 0};
  
  //EMF Vector EC
  float emf_ec[nemf] = {0, 0.000720113, 0.00421442, 0.00876603, 0.015445, 0.0250651, 0.0376762, 0.0522082, 0.0855391, 0.144588, 0.625778, 0};
  
  int bin = 0;
  if (emf == emf_max) {
    bin = nemf-2;
  } else if( emf > emf_max) {
    bin = nemf-1;
  } else if( emf < emf_min) {
    bin = 0;
  } else {
    bin=(int)(1+((emf-emf_min)/(emf_max-emf_min))*(double)(nemf-2));
  }
  
  if(region==1) return emf_cc[bin];
  if(region==2) return emf_ec[bin];
  return 0;
}

double LhoodCalc::hmxSigLookup(int region, double hmx){

  const int nhmx=52;
  
//HMX7 Vector CC
  float hmx_cc[nhmx] = {0, 0.0657824, 0.186946, 0.181879, 0.139654, 0.102585, 0.0696938, 0.053738, 0.039866, 0.0292044, 0.0217007, 0.0166033, 0.0130745, 0.0104893, 0.0084893, 0.00690675, 0.00569587, 0.00480619, 0.00415649, 0.0036286, 0.00314535, 0.00266978, 0.00222348, 0.00185253, 0.00161328, 0.00152073, 0.00152878, 0.00154926, 0.00151451, 0.00141866, 0.00130453, 0.00121673, 0.00117283, 0.00116588, 0.00117612, 0.00117283, 0.0011293, 0.00103345, 0.000909072, 0.00078396, 0.00068848, 0.000624827, 0.000593366, 0.000573611, 0.00054837, 0.000498983, 0.000429111, 0.000352288, 0.000283879, 0.000227176, 0.00017779, 0};

//HMX7 Vector EC
  float hmx_ec[nhmx] = {0, 0.00339924, 0.056654, 0.157498, 0.18016, 0.151833, 0.115574, 0.0875293, 0.0626902, 0.0430268, 0.0289519, 0.0196238, 0.0137835, 0.0101534, 0.00785929, 0.00636717, 0.00539496, 0.00470019, 0.0040637, 0.00340624, 0.00277908, 0.00233844, 0.00214726, 0.00219622, 0.00233611, 0.00243869, 0.00239672, 0.00222886, 0.00197706, 0.00169729, 0.00141052, 0.00116106, 0.000995526, 0.000939571, 0.00095356, 0.000965217, 0.000913925, 0.000785696, 0.000624827, 0.00050126, 0.000482609, 0.00056654, 0.000687775, 0.000757719, 0.000739067, 0.000664461, 0.000575866, 0.000489603, 0.000373031, 0.000207498, 0, 0};

  int bin = 0;  
  if( hmx > hmx_max) {
    bin = nhmx -1;
  } else if( hmx < hmx_min) {
    bin = 0;
  } else {
    bin=(int)(1+((hmx-hmx_min)/(hmx_max-hmx_min))*(double)(nhmx-2));
  } 
  
  if(region==1) return hmx_cc[bin];
  if(region==2) return hmx_ec[bin];
  return 0;
}
double LhoodCalc::hmxBkgLookup(int region, double hmx){

  const int nhmx=52;

//HMX7 Vector CC
  float hmx_cc[nhmx] = {0, 0.0124809, 0.0404143, 0.0454361, 0.0487669, 0.0495938, 0.0482058, 0.0455592, 0.0426122, 0.0399399, 0.0376727, 0.0356792, 0.0337748, 0.0318606, 0.0299519, 0.0281175, 0.0264212, 0.024892, 0.0235331, 0.0223423, 0.0212807, 0.0202936, 0.0193183, 0.0183604, 0.0174602, 0.0166697, 0.0159904, 0.0153866, 0.014797, 0.0141678, 0.0134851, 0.0127748, 0.0120793, 0.0114167, 0.0107894, 0.010202, 0.00966217, 0.00916165, 0.0086676, 0.00815183, 0.00760555, 0.0070482, 0.00649833, 0.00597133, 0.00546635, 0.00497834, 0.00447825, 0.00390607, 0.00318327, 0.00227099, 0.00122362, 0};

//HMX7 Vector EC
  float hmx_ec[nhmx] = {0, 0.00666541, 0.0479597, 0.0725384, 0.0799328, 0.0733715, 0.0683725, 0.0605913, 0.0533269, 0.0468715, 0.041287, 0.0364889, 0.0324021, 0.0289263, 0.0259744, 0.0234588, 0.0213441, 0.0195835, 0.0180761, 0.0166858, 0.0153381, 0.014089, 0.0130464, 0.0122653, 0.011669, 0.0111279, 0.010539, 0.0099006, 0.00928922, 0.00879806, 0.00847019, 0.00825268, 0.00805082, 0.00777609, 0.00741608, 0.00700078, 0.00657969, 0.00618196, 0.00582152, 0.00550822, 0.00523864, 0.00499841, 0.00474404, 0.00444296, 0.0040778, 0.003673, 0.00324956, 0.00282418, 0.00238574, 0.00192993, 0.00145806, 0};

  int bin = 0;
  if( hmx > hmx_max) {
    bin = nhmx -1;
  } else if( hmx < hmx_min) {
    bin = 0;
  } else {
    bin=(int)(1+((hmx-hmx_min)/(hmx_max-hmx_min))*(double)(nhmx-2));
  }   
  if(region==1) return hmx_cc[bin];
  if(region==2) return hmx_ec[bin];
  return 0;
}

    
double LhoodCalc::dcaSigLookup(int region, double dca){

  const int ndca=42;

//DCA Vector CC
  float dca_cc[ndca] = {0.00559167, 0.000916667, 0.00108416, 0.00124033, 0.00138405, 0.00153004, 0.00170168, 0.00192538, 0.00221396, 0.00255799, 0.00292013, 0.00326341, 0.0035784, 0.00391036, 0.00433022, 0.00491568, 0.00570446, 0.00668865, 0.00779167, 0.0207167, 0.4037, 0.426434, 0.0224583, 0.00825001, 0.00703985, 0.00593684, 0.00501602, 0.00430834, 0.00377757, 0.00337092, 0.00302651, 0.00272096, 0.00242898, 0.0021453, 0.0018605, 0.00159681, 0.00137198, 0.00119393, 0.00103097, 0.000849898, 0.000641667, 0.006875};

//DCA Vector EC
  float dca_ec[ndca] = {0.00370929, 0.00028533, 0.000274762, 0.000299421, 0.000375743, 0.000499034, 0.000649332, 0.000814893, 0.000987501, 0.0011742, 0.00138673, 0.00166854, 0.00207598, 0.00268187, 0.00359305, 0.00500326, 0.00715439, 0.010138, 0.0136958, 0.0499328, 0.374638, 0.40631, 0.0459382, 0.0174051, 0.012834, 0.00886989, 0.00587921, 0.00390069, 0.00270887, 0.00201375, 0.00157695, 0.00127283, 0.00103682, 0.000847771, 0.000697474, 0.000587099, 0.00051195, 0.000469679, 0.000463808, 0.000501383, 0.00057066, 0.00456528};
  
  int bin = 0;
  if( dca > dca_max) {
    bin = ndca-1;
  } else if( dca < dca_min) {
    bin = 0;
  }else {
    bin=(int)(1+((dca-dca_min)/(dca_max-dca_min))*(double)(ndca-2));
  }
  if(region==1) return dca_cc[bin];
  if(region==2) return dca_ec[bin];
  return 0;
}
double LhoodCalc::dcaBkgLookup(int region, double dca){

  const int ndca=42;

//DCA Vector CC
  float dca_cc[ndca] = {0.0739179, 0.00413357, 0.00457084, 0.00497209, 0.00532173, 0.00564679, 0.00599801, 0.00643657, 0.00698105, 0.00762416, 0.00832359, 0.00904919, 0.00978064, 0.0105371, 0.0113806, 0.012409, 0.0137175, 0.0153228, 0.0171248, 0.0315401, 0.231341, 0.242213, 0.0340759, 0.0176111, 0.0155365, 0.0136739, 0.0121671, 0.0110271, 0.010171, 0.00945616, 0.00880032, 0.00817836, 0.00762631, 0.00715716, 0.00675291, 0.00636395, 0.0059424, 0.00545739, 0.00488961, 0.00424249, 0.00354306, 0.0689854};
  //DCA Vector EC
  float dca_ec[ndca] = {0.0725028, 0.00265254, 0.00317629, 0.0036386, 0.00401295, 0.00431282, 0.00459106, 0.00490848, 0.00533805, 0.005929, 0.00671472, 0.00768687, 0.0088677, 0.0103456, 0.0123473, 0.0151535, 0.0189854, 0.0238003, 0.02923, 0.0603323, 0.192127, 0.196236, 0.0590841, 0.0324026, 0.0256239, 0.0196951, 0.0151101, 0.0118869, 0.00975766, 0.00835252, 0.00737609, 0.0066336, 0.00603751, 0.00554395, 0.00511095, 0.00470065, 0.00430126, 0.00394553, 0.00368355, 0.00354143, 0.00348471, 0.0708385};

  int bin = 0;
  if( dca > dca_max) {
    bin = ndca-1;
  } else if( dca < dca_min) {
    bin = 0;
  }else {
    bin=(int)(1+((dca-dca_min)/(dca_max-dca_min))*(double)(ndca-2));
  }
  
  if(region==1) return dca_cc[bin];
  if(region==2) return dca_ec[bin];
  return 0;
}

double LhoodCalc::etptSigLookup(int region, double etpt){

  const int netpt=47;

  //Et/Pt Vector CC
  float etpt_cc[netpt] = {0, 0.00150322, 0.00283578, 0.00420436, 0.00589244, 0.00853391, 0.0128449, 0.0191056, 0.0267043, 0.068883, 0.160756, 0.204704, 0.146962, 0.0775486, 0.0532318, 0.0415939, 0.0314727, 0.0237051, 0.0182614, 0.0145664, 0.0119774, 0.0100299, 0.008412, 0.00692115, 0.00545214, 0.00407555, 0.00293258, 0.00214439, 0.00168407, 0.00144864, 0.00131546, 0.0012223, 0.00113024, 0.00102216, 0.000881337, 0.000726685, 0.000587679, 0.00049598, 0.000447583, 0.000424294, 0.000402824, 0.000383538, 0.000373713, 0.000383538, 0.000408647, 0.000442124, 0.0109647};

  //Et/Pt Vector EC
  float etpt_ec[netpt] = {0, 0.00337928, 0.00572715, 0.00817933, 0.010927, 0.0143851, 0.0192014, 0.0260816, 0.0352251, 0.0459018, 0.0673039, 0.0946197, 0.0991254, 0.0782865, 0.0580109, 0.0447754, 0.037938, 0.0315028, 0.0258834, 0.021358, 0.017966, 0.0154895, 0.0136295, 0.0121357, 0.0108563, 0.00979132, 0.00897895, 0.00846905, 0.00813877, 0.00784673, 0.00742606, 0.00689182, 0.00627646, 0.00569007, 0.00513381, 0.00462738, 0.00413834, 0.00371071, 0.00338043, 0.00318227, 0.00308376, 0.00302234, 0.00293659, 0.00278246, 0.0025472, 0.00225285, 0.0918036};

  int bin = 0;
  if( etpt > etpt_max) {
    bin = netpt-1;
  } else if( etpt < etpt_min) {
    bin = 0;
  } else {
    bin=(int)(1+((etpt-etpt_min)/(etpt_max-etpt_min))*(double)(netpt-2));
  }
  
  if(region==1) return etpt_cc[bin];
  if(region==2) return etpt_ec[bin];
  return 0;
}
double LhoodCalc::etptBkgLookup(int region, double etpt){

  const int netpt=47;

//Et/Pt Vector CC
  float etpt_cc[netpt] = {0, 0.00429486, 0.00484573, 0.00558155, 0.00667961, 0.00843588, 0.0114313, 0.016493, 0.0242008, 0.034188, 0.0448663, 0.054039, 0.0601082, 0.0628507, 0.0632102, 0.0623801, 0.0610742, 0.0594349, 0.0572936, 0.0544858, 0.0508678, 0.0463692, 0.0409208, 0.0346559, 0.0279581, 0.0215375, 0.016087, 0.0120109, 0.00922556, 0.00737414, 0.00606568, 0.00505895, 0.00423406, 0.0035384, 0.0029345, 0.00240371, 0.00194248, 0.00156071, 0.00126166, 0.00102907, 0.000837479, 0.000674027, 0.000559215, 0.00052627, 0.000590746, 0.000721536, 0.00316101};

//Et/Pt Vector EC
  float etpt_ec[netpt] = {0, 0.00639452, 0.0068423, 0.00731249, 0.00786532, 0.00862497, 0.00977348, 0.0115026, 0.013891, 0.0167939, 0.0198173, 0.0225027, 0.0245294, 0.0258743, 0.0266985, 0.0272407, 0.0276417, 0.0279842, 0.028258, 0.0284577, 0.0285525, 0.0285401, 0.0283861, 0.0280293, 0.0273823, 0.0263916, 0.0250909, 0.0235762, 0.0219781, 0.0203825, 0.0188443, 0.0173541, 0.0159038, 0.0144686, 0.0130789, 0.0117656, 0.0105792, 0.00952773, 0.00861954, 0.00784988, 0.00722115, 0.0067129, 0.00630317, 0.00597783, 0.00573512, 0.00554896, 0.228194};

  int bin = 0;
  if( etpt > etpt_max) {
    bin = netpt-1;
  } else if( etpt < etpt_min) {
    bin = 0;
  } else {
    bin=(int)(1+((etpt-etpt_min)/(etpt_max-etpt_min))*(double)(netpt-2));
  }

  if(region==1) return etpt_cc[bin];
  if(region==2) return etpt_ec[bin];
  return 0;
}


double LhoodCalc::spchiSigLookup(int region, double spchi){

  //const int nspchi=52;
  const int nspchi=32;

  //Spacial Chi^2 Vector CC
  float spchi_cc[nspchi] = {0.00879553, 0.000362702, 0.000420168, 0.000487708, 0.000568308, 0.000663089, 0.000764959, 0.000873172, 0.000975043, 0.00107169, 0.00116498, 0.00129185, 0.00148402, 0.00176836, 0.00212957, 0.00255422, 0.00305946, 0.00371845, 0.00464647, 0.00595735, 0.00780445, 0.0103915, 0.014032, 0.0189755, 0.0252168, 0.0322805, 0.044431, 0.0704549, 0.12613, 0.214538, 0.392988, 0};
  
  //Spacial Chi^2 Vector EC
  float spchi_ec[nspchi] = {0.0165967, 0.0005723, 0.00053344, 0.000546393, 0.000613515, 0.0006983, 0.000759534, 0.000801927, 0.000852562, 0.00094559, 0.00107041, 0.0012023, 0.00130828, 0.00141309, 0.00155675, 0.00178873, 0.00211963, 0.0025624, 0.00313823, 0.00385184, 0.00473619, 0.0059385, 0.00788267, 0.0110786, 0.0157583, 0.0214612, 0.0300457, 0.0549408, 0.103014, 0.218046, 0.484166, 0};


  if(spchi>0.0){
    double log_spchi=log(spchi);
    int bin = 0;
    if( log_spchi > spchi_max) {
      bin = nspchi-1;
    } else if( log_spchi < spchi_min) {
      bin = 0;
    } else {
      bin=(int)(1+((log_spchi-spchi_min)/(spchi_max-spchi_min))*(double)(nspchi-2));
    }
    if(region==1) return spchi_cc[bin];
    if(region==2) return spchi_ec[bin];
  }
  return 0;
}
double LhoodCalc::spchiBkgLookup(int region, double spchi){

  //const int nspchi=52;
  const int nspchi=32;

//Spacial Chi^2 Vector CC
  float spchi_cc[nspchi] = {0.151561, 0.00789416, 0.00861408, 0.00926904, 0.00984429, 0.0103605, 0.0108716, 0.0113945, 0.0119452, 0.0125061, 0.0130809, 0.0136679, 0.0143076, 0.0150593, 0.0159984, 0.0171684, 0.0185732, 0.0201901, 0.0219822, 0.023941, 0.0260946, 0.0285482, 0.0314485, 0.0349765, 0.0393637, 0.0449415, 0.0521794, 0.0615427, 0.0732829, 0.0871365, 0.102256, 0};

//Spacial Chi^2 Vector EC
  float spchi_ec[nspchi] = {0.301801, 0.00719056, 0.00780301, 0.00829939, 0.00863843, 0.00887391, 0.00909271, 0.00937403, 0.00975767, 0.0102549, 0.0108196, 0.0113602, 0.0117828, 0.0120978, 0.0124108, 0.0128701, 0.0135347, 0.0144105, 0.0154674, 0.016723, 0.0182133, 0.0200261, 0.0222865, 0.0251591, 0.0288458, 0.0337235, 0.040515, 0.050315, 0.064096, 0.0819182, 0.102339, 0};


  if(spchi>0.0){
    double log_spchi=log(spchi);
    
    int bin = 0;
    if( log_spchi > spchi_max) {
      bin = nspchi-1;
    } else if( log_spchi < spchi_min) {
      bin = 0;
    } else {
      bin=(int)(1+((log_spchi-spchi_min)/(spchi_max-spchi_min))*(double)(nspchi-2));
    }
    
    if(region==1) return spchi_cc[bin];
    if(region==2) return spchi_ec[bin];
  }
  return 0;
}

double LhoodCalc::ndr05SigLookup(int region, double ndr05){

  const int nndr05=6;

  //N 0.05 Vector CC
  double ndr05_cc[nndr05] = {0, 0.969176, 0.0245675, 0.00598086, 0.00027604, 0};

  //N 0.05 Vector EC
  double ndr05_ec[nndr05] = {0, 0.955311, 0.0371445, 0.00754498, 0, 0};
  
  int bin = 0;
  if( ndr05 > ndr05_max) {
    bin = nndr05-1;
  } else if( ndr05 < ndr05_min) {
    bin = 0;
  } else {
    bin=(int)(1+((ndr05-ndr05_min)/(ndr05_max-ndr05_min))*(double)(nndr05-2));
  }
  
  if(region==1) return ndr05_cc[bin];
  if(region==2) return ndr05_ec[bin];
  return 0;
}
double LhoodCalc::ndr05BkgLookup(int region, double ndr05){

  const int nndr05=6;

  //N 0.05 Vector CC
  double ndr05_cc[nndr05] = {0, 0.639241, 0.318864, 0.0359831, 0.00532701, 0.000584253};

  //N 0.05 Vector EC
  double ndr05_ec[nndr05] = {0, 0.610383, 0.313193, 0.0615295, 0.0123059, 0.00258794};

  int bin = 0;
  if( ndr05 > ndr05_max) {
    bin = nndr05-1;
  } else if( ndr05 < ndr05_min) {
    bin = 0;
  } else {
    bin=(int)(1+((ndr05-ndr05_min)/(ndr05_max-ndr05_min))*(double)(nndr05-2));
  }
  
  if(region==1) return ndr05_cc[bin];
  if(region==2) return ndr05_ec[bin];
  return 0;
}

double LhoodCalc::tptSigLookup(int region, double tpt){

  const int ntpt=12;

//TrkPT in 0.2 cone Vector CC
  float tpt_cc[ntpt] = {0, 0.686938, 0.105405, 0.0743671, 0.0480433, 0.0289802, 0.0169678, 0.0101295, 0.0063455, 0.00408079, 0.00240729, 0.0163352};

//TrkPT in 0.2 cone Vector EC
  float tpt_ec[ntpt] = {0, 0.720638, 0.0970038, 0.066375, 0.0409545, 0.0233433, 0.0129543, 0.00746791, 0.0045935, 0.00289109, 0.00162123, 0.0221569};

  int bin = 0;
  if( tpt > tpt_max) {
    bin = ntpt-1;
  } else if( tpt < tpt_min) {
    bin = 0;
  } else {
    bin=(int)(1+((tpt-tpt_min)/(tpt_max-tpt_min))*(double)(ntpt-2));
  }
  
  if(region==1) return tpt_cc[bin];
  if(region==2) return tpt_ec[bin];
  return 0;
}
  double LhoodCalc::tptBkgLookup(int region, double tpt){
    
    const int ntpt=12;
    
    //TrkPT in 0.2 cone Vector CC
    float tpt_cc[ntpt] = {0, 0.140013, 0.052951, 0.0514434, 0.0501231, 0.0489616, 0.0477373, 0.0462756, 0.0445764, 0.0427669, 0.0409401, 0.434212};
    
    //TrkPT in 0.2 cone Vector EC
    float tpt_ec[ntpt] = {0, 0.183827, 0.0555695, 0.0521187, 0.0487344, 0.0453987, 0.0420704, 0.0387871, 0.0356573, 0.0327389, 0.0299706, 0.435127};
    
    int bin = 0;
    if( tpt > tpt_max) {
      bin = ntpt-1;
    } else if( tpt < tpt_min) {
      bin = 0;
    } else {
      bin=(int)(1+((tpt-tpt_min)/(tpt_max-tpt_min))*(double)(ntpt-2));
    }
    if(region==1) return tpt_cc[bin];
    if(region==2) return tpt_ec[bin];
    return 0;
  }
  
  double LhoodCalc::cpsBkgLookup(int region, double cps){
    
    const int ncps=32;
    
    //Max # CPS strips
  double cps_cc[ncps] = {0.0616443, 0.0027037, 0.0110685, 0.0184619, 0.0245324, 0.0296939, 0.034624, 0.0397989, 0.045194, 0.0503822, 0.0547949, 0.0578636, 0.0592989, 0.0589618, 0.0570981, 0.053936, 0.0498704, 0.0450527, 0.0398345, 0.0344883, 0.0294748, 0.0249753, 0.0210755, 0.0176363, 0.0145977, 0.0119185, 0.00967101, 0.00795755, 0.00689832, 0.00648442, 0.00648888, 0.0135185};
 
  int bin = 0;
  if( cps > cps_max) {
    bin = ncps-1;
  } else if( cps < cps_min) {
    bin = 0;
  } else{
  bin=(int)(1+((cps-cps_min)/(cps_max-cps_min))*(double)(ncps-2));
  }
  
  if(region==1 && cps != -2) return cps_cc[bin];
  if(region==1 && cps == -2) return 1.0;
  
  return 0;
}

double LhoodCalc::emfSigLookup_not42(int region, double emf){

  const int nemf=22;

  //EMF Vector CC
  float emf_cc[nemf] = {0, 0.000270832, 0.0014032, 0.0029368, 0.00541998, 0.00941002, 0.0150986, 0.0219374, 0.0487498, 0.0831454, 0.405977, 0.176312, 0.0912704, 0.0541664, 0.0371575, 0.0228246, 0.012631, 0.0064565, 0.00315302, 0.00140877, 0.000270832, 0};

  //EMF Vector EC
  float emf_ec[nemf] = {0, 0, 0.000133839, 0.000401517, 0.000887565, 0.00159903, 0.0024725, 0.00342346, 0.00941452, 0.0530637, 0.437347, 0.249057, 0.114686, 0.0556313, 0.0363373, 0.0204492, 0.00971742, 0.00383202, 0.00123273, 0.000313465, 0, 0};

  int bin =0;
  if( emf > emf_max_not42) {
    bin = nemf-1;
  } else if( emf < emf_min) {
    bin = 0;
  } else {
    bin=(int)(1+((emf-emf_min)/(emf_max_not42-emf_min))*(float)(nemf-2));
  }
  
  if(region==1) return emf_cc[bin];
  if(region==2) return emf_ec[bin];
  return 0;
}
double LhoodCalc::emfBkgLookup_not42(int region, double emf){

  const int nemf=22;

  //EMF Vector CC
  float emf_cc[nemf] = {0, 0.00797425, 0.0185885, 0.0298176, 0.0419251, 0.0548187, 0.0682079, 0.0818292, 0.0991937, 0.125352, 0.277459, 0.0690108, 0.0488143, 0.0266057, 0.0194034, 0.0130294, 0.00804785, 0.00461691, 0.0025121, 0.00130282, 0.00052168, 0.000968834};

  //EMF Vector EC
  float emf_ec[nemf] = {0, 0.000285174, 0.00441785, 0.00965659, 0.0169737, 0.0269436, 0.0394332, 0.0534701, 0.0795635, 0.133319, 0.389262, 0.137596, 0.0594587,0.0222436, 0.014308, 0.00782996, 0.00354414, 0.00127917, 0.000349719, 6.57191e-05, 0, 0};

  int bin = 0;
  if( emf > emf_max_not42) {
    bin = nemf-1;
  } else if( emf < emf_min) {
    bin = 0;
  } else {
    bin=(int)(1+((emf-emf_min)/(emf_max_not42-emf_min))*(float)(nemf-2));
  }
  
  if(region==1) return emf_cc[bin];
  if(region==2) return emf_ec[bin];
  return 0;
}

double LhoodCalc::hmxSigLookup_not42(int region, double hmx){

  const int nhmx=52;

  //HMX Vector CC
  float hmx_cc[nhmx] = {0, 0.0511915, 0.159763, 0.1637, 0.130792, 0.085788, 0.059911, 0.0455661, 0.0383502, 0.0319284, 0.0267394, 0.0227402, 0.0196265, 0.0170152, 0.0146898, 0.0125727, 0.010796, 0.00949149, 0.00868702, 0.00819277, 0.00774714, 0.00720311, 0.00657112, 0.00592639, 0.00528051, 0.00461842, 0.00396675, 0.00341346, 0.00301644, 0.00277221, 0.00261247, 0.00249094, 0.00236361, 0.0022386, 0.00212054, 0.00204414, 0.00199668, 0.00194691, 0.00183, 0.00163439, 0.00139826, 0.00120611, 0.00109962, 0.001073, 0.00106953, 0.00104754, 0.000974614, 0.000862336, 0.000736169, 0.000635467, 0.000562544, 0};

  //HMX Vector EC
  float hmx_ec[nhmx] = {0, 0.014841, 0.147537, 0.251424, 0.198171, 0.121347, 0.0689669, 0.0436499, 0.0326387, 0.0233913, 0.0167199, 0.0122974, 0.00930121, 0.00701273, 0.00515177, 0.00370396, 0.00283455, 0.00259026, 0.00284892, 0.00324051, 0.00345607, 0.00337703, 0.00311837, 0.00279503, 0.00243577, 0.002037, 0.00163463, 0.00130411, 0.00104904, 0.000847851, 0.000675406, 0.000560444, 0.000506555, 0.000484999, 0.000456259, 0.00041674, 0.000398777, 0.00041674, 0.000456259, 0.000474222, 0.000456259, 0.000423925, 0.000434703, 0.000524518, 0.000671814, 0.000793962, 0.000801147, 0.000668221, 0.000445481, 0.000211963, 0, 0};

  int bin = 0;
  if( hmx > hmx_max) {
    bin = nhmx-1;
  } else if( hmx < hmx_min) {
    bin = 0;
  } else {
    bin=(int)(1+((hmx-hmx_min)/(hmx_max-hmx_min))*(float)(nhmx-2));
  }
  
  if(region==1) return hmx_cc[bin];
  if(region==2) return hmx_ec[bin];
  return 0;
}

double LhoodCalc::hmxBkgLookup_not42(int region, double hmx){

  const int nhmx=52;
//HMX Vector CC
float hmx_cc[nhmx] = {0, 0.00620536, 0.0291349, 0.0382916, 0.0391997, 0.036135,
0.0334008, 0.0311598, 0.0294323, 0.0281393, 0.0272353, 0.0265987, 0.0260889, 0.0255094, 0.0247866, 0.0239551, 0.0231631, 0.022521, 0.0220716, 0.0217733, 0.0215696, 0.0213765, 0.0211043, 0.0206677, 0.0200673, 0.0193847, 0.0187344, 0.0181714, 0.0176824, 0.0172296, 0.0167893, 0.0163773, 0.0160079, 0.0157059, 0.0154327, 0.0151428, 0.014771, 0.0143397, 0.0139046, 0.0135356, 0.0132382, 0.0129953, 0.0127882, 0.0126187, 0.0124649, 0.0123151, 0.0121731, 0.0120931, 0.012094, 0.0121638, 0.0122594, 0};

//HMX Vector EC
float hmx_ec[nhmx] = {0, 0.0126615, 0.0708161, 0.0917223, 0.0867166, 0.0712578,
0.0597741, 0.0468181, 0.0421136, 0.037687, 0.0336731, 0.0301161, 0.027051, 0.0245851, 0.0227093, 0.0212976, 0.0200701, 0.0188675, 0.01764, 0.0164924, 0.0154715, 0.014587, 0.0137975, 0.0130814, 0.0124022, 0.0117145, 0.0110008, 0.0103138, 0.00972364, 0.00927408, 0.00890753, 0.00855915, 0.00816473, 0.00776546, 0.0073977, 0.00709234, 0.00680334, 0.00651736, 0.00625017, 0.00604721, 0.00590725, 0.00579456, 0.00565884, 0.00546193, 0.0051899, 0.00486697, 0.00454767, 0.00427261, 0.00403329, 0.00379276, 0.00353344, 0};

 int bin = 0;
 if( hmx > hmx_max) {
   bin = nhmx-1;
 } else if ( hmx < hmx_min) {
   bin = 0;
 } else {
   bin=(int)(1+((hmx-hmx_min)/(hmx_max-hmx_min))*(float)(nhmx-2));
 }
 
  if(region==1) return hmx_cc[bin];
  if(region==2) return hmx_ec[bin];
  return 0;
}

} // namespace emid