#include "TCanvas.h"
#include "TPad.h"
#include "TFile.h"
#include "TMultiGraph.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TAxis.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TROOT.h"
#include "TF1.h"
#include <cstdlib> 


double Gamma1=5;
double Gamma2=3;
double Level1=6997.51;
double Level2=6998.36;
double wg1=0.017;
double wg2=0.66;

double a=250000;
double b=300000;
double maxval=0.3;
double max1=0.017;
double max2=0.66;
double relheigh;

double x11=270640;
double x21=270650;
double x31=270660;
double xc1=270670;
double x41=270680;
double x51=270690;
double x61=271000;


double x12=271530;
double x22=271540;
double x32=271550;
double xc2=271560;
double x42=271570;
double x52=271580;
double x62=271590;


double y11;
double y21;
double y31;
double y41;
double y51;
double y61;
double yc1;
double y12;
double y22;
double y32;
double y42;
double y52;
double y62;
double yc2;
double r11;
double r21;
double r31;
double r41;
double r51;
double r61;
double r12;
double r22;
double r32;
double r42;
double r52;
double r62;
double r1m;
double r2m;
double r3m;
double r4m;
double r5m;
double r6m;
double conv=1.e-15;
//double conv=1;


void lorfun() {
	ofstream results;
	results.open("Results.txt",  ios::out);
	TCanvas *c01 = new TCanvas("c01","c01",1280,1024);
	TMultiGraph *mg = new TMultiGraph();
	//TF1 *f1 = new TF1("lorentzian1", "[0]/(2*(TMath::Pi()))*(1/(([0]*[0])/4+(x-[1])*(x-[1])))", a,b);
	//TF1 *f2 = new TF1("lorentzian2", "[0]/(2*(TMath::Pi()))*(1/(([0]*[0])/4+(x-[1])*(x-[1])))", a,b);
	TF1 *f1 = new TF1("level1", "([0]/(2*(TMath::Pi()))*(1/(([0]*[0])/4+(x-[1])*(x-[1]))))*[2]", a,b);
	TF1 *f2 = new TF1("level2", "([0]/(2*(TMath::Pi()))*(1/(([0]*[0])/4+(x-[1])*(x-[1]))))*[2]", a,b);
	c01->SetLogy();
	f1->SetParameter(0,Gamma1); 
	f1->SetParameter(1,xc1);   
	f1->SetParameter(2,wg1);   
	f2->SetParameter(0,Gamma2); 
	f2->SetParameter(1,xc2);   
	f2->SetParameter(2,wg2);   
	f1->SetLineColor(kBlue);
	f2->SetLineColor(kRed);
	f1->SetTitle("");
	f2->SetTitle("");
	f2->GetXaxis()->SetTitle("Ep (eV)");
	f2->GetYaxis()->SetTitle("Height");
	 c01->Modified();
	f1->SetNpx(1000); 
	f2->SetNpx(1000);
	y11 = f1->Eval(x11); 
	y21 = f1->Eval(x21); 
	y31 = f1->Eval(x31); 
	y41 = f1->Eval(x41); 
	y51 = f1->Eval(x51); 
	y61 = f1->Eval(x61); 
	yc1 = f1->Eval(xc1); 
	y12 = f2->Eval(x12); 
	y22 = f2->Eval(x22); 
	y32 = f2->Eval(x32); 
	y42 = f2->Eval(x42); 
	y52 = f2->Eval(x52); 
	y62 = f2->Eval(x62); 
	yc2 = f2->Eval(xc2); 
	r11=y11/yc1*conv;
	r21=y21/yc1*conv;
	r31=y31/yc1*conv;
	r41=y41/yc1*conv;
	r51=y51/yc1*conv;
	r61=y61/yc1*conv;
	r12=y12/yc2*conv;
	r22=y22/yc2*conv;
	r32=y32/yc2*conv;
	r42=y42/yc2*conv;
	r52=y52/yc2*conv;
	r62=y62/yc2*conv;
	gPad->Update();
	f1->SetMaximum(maxval);
	f2->SetMaximum(maxval);
	max1=f1->GetMaximum();
	max2=f2->GetMaximum();
	relheigh=max2/max1;
	f1->Draw();
	f2->Draw("SAME");
	TLegend* leg = new TLegend(0.15, 0.75, .35, .85);
    leg->SetHeader("");
    leg->SetNColumns(1);
    leg->AddEntry(f1, "Level 6997.51 keV", "l");  
    leg->AddEntry(f2, "Level 6998.36 keV", "l");  
    leg->SetBorderSize(0);
    leg->Draw();
	c01->Print("sigma.png");
		results << "**************************************** "<< endl;
   	 		results << "Level 1 (keV)= " << Level1 <<  " " << "Ep1 (eV)= " << xc1 <<  " " <<  "Gamma 1 (eV)= " << Gamma1 <<  endl;
   	 		results << "Level 2(keV)= " << Level2 <<  " " << "Ep2 (eV)= " << xc2 <<  " " <<  "Gamma 2 (eV)= " << Gamma2 <<  endl;
   	 		results << "Maximum 1= " << max1 <<endl;
   	 		results << "Maximum 2= " << max2 <<endl;
   	 		results << "Maximum 2/Maximum 1= " << relheigh <<endl;
			results << "x11= " << x11 << " " << "y11= "  << y11 << " "  << "x12= " << x12 << " "  << "y12= "  << y12 << " "  << "r11= "  << r11 << " "  << "r12= "  << r12   <<endl;
			results << "x21= " << x21 << " " << "y21= "  << y21 << " " << "x22= " << x22  << " " << "y22= "  << y22 << " "  << "r21= "  << r21 << " "  << "r22= "  << r22    <<endl;
			results << "x31= " << x31 << " " << "y31= "  << y31 << " " << "x32= " << x32  << " " << "y32= "  << y32 << " "  << "r31= "  << r31 << " "  << "r32= "  << r32    <<endl;  
			results << "x41= " << x41 << " " << "y41= "  << y41 << " " << "x42= " << x42 << " "  << "y42= "  << y42 << " "  << "r41= "  << r41 << " "  << "r42= "  << r42    <<endl;
			results << "x51= " << x51 << " " << "y51= "  << y51 << " "  << "x52= " << x52 << " " << "y52= "  << y52 << " "  << "r51= "  << r51 << " "  << "r52= "  << r52    <<endl;
			results << "x61= " << x61 << " " << "y61= "  << y61 << " " << "x62= " << x62  << " "  << "y62= "  << y62 << " "  << "r61= "  << r61 << " "  << "r62= "  << r62   <<endl;
			
}