//////////////////////////////////////////////////////////////////////////
//                                                                      //
// GFMath                                                               //
//                                                                      //
// Encapsulate Gero's math routines.                                    //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include "TError.h"
#include "TMath.h"
#include "GFMath.h"
#include "TArrayD.h"
#include "TMatrixD.h"
// #include <iostream>

extern "C" { // FORTRAN function
  Float_t binoml_(Int_t& nRec, Int_t& nGen, Float_t& confLev);
}


ClassImp(GFMath)

//__________________________________________________________________
 TArrayD GFMath::X1ThroughX2CorrErr(Double_t x1, Double_t x2, 
				  Double_t sigmax1, Double_t sigmax2)
{
  // returns an array of length=2 with result[0] = x1/x2
  //                                   result[1] = sigma(x1/x2)
  // assuming that x1 and x2 are 100% correlated

  // cf. logbook 18.12.2002
  TArrayD result(2);
  if(x2 == 0.) {
    ::Error("GFMath::X1ThroughX2CorrErr", "x2 = 0.0, return {0.,0.}");
    return result;
  }
  result[0] = x1/x2;
  result[1] = TMath::Abs((sigmax1 - result[0]*sigmax2) / x2);
  return result;
}

//__________________________________________________________________
 TArrayD GFMath::X1ThroughX2BinomErr(Double_t x1, Double_t x2) 
{
  // returns an array of length=2 with result[0] = x1/x2
  //                                   result[1] = sigma(x1/x2)
  // assuming that x1 is taken with binomial distribution from a sample of size x2

  // cf. logbook 18.12.2002
  TArrayD result(2);
  if(x2 == 0.) {
    ::Error("GFMath::X1ThroughX2BinomErr", "x2 = 0.0, return {0.,0.}");
    return result;
  }
  result[0] = x1/x2;
  result[1] = TMath::Sqrt((result[0] * (1. - result[0])) / x2);
  return result;
}

//__________________________________________________________________
 TArrayD GFMath::NormaliseTo1(const TArrayD& input, const TArrayD* inputErr, 
			     TArrayD* outputErr, TMatrixD* outputVariance)
{
  // return input normalised to its sum (no normalisation if sum == 0.)
  // if inputErr is given, calculate outputErr and/or outputVariance (if given)
  // assuming inputErr to be non-correlated
  // no error calculation if inputErr has wrong length or sum == 0

  const Int_t size = input.GetSize();
  TArrayD result(size);
  const Double_t sum = input.GetSum();
  for(Int_t i = 0; i < size; ++i){
    result[i] = input[i] / (sum ? sum : 1.);
  }

  if(outputErr) {// make valid structure if given - even if not to be filled
    outputErr->Set(size);
    outputErr->Reset();
  }
  if(outputVariance) { // dito
    outputVariance->ResizeTo(size, size);
    outputVariance->Zero();
  }

  if(sum && inputErr && size==inputErr->GetSize() && (outputErr||outputVariance)){
    // The error propagation is done according to:
    // V. Blobel/E. Lohrmann: Statistische und numerische Methoden der Datenanlyse, p. 148
    // y_i = x_i/(sum_j(x_j)), B_ij = dy_i/dx_j
    // cf. log book, May 4th, 2004

    TArrayD inputVar(size); 
    for(Int_t i = 0; i < size; ++i){
      inputVar[i] = inputErr->At(i) * inputErr->At(i);
    }
    const Double_t sumInputVar = inputVar.GetSum();
    const Double_t sumTo2 = sum * sum;
    TArrayD outDiagSigma(size);
    for(Int_t i = 0; i < size; ++i){
      const Double_t first = (sum-input[i]) * (sum-input[i]) * inputVar[i];
      const Double_t second= input[i] * input[i] * (sumInputVar-inputVar[i]);
      outDiagSigma[i] = TMath::Sqrt(first + second)/sumTo2;
    }
    if(outputErr){
      (*outputErr) = outDiagSigma;
    }
    if(outputVariance){
      for(Int_t i = 0; i < size; ++i){
	(*outputVariance)(i,i) = outDiagSigma[i] * outDiagSigma[i];
      }
      for(Int_t i = 0; i < size; ++i){
	for(Int_t j = i+1; j < size; ++j){
	  const Double_t first = input[i] * input[j] * (sumInputVar-inputVar[i]-inputVar[j]);
	  const Double_t second= input[i] * inputVar[j] * (sum - input[j]);
	  const Double_t third = input[j] * inputVar[i] * (sum - input[i]);
	  (*outputVariance)(i,j) = (*outputVariance)(j,i) = 
	    (first - (second + third))/ sumTo2 / sumTo2;
	}
      }
    }// calculate variance
  }// calculate errors

  return result;
}

