#include <iostream>
#include <fstream>
#include <stdio.h>
#include <TFile.h>
#include <TROOT.h>
#include <TMath.h>
#include <TF1.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TGraphAsymmErrors.h>
#include <TLegend.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TColor.h>
#include <TLegend.h>
#include <TLegendEntry.h>

Double_t brkpwlfunction(Double_t* x, Double_t* par);

void plotandfit(double crabunits, double N_0_1TeV, double E_c, double index_alpha1, double index_alpha2, double SS){
   
   TGraphAsymmErrors *grae = new TGraphAsymmErrors(29);
   grae->SetName("Graph0");
   grae->SetTitle("BRKPWL Spectrum");
   grae->SetFillColor(1);
   grae->SetMarkerColor(4);
   grae->SetMarkerStyle(8);
   grae->SetPoint(0,0.1096478,7.625757e-09);
   grae->SetPointError(0,0,0,3.049247e-11,3.049247e-11);
   grae->SetPoint(1,0.1380384,3.818161e-09);
   grae->SetPointError(1,0,0,1.849496e-11,1.849496e-11);
   grae->SetPoint(2,0.1757924,1.853025e-09);
   grae->SetPointError(2,0,0,9.20288e-12,9.20288e-12);
   grae->SetPoint(3,0.2213095,9.256503e-10);
   grae->SetPointError(3,0,0,6.31112e-12,6.31112e-12);
   grae->SetPoint(4,0.2754229,4.812964e-10);
   grae->SetPointError(4,0,0,3.263162e-12,3.263162e-12);
   grae->SetPoint(5,0.3507519,2.333708e-10);
   grae->SetPointError(5,0,0,1.873402e-12,1.873402e-12);
   grae->SetPoint(6,0.4415704,1.171875e-10);
   grae->SetPointError(6,0,0,1.13263e-12,1.13263e-12);
   grae->SetPoint(7,0.5559043,5.903162e-11);
   grae->SetPointError(7,0,0,6.045369e-13,6.045369e-13);
   grae->SetPoint(8,0.699842,2.913175e-11);
   grae->SetPointError(8,0,0,3.785573e-13,3.785573e-13);
   grae->SetPoint(9,0.8709636,1.517041e-11);
   grae->SetPointError(9,0,0,1.987272e-13,1.987272e-13);
   grae->SetPoint(10,1.109175,8.445677e-12);
   grae->SetPointError(10,0,0,1.151211e-13,1.151211e-13);
   grae->SetPoint(11,1.396368,5.992819e-12);
   grae->SetPointError(11,0,0,8.371367e-14,8.371367e-14);
   grae->SetPoint(12,1.737801,4.302218e-12);
   grae->SetPointError(12,0,0,5.519315e-14,5.519315e-14);
   grae->SetPoint(13,2.187762,3.042564e-12);
   grae->SetPointError(13,0,0,3.934229e-14,3.934229e-14);
   grae->SetPoint(14,2.786121,2.156292e-12);
   grae->SetPointError(14,0,0,2.512435e-14,2.512462e-14);
   grae->SetPoint(15,3.548134,1.4709e-12);
   grae->SetPointError(15,0,0,1.747795e-14,1.747811e-14);
   grae->SetPoint(16,4.466836,1.044279e-12);
   grae->SetPointError(16,0,0,1.237983e-14,1.237994e-14);
   grae->SetPoint(17,5.559043,7.707763e-13);
   grae->SetPointError(17,0,0,8.917575e-15,8.917651e-15);
   grae->SetPoint(18,6.91831,5.444112e-13);
   grae->SetPointError(18,0,0,6.062264e-15,6.062303e-15);
   grae->SetPoint(19,8.810489,3.835841e-13);
   grae->SetPointError(19,0,0,3.980762e-15,3.980777e-15);
   grae->SetPoint(20,11.22018,2.662917e-13);
   grae->SetPointError(20,0,0,2.985446e-15,2.985459e-15);
   grae->SetPoint(21,14.12538,1.885494e-13);
   grae->SetPointError(21,0,0,2.148793e-15,2.1488e-15);
   grae->SetPoint(22,17.78279,1.322823e-13);
   grae->SetPointError(22,0,0,1.564582e-15,1.564586e-15);
   grae->SetPoint(23,22.38721,9.408806e-14);
   grae->SetPointError(23,0,0,1.148963e-15,1.148968e-15);
   grae->SetPoint(24,28.18383,6.682623e-14);
   grae->SetPointError(24,0,0,8.507033e-16,8.507069e-16);
   grae->SetPoint(25,35.48134,4.658165e-14);
   grae->SetPointError(25,0,0,6.552391e-16,6.552417e-16);
   grae->SetPoint(26,44.66836,3.316242e-14);
   grae->SetPointError(26,0,0,4.822833e-16,4.822853e-16);
   grae->SetPoint(27,56.23413,2.342079e-14);
   grae->SetPointError(27,0,0,3.549778e-16,3.549795e-16);
   grae->SetPoint(28,70.79458,1.681863e-14);
   grae->SetPointError(28,0,0,2.710382e-16,2.710396e-16);
   
   double  N_0 = N_0_1TeV*TMath::Power(E_c/1., -index_alpha1);
   
   TF1 *f_bpwl = new TF1("f_bpwl", brkpwlfunction, 0.1, 100, 5);
   f_bpwl->SetLineColor(4);
   f_bpwl->SetLineStyle(2);
   f_bpwl->SetParameters(crabunits*N_0, E_c, index_alpha1, index_alpha2, SS);  
   f_bpwl->SetParLimits(0, crabunits*N_0/10, crabunits*N_0*10);  // N_0*crabunits --> flux strength 
   f_bpwl->SetParLimits(1, E_c/10, E_c*10); // Ec --> break energy
   f_bpwl->SetParLimits(2, index_alpha1-2, index_alpha1+2);  // alpha1 --> index of the first power-law
   f_bpwl->SetParLimits(3, index_alpha2-2, index_alpha2+2);   // alpha2 --> index of the second power-law
   f_bpwl->SetParLimits(4, 1.e-4, 1.e2);  // SS --> for the width of the transition region

   grae->Fit("f_bpwl", "RM+");

   TCanvas *spectrum = new TCanvas("spectrum");
   gStyle->SetOptStat(0);
   spectrum->cd(); 
   TAxis *axisX = grae->GetXaxis(); 
   axisX->SetLimits(0.05,200);
   grae->GetHistogram()->SetMaximum(1.e-7);       
   grae->GetHistogram()->SetMinimum(1.e-15);
   grae->Draw("ap");
   gPad->SetLogx();
   gPad->SetLogy();
   gPad->SetGridy();
 
}

// **************************************************************************************************************
// Broken Power Law
// **************************************************************************************************************
Double_t brkpwlfunction(Double_t* x, Double_t* par)
{
  Double_t xx = x[0];
  double scale = xx/par[1];
  double expo1 = TMath::Power(scale, 1./par[4]);
  double expo2 = par[4]*(par[2]-par[3]);
  Double_t f = par[0]*TMath::Power(scale, -par[2])*TMath::Power( 1+expo1, expo2);
  return f;
}