#include <cstdio>

#include "Stntuple/photon/TCesCpr.hh"

 ClassImp(TCesCpr)

//_____________________________________________________________________________
TCesCpr::TCesCpr()
{
  Clear();
}

//_____________________________________________________________________________
TCesCpr::TCesCpr(const char* name)
{
  Clear();
  SetName(name);
}

//_____________________________________________________________________________
TCesCpr::~TCesCpr() {
}


//_____________________________________________________________________________
Float_t TCesCpr::AddPhoton(TStnEvent* event, TStnPhoton* pho) {

  SetEvent(event);
  for(int i=0; i < 8; i++)fCesCprSys[i]=fCprSys[i]=fCesSys[i]=0;

  fCesFlag = OkForCes(pho);
  fCprFlag = OkForCpr(pho);
  fCesWt = CesWeight(pho, fCesProbS, fCesProbB);
  fCprWt = CprWeight(pho, fCprProbS, fCprProbB);
  fCesHit = (fCesWt>0.0);
  fCprHit = (fCprWt<0.0);

  float fWeight = 0.0;
  if(fCesFlag==0) {
    fWeight = fCesWt;
  } else if (fCprFlag==0) {
    fWeight = fCprWt;
  }

  if(fabs(fWeight)>1e-6){
    fSum += fWeight;
    fSum2 += pow(fWeight,2);
    fSum2Bg += pow(1.0-fWeight,2);
    
    // sum up modified weights due to systematics
    if(fCesFlag==0)
      {
	// CES has systematics in element 0-2, empty from 3-8
	for(int is=0; is < 8; is++)
	  {
	    fSumSys[is]   += CesSys(is);
	    fCesCprSys[is] = CesSys(is);
	  }
      }
    else if(fCprFlag==0)
      {
	// CPR has systematics in element 0-3, and 7
	// CP2 has systematics in element 0-7
	for(int is=0; is < 8; is++)
	  {
	    fSumSys[is]   += CprSys(is); 	    
	    fCesCprSys[is] = CprSys(is);
	  }
      }

    fNUsed++;
  } else {
    fNCut++;
  }
  fNEvent++;
  return fWeight;
}

//_____________________________________________________________________________
Float_t TCesCpr::SignalFraction() { return fSum/NUsed();}
Float_t TCesCpr::SignalFractionErr(){ return sqrt(fSum2)/NUsed();}
Float_t TCesCpr::BGFraction(){ return 1.0-SignalFraction();}
Float_t TCesCpr::BGFractionErr(){ return sqrt(fSum2Bg)/NUsed();}
Float_t TCesCpr::Signal(Bool_t scale){
  return fSum*(scale?float(NEvent())/float(NUsed()):1.0);
}
Float_t TCesCpr::SignalErr(Bool_t scale){
  return sqrt(fSum2)*(scale?float(NEvent())/float(NUsed()):1.0);
}
Float_t TCesCpr::BG(Bool_t scale){
  return (fNUsed - fSum)*(scale?float(NEvent())/float(NUsed()):1.0);
}
Float_t TCesCpr::BGErr(Bool_t scale){
  return sqrt(fSum2Bg)*(scale?float(NEvent())/float(NUsed()):1.0);
}


//_____________________________________________________________________________
Float_t TCesCpr::SignalFractionSysErr()
{
  Float_t sysFrac[8];
  for(int is=0; is<8; is++)sysFrac[is]=fSumSys[is]/NUsed();
  
  Float_t ErrSum=0;
  for(int is=0; is<7; is++)
    {
      ErrSum += pow(sysFrac[is]-SignalFraction(),2);
    }

  ErrSum += pow(0.5*(sysFrac[7]-SignalFraction()),2);
  
  ErrSum = sqrt(ErrSum);

  return ErrSum;

}

Float_t TCesCpr::BGFractionSysErr()
{
  return SignalFractionSysErr();

}

Float_t TCesCpr::SignalSysErr(Bool_t scale)
{  
  Float_t ErrSum=0;
  for(int is=0; is<7; is++)
    {
      ErrSum += pow(fSumSys[is]-fSum,2);
    }

  ErrSum += pow(0.5*(fSumSys[7]-fSum),2);
  
  ErrSum = sqrt(ErrSum);

  return ErrSum*(scale?float(NEvent())/float(NUsed()):1.0);
}

Float_t TCesCpr::BGSysErr(Bool_t scale)
{
  return SignalSysErr(scale);
}

//_____________________________________________________________________________
void TCesCpr::Print() {
  printf("Results for Ces/CPR method, %s\n",GetName());
  TCesCprBase::Print();
  printf("Signal fraction = %6.3f +-%6.3f +- %6.3f \n",SignalFraction(),
	 SignalFractionErr(), SignalFractionSysErr());
  printf("BG fraction     = %6.3f +-%6.3f +- %6.3f \n",BGFraction(),
	 BGFractionErr(), BGFractionSysErr());
  printf("Result for useable photons:\n");
  printf("Signal     = %7.1f +-%7.1f +-%7.1f \n",Signal(false),
	 SignalErr(false), SignalSysErr(false));
  printf("Background = %7.1f +-%7.1f +-%7.1f \n",BG(false),
	 BGErr(false), BGSysErr(false));
  printf("Scaled to total numberof phtons:\n");
  printf("Signal     = %7.1f +-%7.1f +-%7.1f \n",Signal(),
	 SignalErr(), SignalSysErr());
  printf("Background = %7.1f +-%7.1f +-%7.1f \n",BG(),
	 BGErr(), BGSysErr());
}

//_____________________________________________________________________________
void TCesCpr::Clear() {
  TCesCprBase::Clear();
  fSum = 0.0;
  fSum2 = 0.0;
  fSum2Bg = 0.0;
  for(int i=0; i<8; i++)
    fSumSys[i]=0;
}