#include <iostream>
#include "TMinuit.h"

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TH1.h>
#include <TF1.h>
#include <TBenchmark.h>
#include <stdio.h>
#endif

static Double_t ConvX, Norm;
static TF1 *func;
static Double_t fpar[4], par[4]; //, errpar[4];

Double_t glconv(Double_t *x, Double_t *par);
Double_t gslan2(Double_t *y, Double_t *par);
void histgl(TH1F **h, const char *sh, Double_t width, Double_t sigma);
void glfitmake(TH1F *h);
void PEAKVAR(Double_t *mean, Double_t *width, Double_t *leftsigma, Double_t *fpar);
void glerror2(TH1F* h);

int glfit(TH1F *h,Double_t norm=5000,Double_t peak=-1,Double_t lwidth=-1, Double_t gwidth=-1)
{

  par[0]=norm;
  if (peak == -1)
    par[1]=h->GetBinCenter(h->GetMaximumBin());
  else
    par[1]=peak;
  if (lwidth < 0) {
    h->Fit("gaus","Q");
    TF1 *gfit = (TF1 *)(h->GetFunction("gaus"));
    par[2]=gfit->GetParameter(2);
  } else 
    par[2]=lwidth;
  if (gwidth < 0) 
    par[3]=par[2];
  else
    par[3]=gwidth;
  glfitmake(h); // makes the fit

  double covar[4][4];
  gMinuit->mnemat(&covar[0][0],4); // stores the covariant matrix

  double par[4],errpar[4];
  // Double_t errpar[4];
  
  for(int i=0;i<4;i++){
    gMinuit->GetParameter(i,par[i],errpar[i]); // get the fit parameters
  }
  glerror2(h);  // calculates MOP, FWHM and their errors
  return 0;
}

void histgl(TH1F **h, const char *sh, Double_t width, Double_t sigma)
{
  Double_t xx[1],par[4],weight;
  (*h) = new TH1F(sh,"G*L",100,-5,35);
  func = new TF1("func",gslan2,0,1,4);
  
  par[0]=70.4e+06; //1
  par[1]=0.19518e-02; //0
  par[2]=width; 
  par[3]=sigma;
  gBenchmark->Reset();
  gBenchmark->Start("gldraw");
  for(Int_t i=0; i<100; i++){
    xx[0]=(((Double_t) i)+0.5)*40.0/100.0-5.0;
    weight=glconv(xx,par);
    (*h)->Fill(xx[0],weight);
  }
  gBenchmark->Show("gldraw");
}

void glfitmake(TH1F *h)
{
  //  Double_t par[4] = {1000.,0.001,0.0004,0.0005};
  func = new TF1("gslan2",gslan2,0,1,4);
  TF1 *funfit = new TF1("glconv",glconv,0,1,4);

  Norm = (((Double_t) h->GetBinLowEdge(h->GetNbinsX()+1))-
	  ((Double_t) h->GetBinLowEdge(1)))/h->GetNbinsX();
  gBenchmark->Reset();
  gBenchmark->Start("glfit");
  funfit->SetParameters(par[0],par[1],par[2],par[3]);
  funfit->SetParNames("Constant","LandauMOP","LandauWidth","GaussSigma");
  printf("Starting the fit with parameter values: %f %f %f %f\n",par[0],par[1],par[2],par[3]);
  h->Fit("glconv");
  gBenchmark->Show("glfit");
}

Double_t glconv(Double_t *x, Double_t *par)
{
  Double_t LandauWidth=par[2]*0.58600;
  Double_t GaussWidth=par[3];
  ConvX=x[0];
  return((par[0]*Norm*func->Integral(-5.0*GaussWidth,5.0*GaussWidth,par,1.0e-8))/(LandauWidth*2.5066282746*GaussWidth)); 
}

Double_t gslan2(Double_t *y, Double_t *par)
{
  return((TMath::Landau(ConvX-y[0],par[1],par[2]*0.58600))
	 *exp(-(y[0]*y[0])/(2*par[3]*par[3])));
}

