// File:  fitPar.cc
// Glen Cowan
// RHUL Physics

#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <cstdlib>
#include <cmath>
#include <TF1.h>
#include "Math/WrappedTF1.h"
#include "Math/BrentRootFinder.h"
#include "fitPar.h"

vector<double> n;         // needs to be global
vector<double> s;
vector<double> b;

double fitPar(vector<double> nArg, vector<double> sArg, vector<double> bArg){

  const double epsilon = 1.e-6;

  n = nArg;
  s = sArg;
  b = bArg;
  int Nbin = n.size();

  double muMin = -100.;
  for (int i=0; i<Nbin; i++){
    if ( b[i] > 0 ) {
      muMin = max(muMin, -s[i]/b[i]);
    }
  }
  muMin += epsilon;
  //  cout << "fitPar:  muMin = " << muMin << endl;

  double muMax = 1000.;
  int npar = 0;
  TF1 f("f", lnLderiv, muMin, muMax, npar);
  ROOT::Math::WrappedTF1 wf(f);
  ROOT::Math::BrentRootFinder brf;
  brf.SetFunction(wf, muMin, muMax);
  brf.Solve();
  double muHat = brf.Root();
  return muHat;

}

double lnLderiv(double* xArray, double* par){

  // derivative of lnL with respect to mu

  double mu = xArray[0];
  int Nbin = n.size();
  double f = 0.;
  for (int i=0; i<Nbin; i++){
    f += n[i]*s[i] / (mu*s[i] + b[i])  - s[i];
  }

  return f;

}