//Run as root -l klein_nishina.cpp+

#include "TMath.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TLine.h"
#include "TGaxis.h"
#include "TROOT.h"
#include "TLegend.h"
#include "TH1.h"
#include "TStyle.h"
#include <iostream>
#include <iomanip>
#include <sstream>
#include <vector>

const UInt_t kWidth=1600;//Canvas size for png output resolution
const UInt_t kHeight=kWidth*0.75;

using namespace TMath;
using std::vector;
using std::cout; using std::endl; using std::setprecision;

static const UInt_t CColorTable[11] = {kBlack, kRed, kBlue, kGreen+3, kOrange+7, kMagenta+2, kViolet+7, kCyan-2, kRed+2, kBlue-9,  kYellow+2};

void ResizeCanvas(TCanvas* const c, const UInt_t kWidth, const UInt_t kHeight)
{	
	if(!c) return;
	if(gROOT->IsBatch())
	{
		c->SetCanvasSize(kWidth,kHeight);
	}
	else
	{
		c->SetWindowSize(kWidth,kHeight);
	}
}

TString d2s(const Double_t x, const UShort_t precision)
{
	stringstream s;
	s.setf(ios::fixed);
	s << setprecision(precision) << x;
	return s.str();
}

TCanvas* CreateCanvas(const TString name)
{
	TCanvas* const c = new TCanvas(name,name,100 ,100, kWidth, kHeight);
	ResizeCanvas(c,kWidth,kHeight);

	c->SetFillColor(0);
	c->SetFillStyle(0);
	c->SetFrameFillColor(0);
	c->SetFrameFillStyle(0);
	c->Draw();
	
	return c;
}

Double_t f_theta(Double_t *x, Double_t *par)
{
	//See PhD thesis
	//constants are cancelled when performing the integral: 2pi/sigmaT dsigma_domega = 2pi/(pi a2 r2 / n3) / (2n...) * a2 * r2 / 2 * (P3...) = n3/(2n...)*(P3...)
	
	const Double_t theta = x[0]*TMath::Pi()/180.;
	const Double_t E = par[0];
	const Double_t me = 0.511;//MeV
	const Double_t eta = E/me;
	const Double_t P = 1./(1.+eta*(1.-Cos(theta)));
	
	//~ st(E)=1.0/g(E)**3*(2.0*g(E)*(2.0+g(E)*(1.0+g(E))*(8.0+g(E)))/(1.0+2.0*g(E))**2+(g(E)*(g(E)-2.)-2.)*log(1.+2.*g(E)))
	//~ s(x,E)=(P(x,E)**3+P(x,E)-P(x,E)**2*sin(x)*sin(x))*sin(x)/st(E)
	
	const Double_t sigmaT = 1./Power(eta,3)*(2.*eta*(2.+eta*(1.+eta)*(8.+eta))/Power(1.+2.*eta,2)+(eta*(eta-2.)-2.)*Log(1.+2.*eta));
	const Double_t dsigma_domega = (Power(P,3)+P-Power(P*Sin(theta),2))*Sin(theta);
	
	return dsigma_domega/sigmaT;
}

Double_t E_theta(Double_t *x, Double_t *par)
{
	//See PhD thesis
	
	const Double_t theta = x[0]*TMath::Pi()/180.;
	const Double_t E = par[0];
	const Double_t me = 0.511;//MeV
	const Double_t eta = E/me;
	const Double_t P = 1./(1.+eta*(1.-Cos(theta)));
	
	return P*E;
}


Double_t T_theta(Double_t *x, Double_t *par)
{
	//See PhD thesis
	
	const Double_t theta = x[0]*TMath::Pi()/180.;
	const Double_t E = par[0];
	const Double_t me = 0.511;//MeV
	const Double_t eta = E/me;
	const Double_t P = 1./(1.+eta*(1.-Cos(theta)));
	
	return E-P*E;
}