//__________________________________________________________________
 Float_t GFMath::BinomialError(Int_t nRec, Int_t nGen, Bool_t upper)
{
  // upper or lower 'error' of nRec/nGen from binomial according to V. Blobel

// *************************************************************
// *
// *           NND=...                           ! n  
// *           NNN=...                           ! k
// *           RATIO=NNN/NND                     ! p = k/n
// *           RTLEFT= BINOML(NNN,NND,0.16)
// *           IF(RTLEFT.GE.RATIO) RTLEFT=0.0
// *           RTRIGT= BINOML(NNN,NND,0.84)  
// *           IF(RTRIGT.LE.RATIO) RTRIGT=1.0
// *           RL=RATIO-RTLEFT
// *           IF(RL.LT.0.0) RL=RATIO            ! left "error" 
// *           RR=RTRIGT-RATIO
// *           IF(RR.LT.0.0) RR=1.0-RATIO        ! right "error"
// *
// *************************************************************

  const Float_t ratio = nRec/static_cast<Float_t>(nGen);
  Float_t result = 0.;

  if(upper){
    Float_t limit = 0.84;
    Float_t rtright = binoml_(nRec, nGen, limit);
    if(rtright <= ratio) rtright = 1.;
    result = rtright - ratio;
    if(result < 0.0) result = 1.0 - ratio;
  } else {
    Float_t limit = 0.16;
    Float_t rtleft = binoml_(nRec, nGen,  limit);
    if(rtleft >= ratio) rtleft = 0.;
    result = ratio - rtleft;
    if(result < 0.0) result = ratio;
  }

  return result;
}

//____________________________________
 TArrayD GFMath::Efficiency(Double_t nCount, Double_t nDenom, Double_t nEffDenom, 
			   Double_t nFillDenom)
{
  // returning 'efficiency = nCount/nDenom' and its error(s): 
  // nCount: surviving, nDenom: reference, nEffDenom: effective number entries for nDenom
  // nFillDenom: hist fill actions for nDenom
  // if(nEffDenom/nDenom >= 0.7): do asymm. binomial error GFMath::BinomialError
  // then result[2/3] are the asymmetric errors
  // result[1]: symmetric easy binomial error, but if asymmetric errors are possible, 
  // 'scaled to nEffDenom'
  
  const Double_t minNeffRatio = 0.7; // cf. GFHistManip::CreateRatioBinomErrors

  TArrayD eff(4);
  eff[0] = nDenom ? nCount/nDenom : 0.;

  if(!nDenom || nEffDenom/nFillDenom < minNeffRatio){
    ::Warning("GFMath::Efficiency", "N_eff^rec %.4g, N^rec %.4g", nEffDenom, nFillDenom);
    eff.Set(2);
    eff[1] = nDenom ? TMath::Sqrt(eff[0] * (1.- eff[0])/nDenom) : 0.;
  } else {
    eff[1] = nDenom ? TMath::Sqrt(eff[0] * (1.- eff[0])/nEffDenom) : 0.;
    const Int_t scaledNcount = TMath::Nint(nCount*nEffDenom/nDenom);//scale to effective number
    eff[2]=  GFMath::BinomialError(scaledNcount, TMath::Nint(nEffDenom), kTRUE);
    eff[3]= -GFMath::BinomialError(scaledNcount, TMath::Nint(nEffDenom), kFALSE);//lower error
  }

  return eff;
}

