#include "GetRunList.C"
#include "GetCharge.C"
#include "TCanvas.h"
#include "TChain.h"
#include "TProof.h"
#include "TList.h"
#include "TH1.h"
#include "TF1.h"
#include "TMath.h"

void SetProofInputList(TProof * p,TList *input_list);
double GetEventNo(TH1 * hist,double tof_min,double tof_max,TString output_name);
double Normalize_by_Z_factor(TH1* z,float min,float max);

int main() {
//int GetEvents() {

// define variables

	TString HRS = "R";
	float p = 1.073;
	int bin_number = 100;

	double tof_min, tof_max, th_min, th_max, ph_min, ph_max;
	double p_min, p_max, zHe3_min, zHe3_max, zC12_min, zC12_max;
	ifstream file("cut_list");   
	char ch[6];             

// read cut definitions from file and add dependant variables

	file >> ch >> tof_min  >> tof_max;
	file >> ch >> th_min   >> th_max;
	file >> ch >> ph_min   >> ph_max;
	file >> ch >> p_min    >> p_max;
	file >> ch >> zHe3_min >> zHe3_max;
	file >> ch >> zC12_min >> zC12_max;
	file.close();

	double C12_density = 4.205964e22;
	double He3_density = 3.168208e23*(zHe3_max-zHe3_min)/0.2;

// define PROOF inputs - will be used for all Process() commands

	TList * input_list = new TList();

	input_list->Add(new TNamed("input_HRS", HRS.Data()            ));
	input_list->Add(new TNamed("input_p",   Form("%g",p)          ));
	input_list->Add(new TNamed("bin_number",Form("%d",bin_number) ));
	input_list->Add(new TNamed("p_min",     Form("%g",p_min)      ));
	input_list->Add(new TNamed("p_max",     Form("%g",p_max)      ));
	input_list->Add(new TNamed("zHe3_min",  Form("%g",zHe3_min)   ));
	input_list->Add(new TNamed("zHe3_max",  Form("%g",zHe3_max)   ));
	input_list->Add(new TNamed("zC12_min",  Form("%g",zC12_min)   ));
	input_list->Add(new TNamed("zC12_max",  Form("%g",zC12_max)   ));
	input_list->Add(new TNamed("ph_min",    Form("%g",ph_min)     ));
	input_list->Add(new TNamed("ph_max",    Form("%g",ph_max)     ));
	input_list->Add(new TNamed("th_min",    Form("%g",th_min)     ));
	input_list->Add(new TNamed("th_max",    Form("%g",th_max)     ));
	input_list->Add(new TNamed("tof_min",   Form("%g",tof_min)    ));
	input_list->Add(new TNamed("tof_max",   Form("%g",tof_max)    ));
	
// Open PROOF session and add inputs
	input_list->Add(new TNamed("input_target","C12"));	
	TProof * proof_instance = TProof::Open("");
	SetProofInputList(proof_instance,input_list);

// Process C12 in runs

	TChain * chain = GetRunList("run_list",p,"C12","in");
	chain->SetProof(); chain->Process("T.C+");		
	TH1 * C12_in_tof_hist = (TH1 *)proof_instance->GetOutputList()->FindObject("thist");

// Reopen PROOF session and add inputs

	proof_instance->Close();
	proof_instance  = TProof::Open("");
	SetProofInputList(proof_instance,input_list);

// Process C12 out runs

	chain = GetRunList("run_list",p,"C12","out");
	chain->SetProof(); chain->Process("T.C+");		
	TH1 * C12_out_tof_hist = (TH1 *)proof_instance->GetOutputList()->FindObject("thist");

// Reopen PROOF session and add inputs
	input_list->Remove(input_list->FindObject("input_target"));
	input_list->Add(new TNamed("input_target","He3"));	
	proof_instance->Close();
	proof_instance  = TProof::Open("");
	SetProofInputList(proof_instance,input_list);

// Process He3 in runs

	chain = GetRunList("run_list",p,"He3","in");
	chain->SetProof(); chain->Process("T.C+");		
	TH1 * He3_in_tof_hist = (TH1 *)proof_instance->GetOutputList()->FindObject("thist");
	TH1 * He3_in_z_hist   = (TH1 *)proof_instance->GetOutputList()->FindObject("zhist");
// Reopen PROOF session and add inputs

	proof_instance->Close();
	proof_instance  = TProof::Open("");
	SetProofInputList(proof_instance,input_list);

// Process He3 in runs

	chain = GetRunList("run_list",p,"He3","out");
	chain->SetProof(); chain->Process("T.C+");		
	TH1 * He3_out_tof_hist = (TH1 *)proof_instance->GetOutputList()->FindObject("thist");
	TH1 * He3_out_z_hist   = (TH1 *)proof_instance->GetOutputList()->FindObject("zhist");

// Fit and plot

	float z_in_factor  = Normalize_by_Z_factor(He3_in_z_hist, zHe3_min,zHe3_max);
	float z_out_factor = Normalize_by_Z_factor(He3_out_z_hist,zHe3_min,zHe3_max);

	double Q_C12_in  = GetCharge("run_data.csv",p,"C12","in" );
	double Q_C12_out = GetCharge("run_data.csv",p,"C12","out");
	double Q_He3_in  = GetCharge("run_data.csv",p,"He3","in" );
	double Q_He3_out = GetCharge("run_data.csv",p,"He3","out");
  
	double div = (tof_max-tof_min)/bin_number*(th_max-th_min)*(ph_max-ph_min);
	float C12_in  = GetEventNo(C12_in_tof_hist ,tof_min,tof_max,"ci.ps")/div/Q_C12_in /C12_density;
	float C12_out = GetEventNo(C12_out_tof_hist,tof_min,tof_max,"co.ps")/div/Q_C12_out/C12_density;
	float He3_in  = GetEventNo(He3_in_tof_hist ,tof_min,tof_max,"hi.ps")/div*z_in_factor /Q_He3_in /He3_density;
	float He3_out = GetEventNo(He3_out_tof_hist,tof_min,tof_max,"ho.ps")/div*z_out_factor/Q_He3_out/He3_density;

   printf("He3 Q: %g %g\n",Q_He3_in,Q_He3_out);
   printf("C12 Q: %g %g\n",Q_C12_in,Q_C12_out);
	printf("Normalized Yields: %g %g %g %g\n",C12_in,C12_out,He3_in,He3_out);
	printf("He3 PP,EP: %g %g\n",He3_in-He3_out,He3_out);
	printf("C12 PP,EP: %g %g\n",C12_in-C12_out,C12_out);
	printf("Ratios: %g %g\n",C12_out/He3_out,(C12_in-C12_out)/(He3_in-He3_out));	
	return 0;
}

