///////////////////////////////////////////////////////////////////
// 
// This class computes weights used for finding backgrounds for photons.
// Run 2 b    - Shin-Shan Yu (Eiko)
// Run 2 port - Ray Culbertson
// 
//  12/21/07- added corrections to CES, CPR weights and model CP2 
//  2/13/03 - added Run II Cpr version fixed Run I to final numbers
//  2/17/03 - fix bug in calibration numbers, Run I and Run II 
//  3/21/03 - remove calibration corrections, not needed for Run II
//  3/ 8/04 - for Run II add calibration const back in, add 5% correction 
//            to CES this is "Yanwen's tune" from Spring 2003
///////////////////////////////////////////////////////////////////

#include <cmath>
#include "HighLevelObjects/PhotonBackgroundComputer.hh"


/////////////////////////////////////////////
double PhotonBackgroundComputer::cesWeight(const double et, 
					   const double chisq)
{
  double sys[8]={0,0,0,0,0,0,0,0};
  double phoeff=0, baceff=0;
  return getCesWeight(et, chisq, phoeff, baceff, sys);

}

double PhotonBackgroundComputer::getCesWeight(const double et, 
					      const double chisq, 
					      double& photoneff,
					      double& backgreff,
					      double cessys[8]) {
  // same as Run 1b except for this 5% correction
  photoneff = phoEff(et) - 0.05;
  backgreff = backEf(et) - 0.05;

  double phocorr = 0.0272;        // from W, Z gamma
  double bakcorr = et < 22.5? -0.06: 0.03-0.004*et;  // from isolation correction
  double cesw = 0;

  // my correction
  photoneff += phocorr;
  backgreff += bakcorr;

  if(chisq>20.0 || chisq<0.0) cesw = 0;
  if( photoneff==backgreff ) cesw = 0.5;
  else{
    if( chisq<4.0 ) {
      cesw = (1.0 - backgreff)/(photoneff-backgreff);
    } else {
      cesw = -backgreff/(photoneff-backgreff);
    }
  }

  // now look at systematic uncertainties
  double signal_syseff = photoneff;
  double backgr_syseff = backgreff;


  // #0: uncertainty on correction from W,Z gamma 
  signal_syseff = photoneff - 0.009;
  backgr_syseff = backgreff;

  if(chisq>20.0 || chisq<0.0) cessys[0] = 0;
  if( signal_syseff==backgr_syseff ) cessys[0] = 0.5;
  else {
    if( chisq<4.0 ) {
      cessys[0] = (1.0 - backgr_syseff)/(signal_syseff-backgr_syseff);
    } else {
      cessys[0] = -backgr_syseff/(signal_syseff-backgr_syseff);
    }
  }

  // #1: uncertainty on isolation formula
  signal_syseff = photoneff;
  backgr_syseff = backgreff + 0.04 - 0.001*et; 

  if(chisq>20.0 || chisq<0.0) cessys[1] = 0;
  if( signal_syseff==backgr_syseff ) cessys[1] = 0.5;
  else{
    if( chisq<4.0 ) {
      cessys[1] = (1.0 - backgr_syseff)/(signal_syseff-backgr_syseff);
    } else {
      cessys[1] = -backgr_syseff/(signal_syseff-backgr_syseff);
    }
  }
  // #2: gas saturation (guess from run I's paper)
  signal_syseff = photoneff + 0.045*et/30 - 0.015;
  backgr_syseff = backgreff + 0.015;


  if(chisq>20.0 || chisq<0.0) cessys[2] = 0;
  if( signal_syseff==backgr_syseff ) cessys[2] = 0.5;
  else {
    if( chisq<4.0 ) {
      cessys[2] = (1.0 - backgr_syseff)/(signal_syseff-backgr_syseff);
    } else {
      cessys[2] = -backgr_syseff/(signal_syseff-backgr_syseff);
    }
  }

  // unmodifed   CES weights
  for(int is=3; is<8; is++)
    cessys[is] = cesw;

  return cesw;
}

/////////////////////////////////////////////
double PhotonBackgroundComputer::cprWeight(
    const double cprph, const double pt, const double sth, 
    const int nVert, const double sumEtNojets) {

  double phoeff=0, baceff=1, nphobac=-1;
  double cprsys[8]={0,0,0,0,0,0,0,0};
  return  getCprWeight( cprph, pt, sth, phoeff, baceff, nphobac,
			cprsys,
			-1.0, nVert, sumEtNojets);

}