//__________________________________________________________________
 TArrayD GFMath::EfficiencyGassnerError(Double_t nCount, Double_t nDenom, 
					Double_t errCount, Double_t errDenom)
{
  // efficiency nCount/nDenom with 'Gassner' errors
  TArrayD effPlusErr(2);
  effPlusErr[0] = (nDenom ? nCount/nDenom : 0.);

  Double_t relErrEff2 = 0.;
  if(nCount && nDenom) {
    relErrEff2 = errDenom * errDenom / nDenom / nDenom 
      + (1. - 2.*effPlusErr[0]) * errCount * errCount / nCount / nCount;
  }
  if(relErrEff2 >= 0.) {
    effPlusErr[1] = TMath::Sqrt(relErrEff2) * effPlusErr[0];
  } else {
    ::Error("GFMath::EfficiencyGassnerError", "square of relative error is %.2g", relErrEff2);
  }
  
  return effPlusErr;
}

//__________________________________________________________________
 TArrayD GFMath::EfficiencyBlobelError(Double_t nCount, Double_t nDenom,
				      Double_t errCount, Double_t errDenom)
{
  // according to Blobel, 2004-09-08
  TArrayD effPlusErr(2);
  effPlusErr[0] = (nDenom ? nCount/nDenom : 0.);

  // error from correlation, approximating rho by 
  // common events/all events ~ errCount^2/errDenom^2
  const Double_t errRelDenom = (nDenom ? errDenom / nDenom : 0.);
  const Double_t errRelCount = (nCount ? errCount / nCount : 0.);
  const Double_t rho = (errDenom ? errCount * errCount / (errDenom * errDenom) : 0.);
  const Double_t relErrEff2 = 
    errRelDenom * errRelDenom + errRelCount * errRelCount - 2.*rho * errRelDenom * errRelCount;

  if(relErrEff2 >= 0.) {
    effPlusErr[1] = TMath::Sqrt(relErrEff2) * effPlusErr[0];
  } else {
    ::Error("GFMath::EfficiencyBlobelError", "square of relative error is %.2g", relErrEff2);
  }
  // now unsymmetric error from correct binomial, symmetrise by averaging squares
  TArrayD unsymErrs(4);
  // if no denominator, the next would put an error... But unsymm. Binomial errors are senseless
  if(nDenom) unsymErrs = GFMath::Efficiency(nCount, nDenom, nDenom, nDenom);
  const Double_t unsymmErr2 = (unsymErrs[2]*unsymErrs[2] + unsymErrs[3]*unsymErrs[3])/2.;
  // add them quadratically:
  effPlusErr[1] = TMath::Sqrt(effPlusErr[1]*effPlusErr[1] + unsymmErr2);

  return effPlusErr;
}

