///////////////////////////////////////////////////////////////////
// 
// This class computes weights used for finding backgrounds for photons.
// Run 2 port - Ray Culbertson
// 
//  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) {
  // same as Run 1b except for this 5% correction
  double photoneff = phoEff(et) - 0.05;
  double backgreff = backEf(et) - 0.05;

  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);
  }
}

/////////////////////////////////////////////
double PhotonBackgroundComputer::cprWeight2b(
       const double cprph, const double pt, const double sth, const int nVert, const bool All4Pads) {

  double phoeff, baceff, nphobac;
  double cprsys1, cprsys2, cprsys3, cprsys4;
  return  getCprWeight2b( cprph, pt, sth, phoeff, baceff, nphobac,
			cprsys1, cprsys2, cprsys3, cprsys4, nVert, All4Pads);

}

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

  double phoeff, baceff, nphobac;
  double cprsys1, cprsys2, cprsys3, cprsys4;
  return  getCprWeight( cprph, pt, sth, phoeff, baceff, nphobac,
			cprsys1, cprsys2, cprsys3, cprsys4,
			-1.0, nVert, sumEtNojets);

}

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

  double phoeff, baceff, nphobac;
  double cprsys1, cprsys2, cprsys3, cprsys4;

  return  getCprWeight1b( cprph, pt, sth, phoeff, baceff, nphobac,
			  cprsys1, cprsys2, cprsys3, cprsys4);
}


           /////////////////////////////////////////////
double PhotonBackgroundComputer::getCprWeight2b(
    const double cprph, const double pt, const double sth, 
    double& phoeff, double& baceff, double& nphobac,
    double& cprsys1, double& cprsys2, double& cprsys3, double& cprsys4,
    const int nVert, const bool All4Pads) {

  double epho,backprob,cprweight;
  double truepp, ratpp, exppp = 1.4564199;
  double mat = 1.096;
  double ue0 = .934;
  double allded = .998;
  double cprcut = 125./sth;
  double pef2, bef2;
  
  //  get photon efficiency
  epho = pt/sth;
  truepp = pairCrossSection( 13., epho );
  ratpp = truepp/exppp;
  phoeff = 1. - exp((-7.*mat*ratpp)/(9.*sth));

  // The background Nphoton fits were done for a mixture of three components: Pizero, Eta, Kshort.
  // The official mixture is pi0:eta:Kshort = 1.0:1.0:0.4
  // Fit results using the 4 pads algorithm: add energy from the 4 CPR pads in the photon stntuple block;
  // Fit results using the 1+1 pads algorithm: count the energy in only one pad but if  within 2 cm of  
  // the pad margin add the second pad too.
  // 
  float back_4_Etmax = 100.0;
  float back_4_nphMax  = 1.815; 
  float back_2_Etmax = 95.0;
  float back_2_nphMax  = 1.791; 
  float back_4_pad[7] = {1.07140e+00, 2.18902e-02, -2.32828e-04, 1.05743e-06, -1.73011e-09};
  float back_2_pad[7] = {9.00140e-01, 2.60009e-02, -2.76363e-04, 1.26754e-06, -2.11909e-09 };
  
  nphobac = back_2_pad[0] + back_2_pad[1]*pt + back_2_pad[2]*pow(pt,2) +  back_2_pad[3]*pow(pt,3) +
                     back_2_pad[4]*pow(pt,4);
  if (All4Pads) 
        nphobac = back_4_pad[0] + back_4_pad[1]*pt + back_4_pad[2]*pow(pt,2) +  back_4_pad[3]*pow(pt,3) +
                     back_4_pad[4]*pow(pt,4);

  if (All4Pads && pt>back_4_Etmax) 
        nphobac = back_4_nphMax;
  if (!All4Pads && pt>back_2_Etmax) 
        nphobac = back_2_nphMax;
  
  //  get background efficiency 
  baceff = (1. - exp((-7.*mat*nphobac)/(9.*sth)));

  // use requested underlying event 
  double ue;
  if (nVert>=0) {
    ue = 1.0 - ( 0.0712 + 0.00141*nVert ) ;
    if (All4Pads) ue = 1.0 - ( 0.0931 + 0.0180*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;

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

  //   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) {
    cprsys1 =  (1. - bef2)/(pef2 - bef2);
  } else {
    cprsys1 = - 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) {
    cprsys2 =  (1. - bef2)/(pef2 - bef2);
  } else {
    cprsys2 = - bef2/(pef2 - bef2);
  }

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

  //   cpr photon weight
  if( cprph > cprcut) {
    cprsys3 =  (1. - bef2)/(phoeff - bef2);
  } else {
    cprsys3 = - bef2/(phoeff - bef2);
  }
 
  // systematic #4 for eta/pi0 and ksh/pi0 and else
  bef2 = baceff + 0.008; 
  
  //  cpr photon weight
  if( cprph > cprcut) {
    cprsys4 =  (1. - bef2)/(phoeff - bef2);
  } else {
    cprsys4 = - bef2/(phoeff - bef2);
  }
  
  return cprweight;
}

/////////////////////////////////////////////
double PhotonBackgroundComputer::getCprWeight(
    const double cprph, const double pt, const double sth, 
    double& phoeff, double& baceff, double& nphobac,
    double& cprsys1, double& cprsys2, double& cprsys3, double& cprsys4,
    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;

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

  //   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) {
    cprsys1 =  (1. - bef2)/(pef2 - bef2);
  } else {
    cprsys1 = - 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) {
    cprsys2 =  (1. - bef2)/(pef2 - bef2);
  } else {
    cprsys2 = - bef2/(pef2 - bef2);
  }


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

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

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

  return cprweight;
}

           /////////////////////////////////////////////
double PhotonBackgroundComputer::getCprWeight1b(
    const double cprph, const double pt, const double sth, 
    double& phoeff, double& baceff, double& nphobac,
    double& cprsys1, double& cprsys2, double& cprsys3, double& cprsys4,
    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) {
    cprsys1 =  (1. - bef2)/(pef2 - bef2);
  } else {
    cprsys1 = - 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) {
    cprsys2 =  (1. - bef2)/(pef2 - bef2);
  } else {
    cprsys2 = - bef2/(pef2 - bef2);
  }


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

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

  // systematic #4 for eta/pi0 and ksh/pi0 and else
  bef2 = baceff + 0.008;
  
  //  cpr photon weight
  if( cprph > cprcut) {
    cprsys4 =  (1. - bef2)/(phoeff - bef2);
  } else {
    cprsys4 = - 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);
  }

}

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