/////////////////////////////////////////////
double PhotonBackgroundComputer::getCprWeight(
    const double cprph, const double pt, const double sth, 
    double& phoeff, double& baceff, double& nphobac,
    double cprsys[8],
    const double nPhoBac, const int nVert, 
    const double sumEtNojets) {

  double epho,backprob,cprweight;
  double cprcut = 500.0;
  double pef2, bef2;
  double truepp, ratpp, exppp = 1.4564199;

  //double mat = 1.0748; // Run I number
  double mat = 1.072;

  //double crack = .9812, chimney = .9896, dead = .9941;
  //double ue = .937;  // Run I
  double ue0 = .934;

  //double allded = .96527;  // Run I
  double allded = .978;

  double nphoa = .66641, nphob = .12889, nphoc = -.48781e-2;
  double nphod = .95325e-4, nphoe = -.97846e-6, nphof = .50027e-8;
  double nphog = -.10038e-10;


  //  get photon efficiency
  epho = pt/sth;
  truepp = pairCrossSection( 13., epho );
  ratpp = truepp/exppp;
  phoeff = 1. - exp((-7.*mat*ratpp)/(9.*sth));

  //  get background efficiency

    // set a constant
  if(nPhoBac>0.0001) {
    nphobac = nPhoBac;

    // -2 means pi0 fit from Run I
  } else if( fabs(nPhoBac+2.0)<0.001 ) {
    nphobac = 0.82168 + 0.11238*pt - 0.42197e-2*pow(pt,2) + 
	0.78249e-4*pow(pt,3) - 0.75879e-6*pow(pt,4) + 
	0.36828e-8*pow(pt,5) - 0.70565e-11*pow(pt,6);
    if(pt>125.0) nphobac = 1.98671;

    // -1 means background fit from Run I
  } else {
    nphobac = nphoa + nphob*pt + nphoc*pow(pt,2) + nphod*pow(pt,3)
      + nphoe*pow(pt,4) + nphof*pow(pt,5) + nphog*pow(pt,6);
    if(pt>125.0) nphobac = 2.24;
  }

  baceff = (1. - exp((-7.*mat*nphobac)/(9.*sth)));

  // use requested underlying event 
  double ue;
  if(sumEtNojets>0.0) {
    ue = 1.0 - ( 0.0173 + 0.000889*sumEtNojets ) ;
    if(ue<0.1) ue = 0.1;
  } else if (nVert>=0) {
    ue = 1.0 - ( 0.035 + 0.0259*nVert ) ;
    if(ue<0.1) ue = 0.1;
  } else {
    ue = ue0;
  }

  // correct photon and background effs for cracks/ue
  phoeff = (1. - ((1.-phoeff)*ue))*allded;
  baceff = (1. - ((1.-baceff)*ue))*allded;

  // correct for backscattered showers
  backprob = 7.4e-4*pt/sth;
  phoeff = ((1. - phoeff)*backprob) + phoeff;
  baceff = ((1. - baceff)*backprob) + baceff;

  // Dana's calibration correction
  
  phoeff = phoeff + 0.010; 
  baceff = baceff + 0.033; 

  // Eiko's calibration from W and Z gamma

  phoeff = phoeff + 0.0174;
  baceff = baceff + 0.0174;


  //   cpr photon weight
  if( cprph > cprcut) {
    cprweight =  (1. - baceff)/(phoeff - baceff);
  } else {
    cprweight = - baceff/(phoeff - baceff);
  }


// #0: W gamma, Z gamma uncertainty
  pef2 = phoeff + .0221;
  bef2 = baceff + .0221;

  //   cpr photon weight
  if( cprph > cprcut) {
    cprsys[0] =  (1. - bef2)/(pef2 - bef2);
  } else {
    cprsys[0] = - bef2/(pef2 - bef2);
  }

  // systematic #1: for backscattered showers
  backprob = -7.4e-4*pt/sth;
  pef2 = ((1. - phoeff)*backprob) + phoeff;
  bef2 = ((1. - baceff)*backprob) + baceff;

  //   cpr photon weight
  if( cprph > cprcut) {
    cprsys[1] =  (1. - bef2)/(pef2 - bef2);
  } else {
    cprsys[1] = - bef2/(pef2 - bef2);
  }


  // systematic #2 for eta/pi0 and ksh/pi0 and else
  bef2 = baceff - 0.008;

  //   cpr photon weight
  if( cprph > cprcut) {
    cprsys[2] =  (1. - bef2)/(phoeff - bef2);
  } else {
    cprsys[2] = - bef2/(phoeff - bef2);
  }

  // systematic #3 for hit rate error
  pef2 = phoeff + .0078;
  bef2 = baceff + .006;

  //  cpr photon weight
  if( cprph > cprcut) {
    cprsys[3] =  (1. - bef2)/(pef2 - bef2);
  } else {
    cprsys[3] = - bef2/(pef2 - bef2);
  }


  // systematic #7 for isolation correction
  // need to consider twice how to include this uncertainty
  double isocorr = 0;
  if(pt > 35 && pt < 50)isocorr=-0.03;
  else if(pt>50)isocorr =  -0.075;
  bef2 = baceff + isocorr;
  
  //  cpr photon weight
  if( cprph > cprcut) {
    cprsys[7] =  (1. - bef2)/(phoeff - bef2);
  } else {
    cprsys[7] = - bef2/(phoeff - bef2);
  }

  for(int is=4; is<7; is++)
    cprsys[is] = cprweight;

  return cprweight;
}


/////////////////////////////////////////////
double PhotonBackgroundComputer::isoWeight(const double iso, const double pt) {
  double w;
  double phoeff = isoPhoProb(pt);
  double baceff = isoBkgProb(pt);
  if( iso < 2.0 ) {
    w =  (1. - baceff)/(phoeff - baceff);
  } else {
    w = - baceff/(phoeff - baceff);
  }
  return w;
}

/////////////////////////////////////////////
double PhotonBackgroundComputer::cprWeight1b(
    const double cprph, const double pt, const double sth) {

  double phoeff=0, baceff=1, nphobac=-1;
  double cprsys[4]={0,0,0,0};

  return  getCprWeight1b( cprph, pt, sth, phoeff, baceff, nphobac,
			  cprsys);
}




/////////////////////////////////////////////
double PhotonBackgroundComputer::getCprWeight1b(
    const double cprph, const double pt, const double sth, 
    double& phoeff, double& baceff, double& nphobac,
    double cprsys[4],
    const double nPhoBac) {

  double epho,backprob,cprweight;
  double cprcut = 500.0;
  double pef2, bef2;
  double truepp, ratpp, exppp = 1.4564199;

  double mat = 1.0748;
  //double crack = .9812, chimney = .9896, dead = .9941;
  double ue = .937;
  double allded = .96527;

  double nphoa = .66641, nphob = .12889, nphoc = -.48781e-2;
  double nphod = .95325e-4, nphoe = -.97846e-6, nphof = .50027e-8;
  double nphog = -.10038e-10;

  //     save cprcut,exppp,mat,crack,chimney,dead,ue,allded
  //      save nphoa,nphob,nphoc,nphod,nphoe,nphof,nphog,epho,backprob

  //  get photon efficiency
  epho = pt/sth;
  //truepp = qprsgg(13.,epho);
  //std::cout << "epho " << epho << std::endl;
  truepp = pairCrossSection( 13., epho );
  ratpp = truepp/exppp;
  phoeff = 1. - exp((-7.*mat*ratpp)/(9.*sth));

  //  get background efficiency

    // set a constant
  if(nPhoBac>0.0001) {
    nphobac = nPhoBac;

    // -2 means pi0 fit from Run I
  } else if( fabs(nPhoBac+2.0)<0.001 ) {
    nphobac = 0.82168 + 0.11238*pt - 0.42197e-2*pow(pt,2) + 
	0.78249e-4*pow(pt,3) - 0.75879e-6*pow(pt,4) + 
	0.36828e-8*pow(pt,5) - 0.70565e-11*pow(pt,6);
    if(pt>125.0) nphobac = 1.98671;
    // -1 means background fit from Run I
  } else {
    nphobac = nphoa + nphob*pt + nphoc*pow(pt,2) + nphod*pow(pt,3)
      + nphoe*pow(pt,4) + nphof*pow(pt,5) + nphog*pow(pt,6);
    if(pt>125.0) nphobac = 2.24;
  }

  baceff = (1. - exp((-7.*mat*nphobac)/(9.*sth)));

  //std::cout <<"true "<< truepp << " " << phoeff << " " << baceff << std::endl;

  // correct photon and background effs for cracks/ue
  phoeff = (1. - ((1.-phoeff)*ue))*allded;
  baceff = (1. - ((1.-baceff)*ue))*allded;

  // correct for backscattered showers
  backprob = 7.4e-4*pt/sth;
  phoeff = ((1. - phoeff)*backprob) + phoeff;
  baceff = ((1. - baceff)*backprob) + baceff;

  // calibration correction
  phoeff = phoeff + 0.010;
  baceff = baceff + 0.033;

  //  std::cout <<"true " << phoeff << " " << baceff <<  std::endl;
  //   cpr photon weight
  if( cprph > cprcut) {
    cprweight =  (1. - baceff)/(phoeff - baceff);
  } else {
    cprweight = - baceff/(phoeff - baceff);
  }


  // sys for hit rate error
  pef2 = phoeff + .0078;
  bef2 = baceff + .006;
  //   cpr photon weight
  if( cprph > cprcut) {
    cprsys[0] =  (1. - bef2)/(pef2 - bef2);
  } else {
    cprsys[0] = - bef2/(pef2 - bef2);
  }

  // systematic for backscattered showers
  backprob = -7.4e-4*pt/sth;
  pef2 = ((1. - phoeff)*backprob) + phoeff;
  bef2 = ((1. - baceff)*backprob) + baceff;


  //   cpr photon weight
  if( cprph > cprcut) {
    cprsys[1] =  (1. - bef2)/(pef2 - bef2);
  } else {
    cprsys[1] = - bef2/(pef2 - bef2);
  }


  // systematic #3 for eta/pi0 and ksh/pi0 and else
  bef2 = baceff - 0.008;

  //   cpr photon weight
  if( cprph > cprcut) {
    cprsys[2] =  (1. - bef2)/(phoeff - bef2);
  } else {
    cprsys[2] = - bef2/(phoeff - bef2);
  }

  // systematic #4 for eta/pi0 and ksh/pi0 and else
  bef2 = baceff + 0.008;
  
  //  cpr photon weight
  if( cprph > cprcut) {
    cprsys[3] =  (1. - bef2)/(phoeff - bef2);
  } else {
    cprsys[3] = - bef2/(phoeff - bef2);
  }

  return cprweight;
}

