#include "TH1.h"
#include "TF1.h"
#include "iostream.h"
#include "fstream.h"
#include "iomanip"
#include <sstream>
#include "stdio.h"
#include <stdio.h>
#include <stdarg.h>
#include "math.h"
#include <string>
#include <exception>

using namespace std;
using std::ofstream;

const int size = 8192;

float eff1(float x) //x is energy
{
	float a0=6.;
	float a1=-0.471494;
	float a2=-0.126186;
	float a3=-0.261277; 
	float lne, eff;
	if (x<=1022) lne=a0+a1*log(x);
	if (x>1022) lne=a0+a1*log(x)+a2*log(x/1022)+a3*log(x/1022)*log(x/1022);
	eff=exp(lne);
	return eff;
}

float eff2(float x) //x is energy
{
	float a0=6.;
	float a1=-0.521824;
	float a2=-0.077314;
	float a3=-0.240241;
	float lne, eff;
	if (x<=1022) lne=a0+a1*log(x);
	if (x>1022) lne=a0+a1*log(x)+a2*log(x/1022)+a3*log(x/1022)*log(x/1022);
	eff=exp(lne);
	return eff;
}

float CtE(float C)
{
	float a = 1.0006;
	float b = 18.5;
	return a*C+b;
}

void read()
{
	ifstream if1; 
	if1.open("Yb.dat");
	int E1, E2;
	TH2F *h2 = new TH2F("h2","",size,0,size+.5,size,0,size+.5);
	while (if1 >> E1 >> E2)
	{
	h2->Fill(E1,E2);
	}
	if1.close();
	
	cout<<"Read file to histogram is done!"<<endl;
	cout<<"Efficiency calibration processing, it may take 30 minutes"<<endl;
	Double_t Nbefore = h2->Integral();
	cout<<"Nbefore = "<<Nbefore<<endl;
	for (int i=1; i<=size; i++)
	{
		for (int j=1; j<=size; j++)
		{
			float value = h2->GetBinContent(i,j);
			value = value / (eff1(CtE(i))*eff2(CtE(j)));
			h2->SetBinContent(i,j,value);
		}
	}
	
	cout<<"Normalizing ..."<<endl;
	Double_t Nafter = h2->Integral();
	cout<<"Nafter ="<<Nafter<<endl;
	for (int i=1; i<=size; i++)
	for (int j=1; j<=size; j++)
	{
		float value = h2->GetBinContent(i,j);
		h2->SetBinContent(i,j,value*Nbefore/Nafter);
	}
	
	cout<<"h2 entries = "<<h2->Integral()<<endl;
	
	cout<<"Save histogram"<<endl;
	TFile *f = new TFile("Yb.root","new");
	h2->Write();
	f->Close();
}

void det(string spename)
{
	TFile f("Yb.root");
	TH2F *h = (TH2F*)f.Get("h2");
	h->SetDirectory(gROOT);
	f.Close();

	TH1D *px = h->ProjectionX();
	TH1D *py = h->ProjectionY();
	//px->Draw();
	
	string s1, s2; 
	s1 = spename+"_det1.asc";
	s2 = spename+"_det2.asc";
	ofstream of1, of2;
	of1.open(s1.c_str());
	of2.open(s2.c_str());
	for (int i =1; i<=8192; i++)
		{
			of1 << i <<" "<< px->GetBinContent(i)<<endl;
			of2 << i <<" "<< py->GetBinContent(i)<<endl;
		}
	of1.close();
	of2.close();
}


void coin(string name, int BL1, int BR1, int L, int R, int BL2, int BR2)
{
	TFile f("Yb.root");
	TH2F *h = (TH2F*)f.Get("h2");
	h->SetDirectory(gROOT);
	f.Close();
	
	string s1, s2, s3; 
	s1 = name + "d.coin";
	s2 = name + "b.coin";
	TH1D *b11 = h->ProjectionX("b11",BL1,BR1);
	TH1D *b12 = h->ProjectionY("b12",BL1,BR1);
	TH1D *b21 = h->ProjectionX("b21",BL2,BR2);
	TH1D *b22 = h->ProjectionY("b22",BL2,BR2);
	TH1D *d1  = h->ProjectionX("d1",L,R);
	TH1D *d2  = h->ProjectionY("d2",L,R);
	
	TH1D *d = d1;
	TH1D *b = b11;
	d->Add(d2);
	b->Add(b12);
	b->Add(b21);
	b->Add(b22);
	
	// d->Draw();
	// b->SetLineColor(kRed);
	// b->Draw("SAME");
	
	ofstream of1, of2;
	of1.open(s1.c_str()); 
	of2.open(s2.c_str());
	
	for (int i=1; i<=size; i++)
	{
		of1<<i<<" "<<d->GetBinContent(i)<<endl;
		of2<<i<<" "<<b->GetBinContent(i)<<endl;
	}
		
	of1.close();
	of2.close();	
}

void coinyb()
{
// coin("6976",6939,6945,6946,6960,6961,6968);
// coin("6901",6865,6871,6872,6885,6886,6892);
// coin("6864",6812,6823,6824,6847,6848,6859);
// coin("6820",6779,6787,6788,6805,6806,6814);
// coin("6614",6574,6582,6583,6600,6601,6609);
// coin("6542",6504,6511,6512,6527,6528,6535);
// coin("6225",6190,6195,6196,6208,6209,6215);
// coin("6169",6133,6139,6140,6153,6154,6160);
// coin("6124",0,0,6096,6107,6108,6119);
// coin("6100",6055,6071,6072,6088,0,0);
// coin("6009",5971,5978,5979,5994,5995,6002);
// coin("5824",5785,5793,5794,5810,5811,5818);
// coin("5690",5655,5661,5662,5675,5676,5682);
// coin("5677",5642,5650,5651,5659,0,0);
// coin("5643",5602,5613,5614,5625,0,0);
// coin("5538",5501,5508,5509,5523,5524,5530);
// coin("5495",5464,5469,5470,5481,5482,5487);
// coin("5394",5358,5363,5364,5374,5375,5379);
// coin("5313",5280,5286,5287,5299,5300,5305);
coin("5271",5234,5240,5241,5254,5255,5261);
}