Double_t f_E(Double_t *x, Double_t *par)
{
	//See PhD thesis
	//constants are cancelled when performing the integral: 2pi/sigmaT dsigma_domega = 2pi/(pi a2 r2 / n3) / (2n...) * a2 * r2 / 2 * (P3...) = n3/(2n...)*(P3...)
	
	const Double_t E_prima = x[0];//MeV
	const Double_t me = 0.511;//MeV
	const Double_t E = par[0];
	const Double_t theta = ACos(1.+me*(1./E-1./E_prima));
	
	const Double_t eta = E/me;
	const Double_t P = 1./(1.+eta*(1.-Cos(theta)));
	
	//~ st(E)=1.0/g(E)**3*(2.0*g(E)*(2.0+g(E)*(1.0+g(E))*(8.0+g(E)))/(1.0+2.0*g(E))**2+(g(E)*(g(E)-2.)-2.)*log(1.+2.*g(E)))
	//~ s(x,E)=(P(x,E)**3+P(x,E)-P(x,E)**2*sin(x)*sin(x))*sin(x)/st(E)
	
	//We apply change of variable theta(E') -->  dtheta = dE' * dtheta/dE' = dE (-)/|sin(theta)|*me/E'^2 (derivative of arccosine)
	//Note: the - is cancelled when interchanging the limits, the module of the sinus is not necessary as theta is always positive, 0<theta< pi
	
	const Double_t sigmaT = 1./Power(eta,3)*(2.*eta*(2.+eta*(1.+eta)*(8.+eta))/Power(1.+2.*eta,2)+(eta*(eta-2.)-2.)*Log(1.+2.*eta))/(me/Power(E_prima,2));
	const Double_t dsigma_domega = (Power(P,3)+P-Power(P*Sin(theta),2));//**Sin(theta);
	
	return dsigma_domega/sigmaT;
}

Double_t f_E2(Double_t *x, Double_t *par)
{
	//scale to maximum
	const Double_t minEprima = 1./(2./0.511+1./par[0]);
	const Double_t aux = x[0];
	x[0]=minEprima;
	const Double_t max = f_E(x,par);
	x[0] = aux;
	return f_E(x,par)/max;
}

Double_t f_T(Double_t *x, Double_t *par)
{
	x[0]=par[0]-x[0];
	return f_E(x,par);
}