/////////////////////////////////////////////
double PhotonBackgroundComputer::
cesWeight1b( const double et, const double chisq) {
  double photoneff = phoEff(et);
  double backgreff = backEf(et);

  if(chisq>20.0 || chisq<0.0) return 0.0;
  if( photoneff==backgreff ) return 0.5;

  if( chisq<4.0 ) {
     return  (1.0 - backgreff)/(photoneff-backgreff);
  } else {
     return -backgreff/(photoneff-backgreff);
  }

}


//_____________________________________________________________________________
// 
//  Core function
//


//_____________________________________________________________________________
double PhotonBackgroundComputer::cprWeight2b(
       const double cprph, const double corre, const double sth, 
       const int nVert, const bool UseDifferentMat) {

  double phoeff=0, baceff=1, nphobac=-1;
  double cprsys[8]={0,0,0,0,0,0,0,0};
  return  getCprWeight2b( cprph, corre, sth, phoeff, baceff, nphobac,
			  cprsys, nVert, UseDifferentMat);
}



//_____________________________________________________________________________
double PhotonBackgroundComputer::getCprWeight2b(
    const double cprph, const double epho, const double sth, 
    double& phoeff, double& baceff, double& nphobac,
    double cprsys[8],
    const int nVert, const bool UseDifferentMat,		
    const bool All4Pads) {


  double backprob,cprweight;
  double mat    = UseDifferentMat? 1.125: 1.105;

  double cprcut = 125.;
  double sef2, bef2;
  double pt     = epho*sth;
  //  get photon efficiency
  phoeff = conversionRate(1.0, epho, sth, mat);

  
  //  get background efficiency from MC
  baceff = CP2BackEff(epho, sth, mat, nphobac, 
   		      TVector3(0,0,0), All4Pads);


  // use requested underlying event 
  double ue = CP2UnderlyingEvent(nVert, All4Pads);
  
  // correct photon and background effs for cracks/ue
  phoeff = (1. - ((1.-phoeff)*(1.-ue)));
  baceff = (1. - ((1.-baceff)*(1.-ue)));
  
  
  double mycorr = 0.0316; // from W and Z gamma

  phoeff = phoeff + mycorr;
  baceff = baceff + mycorr;
  

  const int ndim=4;

  double par[ndim]={
    -1.39186e-02,
    2.79093e-03,
    -7.03511e-05,
    2.97155e-07
  };

  double parerr[ndim]={
    6.16026e-02,
    3.63753e-03,
    6.58074e-05,
    3.68036e-07
  };

  double isocorr = 0;
  double temppt= pt;
  if(pt < 27.5)temppt=27.5;
  else if(pt > 42.5)temppt=42.5;
  for(int ip=0; ip<ndim; ip++)isocorr += par[ip]*pow(temppt,ip);

  baceff = baceff + isocorr;

  //   cpr photon weight
  if( cprph > cprcut) {
    cprweight =  (1. - baceff)/(phoeff - baceff);
  } else {
    cprweight = - baceff/(phoeff - baceff);
  }

  // #0 systematics due to the uncertainty of measured photon efficiency
  sef2   = phoeff + 0.0106;
  bef2   = baceff + 0.0106;

  //   cpr photon weight #0
  if( cprph > cprcut) {
    cprsys[0] =  (1. - bef2)/(sef2 - bef2);
  } else {
    cprsys[0] = - bef2/(sef2 - bef2);
  }

  // #1 systematics: amount of material
  sef2 = conversionRate(1.0, epho, sth, mat-0.018);
  bef2 = CP2BackEff(epho, sth, mat-0.018, nphobac, 
   		    TVector3(0,0,0), All4Pads);

  sef2 = (1. - ((1.-sef2)*(1.-ue)));
  bef2 = (1. - ((1.-bef2)*(1.-ue)));

  sef2 += mycorr;
  bef2 += (mycorr + isocorr);

  if( cprph > cprcut) {
    cprsys[1] =  (1. - bef2)/(sef2 - bef2);
  } else {
    cprsys[1] = - bef2/(sef2 - bef2);
  }

  // #2 systematic for backscattered showers

  backprob = -7.4e-4*pt/sth;
  sef2 = ((1. - phoeff)*backprob) + phoeff;
  bef2 = ((1. - baceff)*backprob) + baceff;


  if( cprph > cprcut) {
    cprsys[2] =  (1. - bef2)/(sef2 - bef2);
  } else {
    cprsys[2] = - bef2/(sef2 - bef2);
  }


  // #3 systematic for underlying event

  sef2 = phoeff;
  bef2 = baceff;

  sef2 = (1. - ((1.-sef2)*(1.-ue)));
  bef2 = (1. - ((1.-bef2)*(1.-ue)));


  if( cprph > cprcut) {
    cprsys[3] =  (1. - bef2)/(sef2 - bef2);
  } else {
    cprsys[3] = - bef2/(sef2 - bef2);
  }


  // systematic #4 for Nphoton (how function describes MC)


  sef2 = phoeff;

  bef2 = (1- 0.0104)*CP2BackEff(epho, sth, mat, nphobac, 
				TVector3(0,0,0), All4Pads);
  bef2 = (1. - ((1.-bef2)*(1.-ue)));
  bef2 += (mycorr + isocorr);


  if( cprph > cprcut) {
    cprsys[4] =  (1. - bef2)/(sef2 - bef2);
  } else {
    cprsys[4] = - bef2/(sef2 - bef2);
  }

  // #5: remove eta

  sef2 = phoeff;

  bef2 = CP2BackEff(epho, sth, mat, nphobac, 
 		    TVector3(0,0,0), All4Pads, +3);
  bef2 = (1. - ((1.-bef2)*(1.-ue)));
  bef2 += (mycorr + isocorr);



  if( cprph > cprcut) {
    cprsys[5] =  (1. - bef2)/(sef2 - bef2);
  } else {
    cprsys[5] = - bef2/(sef2 - bef2);
  }


  // #6: only pi0 

  sef2 = phoeff;

  bef2 = CP2BackEff(epho, sth, mat, nphobac, 
 		    TVector3(1,0,0), All4Pads);
  bef2 = (1. - ((1.-bef2)*(1.-ue)));
  bef2 += (mycorr+isocorr);



  if( cprph > cprcut) {
    cprsys[6] =  (1. - bef2)/(sef2 - bef2);
  } else {
    cprsys[6] = - bef2/(sef2 - bef2);
  }


//   // systematic #3 for eta/pi0 and ksh/pi0 and smearing of generator level 
//   very small and can be ignored
//   bef2 = CP2BackEff(epho, sth, mat, nphobac, 
// 		    TVector3(0,0,0), All4Pads,+1);
//   bef2 = CP2BackEff(epho, sth, mat, nphobac, 
// 		    TVector3(0,0,0), All4Pads,-1);



  // #7: isolation correction
  // need to consider twice if you want to include this uncertainty
  sef2 = phoeff;

  double isocorr_highpt = 0;
  temppt = pt;
  if(pt > 100)temppt=100;
  if(pt > 42.5)
    for(int ip=0; ip<ndim; ip++)isocorr_highpt += par[ip]*pow(temppt,ip);
  else isocorr_highpt = 0;
  bef2 = baceff + isocorr_highpt;

  if( cprph > cprcut) {
    cprsys[7] =  (1. - bef2)/(sef2 - bef2);
  } else {
    cprsys[7] = - bef2/(sef2 - bef2);
  }
  
  return cprweight;
}




