#include "TF1.h"
#include "TMath.h"
#include "TCanvas.h"
#include "Math.h"
gSystem->Load("libMathMore");
using namespace ROOT::Math; 
using namespace std;

TF1 *f1, *f2, *f_12;

double Func_1(double *, double *);
double Func_2(double *, double *);

double Func_12( double *, double * );
double Conv_12( double *, double * );

void main()
{
  TCanvas *c1 = new TCanvas("c1", "", 750, 750);
  c1->Divide(2, 2);

  f1 = new TF1("Differential Cross Section of Compton Scattering for 5MeV Photons",myfunction,0,10,3);
  f1->SetParameters(0.0794, 9.769, 4.9);
  f2 = new TF1("f2", Func_2, -20, 20, 2);
  f2->SetParameters(2.5, 0);

  f1->SetNpx(1000);
  f2->SetNpx(1000);
  
  c1->cd(1);
  gPad->SetLogy();
  f1->SetTitle("Differential Cross Section of Compton Scattering for 5MeV Photons");
  f1->Draw();
  f1->GetHistogram()->GetYaxis()->SetTitle("d#sigma/dE_{k}(barns}/MeV)");
  f1->GetHistogram()->GetXaxis()->SetTitle("Kinetic Energy of Ejected Electron = E_{k}(MeV)");  

  c1->cd(2);
  f2->SetTitle("Gaussian Distribution");
  f2->Draw();
  

  f_12 = new TF1("f_12", Func_12, 0, 30, 1);
  f_12->SetNpx(1000);
  f_12->SetParameter(0, -20);

  c1->cd(3);
  f_12->SetTitle("Product of Differential Cross Section of Compton Scattering for 5MeV Photons and Gaussian Distribution");
  f_12->Draw();
  f_12->GetHistogram()->GetYaxis()->SetTitle("d#sigma/dE_{k}(barns/MeV)");
  f_12->GetHistogram()->GetXaxis()->SetTitle("Kinetic Energy of Ejected Electron = E_{k}(MeV)");
  
  
  //    TF1 for convolution function
  TF1 *f_Conv = new TF1("f_Conv", Conv_12, 0, 30, 0 );
  f_Conv->SetNpx(1000);
  c1->cd(4);
  f_Conv->SetTitle("Convolution of Differential Cross Section of Compton Scattering for 5MeV Photons and Gaussian Distribution");
  f_Conv->Draw();
  f_Conv->GetHistogram()->GetYaxis()->SetTitle("d#sigma/dE_{k}(barns/MeV)");
  f_Conv->GetHistogram()->GetXaxis()->SetTitle("Kinetic Energy of Ejected Electron = E_{k}(MeV)");
    
}

Double_t myfunction(Double_t *x, Double_t *par)
{
    const Double_t pi = TMath::Pi(); 
    Float_t y = x[0];
    Double_t f = (pi*par[0])/(par[1]*par[2])*(2-(2*y)/(par[1]*(par[2]-y))+(y*y)/((par[1]*par[1])*(par[2]-y)*(par[2]-y))+(y*y/(par[2]*(par[2]-y))));
  
    return f;
}

double Func_2(double *x, double *par)
{
  const double pi = TMath::Pi();
  const double e = TMath::E();
  double result = 1/(par[0]*sqrt(2*pi))*pow(e, -pow((x[0]-par[1]),2)/(2*pow(par[0],2)));
  
  return result;
}

double Func_12( double *x, double *par )
{
  return f1->Eval(x[0])*f2->Eval(par[0]-x[0]);
}

double Conv_12( double *x, double *par )
{
  f_12->SetParameter(0, x[0]);
  return f_12->Integral(0, 4.9); //go beyond
}