#ifndef PHOTONBACKGROUNDCOMPUTER_HH
#define PHOTONBACKGROUNDCOMPUTER_HH 1

///////////////////////////////////////////////////////////////////
// 
// This class computes background weights for photons.  The CES method
// weights events with CES average chi^2 less than 4 positively and 
// events with chi^2 > 4 negatively, the sum is the number of real photons.
// The CPR method weights events with a significant CPR hit postively and
// with a hit negatively, the sum is the number of real photons.
// 
// Typical use:
// 
// PhotonBackgroundComputer phoBg;
// ces weight = phoBg.cesWeight(et, chiSq);
// or
// cpr weight = phoBg.cprWeight(cpr_pulse_height, et, sin_theta);
// then sum the weights over the photon sample
// 
// Note: currently the Run 2 methods return Run 1 results.
// 
// Run 1b routines from Rob Harris, Steve Kuhlmann, Bob Blair, et al.
// Run 2 port - Ray Culbertson
// 
///////////////////////////////////////////////////////////////////

#include <iostream>
#include <cmath>

class PhotonBackgroundComputer {
public:
  // Run 2 CES weight from the photon Et and CES chisq
  static double cesWeight(const double et, const double chisq);

  // Run 1b CES weight from the photon Et and CES chisq
  static double cesWeight1b(const double et, const double chisq);


  // Run 2 CPR2 weight from the photon Et and sin(theta).
  // nVert is # of Class12 ZVertexColl vertices should drive  
  // corrections to the underlying event: not implemented yet.
  static double cprWeight2b(const double cprph,
			    const double pt, 
			    const double sth,
			    const int nVert=-1,
			    const bool All4Pads = false);

  // Run 2 CPR weight from CPR pulse height, the photon Et and sin(theta)
  // nVert and sumEtNojets drive corrections to the underlying event
  // hit rate, using one is recommended, but not critical.
  // if both are specified, then sumEtNojets is used
  // nVert is # of Class12 ZVertexColl vertices
  // sumEtNojets is Metblock raw sumet minus all jets with raw Et>10GeV.
  static double cprWeight(const double cprph, const double pt, 
			  const double sth,
			  const int nVert=-1, 
			  const double sumEtNojets=-1.0);

  // Run 1b CPR weight CPR pulse height, the photon Et and sin(theta)
  static double cprWeight1b(const double cprph, const double pt, 
			    const double sth);


  // the following are expert interfaces, not typically used


  // CPR2 function: nvertices uev corrections and systematics not implemented.
  // make 4Pads true if want to use 4 pads instead of 1+1 pads 
  static double 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=-1, const bool All4Pads = false);

  // same function, more general interface that return extra information
  // nPhoBac: -1 = Run I jet function, -2 = Run I pi0 function
  //              >0 = set to this value
  static double 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=-1.0,
    const int NVert=-1, const double sumEtNojets=-1.0);

  static double 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=-1.0);

  // probability of a photon with chisq<20 to have chisq<4
  static double phoEff(const double pt);

  // probability of background pi0/eta/Ks with chisq<20 to have chisq<4
  static double backEf(const double pt);

  // pair production cross section in barn/atom
  static double pairCrossSection( const double z, const double energy);
};

#endif