/////////////////////////////////////////////
double PhotonBackgroundComputer::phoEff( const double pt  ) {
  //    port of Run 1B fortran routine - 2/2001
  //  Description:
  //  ------------
  //  This function returns the efficiency for a photon within
  //  chisquare<20 to have chisquare<4 (i.e. N(Chisq<4)/N(Chisq<20) ).  The
  //  efficiency was first calculated using QFL.  QFL uses 1990 testbeam
  //  electron data and parameterizes electron/photon differences based
  //  on a differnce in shower maximum of 0.6 radiation lengths.
  //  It estimates the affects of gas saturation using testbeam data
  //  taken at two HV for 10 GeV and 50 GeV electrons.
  //  The resulting efficiencies are smoothed at level 2 with SMCTRL in
  //  C$JET:JETLIB (topdrawer smoothing function).  This routine simply
  //  linearly interpolates on the smoothed values.
  //
  //  Isolation: Less than 15% of the Et of EM cluster within cone of 0.7
  //  Track:     No track as defined by ELES (3D tracks)
  //  Fiducial:  14<|Z|<217, and |X|<17.5, for lead strip and wire clust resp
  //  Z_Vertex:  |Z_Vertex|<50 (Z vertex dist had mean of -2, sigma of 31)
  //  Eta:       |Eta|<0.9, This is physics eta of lead strip cluster
  //  Chisquare: Strip Chisquare<20, just like in electron stream.
  //  2nd Clust: 2nd Strip Clust absolute Energy<1.00 GeV
  //             2nd Wire Clust<1.00 GeV or |2nd X - 1st X|<7 cm. (diff view)
  //
  //  The efficiencies are:
  //  ---------------------
  //
  //  Mean Pt   Sta Err    Efficiency     Stat Err.      Smoothed Efficiency
  //  -------   -------    ----------     ---------      -------------------
  //   6.0031919   0.012     0.7592181     0.0090117      0.7592181    
  //   7.9860773   0.012     0.8117965     0.0082940      0.7791888    
  //   9.9974022   0.012     0.7876067     0.0086669      0.7916738    
  //  12.0166178   0.012     0.7940125     0.0085488      0.7989987    
  //  14.0012236   0.011     0.8062965     0.0082639      0.8029454    
  //  16.0219593   0.012     0.8137615     0.0083379      0.8041915    
  //  18.0043316   0.012     0.7963701     0.0084727      0.8045152    
  //  20.0179405   0.012     0.8015803     0.0083558      0.8043801    
  //  21.9937038   0.012     0.8019982     0.0084921      0.8035865    
  //  23.9993629   0.012     0.8087902     0.0083709      0.8020113    
  //  25.9964428   0.012     0.7949674     0.0087151      0.7999562    
  //  27.9866600   0.012     0.8000930     0.0086231      0.7976870    
  //  29.9861546   0.012     0.7949426     0.0086572      0.7953099    
  //  32.0117874   0.012     0.7903603     0.0088054      0.7931799    
  //  33.9873238   0.012     0.7809121     0.0089686      0.7915941    
  //  35.9735260   0.012     0.7987249     0.0085561      0.7901247    
  //  38.0142784   0.012     0.7853081     0.0089390      0.7878305    
  //  39.9864922   0.012     0.7978967     0.0085869      0.7843165    
  //  41.9912796   0.012     0.7673196     0.0092360      0.7799805    
  //  43.9938583   0.012     0.7801191     0.0088644      0.7749130    
  //  45.9858208   0.012     0.7669209     0.0092305      0.7689824    
  //  48.0016518   0.012     0.7645643     0.0093095      0.7623571    
  //  49.9860001   0.012     0.7572547     0.0092017      0.7551474    
  //  52.0127144   0.012     0.7424670     0.0094881      0.7471092    
  //  54.0094490   0.012     0.7247002     0.0097820      0.7385852    
  //  56.0148087   0.012     0.7451991     0.0094306      0.7302781    
  //  58.0150146   0.012     0.7203837     0.0098290      0.7220460    
  //  60.0038528   0.012     0.7173268     0.0100190      0.7139555    
  //  61.9903336   0.012     0.7083333     0.0100635      0.7069032    
  //  64.0139008   0.013     0.6923465     0.0103563      0.7007594    
  //  65.9923553   0.012     0.6945933     0.0101651      0.6938555    
  //  68.0141602   0.012     0.6831055     0.0102810      0.6863536    
  //  69.9921951   0.012     0.6977777     0.0102049      0.6787895    
  //  72.0039520   0.012     0.6563408     0.0104487      0.6713139    
  //  73.9832458   0.012     0.6641260     0.0106464      0.6641260    

  const int num_val = 35;
  int i;
  double pt_lo,phoeff;
  double pt_max = 74.;
  double pt_min = 6.0;
  double pt_bin = 2.0;

  double eff[num_val] = {
    0.7592181,0.7791888,0.7916738,0.7989987,
    0.8029454,0.8041915,0.8045152,0.8043801,
    0.8035865,0.8020113,0.7999562,0.7976870,
    0.7953099,0.7931799,0.7915941,0.7901247,
    0.7878305,0.7843165,0.7799805,0.7749130,
    0.7689824,0.7623571,0.7551474,0.7471092,
    0.7385852,0.7302781,0.7220460,0.7139555,
    0.7069032,0.7007594,0.6938555,0.6863536,
    0.6787895,0.6713139,0.6641260};

  if( (pt >= pt_min) && (pt < pt_max) ) {
    // pt within our values, interpolate to find efficiency.
    i = int( (pt-pt_min)/pt_bin );
    pt_lo = pt_min + i*pt_bin;
    //pt_hi = pt_lo + pt_bin;
    phoeff = eff[i]+(eff[i+1]-eff[i])*(pt-pt_lo)/pt_bin;
  } else if(pt < pt_min) {
    phoeff = eff[0]+(eff[1]-eff[0])*(pt-pt_min)/pt_bin;
  } else {
    //	  pt above our values, extrapolate to find efficiency.
    //	  use 5 bins for extra lever arm.
    phoeff = eff[num_val-1]+(eff[num_val-1]-eff[num_val-6])*
                                         (pt-pt_max)/(pt_bin*5);
  }     // end if 

  if( phoeff < 0.0 ) phoeff=0.0;
  if( phoeff > 1.0 ) phoeff=1.0;
  return phoeff;
}