void SetProofInputList(TProof * p,TList *input_list) {

	int k = 0;
	while(input_list->At(k)) {
		p->AddInput(input_list->At(k));
		k++;
	}

}

double Normalize_by_Z_factor(TH1* z,float min,float max) {

	TAxis *axis = z->GetXaxis();
	int bmin = axis->FindBin(min);
	int bmax = axis->FindBin(max); 
	double cur_int =  z->Integral(bmin,bmax);
	double max_int  = (bmax-bmin+1)*z->GetMaximum();
	printf("z correction: %g\n",max_int/cur_int);
	TCanvas * canvas = new TCanvas("canvas","test canvas",100,100,500,400);
	z->Draw();
	canvas->SaveAs("z.ps");
	return max_int/cur_int;	

}

double GetEventNo(TH1 * hist,double tof_min,double tof_max,TString output_name = "") {

	TCanvas * canvas = new TCanvas("canvas","test canvas",100,100,500,400);
	TF1 * gaussian = new TF1("gaussian","gaus",tof_min,tof_max);
	hist->Draw();
	hist->Fit("gaussian","QW","",tof_min,tof_max);
	if (output_name != "") canvas->SaveAs(output_name);
	delete canvas;
	return gaussian->GetParameter(0)*gaussian->GetParameter(2)*TMath::Sqrt(2*TMath::Pi());

}