void klein_nishina(bool texOutput=false)
{
	gStyle->SetStripDecimals(kFALSE);
	gStyle->SetOptStat(0);
	gStyle->SetOptTitle(0);
	
	TString descriptor = "klein_nishina_theta";
	TCanvas * c = CreateCanvas(descriptor);
	c->SetGridx();

	TF1* f = nullptr;
	
	TLegend* aLegend = new TLegend(0.6,0.5,0.8,0.8);
	if(texOutput)
	{
		aLegend->SetHeader("E_{#gamma}#,/#,#text{MeV}");
		aLegend->SetTextSize(0.04);
	}
	else
	{
		aLegend->SetHeader("E_{#gamma} / MeV");
	}
			
	aLegend->SetFillColor(0);
	//~ title->SetTextAlign(11);
	//title->SetTextSize(1);
	aLegend->SetShadowColor(0);
	//~ title->SetFillStyle(4000);
	//~ title->SetBorderSize(0);
	
	Double_t E;
	
	vector<Double_t> energies = {0.511,1.,2., 4., 6.};
	
	for(UInt_t i=0;i<energies.size();++i)
	{
		E = energies[i];//MeV
		f = new TF1(d2s(E,3),f_theta,0.,180.,1);
		f->SetParameter(0,E);
		f->SetNpx(1000);
		f->SetLineColor(CColorTable[i%11]);
		if(i==0)
		{
			//~ f->GetHistogram()->GetXaxis()->SetTitle("#theta (#circ)");
			f->GetHistogram()->GetXaxis()->SetTitle("#theta");
			f->GetHistogram()->GetYaxis()->SetTitle("f_{#theta}(E_{#gamma})");
			//~ f->SetLineStyle(7);
			f->GetHistogram()->SetMaximum(1.);
			
			//~ for(UInt_t j=1;j<f->GetHistogram()->GetXaxis()->GetNbins();j+=100)
			//~ {
				//~ f->GetHistogram()->GetXaxis()->SetBinLabel(j,"a");
				//~ break;
			//~ }
			TAxis* a = f->GetHistogram()->GetXaxis();
			a->SetBit(TAxis::kLabelsHori);
			if(texOutput) a->SetBinLabel(1.,"$0$");
			else a->SetBinLabel(1.,"0");
			a->SetBinLabel(1000/4.,"#pi/4");
			a->SetBinLabel(1000/2.,"#pi/2");
			a->SetBinLabel(1000*3./4.,"3#pi/4");
			a->SetBinLabel(1000,"#pi");
			a->SetLabelSize(0.075);
			a->SetNdivisions(0);
			a->SetTickLength(0);
			f->Draw();
			//~ https://root.cern.ch/phpBB3/viewtopic.php?f=3&t=5063
			//~ https://root.cern.ch/phpBB3/viewtopic.php?f=3&t=19665&sid=e12692c6c06cfb9bbd0cbc7e196721a8&p=84335#p84335
			gPad->Update();

			TGaxis* ax = new TGaxis(0., gPad->GetUymin(),
								180., gPad->GetUymin(),
							   0, 1, 4+100*3, "UNW");
			//We need to repaint the gridx as we have overwritten the x-axis... ups
			ax->SetGridLength(0.7);
			//~ ax->SetLabelSize(0.);
			ax->Draw();
			
			//We add some extra grid on secondary tick marks
			TLine* t;
			
			vector<Double_t> extraGrid = {15.,30.,60.,75.,105.,120.,150.,165.};
			for(Double_t xg : extraGrid)
			{
				t = new TLine(xg,gPad->GetUymin(),xg,gPad->GetUymax());
				t->SetLineColor(1);
				t->SetLineStyle(3);
				t->SetLineWidth(1);
				t->Draw();
			}
			
			gPad->Update();
			
		}
		else
		{
			f->Draw("same");
		}
		//~ cout << f->Integral(0.,180.)/180.*TMath::Pi();
		cout << E << "\t" << f->Integral(0.,90.)/180*TMath::Pi() << endl;
		aLegend->AddEntry(f,f->GetName(),"l");
	}


	aLegend->Draw("same");
	
	c->Modified(); c->Update();
	if(!texOutput) c->SaveAs(descriptor+".png");
	else c->SaveAs(descriptor+".tex");

	
	descriptor = "compton_E_theta";
	c = CreateCanvas(descriptor);
	c->SetGridx();
	c->SetGridy();
	
	aLegend = new TLegend(0.6,0.5,0.8,0.8);
	if(texOutput)
	{
		aLegend->SetHeader("E_{#gamma}#,/#,#text{MeV}");
		aLegend->SetTextSize(0.04);
	}
	else
	{
		aLegend->SetHeader("E_{#gamma} / MeV");
	}
	
	aLegend->SetFillColor(0);
	//~ title->SetTextAlign(11);
	//title->SetTextSize(1);
	aLegend->SetShadowColor(0);
	//~ title->SetFillStyle(4000);
	//~ title->SetBorderSize(0);
	
	for(UInt_t i=0;i<energies.size();++i)
	{
		E = energies[i];//MeV
		f = new TF1(d2s(E,3),E_theta,0.,180.,1);
		f->SetParameter(0,E);
		const UInt_t Npx = 1000;
		f->SetNpx(Npx);
		f->SetLineColor(CColorTable[i%11]);
		//~ f->SetLineWidth(5);
		if(i==0)
		{
			f->GetHistogram()->GetXaxis()->SetTitle("#theta (#circ)");
			f->GetHistogram()->GetXaxis()->SetTitle("#theta");
			if(texOutput) f->GetHistogram()->GetYaxis()->SetTitle("{E_{#gamma}}'#,/#,#text{MeV}");
			else f->GetHistogram()->GetYaxis()->SetTitle("E_{#gamma}' / MeV");
			//~ f->SetLineStyle(7);
			f->GetHistogram()->SetMinimum(0.);
			f->GetHistogram()->SetMaximum(energies.back());
			
			TAxis* a = f->GetHistogram()->GetXaxis();
			a->SetBit(TAxis::kLabelsHori);
			if(texOutput) a->SetBinLabel(1.,"$0$");
			else a->SetBinLabel(1.,"0");
			a->SetBinLabel(Npx/4.,"#pi/4");
			a->SetBinLabel(Npx/2.,"#pi/2");
			a->SetBinLabel(Npx*3./4.,"3#pi/4");
			a->SetBinLabel(Npx,"#pi");
			a->SetLabelSize(0.075);
			a->SetNdivisions(0);
			a->SetTickLength(0);
			f->Draw();
			//~ https://root.cern.ch/phpBB3/viewtopic.php?f=3&t=5063
			//~ https://root.cern.ch/phpBB3/viewtopic.php?f=3&t=19665&sid=e12692c6c06cfb9bbd0cbc7e196721a8&p=84335#p84335
			gPad->Update();

			TGaxis* ax = new TGaxis(0., gPad->GetUymin(),
								180., gPad->GetUymin(),
							   0, 1, 4+100*3, "UNW");
			//~ ax->SetLabelSize(0.);
			//We need to repaint the gridx as we have overwritten the x-axis... ups
			ax->SetGridLength(0.7);
			//~ ax->SetLabelSize(0.);
			ax->Draw();
			
			//We add some extra grid on secondary tick marks
			TLine* t;
			
			vector<Double_t> extraGrid = {15.,30.,60.,75.,105.,120.,150.,165.};
			for(Double_t xg : extraGrid)
			{
				t = new TLine(xg,gPad->GetUymin(),xg,gPad->GetUymax());
				t->SetLineColor(1);
				t->SetLineStyle(3);
				t->SetLineWidth(1);
				t->Draw();
			}
			
			gPad->Update();
		}
		else
		{
			f->Draw("same");
		}
		
		aLegend->AddEntry(f,f->GetName(),"l");
	}


	aLegend->Draw("same");

	c->Modified(); c->Update();
	if(!texOutput) c->SaveAs(descriptor+".png");
	else c->SaveAs(descriptor+".tex");
	
	descriptor = "compton_T_theta";
	c = CreateCanvas(descriptor);
	c->SetGridx();
	c->SetGridy();
	
	aLegend = new TLegend(0.6,0.5,0.8,0.8);
	if(texOutput)
	{
		aLegend->SetHeader("E_{#gamma}#,/#,#text{MeV}");
		aLegend->SetTextSize(0.04);
	}
	else
	{
		aLegend->SetHeader("E_{#gamma} / MeV");
	}
	aLegend->SetFillColor(0);
	//~ title->SetTextAlign(11);
	//title->SetTextSize(1);
	aLegend->SetShadowColor(0);
	//~ title->SetFillStyle(4000);
	//~ title->SetBorderSize(0);
	
	for(UInt_t i=0;i<energies.size();++i)
	{
		E = energies[i];//MeV
		f = new TF1(d2s(E,3),T_theta,0.,180.,1);
		f->SetParameter(0,E);
		f->SetNpx(1000);
		f->SetLineColor(CColorTable[i%11]);
		if(i==0)
		{
			//~ f->GetHistogram()->GetXaxis()->SetTitle("#theta (#circ)");
			f->GetHistogram()->GetXaxis()->SetTitle("#theta");
			if(texOutput) f->GetHistogram()->GetYaxis()->SetTitle("{T_#mathrm{e}}'#,/#,#text{MeV}");
			else f->GetHistogram()->GetYaxis()->SetTitle("T_{e}' / MeV");
			//~ f->SetLineStyle(7);
			f->GetHistogram()->SetMaximum(energies.back());
			
			TAxis* a = f->GetHistogram()->GetXaxis();
			a->SetBit(TAxis::kLabelsHori);
			if(texOutput) a->SetBinLabel(1.,"$0$");
			else a->SetBinLabel(1.,"0");
			a->SetBinLabel(1000/4.,"#pi/4");
			a->SetBinLabel(1000/2.,"#pi/2");
			a->SetBinLabel(1000*3./4.,"3#pi/4");
			a->SetBinLabel(1000,"#pi");
			a->SetLabelSize(0.075);
			a->SetNdivisions(0);
			a->SetTickLength(0);
			f->Draw();
			//~ https://root.cern.ch/phpBB3/viewtopic.php?f=3&t=5063
			//~ https://root.cern.ch/phpBB3/viewtopic.php?f=3&t=19665&sid=e12692c6c06cfb9bbd0cbc7e196721a8&p=84335#p84335
			gPad->Update();

			TGaxis* ax = new TGaxis(0., gPad->GetUymin(),
								180., gPad->GetUymin(),
							   0, 1, 4+100*3, "UNW");
			//~ ax->SetLabelSize(0.);
			//We need to repaint the gridx as we have overwritten the x-axis... ups
			ax->SetGridLength(0.7);
			//~ ax->SetLabelSize(0.);
			ax->Draw();
			
			//We add some extra grid on secondary tick marks
			TLine* t;
			
			vector<Double_t> extraGrid = {15.,30.,60.,75.,105.,120.,150.,165.};
			for(Double_t xg : extraGrid)
			{
				t = new TLine(xg,gPad->GetUymin(),xg,gPad->GetUymax());
				t->SetLineColor(1);
				t->SetLineStyle(3);
				t->SetLineWidth(1);
				t->Draw();
			}
			
			gPad->Update();
		}
		else
		{
			f->Draw("same");
		}
		
		aLegend->AddEntry(f,f->GetName(),"l");
	}


	aLegend->Draw("same");

	c->Modified(); c->Update();
	if(!texOutput) c->SaveAs(descriptor+".png");
	else c->SaveAs(descriptor+".tex");
	
	descriptor = "f_E";
	c = CreateCanvas(descriptor);
	c->SetGridx();
	//~ c->SetGridy();
	//~ c->SetLogx();
	
	aLegend = new TLegend(0.6,0.5,0.8,0.8);
	if(texOutput)
	{
		aLegend->SetHeader("E_{#gamma}#,/#,#text{MeV}");
		aLegend->SetTextSize(0.04);
	}
	else
	{
		aLegend->SetHeader("E_{#gamma} / MeV");
	}
	
	aLegend->SetFillColor(0);
	//~ title->SetTextAlign(11);
	//title->SetTextSize(1);
	aLegend->SetShadowColor(0);
	//~ title->SetFillStyle(4000);
	//~ title->SetBorderSize(0);
	
	
	
	for(UInt_t i=0;i<energies.size();++i)
	{
		E = energies[i];//MeV
		
		//~ f = new TF1(d2s(E,3),f_E2,0.,energies.back(),1);
		f = new TF1(d2s(E,3),f_E,0.,energies.back(),1);
		f->SetParameter(0,E);
		f->SetNpx(1000);
		f->SetLineColor(CColorTable[i%11]);
		const Double_t minEprima = 1./(2./0.511+1./E);
		const Double_t maxEprima = E;
		f->SetRange(minEprima,E);
		//~ cout << E << " " << minEprima << endl;
		
		if(i==0)
		{
			TH1* h = new TH1F("Skeleton","",1,0.,energies.back());
			h->Draw();
			
			if(texOutput) h->GetXaxis()->SetTitle("{E_{#gamma}}'#,/#,#text{MeV}");
			else h->GetXaxis()->SetTitle("E_{#gamma}' / MeV");
			if(texOutput) h->GetYaxis()->SetTitle("f_{{E_{#gamma}}'}(E_{#gamma})");
			else h->GetYaxis()->SetTitle("f_{E_{#gamma}'}(E_{#gamma})");
			h->GetYaxis()->SetRangeUser(0,3.2);
			h->Draw();
			
			//~ f->GetHistogram()->GetXaxis()->SetTitle("E_{#gamma}' (MeV)");
			//~ f->GetHistogram()->GetYaxis()->SetTitle("f_{E_{#gamma}'}(E_{#gamma})");
			//~ f->SetLineStyle(7);
			//~ f->GetHistogram()->SetMaximum(1.);
			if(i!=0&&i!=4) f->Draw("same");
		}
		else
		{
			if(i!=0&&i!=4) f->Draw("same");
		}
		
		TLine* t = new TLine(minEprima,0.,minEprima,f->Eval(minEprima));
		t->SetLineColor(CColorTable[i%11]);
		t->SetLineStyle(2);
		if(i!=0&&i!=4) t->Draw();
		
		t = new TLine(maxEprima,0.,maxEprima,f->Eval(maxEprima));
		t->SetLineColor(CColorTable[i%11]);
		t->SetLineStyle(2);
		if(i!=0&&i!=4) t->Draw();
		
		if(i!=0&&i!=4) aLegend->AddEntry(f,f->GetName(),"l");
		else aLegend->AddEntry(f,"","");
		
		//~ cout << f->Integral(minEprima,maxEprima) << endl;
		
	}
	aLegend->Draw("same");

	c->Modified(); c->Update();
	if(!texOutput) c->SaveAs(descriptor+".png");
	else c->SaveAs(descriptor+".tex");
	
	descriptor = "f_T";
	c = CreateCanvas(descriptor);
	c->SetGridx();
	//~ c->SetGridy();
	//~ c->SetLogx();
	
	aLegend = new TLegend(0.6,0.5,0.8,0.8);
	if(texOutput)
	{
		aLegend->SetHeader("E_{#gamma}#,/#,#text{MeV}");
		aLegend->SetTextSize(0.04);
	}
	else
	{
		aLegend->SetHeader("E_{#gamma} / MeV");
	}
	aLegend->SetFillColor(0);
	//~ title->SetTextAlign(11);
	//title->SetTextSize(1);
	aLegend->SetShadowColor(0);
	//~ title->SetFillStyle(4000);
	//~ title->SetBorderSize(0);
	
	
	
	for(UInt_t i=0;i<energies.size();++i)
	{
		E = energies[i];//MeV
		
		//~ f = new TF1(d2s(E,3),f_E2,0.,energies.back(),1);
		f = new TF1(d2s(E,3),f_T,0.,energies.back(),1);
		f->SetParameter(0,E);
		f->SetNpx(1000);
		f->SetLineColor(CColorTable[i%11]);
		const Double_t minEprima = 1./(2./0.511+1./E);
		const Double_t maxEprima = E;
		f->SetRange(0,E-minEprima);
		//~ cout << E << " " << minEprima << endl;
		
		if(i==0)
		{
			TH1* h = new TH1F("Skeleton2","",1,0.,energies.back());
			h->Draw();
			
			if(texOutput) h->GetXaxis()->SetTitle("{T_#mathrm{e}}'#,/#,#text{MeV}");
			else h->GetXaxis()->SetTitle("T_{e}' / MeV");
			
			if(texOutput) h->GetYaxis()->SetTitle("f_{{T_#mathrm{e}}'}(E_{#gamma})");
			else h->GetYaxis()->SetTitle("f_{T_{e}'}(E_{#gamma})");
			//~ h->GetYaxis()->SetRangeUser(0,1.);
			h->GetYaxis()->SetRangeUser(0,3.2);
			h->Draw();
			
			//~ f->GetHistogram()->GetXaxis()->SetTitle("E_{#gamma}' (MeV)");
			//~ f->GetHistogram()->GetYaxis()->SetTitle("f_{{E_{#gamma}}'}(E_{#gamma})");
			//~ f->SetLineStyle(7);
			//~ f->GetHistogram()->SetMaximum(1.);
			if(i!=0&&i!=4) f->Draw("same");
		}
		else
		{
			if(i!=0&&i!=4) f->Draw("same");
		}
		
		TLine* t = new TLine(E-minEprima,0.,E-minEprima,f->Eval(E-minEprima));
		t->SetLineColor(CColorTable[i%11]);
		t->SetLineStyle(2);
		if(i!=0&&i!=4) t->Draw();
		
		t = new TLine(E-maxEprima,0.,E-maxEprima,f->Eval(E-maxEprima));
		t->SetLineColor(CColorTable[i%11]);
		t->SetLineStyle(2);
		//~ if(i!=0&&i!=4) t->Draw();
		
		if(i!=0&&i!=4) aLegend->AddEntry(f,f->GetName(),"l");
		else aLegend->AddEntry(f,"","");
		
		//~ cout << f->Integral(minEprima,maxEprima) << endl;
		
	}

	aLegend->Draw("same");

	c->Modified(); c->Update();
	if(!texOutput) c->SaveAs(descriptor+".png");
	else c->SaveAs(descriptor+".tex");

}