/////////////////////////////////////////////
double PhotonBackgroundComputer::backEf( const double pt  ) {
  //    port of Run 1B fortran routine 2/2001
  //	Description:
  //	------------
  //	This function returns the efficiency for the photon background within
  //	chisquare<20 to have chisquare<4 (i.e. N(Chisq<4)/N(Chisq<20) ).  The
  //	efficiency was first calculated using QFL for a background that
  //	contained Pi0s, etas, and kshorts with production ratio 1,1.02,0.4.
  //	The cuts were:
  //
  //   Isolation: Less than 15% of the Et of EM cluster within cone of 0.7
  //   Track:     No track as defined by ELES (3D tracks)
  //   Fiducial:  14<|Z|<217, and |X|<17.5, for lead strip and wire clust resp.
  //   Z_Vertex:  |Z_Vertex|<50 (Z vertex dist had mean of -2, sigma of 31)
  //   Eta:       |Eta|<0.9, This is physics eta of lead strip cluster
  //   2nd Clust: 2nd Strip Clust absolute Energy<1.00 GeV
  //              2nd Wire Clust<1.00 GeV or |2nd X - 1st X|<7 cm. (diff view)
  //	Chisquare: Strip Chisquare<20, just like in electron stream.
  //
  //	The resulting efficiency was then smoothed using SMCTRL level 2.
  //	SMCTRL is a subroutine in JETLIB used by topdrawer to 
  //	smooth curves, to remove statistical fluctuations. 
  //
  //	The efficiencies are:
  //	---------------------
  //
  //	Mean Pt         Err      Efficiency    Stat Err.     Smoothed Efficiny
  //	-------      ----------  ---------     ---------     -----------------
  //   5.9978805       0.011    0.5222654     0.0103861     0.5222654    
  //   7.9951830       0.011    0.4387546     0.0100439     0.4693221    
  //  10.0165625       0.011    0.4205046     0.0098787     0.4315549    
  //  12.0148392       0.011    0.4056488     0.0099342     0.4049621    
  //  14.0056887       0.011    0.3953298     0.0098959     0.3876600    
  //  15.9756346       0.011    0.3566462     0.0096873     0.3811498    
  //  17.9896278       0.011    0.3877791     0.0096432     0.3837475    
  //  19.9920063       0.011    0.3913043     0.0095311     0.3938328    
  //  21.9921265       0.011    0.4197871     0.0097982     0.4101248    
  //  24.0194168       0.011    0.4356976     0.0097886     0.4325328    
  //  25.9849548       0.011    0.4525714     0.0097150     0.4586802    
  //  28.0096378       0.011    0.4847594     0.0095773     0.4841227    
  //  29.9945927       0.011    0.5045771     0.0095673     0.5074323    
  //  31.9946651       0.011    0.5421013     0.0096169     0.5289757    
  //  33.9883194       0.011    0.5446623     0.0094896     0.5486055    
  //  36.0095558       0.010    0.5789855     0.0093978     0.5665676    
  //  38.0097542       0.011    0.5649161     0.0094677     0.5837817    
  //  40.0143051       0.011    0.5925241     0.0094528     0.6001421    
  //  42.0038185       0.011    0.6290323     0.0092487     0.6137991    
  //  44.0104103       0.011    0.6238430     0.0093209     0.6233943    
  //  46.0075226       0.011    0.6268436     0.0092871     0.6288191    
  //  47.9859047       0.010    0.6404770     0.0089870     0.6306612    
  //  49.9962082       0.010    0.6207635     0.0091645     0.6309066    
  //  52.0044861       0.010    0.6244898     0.0089310     0.6316956    
  //  53.9942513       0.010    0.6357143     0.0090944     0.6335500    
  //  56.0117798       0.010    0.6450124     0.0090157     0.6355751    
  //  57.9979095       0.010    0.6392046     0.0090497     0.6373885    
  //  60.0031509       0.010    0.6310105     0.0090071     0.6394754    
  //  61.9920464       0.010    0.6473073     0.0089643     0.6420847    
  //  64.0003433       0.010    0.6375136     0.0091387     0.6443589    
  //  65.9967422       0.010    0.6473570     0.0090913     0.6447049    
  //  67.9907761       0.010    0.6472789     0.0088123     0.6420603    
  //  70.0068817       0.010    0.6449115     0.0091891     0.6363189    
  //  71.9785690       0.010    0.6186770     0.0091351     0.6279040    
  //  74.0104980       0.010    0.6180727     0.0091282     0.6180727    

    const int num_val = 35;
    int i;
    double pt_lo, backef;
    double pt_max = 74.;
    double pt_min = 6.0;
    double pt_bin = 2.0;
    double eff[num_val] = {
      0.5222654,0.4693221,0.4315549,0.4049621,
      0.3876600,0.3811498,0.3837475,0.3938328,
      0.4101248,0.4325328,0.4586802,0.4841227,
      0.5074323,0.5289757,0.5486055,0.5665676,
      0.5837817,0.6001421,0.6137991,0.6233943,
      0.6288191,0.6306612,0.6309066,0.6316956,
      0.6335500,0.6355751,0.6373885,0.6394754,
      0.6420847,0.6443589,0.6447049,0.6420603,
      0.6363189,0.6279040,0.6180727 };

    //	branch on interpolation or exptrapolation
    if( (pt >= pt_min) && (pt < pt_max) ) {
     //	  pt within our values, interpolate to find efficiency.
      i = int( (pt-pt_min)/pt_bin );
      pt_lo = pt_min + i*pt_bin;
      //pt_hi = pt_lo + pt_bin;
      backef=eff[i]+(eff[i+1]-eff[i])*(pt-pt_lo)/pt_bin;
    } else if(pt < pt_min) {
      //	  pt beneath our values, extrapolate to find efficiency.
      backef = eff[0]+(eff[1]-eff[0])*(pt-pt_min)/pt_bin;
    } else {
      // pt above our values, extrapolate to find efficiency.
      // use 5 bins for extra lever arm.  the efficiency is pretty flat.
	  backef = eff[num_val-1]+(eff[num_val-1]-eff[num_val-6])*
          (pt-pt_max)/(pt_bin*5);
    }
    if( backef < 0.0 ) backef = 0.0;
    if( backef > 1.0 ) backef = 1.0;
    return backef;
}