void glerror2(TH1F *h) {
  Double_t covar[4][4];
  Double_t par[4],errpar[4], dx, sigleft, sigright;
  // SLH 040324
  Double_t peak, width, peak1, width1, peak3, width3, errp, errw;

  gBenchmark->Reset();
  gBenchmark->Start("glerror");

  Norm = (((Double_t) h->GetBinLowEdge(h->GetNbinsX()+1))-
	  ((Double_t) h->GetBinLowEdge(1)))/h->GetNbinsX();

  /* First obtain the parameter values and covariance matrix */

  for(int i=0; i<4; i++) {
    gMinuit->GetParameter(i,par[i],errpar[i]);
  }
  gMinuit->mnemat(&covar[0][0],4);


  /* Fit the Peak position and FWHM of the obtained Landau*Gauss
     convolution: 
  */

  PEAKVAR(&peak,&width,&sigleft, par);

  sigright=width/(sqrt(2.*log(2.)))-sigleft;

  fpar[0] = par[0];
  fpar[1] = peak; 
  fpar[2] = par[2];
  fpar[3] = par[3];
  Double_t dummy, dPPd2, dFWd2, dPPd3, dFWd3;

  /* Computing the partial derivatives */
  /* -: d(PeakPosition)/d(fpar[2]) = dPPd2, d(FWHM)/d(fpar[2]) = dFWd2 */

  dx=sqrt(covar[2][2])/5.;  
  fpar[2]=par[2]-dx;
  PEAKVAR(&peak1,&width1,&dummy, fpar);
  fpar[2]=par[2]+dx;
  PEAKVAR(&peak3,&width3,&dummy, fpar);
  dPPd2=(peak3-peak1)/(2.*dx);
  dFWd2=(width3-width1)/(2.*dx);
  
  /* -: d(PeakPosition)/d(fpar[3]) = dPPd3, d(FWHM)/d(fpar[3]) = dFWd3 */

  dx=sqrt(covar[3][3])/5.;
  fpar[2]=par[2];
  fpar[3]=par[3]-dx;
  PEAKVAR(&peak1,&width1,&dummy, fpar);
  fpar[3]=par[3]+dx;
  PEAKVAR(&peak3,&width3,&dummy, fpar);
  dPPd3=(peak3-peak1)/(2.*dx);
  dFWd3=(width3-width1)/(2.*dx);

  /* Getting the right errors of peak and width: */
 
  errp = sqrt(covar[1][1] + covar[2][2]*dPPd2*dPPd2 +
	      covar[3][3]*dPPd3*dPPd3 + 2.*covar[1][2]*dPPd2 +
	      2.*covar[1][3]*dPPd3       + 2.*covar[2][3]*dPPd2*dPPd3);

  errw = sqrt(covar[2][2]*dFWd2*dFWd2 + covar[3][3]*dFWd3*dFWd3 +
  	      2.*covar[2][3]*dFWd2*dFWd3); 

  //  cout << endl;
  printf(" Peak position = %11.5e +/- %9.3e \n",peak,errp);
  printf(" FWHM          = %11.5e +/- %9.3e \n",width,errw);
  printf(" Left sigma    = %11.5e Right sigma = %11.5e\n",sigleft,sigright);
  gBenchmark->Show("glerror");
}

void PEAKVAR(Double_t* mean, Double_t* width, Double_t* leftsigma, Double_t* fpar) {
  Double_t delta, xmin, xmax, step, eps, r, dummy;
  Double_t f1, f2, fxmin, fxmax;
  Double_t yHalf, xHalfLow, xHalfHigh;

  int maxcall=1000;
  delta = sqrt(fpar[2]*fpar[2]+fpar[3]*fpar[3]);
  xmin  = fpar[1]-delta;
  xmax  = fpar[1]+3.*delta;

  Double_t xtop, ytop;
  par[0] = fpar[0];
  par[1] = fpar[1];
  par[2] = fpar[2];
  par[3] = fpar[3];
  ConvX = fpar[0];

  TF1 *funfit2 = new TF1("glconv",glconv,0,1,4);
  funfit2->SetParameters(par);

  ytop = funfit2->GetMaximum(xmin, xmax);
  xtop = funfit2->GetMaximumX(xmin, xmax);

  /* Find halfmaximum */

  yHalf = ytop / 2.;
  xHalfLow = funfit2->GetX(yHalf, xmin, par[1]);
  xHalfHigh = funfit2->GetX(yHalf, par[1], xmax);

  *mean = xtop;
  *leftsigma = (xtop - xHalfLow)/sqrt(2.*log(2.));
  *width = xHalfHigh - xHalfLow;
}