#include <iostream>
#include <iomanip>
using namespace std;
#include "TH1D.h"
#include "TF1.h"

#include "fitTools.h"


const int nbin=32;
int mytest(){
  
  double bins[nbin+1]={ 0,  2.90158,  5.7519,  8.78677,  11.6986,  14.5387,  17.5428,  20.4444,  23.3152,  26.3193,  29.2209,  32.1122,  35.1061,  38.0179,  40.899,  43.8723,  46.7739,  49.6755,  52.8437,
			55.6837,  58.5033,  61.4561,  64.3885,  67.2798,  70.2224,  73.1445,  75.9538,  78.6606,  81.3879,  84.4125,  87.109,  89.785,  92.4508};
  TH1D *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 = EstimateR0(hd);
  double E0 = pow(R0/alpha,1/p)/mp;
  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;

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


  hd->Fit("Fitted Light Output","IRQ0"); //updates roughly 5x faster without I option		
  fphotons->SetRange(0, 92.4508);
  R0 = fphotons->GetParameter(0);      
  sigma = fphotons->GetParameter(1);
  E0 = pow(R0/alpha,1/p)/mp;

  cout<<"R0 = "<<R0<<endl;
  cout<<"sigma = "<<sigma<<endl;
  cout<<"E0 = "<<E0<<endl;
  cout<<"phi = "<<fphotons->GetParameter(2)<<endl;
  cout<<"?? = "<<fphotons->GetParameter(3)<<endl;  
  cout<<"kB = "<<fphotons->GetParameter(4)<<endl;  
  
  //  hd->Draw("");

  return 0;
}