/////////////////////////////////////////////
double PhotonBackgroundComputer::pairCrossSection( const double z, const double energy  ) {
  // by Bob Blair
  // ported to c++ 2/2001
  // The tables below are basically a transcription of table III.5
  //    in Rev. Mod. Phys., Vol. 46, No. 4 p815 (Oct. 1974)
  //    which summarize the pair production cross section
  //    The article was written by Yung-Su Tsai.
  const double tsaidel[14][9] = {
      0.011,0.028,0.039,0.079,0.126,0.174,0.222,0.323,0.441
     ,0.012,0.023,0.030,0.058,0.091,0.128,0.166,0.253,0.367
     ,0.004,0.024,0.034,0.073,0.113,0.154,0.195,0.283,0.391
     ,0.003,0.020,0.029,0.064,0.101,0.139,0.178,0.263,0.370
     ,0.002,0.016,0.023,0.053,0.087,0.122,0.158,0.238,0.343
     ,0.002,0.015,0.022,0.050,0.082,0.116,0.151,0.230,0.334
     ,0.002,0.012,0.019,0.044,0.073,0.105,0.137,0.213,0.315
     ,0.002,0.011,0.017,0.040,0.068,0.098,0.129,0.202,0.302
     ,0.001,0.009,0.014,0.033,0.057,0.084,0.112,0.180,0.275
     ,0.001,0.009,0.013,0.032,0.056,0.082,0.110,0.177,0.272
     ,0.001,0.008,0.012,0.029,0.051,0.075,0.102,0.166,0.259
     ,0.001,0.007,0.011,0.028,0.049,0.073,0.099,0.162,0.256
     ,0.001,0.007,0.011,0.028,0.049,0.072,0.098,0.162,0.257
     ,0.001,0.007,0.011,0.028,0.048,0.072,0.098,0.162,0.258};
  const double tsaiz[14]= { 1, 2, 3, 4, 6, 7,10,13,26,29,50,74,82,92};
  const double tsaisig[14] ={20.73,55.06,108.8,179.4,361.5,473.8,896.1,
                             1443.,5182.,6343.,17276,34869,41720,50870};
  const double tsaie[9] ={ 100.,10.,6.,2.,1.,0.6,0.4,0.2,0.1 };
  const double emass = 0.000511;
  const double rad = 1.0093;

  if(energy<=2.0*emass) return 0.0;
  double sig = 0.0;

  // find energy bin
  int je = 0;
  while(energy<tsaie[je] && je<8) je++;

  // find nearest z
  double zdif = 1.0e10;
  int jz;
  for(int j=0; j<14; j++) {
    if(fabs(tsaiz[j]-z) < fabs(zdif) ) {
      zdif = z-tsaiz[j];
      jz = j;
    }
  }

  // Z<1 HUH?
  if(jz<1 && z<tsaiz[0]) jz=1;
  // Z>92 well maybe.
  if(jz>12 && z>tsaiz[13]) jz=12;

  // interpolate z if needed
  double siginf;
  if( zdif > 0.99 ) {
    siginf = tsaisig[jz] + (z*z - tsaiz[jz]*tsaiz[jz])*(tsaisig[jz+1] -
	     tsaisig[jz])/(tsaiz[jz+1]*tsaiz[jz+1] - tsaiz[jz]*tsaiz[jz]);
  } else if( zdif < -.99 ) {
    siginf = tsaisig[jz] + (z*z-tsaiz[jz]*tsaiz[jz])*(tsaisig[jz-1] -
      tsaisig[jz])/(tsaiz[jz-1]*tsaiz[jz-1] - tsaiz[jz]*tsaiz[jz]);
  } else {
    siginf = tsaisig[jz];
  }

  // interpolate in 1/log(e)  (at low energy vs. log(e))
  double sigpl,sigmn,xpl,xmn,xx;
  if(je == 9 ) {
    // at the low end interpolate between threshold and the first point
    sigpl=tsaidel[jz][8];
    sigmn=1.0;
    xpl=log(tsaie[8]/emass);
    xmn=log(2.);
    xx=log(energy/emass);
  } else if( je == 0 ) {
    // in the high energy end interpolate to infinity via 1/log(e)
    sigpl=0.0;
    sigmn=tsaidel[jz][0];
    xmn=1.0/log(tsaie[0]/emass);
    xpl=0.0;
    xx=1.0/log(energy/emass);
  } else {
    sigpl=tsaidel[jz][je-1];
    sigmn=tsaidel[jz][je];
    xmn=1.0/log(tsaie[je]/emass);
    xpl=1.0/log(tsaie[je-1]/emass);
    xx=1.0/log(energy/emass);
  }

  // calculate adjustment to value at infinity
  sig = sigmn + (xx-xmn)*(sigpl-sigmn)/(xpl-xmn);
  // include cross section at infinity and radiative correction (rad)
  sig=(1.0-sig)*siginf*rad;
  return sig/1000.0;
    
} // pair cross section


// Effective number of photons
//_____________________________________________________________________________
double PhotonBackgroundComputer::NPho_Pi0(const double e, const double sintheta, const bool all4pads)
{
  double ne = 0;
  if(all4pads)
    {
      double constant = 1.97793e+00*1.01031e+00;
      double slope = 2.54095e-02 + exp(-2.72722e+01+2.33353e+01*sintheta);
      double p2       = -1.28708e+00;
      ne      = constant*(1-exp(-slope*e+p2));
    }
  else
    {
      double constant = 1.97566e+00*1.00177e+00;
      double slope    = 3.35977e-02 + exp(-1.69838e+01+1.36113e+01*sintheta);
      double p2       = -7.10228e-01;
      ne      = constant*(1-exp(-slope*e+p2));


    } // if 1+2+4 pads

      return ne;
}