//__________________________________________________________________
 TArrayD GFMath::EfficiencyTestError(Double_t nCount, Double_t nDenom,
				      Double_t errCount, Double_t errDenom)
{
  // similar ansatz as Blobel, but rho = nCount/nDenom
  // EXPERIMENTAL!!!
  TArrayD effPlusErr(2);
  effPlusErr[0] = (nDenom ? nCount/nDenom : 0.);

  // error from correlation, approximating rho by 
  // common events/all events ~ errCount^2/errDenom^2
  const Double_t errRelDenom = (nDenom ? errDenom / nDenom : 0.);
  const Double_t errRelCount = (nCount ? errCount / nCount : 0.);
  const Double_t rho = (nDenom ? nCount/nDenom : 0.);
  const Double_t relErrEff2 = 
    errRelDenom * errRelDenom + errRelCount * errRelCount - 2.*rho * errRelDenom * errRelCount;

  if(relErrEff2 >= 0.) {
    effPlusErr[1] = TMath::Sqrt(relErrEff2) * effPlusErr[0];
  } else {
    ::Error("GFMath::EfficiencyTestError", "square of relative error is %.2g", relErrEff2);
  }
  // now unsymmetric error from correct binomial, symmetrise by averaging squares
  TArrayD unsymErrs(4);
  // if no denominator, the next would put an error... But unsymm. Binomial errors are senseless
  if(nDenom) unsymErrs = GFMath::Efficiency(nCount, nDenom, nDenom, nDenom);
  const Double_t unsymmErr2 = (unsymErrs[2]*unsymErrs[2] + unsymErrs[3]*unsymErrs[3])/2.;
  // add them quadratically:
  effPlusErr[1] = TMath::Sqrt(effPlusErr[1]*effPlusErr[1] + unsymmErr2);

  return effPlusErr;
}

//__________________________________________________________________
 Bool_t GFMath::EquidistLogBins(Double_t* bins, Int_t nBins,Double_t first,Double_t last)
{
  // Filling 'bins' with borders of 'nBins' bins between 'first' and 'last' that are
  // equidistant when viewed in log scale, so 'bins' must have length nBins+1;
  // If 'first', 'last' or 'nBins' are not positive, failure is reported.
  // 
  if(nBins < 1 || first <= 0. || last <= 0.) return kFALSE;
  
  bins[0] = first;
  bins[nBins] = last;
  const Double_t firstLog = TMath::Log10(bins[0]);
  const Double_t lastLog = TMath::Log10(bins[nBins]);
  for (Int_t i = 1; i < nBins; ++i){
    bins[i] = TMath::Power(10., firstLog + i*(lastLog-firstLog)/(nBins));
  }
  
  return kTRUE;
}

//__________________________________________________________________
 Bool_t GFMath::EquidistBins(Double_t* bins, Int_t nBins,Double_t first,Double_t last)
{
  // Filling 'bins' with borders of 'nBins' bins between 'first' and 'last' that are
  // equidistant, so 'bins' must have length nBins+1;
  // If 'first' >= 'last' or 'nBins' < 1, failure is reported.
  // 
  if(nBins < 1 || first >= last) return kFALSE;
  
  bins[0] = first;
  bins[nBins] = last;
  for (Int_t i = 1; i < nBins; ++i){
    bins[i] = first + i*(last-first)/(nBins);
  }
  
  return kTRUE;
}

//__________________________________________________________________
 Bool_t GFMath::IsEqualFloat(Double_t x1, Double_t x2)
{
  // true if x1 , x2 are equal besides floating point surprises:
  // difference < 1.e-6 and 
  // (if x1 != 0 && x2 != 0) relative deviation < 1.e-4 
  // (for the latter 0 means |x| < 1.e-6

  const Double_t posZeroDiff = 1.e-6;
  const Double_t posZeroQuot = 1.e-4;
  
  if(TMath::Abs(x1 - x2) > posZeroDiff
     || (TMath::Abs(x2) > posZeroDiff && TMath::Abs((x1/x2) -1.) > posZeroQuot)
     || (TMath::Abs(x1) > posZeroDiff && TMath::Abs((x2/x1) -1.) > posZeroQuot)
     ){
//     std::cout << "IsEqualFloat " << " x1 x2: " << x1 << "  " << x2;
//     if(x2) std::cout	<< ", (x1/x2)-1 = " << (x1/x2) -1.;
//     if(x1) std::cout << ", (x2/x1)-1 = " << (x2/x1) -1.;
//     std::cout << ", (x2-x1) = " << x2-x1  
// 	 << std::endl;
    return kFALSE;
  }
  return kTRUE;
}