#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
// 
// 
// Run 1b routines from Rob Harris, Steve Kuhlmann, Bob Blair, et al.
// Run 2 port - Ray Culbertson
// Run 2b update of CES, CPR and CP2 weights - Shin-Shan Yu (Eiko)
// 
///////////////////////////////////////////////////////////////////

#include <iostream>
#include <cmath>
#include "TVector3.h"

class PhotonBackgroundComputer {
public:

  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...

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

  // same function but with more general interface
  //  input: photon Et, CES wire and strip average chi2   
  // output: photon and background CES chi2 < 4 efficiencies, 
  //  3 systematically modified CES weights (need to take the difference from 
  //                                 central value to get the uncertainty)
  static double getCesWeight(const double et, const double chisq,
			     double& phoeff,
			     double& baceff,
			     double cessys[8]);

  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...

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

  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...

  // Run 2a 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);


  // 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
  //  input: CPR ADC counts, photon Et, sintheta, nvertex
  // output: CPR ADC > 500 count rates of photon and background, 
  //         effective ngamma, 
  //  5 systematically modified CPR weights (need to take the difference from 
  //                                 central value to get the uncertainty)
  //   cprsys[0], cprsys[1], cprsys[2], cprsys[3] and cprsys[7]
  // recommend to take 50% of cprsys[7] (systematics due to isolation corr)  

  static double getCprWeight(
    const double cprph, const double pt, const double sth, 
    double& phoeff, double& baceff, double& nphobac,
    double cprsys[8],
    const double nPhoBac=-1.0,
    const int NVert=-1, const double sumEtNojets=-1.0);


  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...

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


  // same function, more general interface that return extra information
  //  input: CPR ADC counts, photon Et, sintheta
  // output: CPR ADC > 500 count rates of photon and background, 
  //         effective ngamma, 
  //   4 systematic modified CPR weights (need to take the difference from 
  //                                 central value to get the uncertainty)
  //                           

  static double getCprWeight1b(
    const double cprph, const double pt, const double sth, 
    double& phoeff, double& baceff, double& nphobac,
    double cprsys[4],
    const double nPhoBac=-1.0);



  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
  // Run 2 CP2 weight 
  // CP2 function from the CP2 ADC counts, corrected photon Energy, sin(theta),
  //                       nvertex, whether extrapolated TOF z position of 
  //                       photon is larger than TOF PMT bar Z position
  // calibration is mainly based on 1+2+4 pad algorithm
  // 4-pad algorithm is discouraged

  static double cprWeight2b(const double cprph,
			    const double corre, 
			    const double sth,
			    const int nVert,
			    const bool ReachTOFPMT);


  //  input: 1+2+4 pad CP2 ADC counts, corrected photon energy, sintheta,
  //  nvertex, whether extrapolated TOF z position of photon is larger 
  //  than TOF PMT bar Z position
  //  4-pad algorithm is discouraged. 
  // output: 1+2+4 pad CP2 ADC > 125 count rates of photon and background, 
  //         effective ngamma, 
  //  8 systematically modified CPR weights (need to take the difference from 
  //                                 central value to get the uncertainty)
  //  recommend to take 50% of cprsys[7] (systematics due to isolation corr)  
   

  static double getCprWeight2b(
    const double cprph, const double corre, const double sth, 
    double& phoeff, double& baceff, double& nphobac,
    double cprsys[8],
    const int NVert, 
    const bool ReachTOFPMT,
    const bool All4Pads = false);


  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...

  // weight based on fully corrected isolation
  static double isoWeight(const double iso, const double pt);


  // the following are expert interfaces, not typically used

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

  // effective number of photons for pi0
  static double NPho_Pi0(const double e, const double sintheta, const bool all4Pads=false);

  // effective number of photons for inclusive eta
  static double NPho_InclusiveEta(const double e, const double sintheta, const bool all4Pads=false);

  // effective number of photons for Ks
  static double NPho_Ks(const double e, const double sintheta, const bool all4Pads=false);

  // combined efficiency
  static double CP2BackEff(const double e, const double sintheta, 
			   const double mat, double& nbackpho, 
			   const TVector3 externpfrac, 
			   const bool all4Pads=false,
			   const int sys =0);

  static double conversionRate(const double npho, const double epho,
			       const double sth,
			       const double mat);

  // particle fraction in a jet 
  // pi0, eta, ks
  static void ParticleFraction_InclusiveEta(const double e, double* pfrac, const int sys=0);

  static double CP2UnderlyingEvent(int nvertex, const bool all4Pads=false);

  // the signal and background probabilities from the iso method
  static double isoPhoProb(const double pt);
  static double isoBkgProb(const double pt);

};

#endif