double PhotonBackgroundComputer::NPho_InclusiveEta(const double e, const double sintheta, const bool all4pads)
{
  double ne=0;
  if(all4pads)
    {
      double p0= 5.18863e-01+exp(-1.11620e+01+9.65079e+00*sintheta);
      double p1=-4.19489e-02;
      double p2=exp(6.29109e+00-2.99021e+00*sintheta);
      double p3=1.79321e+00;
      ne = p0*tanh(p1*(p2-e))+p3;
      ne *= 9.94298e-01;     
    }
  else
    {
       double p0 = 5.42730e-01+exp(-1.01235e+01+8.80641e+00*sintheta);
       double p1 =-3.06012e-02;
       double p2 = exp(6.47690e+00-2.98944e+00*sintheta);
       double p3 =1.69662e+00;
       ne = p0*tanh(p1*(p2-e))+p3;
       ne *= 9.96172e-01;

    } // 1+2+4 pad
  return ne;
}


double PhotonBackgroundComputer::NPho_Ks(const double e, const double sintheta,
 const bool all4pads)
{
  double ne=0;
  if(all4pads)
    {
      double constant = 3.90376e+00*9.88983e-01;
      double slope    = 3.44421e-02+exp(-1.24148e+01+8.86053e+00*sintheta);
      double p2       = -3.50599e-01;
      ne       = constant*(1-exp(-slope*e+p2));	
    }
  else
    {
       double constant = 3.88374e+00*9.95271e-01;
       double slope    = 2.83315e-02 + exp(-1.24119e+01+8.54268e+00*sintheta);
       double p2       = -1.97517e-01;
       ne       = constant*(1-exp(-slope*e+p2));	

    } // 1+2+4 pad
  return ne;
}


// particle fraction (returns an array of particle fraction given an energy)
//_____________________________________________________________________________

void PhotonBackgroundComputer::ParticleFraction_InclusiveEta(const double e, double* pfrac, const int sys)
{
  pfrac[0]=pfrac[1]=pfrac[2]=0;
  const int NPART=3;
  const int NBIN=20;
  // pi0, eta, ks
  double frac[NPART][NBIN] ;
  double mother[7][NPART][NBIN]=
    {
      // -1 sigma Ks/pi0
      {
 	{0.745778, 0.740077, 0.744917, 0.749959, 0.746884, 0.736854, 0.722492, 0.706172, 0.68903, 0.672289, 0.658187, 0.646876, 0.63816, 0.632017, 0.626632, 0.622396, 0.620226, 0.61859, 0.619049, 0.623774},
 	{0.204881, 0.207287, 0.204812, 0.204572, 0.21136, 0.224485, 0.242237, 0.262063, 0.282352, 0.301386, 0.317363, 0.3307, 0.341205, 0.3487, 0.355335, 0.360823, 0.364093, 0.366588, 0.366777, 0.362547},
 	{0.0493411, 0.0526357, 0.0502708, 0.0454694, 0.041756, 0.0386614, 0.0352705, 0.0317649, 0.0286187, 0.0263251, 0.02445, 0.022424, 0.0206346, 0.0192829, 0.0180331, 0.0167805, 0.0156807, 0.0148222, 0.0141741, 0.0136792}
      },

      // - 1 sigma eta/pi0
      {
	{0.710431, 0.70481, 0.709693, 0.714177, 0.709733, 0.697628, 0.680704, 0.661752, 0.642186, 0.623453, 0.60783, 0.595284, 0.58563, 0.57884, 0.572895, 0.568182, 0.565683, 0.563795, 0.564151, 0.56902},
	{0.247789, 0.250632, 0.247735, 0.247334, 0.254997, 0.269836, 0.289758, 0.311789, 0.334105, 0.354847, 0.372099, 0.386373, 0.397538, 0.405462, 0.41245, 0.418201, 0.421604, 0.424196, 0.424367, 0.419888}, 
	{0.04178, 0.0445577, 0.0425722, 0.0384889, 0.0352702, 0.0325362, 0.0295382, 0.0264594, 0.0237094, 0.0217003, 0.0200705, 0.0183427, 0.016832, 0.0156982, 0.0146548, 0.0136167, 0.0127126, 0.0120082, 0.0114819, 0.0110919}
      },

      // central value
      {
	{0.749889, 0.744431, 0.749101, 0.753767, 0.750365, 0.740033, 0.725335, 0.708673, 0.691228, 0.674261, 0.65998, 0.648492, 0.639627, 0.633374, 0.62789, 0.623559, 0.621309, 0.61961, 0.620025, 0.624724}, 
	{0.20601, 0.208506, 0.205962, 0.205611, 0.212346, 0.225453, 0.24319, 0.262992, 0.283252, 0.30227, 0.318227, 0.331526, 0.341989, 0.349448, 0.356049, 0.361497, 0.364729, 0.367193, 0.367356, 0.363099}, 
	{0.0441005, 0.0470625, 0.0449362, 0.0406225, 0.0372895, 0.0345139, 0.0314749, 0.0283355, 0.02552, 0.0234687, 0.0217925, 0.0199822, 0.018384, 0.0171772, 0.0160616, 0.0149439, 0.0139627, 0.013197, 0.0126191, 0.0121778}
      },
      // +1 sigma eta/pi0
      {
	{0.793989, 0.788772, 0.793144, 0.798003, 0.795932, 0.787926, 0.776229, 0.762756, 0.748379, 0.734085, 0.721918, 0.712145, 0.704592, 0.699254, 0.694563, 0.690895, 0.689067, 0.68769, 0.688184, 0.692517},
	{0.159317, 0.161362, 0.159278, 0.15899, 0.164514, 0.175326, 0.190088, 0.206746, 0.223991, 0.240364, 0.254244, 0.265912, 0.275157, 0.281782, 0.28767, 0.292547, 0.295447, 0.297663, 0.297809, 0.293984}, 
	{0.046694, 0.0498657, 0.0475781, 0.0430065, 0.0395539, 0.0367476, 0.0336834, 0.0304979, 0.02763, 0.025551, 0.0238377, 0.0219436, 0.0202512, 0.0189639, 0.0177671, 0.0165576, 0.0154854, 0.014647, 0.0140063, 0.0134993}
      },

      // +1 sigma Ks/pi0
      {
	{0.754046, 0.748836, 0.753333, 0.757614, 0.753879, 0.743239, 0.7282, 0.711192, 0.69344, 0.676245, 0.661783, 0.650116, 0.6411, 0.634737, 0.629153, 0.624726, 0.622395, 0.620634, 0.621005, 0.625676}, 
	{0.207152, 0.20974, 0.207126, 0.20666, 0.21334, 0.22643, 0.244151, 0.263926, 0.284159, 0.30316, 0.319097, 0.332356, 0.342777, 0.3502, 0.356765, 0.362174, 0.365366, 0.3678, 0.367936, 0.363652}, 
	{0.0388019, 0.0414234, 0.0395413, 0.0357261, 0.0327811, 0.0303305, 0.0276493, 0.0248817, 0.0224014, 0.0205956, 0.0191205, 0.0175282, 0.016123, 0.0150624, 0.0140822, 0.0131004, 0.0122387, 0.0115664, 0.0110592, 0.0106718} 
      },

      // remove eta
      {
	{0.944457, 0.94054, 0.943408, 0.948863, 0.952658, 0.95544, 0.958411, 0.961553, 0.964395, 0.966364, 0.968036, 0.970108, 0.972061, 0.973596, 0.975058, 0.976596, 0.978021, 0.979145, 0.980053, 0.98088},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0.0555429, 0.0594604, 0.056592, 0.0511368, 0.0473424, 0.0445601, 0.0415889, 0.0384466, 0.0356052, 0.0336358, 0.0319645, 0.0298923, 0.0279387, 0.0264041, 0.0249422, 0.0234045, 0.0219791, 0.0208547, 0.0199466, 0.0191203}
      },
      
      // remove Ks

      {
	{0.784486, 0.781196, 0.784347, 0.785683, 0.77943, 0.766487, 0.748906, 0.729339, 0.70933, 0.690465, 0.674683, 0.661714, 0.651606, 0.644444, 0.638139, 0.633019, 0.630107, 0.627896, 0.62795, 0.632425},
	{0.215514, 0.218804, 0.215653, 0.214317, 0.22057, 0.233513, 0.251094, 0.270661, 0.29067, 0.309535, 0.325317, 0.338286, 0.348394, 0.355556, 0.361861, 0.366981, 0.369893, 0.372104, 0.37205, 0.367575},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
      }
    }
  ;


  for(int i=0; i < NPART; i++)
    for(int j=0; j< NBIN; j++)
      frac[i][j] = mother[sys+2][i][j];

  const double elo=5.0;
  const double ehi=205;
  const double binwidth = (ehi-elo)/NBIN;
  int bin = int((e-elo)/binwidth);
  int lobin=0;
  int hibin=NBIN-1;

  if(e < elo)
    {
      for(int ipart=0; ipart < NPART; ipart ++)
	pfrac[ipart] = frac[ipart][lobin];
    }
  else if( e > ehi)
    {
      for(int ipart=0; ipart < NPART; ipart ++)
	pfrac[ipart] = frac[ipart][hibin];

    }
  // within 5-205 GeV
  else{
    // do the extrapolation
    double remain = (e-elo)-bin*binwidth;


    lobin = bin-1 + (int)(remain/0.5/binwidth);
    hibin = bin   + (int)(remain/0.5/binwidth);

    if(lobin < 0){lobin=0; hibin=1;}
    if(hibin >= NBIN){lobin= NBIN-2; hibin=NBIN-1;}

    double lov = elo + (lobin+0.5)*binwidth;
    double hiv = elo + (hibin+0.5)*binwidth;

    for(int ipart=0; ipart < NPART; ipart ++)
     pfrac[ipart] = ((e-lov)*frac[ipart][hibin]+(hiv-e)*frac[ipart][lobin])/binwidth;

  } 

  return;
  
}


