#include <iostream>
#include <iomanip>
using namespace std;
#include "TH1D.h"
#include "TF1.h"
#include <TSystem.h>
#include <TCanvas.h>
#include <TApplication.h>
#include <TF1.h>
#include <TH1.h>
#include <TAxis.h>
#include <TLegend.h>
#include <TLine.h>
#include <TStyle.h>
#include <TPaveText.h>
#include "Math/IntegratorOptions.h"
#include "Math/MinimizerOptions.h"
#include "TMinuit.h"
#include "TRandom.h"
#include "fitTools.h"
#include "my_fitTools.h"
#include "fw_fitTools.h"

//TH1D *hd ;
//int mytest(){
int main(){
  
  hd = new TH1D("hd","hd", nbin, bins);
  hd->SetBinContent(0  , 0          ); 
  hd->SetBinContent(1  , 4.53902    ); 
  hd->SetBinContent(2  , 4.70824    ); 
  hd->SetBinContent(3  , 4.73104    ); 
  hd->SetBinContent(4  , 4.8716     ); 
  hd->SetBinContent(5  , 4.94311    ); 
  hd->SetBinContent(6  , 4.92585    ); 
  hd->SetBinContent(7  , 5.0351     ); 
  hd->SetBinContent(8  , 5.25363    ); 
  hd->SetBinContent(9  , 5.30603    ); 
  hd->SetBinContent(10 ,  5.44443   ); 
  hd->SetBinContent(11 ,  5.53738   ); 
  hd->SetBinContent(12 ,  5.54719   ); 
  hd->SetBinContent(13 ,  5.76489   ); 
  hd->SetBinContent(14 ,  5.94505   ); 
  hd->SetBinContent(15 ,  6.15808   ); 
  hd->SetBinContent(16 ,  6.43717   ); 
  hd->SetBinContent(17 ,  6.73969   ); 
  hd->SetBinContent(18 ,  6.9798    ); 
  hd->SetBinContent(19 ,  7.25404   ); 
  hd->SetBinContent(20 ,  7.69753   ); 
  hd->SetBinContent(21 ,  8.0544    ); 
  hd->SetBinContent(22 ,  8.73745   ); 
  hd->SetBinContent(23 ,  9.82025   ); 
  hd->SetBinContent(24 ,  11.7717   ); 
  hd->SetBinContent(25 ,  13.485    ); 
  hd->SetBinContent(26 ,  7.77343   ); 
  hd->SetBinContent(27 ,  1.36615   ); 
  hd->SetBinContent(28 ,  1.25825   ); 
  hd->SetBinContent(29 ,  0.556537  );  
  hd->SetBinContent(30 ,  0.316932  );  
  hd->SetBinContent(31 ,  0.348092  );  
  hd->SetBinContent(32 ,  0.472119  );  
  hd->SetBinContent(33 ,  0  	    ); 

  hd->SetBinError(0,   0            );
  hd->SetBinError(1,   0.41905      );
  hd->SetBinError(2,   0.480005     );
  hd->SetBinError(3,   0.256493     );
  hd->SetBinError(4,   0.334958     );
  hd->SetBinError(5,   0.368144     );
  hd->SetBinError(6,   0.214421     );
  hd->SetBinError(7,   0.272643     );
  hd->SetBinError(8,   0.353102     );
  hd->SetBinError(9,   0.34398      );
  hd->SetBinError(10,   0.345476    );
  hd->SetBinError(11,   0.289041    );
  hd->SetBinError(12,   0.26208     );
  hd->SetBinError(13,   0.212691    );
  hd->SetBinError(14,   0.203182    );
  hd->SetBinError(15,   0.310943    );
  hd->SetBinError(16,   0.360032    );
  hd->SetBinError(17,   0.489854    );
  hd->SetBinError(18,   0.461181    );
  hd->SetBinError(19,   0.639442    );
  hd->SetBinError(20,   0.510431    );
  hd->SetBinError(21,   0.426875    );
  hd->SetBinError(22,   0.232739    );
  hd->SetBinError(23,   0.310051    );
  hd->SetBinError(24,   0.274143    );
  hd->SetBinError(25,   0.569625    );
  hd->SetBinError(26,   0.410178    );
  hd->SetBinError(27,   0.215306    );
  hd->SetBinError(28,   0.229667    );
  hd->SetBinError(29,   0.107403    );
  hd->SetBinError(30,   0.0644741   );
  hd->SetBinError(31,   0.0677504   );
  hd->SetBinError(32,   0.0930139   );
  hd->SetBinError(33,   0  	    );

  QuenchedBragg qb; //TF1 Functors
  //  double R0 =hd->GetBinCenter(hd->GetMaximumBin())+10; //mod aa.... I see from back is important.  not  good for these dependence


  //  double R0 = EstimateR0(hd);
  double R0 =hd->GetBinCenter(hd->GetMaximumBin())+10; //mod aa.... I see from back is important.  not  good for these dependence
  double E0 = pow(R0/alpha,1/p)/mp;
  //  double E0 = 0;
  double xOffset = 0.0;
  double sigma_mono = 0.012*pow(R0,0.935);
  double sigma_E0 = 0.002*E0;                                         //in MeV. realistic initial beam energy spread?
  double sigma_0 = pow(sigma_E0*alpha*p,2)*pow(E0,2*p-2);             //in mm, estimated offset due to detector effects and initial beam energy spread
  double sigma = sqrt(pow(sigma_mono,2) + pow(sigma_0,2));
  double phi0 = 1e-7*hd->GetBinContent(hd->GetMaximumBin());  //rough guess
  double fit_start = 0;
  if(R0 < 30) fit_start = (R0-xOffset-7)/1.8;                         // protons hit sensor directly at low beam energies/ranges
  double fit_end = hd->GetXaxis()->GetBinLowEdge(hd->GetMaximumBin()+3); //light back to zero normally 3 sheets after peak
  double kb = 0.07;

  cout<<"init para = "<<R0<<" "<<sigma<<" "<<phi0<<" "<<xOffset<<" "<<kb<<endl;
  TF1* fphotons = new TF1("Fitted Light Output",&qb,0, 78.6606, 5);
  fphotons->SetParameters(R0,sigma,phi0,xOffset,kb);
  fphotons->SetParNames("R_{0}","#sigma_{R}","#Phi_{0}","xOffset","kB");
  fphotons->FixParameter(3,0.0);
  fphotons->SetNpx(100);
  fphotons->SetRange(0, 92.4508);
  E0 = pow(fphotons->GetParameter(0)/alpha,1/p)/mp;
  cout<<"Initial value -----------------------------------"<<endl;
  cout<<"PreCalc R0    = "<<fphotons->GetParameter(0)<<endl;
  cout<<"PreCalc sigma = "<<fphotons->GetParameter(1)<<endl;
  cout<<"PreCalc phi   = "<<fphotons->GetParameter(2)<<endl;
  //  cout<<"PreCalc off   = "<<fphotons->GetParameter(3)<<endl;  
  cout<<"PreCalc kB    = "<<fphotons->GetParameter(4)<<endl;  
  cout<<"PreCalc E0    = "<<E0<<endl;

  hd->Fit("Fitted Light Output","IRQ0"); //updates roughly 5x faster without I option		// out to reduce the test time
  E0 = pow(fphotons->GetParameter(0)/alpha,1/p)/mp;
  cout<<"original fit results -----------------------------------"<<endl;
  cout<<"Orig Fit R0    = "<<fphotons->GetParameter(0)<<endl;
  cout<<"Orig Fit sigma = "<<fphotons->GetParameter(1)<<endl;
  cout<<"Orig Fit phi   = "<<fphotons->GetParameter(2)<<endl;
  //  cout<<"Orig Fit off   = "<<fphotons->GetParameter(3)<<endl;  
  cout<<"Orig Fit kB    = "<<fphotons->GetParameter(4)<<endl;  
  cout<<"Orig  E0       = "<<E0<<endl;

  //HH Full Minuit HH//
  TF1* fphotons2 = new TF1("Fitted Light Output2",&qb,0, 78.6606, 5);
   fphotons2->SetParNames("R_{0}","#sigma_{R}","#Phi_{0}","xOffset","kB");
   fphotons2->FixParameter(3,0.0);
   fphotons2->SetNpx(100);
   fphotons2->SetRange(0, 92.4508);

  const int npar=4;
  int ier=0;
  TMinuit *gMinuit = new TMinuit(npar);
  gMinuit->SetPrintLevel(-1); // -1:no , 0:summary, 1:many
  gMinuit->SetFCN(QuenchedBraggFitFunc);     
  //H set parameter   ( start, step, min, max)H//
   gMinuit->mnparm(0, "R0"   ,    71,  0.01,     0,   100, ier);
   gMinuit->mnparm(1, "sigma",     2, 0.001,   0.0,    10, ier);
   gMinuit->mnparm(2, "phi0" ,10e-07,  1e-9, 1e-07, 20e-07, ier);
   gMinuit->mnparm(3, "birks",  0.07, 0.001,  0.01,   0.5, ier);
  //  gMinuit->FixParameter(1);

  Double_t arg[] = {2000,1.0};
  gMinuit->mnexcm("MIGRAD", arg, npar, ier);

  TString name;
  double value;
  double error;
  double other;
  double eR0, esigma, ephi0, ekb, initR0;
  gMinuit->mnpout(0, name,    R0,    eR0, other, other,  ier);
  gMinuit->mnpout(1, name, sigma, esigma, other, other,  ier);
  gMinuit->mnpout(2, name,  phi0,  ephi0, other, other,  ier);
  gMinuit->mnpout(3, name,    kb,    ekb, other, other,  ier);

  E0 = pow(R0/alpha,1/p)/mp;
  cout<<"my fit results -----------------------------------"<<endl;
  cout<<"MyF_Min R0 = "<<R0<<" +/- "<<eR0<<endl;
  cout<<"MyF_Min sigma = "<<sigma<<" +/- "<<esigma<<endl;
  cout<<"MyF_Min phi0 = "<<phi0<<" +/- "<<ephi0<<endl;
  cout<<"MyF_Min kb = "<<kb<<" +/- "<<ekb<<endl;
  cout<<"MyF_Min E0    = "<<E0<<endl;

  fphotons2->SetParameters(R0,sigma,phi0,xOffset,kb);


  //HHHHH fw test HHHHHH//
  TF1* fphotons3 = new TF1("Fitted Light Output3",&qb,0, 78.6606, 5);
   fphotons3->SetParNames("R_{0}","#sigma_{R}","#Phi_{0}","xOffset","kB");
   fphotons3->FixParameter(3,0.0);
   fphotons3->SetNpx(100);
   fphotons3->SetRange(0, 92.4508);

  double ob[nbin], obe[nbin];
  for(int kk=1; kk<nbin+1; kk++){
        ob[kk-1]=hd->GetBinContent(kk);
        obe[kk-1]=hd->GetBinError(kk);
  }
  initR0=hd->GetBinCenter(hd->GetMaximumBin())+10;

  cout<<"fw fit before -----------------------------------"<<endl;
  cout<<"MyF_Min R0 = "<<R0<<endl;
  cout<<"MyF_Min sigma = "<<sigma<<endl;
  cout<<"MyF_Min phi0 = "<<phi0<<endl;
  cout<<"MyF_Min kb = "<<kb<<endl;
  cout<<"MyF_Min E0    = "<<E0<<endl;
  fw_binned(ob, obe, R0, sigma, phi0, kb, initR0);
  phi0=phi0*1e-7;
  cout<<"fw fit results -----------------------------------"<<endl;
  cout<<"MyF_Min R0 = "<<R0<<endl;
  cout<<"MyF_Min sigma = "<<sigma<<endl;
  cout<<"MyF_Min phi0 = "<<phi0<<endl;
  cout<<"MyF_Min kb = "<<kb<<endl;
  cout<<"MyF_Min E0    = "<<E0<<endl;

  fphotons3->SetParameters(R0,sigma,phi0,xOffset,kb);
  cout<<" Maybe origianl best fit ------------------------------ exmaple"<<endl;
  cout<<"  //fit results"<<endl;
  cout<<"  //R0 = 73.8947"<<endl;
  cout<<"  //sigma = 2.00352"<<endl;
  cout<<"  //phi = 6.12822e-07"<<endl;
  cout<<"  //off = 0"<<endl;
  cout<<"  //kB = 0.0751327"<<endl;
  cout<<"  //E0 = 96.5939"<<endl;


  //HHH  Draw  HHH// ====================================================
  gStyle->SetOptStat("000");
  TCanvas *c1 = new TCanvas("c1","c1",600, 600);
  c1->SetGrid();
  c1->SetGridx();
  c1->SetGridy();  
  hd->SetMaximum(hd->GetMaximum()*1.8);
  hd->SetTitle("");
  hd->SetXTitle("WET (mm)");
  hd->SetYTitle("Charge (pC)");  
  hd->Draw();

  fphotons->SetLineColor(1);
  fphotons->Draw("same");  
  fphotons2->SetLineColor(2);
  fphotons2->Draw("same");  
  fphotons3->SetLineColor(kGreen);
  fphotons3->Draw("same");  
  TLegend* leg = new TLegend(0.10, 0.73, 0.7, 0.88);
  leg->SetBorderSize(0);
  leg->SetFillColor(0);
  leg->SetFillStyle(0);
  leg->AddEntry(fphotons, "origFit");
  leg->AddEntry(fphotons2, "FW+Minui");
  leg->AddEntry(fphotons3, "Full FW");
  leg->Draw();

  c1->Update();
  c1->SaveAs("test.pdf");
  
  return 0;
}