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

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

#endif

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

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

//gerolf 8.10.05
//int glfit(TH1F *h,Double_t norm=-1,Double_t peak=-1,Double_t lwidth=-1, Double_t gwidth=-1)
Double_t glfit(TH1F *h,Double_t norm=-1,Double_t peak=-1,Double_t lwidth=-1, Double_t gwidth=-1)
{
  TF1 *lfit1,*lfit2,*lfit3;

  if(norm < 0) {
    h->Fit("landau","QE");
    lfit1 = (TF1*)(h->GetFunction("landau"));
    par[0] =lfit1->GetParameter(0);
    //printf("norm %f par %f \n",norm, par[0]);
    //printf("1: %f %f %f %f\n", par[0], par[1], par[2], par[3]);
  }
  else
    par[0] = norm;
  if (peak == -1){
    //gerolf 9.10.05
    //par[1]=h->GetBinCenter(h->GetMaximumBin());
    h->Fit("landau","QE");
    lfit2 = (TF1*)(h->GetFunction("landau"));
    par[1] = lfit2->GetParameter(1);
    //printf("2: %f %f %f %f\n", par[0], par[1], par[2], par[3]);
    //gerolf
  }
  else
    par[1]=peak;
  if (lwidth < 0) {
    //gerolf 25.10.05
    //h->Fit("gaus","Q");
    //TF1 *gfit = (TF1 *)(h->GetFunction("gaus"));
    h->Fit("landau","QE");
    lfit3 = (TF1*)(h->GetFunction("landau"));
    par[2]= lfit3->GetParameter(2);
    //printf("3: %f %f %f %f\n", par[0], par[1], par[2], par[3]);
  } else 
    par[2]=lwidth;
  if (gwidth < 0) 
    par[3]=par[2];
  else
    par[3]=gwidth;

  //printf("Suggested: %f %f %f %f\n", par[0], par[1], par[2], par[3]);
  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
  //gerolf 8.10.05
  //return 0;

  
 

  return peakMOP;
}

 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,5);

  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],Norm);
  funfit->FixParameter(4,Norm);
  funfit->SetParNames("Constant","LandauMOP","LandauWidth","GaussSigma","Norm");
  printf("Starting the fit with parameter values: %f %f %f %f\n",par[0],par[1],par[2],par[3]);


  //gerolf
  h->Fit("glconv");
  //h->Fit("glconv","q");
  //end gerolf
  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]*par[4]*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
  //gerolf 8.10.05
  //Double_t peakMOP;
  Double_t width, peak1, width1, peak3, width3, errp, errw;
  char glMop[50],glFWHM[50];

  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(&peakMOP,&width,&sigleft, par);

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

  fpar[0] = par[0];
  fpar[1] = peakMOP; 
  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",peakMOP,errp);
  printf(" FWHM          = %11.5e +/- %9.3e \n",width,errw);
  printf(" Left sigma    = %11.5e Right sigma = %11.5e\n",sigleft,sigright);

  /*
  //add 2 custom entries to statistics display
  gPad->Update(); //important otherwise stats will be NULL pointer
  TPaveStats *stats =(TPaveStats*)h->GetListOfFunctions()->FindObject("stats");
  h->SetBit(TH1::kNoStats);
  sprintf(glMop," lgMOP  %11.3e +/- %9.3e  \n",peakMOP,errp);
  sprintf(glFWHM," FWHM  %11.3e +/- %9.3e  \n",width,errw);
  stats->AddText(glMop);
  stats->AddText(glFWHM); 
  stats->Draw();
  gPad->Update();
  //gPad->Modified();
  //h->DrawCopy();
  h->ResetBit(TH1::kNoStats);
*/


  //create new box for additional entries -- old solution
  gPad->Update(); //important otherwise stats will be NULL pointer
  TPaveStats *stats =(TPaveStats*)h->GetListOfFunctions()->FindObject("stats");
  TPaveStats *newstats = new TPaveStats;
  //stats->Copy((TObject&) newstats); //doesn't work while compiling
  newstats->SetX1NDC(stats->GetX1NDC());
  newstats->SetX2NDC(stats->GetX2NDC());
  newstats->SetY2NDC(stats->GetY1NDC());
  newstats->SetY1NDC(stats->GetY1NDC()-0.1);
  //  newstats->Draw("SAME");
  newstats->SetTextSize(0);
  sprintf(glMop,"  lgMOP  %11.3e +/- %9.3e  \n",peakMOP,errp);
  sprintf(glFWHM,"  FWHM  %11.3e +/- %9.3e  \n",width,errw);
  newstats->AddText(glMop);
  newstats->AddText(glFWHM); 
  newstats->Draw();
  

  gBenchmark->Show("glerror");



}

void PEAKVAR(Double_t* mean, Double_t* width, Double_t* leftsigma, Double_t* fpar) {
  Double_t delta, xmin, xmax;
  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;
 

}