// CP2 rate from underlying event of Inclusive Pho 25 triggers
//_____________________________________________________________________________

double PhotonBackgroundComputer::CP2UnderlyingEvent(int nvertex, bool all4Pads)
{
  
  double UE_4pad[11]={
    0.0205357,
    0.0351439,
    0.0549624,
    0.0748336,
    0.0943047,
    0.113478,
    0.131004,
    0.142239,
    0.144152,
    0.175549,
    0.159574,
  };

  double UE_124pad[11]={
    0.00678571,
    0.0141023,
    0.0220098,
    0.0299559,
    0.0388937,
    0.0443694,
    0.0524017,
    0.0575064,
    0.0617792,
    0.0532915,
    0.0638298
  };



  if(nvertex < 0 )nvertex=0;
  if(nvertex > 10)nvertex=10;
  
return ((all4Pads==true)? UE_4pad[nvertex]: UE_124pad[nvertex]);


}


//_____________________________________________________________________________
double PhotonBackgroundComputer::conversionRate(const double npho, 
						const double epho,
						const double sth,
						const double mat)
{

  double truepp, ratpp, exppp = 1.4564199;
  truepp = pairCrossSection( 13., epho );
  ratpp = truepp/exppp;

  double rate = 1. - exp((-7.*mat*npho*ratpp)/(9.*sth));
  return rate;

}

//_____________________________________________________________________________
// combined efficiency
double PhotonBackgroundComputer::CP2BackEff(
					    const double epho, 
					    const double sth, 
					    const double mat,
					    double& nbackpho,
					    const TVector3 externpfrac, 
					    const bool all4pads,
					    const int sys)

{
  double et = fabs(epho*sth);
  double effcorr_pi0 = (1.00633+exp(-3.075*et-0.044237)-1)*0.5+1;
  double effcorr_eta = (1.01024+exp(- 3.00596*et-0.0335248)-1)*0.5+1;
  double effcorr_ks  = (1.00678-1)*0.5+1;

  double eff_pi0 = conversionRate(NPho_Pi0(epho,sth,all4pads), epho, sth, mat)*effcorr_pi0;
  double eff_eta = conversionRate(NPho_InclusiveEta(epho,sth,all4pads), epho, sth, mat)*effcorr_eta;
  double eff_ks  = conversionRate(NPho_Ks(epho,sth,all4pads), epho, sth, mat)*effcorr_ks;

  double combineff = 0;

  double pfrac[3];
  if(externpfrac.Mag()==0)
    ParticleFraction_InclusiveEta(epho, pfrac, sys);
  else
    {
      for(int ndim=0; ndim < 3; ndim++)
	pfrac[ndim] = externpfrac[ndim];
    }
  combineff = eff_pi0*pfrac[0] + 
    eff_eta*pfrac[1] + eff_ks*pfrac[2];

  double truepp, ratpp, exppp = 1.4564199;
  truepp = pairCrossSection( 13., epho );
  ratpp = truepp/exppp;

  nbackpho = -log(1-combineff)*9.*sth/7./mat/ratpp;
  
  return combineff;
}



/////////////////////////////////////////////
double PhotonBackgroundComputer::isoPhoProb( const double pt  ) {
  return 0.98;
}

/////////////////////////////////////////////
double PhotonBackgroundComputer::isoBkgProb( const double pt  ) {
  if(pt>100.0) return 0.1957;
  return  (0.9621 - 0.007735*pt);
}