/******************************************************************
 * MC plot LCE
 ******************************************************************/
// include own functions
#include "fileoperations.cc"
#include "TPC_Definition.cc"

// include C++ libs
#include <vector>
#include <iostream>
#include <stdio.h>
#include <fstream>
#include <string>
#include <sstream>
#include <algorithm>

#include <stdlib.h>

// include ROOT
#include "TROOT.h"
#include "TChain.h"
#include "TApplication.h"
#include "TSystemDirectory.h"
#include "TSystem.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TFile.h"
#include "TLine.h"
#include "TColor.h"
#include "TEntryList.h"
#include "TGaxis.h"
#include "TError.h"

using namespace std;

//void optPhot_S2(string datafile, string output_dir, string export_format, bool batch);
//void optPhot_S2(string datafile, int bin_z, int bin_r, int bin_rr, string output_dir, string export_format, bool batch);

/*=================================================================*/
void optPhot_S2(string datafile, string output_dir = "", string export_format = "png", bool batch = true) {
	// Some good binnings
	//TPC.Set_Bins(26,50,22) - default
	//TPC.Set_Bins(52,100,44)- nevents > 10000000
	optPhot_S2(datafile,26,50,22,output_dir,export_format,batch);
}

/*=================================================================*/

void optPhot_S2(string datafile, int bin_z, int bin_r, int bin_rr, string output_dir = "", string export_format = "png", bool batch = true) {
	
	//gErrorIgnoreLevel = kPrint, kInfo, kWarning, kError, kBreak, kSysError, kFatal;
	gErrorIgnoreLevel = kFatal;
	
	// read in datafilename and get working directory
	size_t found=datafile.find_last_of("/\\");
	string workingdirectory = datafile.substr(0,found);
	string datafilename = datafile.substr(found+1);
	string rawdatafilename = "";
	if (datafilename == "") {
		rawdatafilename = "optPhot";
	} else {
		size_t lastindex = datafilename.find_last_of("."); 
		rawdatafilename = datafilename.substr(0, lastindex); 
	}
	
	if (output_dir == "") {output_dir = workingdirectory;}
	else if (output_dir.find("\\", output_dir.size())) {output_dir = output_dir.substr(0,output_dir.size()-1);}
	if (fileexists(output_dir) == false) {
		cout << endl;
		cout << "x Error xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" << endl;
		cout << "output directory not found:" << endl;
		cout << "-> " << output_dir << endl;
		cout << "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" << endl;
		cout << endl;
		gApplication->Terminate();
	}
	
	Int_t canvas_x = 850;
	Int_t canvas_y = 800;

	char canvasfile[10000];
	char draw_selection[10000];
	
	TPC_Definition TPC(bin_z, bin_r, bin_rr);
	TChain *file_input_tree = new TChain("events/events");
	TNamed *G4MCname;
	
	char file_outname[10000];
	sprintf(file_outname,"%s/%s_S2.dat", output_dir.c_str(), rawdatafilename.c_str());
	
	ofstream file_outstat;
	file_outstat.open(file_outname);
	file_outstat << "============================================================" << "\n";
	sprintf(file_outname,"%s_S2.root", rawdatafilename.c_str());
	int nevents = 0;
	
	if (fileexists(datafile) == false) {
		cout << endl;
		cout << "x Error xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" << endl;
		cout << "file not found:" << endl;
		cout << "-> " << datafile << endl;
		cout << "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" << endl;
		cout << endl;
		gApplication->Terminate();
	}
	else if ((fileexists(datafile) == true) && (datafilename=="")) {
		file_outstat << "= reading datafiles ==== dir mode ==========================" << "\n";
		cout << "= reading datafiles ==== dir mode ==========================" << endl;
		string ext = ".root";
		TSystemDirectory dir(workingdirectory.c_str(), workingdirectory.c_str());
		TList *files = dir.GetListOfFiles();
		if (files) {
			TSystemFile *file;
			TString fname;
			TIter next(files);
			while ((file=(TSystemFile*)next())) {
				fname = file->GetName();
				if (!file->IsDirectory() && fname.EndsWith(ext.c_str()) && !(fname == file_outname)) {
					char filename[10000];
					sprintf(filename,"%s/%s", workingdirectory.c_str(), fname.Data());
					
					if (file_input_tree->GetEntries() == 0) {
						TFile *f = new TFile(filename,"READ");
						if (f->GetListOfKeys()->Contains("MC_TAG")) {
							f->GetObject("MC_TAG",G4MCname);
						}
						else {
							G4MCname = new TNamed("MC_TAG","Xenon1t");
						}
						TPC.Init(G4MCname);
						f->Close();
					}
					
					file_input_tree->AddFile(filename); 
					nevents = file_input_tree->GetEntries();
					file_outstat << " file: " << fname.Data() << " " << nevents << " events total" << "\n";
					cout << " file: " << fname.Data() << " " << nevents << " events total" << endl;
				}
			}
		}		
	}
	else {
		file_outstat << "= reading datafile ===== single file =======================" << "\n";
		cout << "= reading datafile ===== single file =======================" << endl;
		if (file_input_tree->GetEntries() == 0) {
			TFile *f = new TFile(datafile.c_str(),"READ");
			if (f->GetListOfKeys()->Contains("MC_TAG")) {
				f->GetObject("MC_TAG",G4MCname);
			}
			else {
				G4MCname = new TNamed("MC_TAG","Xenon1t");
			}
			TPC.Init(G4MCname);
			f->Close();
		}
		file_input_tree->AddFile(datafile.c_str()); 
		nevents = file_input_tree->GetEntries();
		file_outstat << " file: " << datafilename << " " << nevents << " events " << "\n";
		cout << " file: " << datafilename << " " << nevents << " events " << endl;
	}
	file_outstat << "============================================================" << "\n";
	cout << "============================================================" << endl;
	
	file_input_tree->SetAlias("rrp_pri","(xp_pri*xp_pri + yp_pri*yp_pri)/10./10.");  
	
	// generate plots  
	sprintf(file_outname,"%s/%s_S2.root", output_dir.c_str(), rawdatafilename.c_str());
	TFile *file_outplot = new TFile(file_outname,"RECREATE");
	
	file_outstat << "= geometry parameters ======================================" << "\n";
	file_outstat << "binning: " << bin_z << " " << bin_r << " " << bin_rr << "\n";
	file_outstat << "chamber_minZ: " << TPC.Get_chamber_minZ() << "\n";
	file_outstat << "chamber_maxZ: " << TPC.Get_chamber_maxZ() << "\n";
	file_outstat << "chamber_minRR: " << TPC.Get_chamber_minRR() << "\n";
	file_outstat << "chamber_maxRR: " << TPC.Get_chamber_maxRR() << "\n";
	file_outstat << "GXe_minZ: " << TPC.Get_LXe_minZ() << "\n";
	file_outstat << "GXe_maxZ: " << TPC.Get_LXe_maxZ() << "\n";
	file_outstat << "GXe_minRR: " << TPC.Get_LXe_minRR() << "\n";
	file_outstat << "GXe_maxRR: " << TPC.Get_LXe_maxRR() << "\n";
	file_outstat << "LCE_min: " << TPC.Get_LCE_min() << "\n";
	file_outstat << "LCE_max: " << TPC.Get_LCE_max() << "\n";
	file_outstat << "============================================================" << "\n";
	
	TGaxis::SetMaxDigits(3);
	//TGaxis::SetExponentOffset(-0.01, 0.01, "y"); // X and Y offset for Y axis
	TGaxis::SetExponentOffset(0.01, -0.0325, "x"); // Y and Y offset for X axis
	
	const Int_t NRGBs = 5;
	const Int_t NCont = 255;
	static Int_t ColPalette[255];
	Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
	Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
	Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
	Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
	
	Int_t FI = TColor::CreateGradientColorTable(NRGBs,stops,red,green,blue,NCont);
    for (int i=0; i<NCont; i++) ColPalette[i] = FI+i;
		
	TStyle *style_1D = new TStyle("1D","1D");
	style_1D->SetCanvasColor(10);
	style_1D->SetTitleFillColor(0);
	style_1D->SetTitleBorderSize(1);
	style_1D->SetOptStat(0);
	style_1D->SetPadLeftMargin(0.105);
	style_1D->SetPadRightMargin(0.09);
	style_1D->SetPadTopMargin(0.075);
	style_1D->SetPadBottomMargin(0.075);
	
	style_1D->SetTitleOffset(1.,"X");
	style_1D->SetTitleOffset(1.45,"Y");
	style_1D->SetTitleOffset(1.35,"Z");
	
	style_1D->SetPalette(NCont,ColPalette);
	style_1D->SetNumberContours(NCont);
	style_1D->cd();
	
	TStyle *style_2D = new TStyle("2D","2D");
	style_2D->SetCanvasColor(10);
	style_2D->SetTitleFillColor(0);
	style_2D->SetTitleBorderSize(1);
	style_2D->SetOptStat(0);
	style_2D->SetPadLeftMargin(0.105);
	style_2D->SetPadRightMargin(0.15);
	style_2D->SetPadTopMargin(0.075);
	style_2D->SetPadBottomMargin(0.075);
	
	style_2D->SetTitleOffset(1.,"X");
	style_2D->SetTitleOffset(1.45,"Y");
	style_2D->SetTitleOffset(1.35,"Z");
	
	style_2D->SetPalette(NCont,ColPalette);
	style_2D->SetNumberContours(NCont);
	style_2D->cd();
	
	gStyle->SetPalette(NCont,ColPalette);
	
	file_outstat << "= analyse ttree ============================================" << "\n";
	file_outstat << "generated events: " << file_input_tree->GetEntries() << "\n";
	
	file_input_tree->Draw(">>elist_top","(ntpmthits > 0)","goff");
	TEntryList *elist_top = (TEntryList*)gDirectory->Get("elist_top");
	file_outstat << "Top PMT hits: " << elist_top->GetEntriesToProcess() << "\n";
	
	file_input_tree->Draw(">>elist_all","(nbpmthits > 0 || ntpmthits > 0)","goff");
	TEntryList *elist_all = (TEntryList*)gDirectory->Get("elist_all");
	file_outstat << "All PMT hits: " << elist_all->GetEntriesToProcess() << "\n";
	file_outstat << "AFT (without QEs and CE): " << (double)elist_top->GetEntriesToProcess()/(double)elist_all->GetEntriesToProcess() << "\n";
	
	file_outstat << "============================================================" << "\n";
	
	gROOT->SetBatch(kTRUE);
	/*=================================================================*/
	// X vs. Y of generated events LXe
	/*=================================================================*/
	style_2D->cd();
	TCanvas *c_xy = new TCanvas("xy","xy",canvas_x,canvas_x);
	TH2F* h_xy = new TH2F("xy_pri", "X vs. Y generated events", TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR(), TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR());
	h_xy->SetXTitle("X [cm]");
	h_xy->GetXaxis()->CenterTitle();
	h_xy->SetYTitle("Y [cm]");
	h_xy->GetYaxis()->CenterTitle();
	h_xy->SetZTitle("events [#]");
	h_xy->GetZaxis()->CenterTitle();
	sprintf(draw_selection,"zp_pri/10<=%f && zp_pri/10>=%f && rrp_pri>=%f && rrp_pri<=%f",TPC.Get_GXe_maxZ(),TPC.Get_GXe_minZ(),TPC.Get_LXe_minRR(),TPC.Get_LXe_maxRR());
	file_input_tree->Draw("xp_pri/10. : yp_pri/10. >> xy_pri", draw_selection, "goff");
	h_xy->Draw("colz");
	if (file_outplot) c_xy->Write();
	
	/*=================================================================*/
	// X vs. Y detected events (TOP + BOTTOM PMTs)
	/*=================================================================*/
	style_2D->cd();
	TCanvas *c_xy_det = new TCanvas("xy_det","xy_det",canvas_x,canvas_x);
	TH2F* h_xy_det = new TH2F("xy_det", "X vs. Y detected events (ALL PMTs)", TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR(), TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR());
	h_xy_det->SetXTitle("X [cm]");
	h_xy_det->GetXaxis()->CenterTitle();
	h_xy_det->SetYTitle("Y [cm]");
	h_xy_det->GetYaxis()->CenterTitle();
	h_xy_det->SetZTitle("detected [#]");
	h_xy_det->GetZaxis()->CenterTitle();
	sprintf(draw_selection,"zp_pri/10<=%f && zp_pri/10>=%f && rrp_pri>=%f && rrp_pri<=%f && (nbpmthits > 0 || ntpmthits > 0)",TPC.Get_GXe_maxZ(),TPC.Get_GXe_minZ(),TPC.Get_LXe_minRR(),TPC.Get_LXe_maxRR());
	file_input_tree->Draw("xp_pri/10. : yp_pri/10. >> xy_det", draw_selection, "goff");
	h_xy_det->Draw("colz");
	if (file_outplot) c_xy_det->Write();
	
	/*=================================================================*/
	// X vs. Y detected events (TOP PMTs)
	/*=================================================================*/
	style_2D->cd();
	TCanvas *c_xy_det_top = new TCanvas("xy_det_top","xy_det_top",canvas_x,canvas_x);
	TH2F* h_xy_det_top = new TH2F("xy_det_top", "X vs. Y detected events (TOP PMTs)", TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR(), TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR());
	h_xy_det_top->SetXTitle("X [cm]");
	h_xy_det_top->GetXaxis()->CenterTitle();
	h_xy_det_top->SetYTitle("Y [cm]");
	h_xy_det_top->GetYaxis()->CenterTitle();
	h_xy_det_top->SetZTitle("TOP detected [#]");
	h_xy_det_top->GetZaxis()->CenterTitle();
	sprintf(draw_selection,"zp_pri/10<=%f && zp_pri/10>=%f && rrp_pri>=%f && rrp_pri<=%f && (ntpmthits > 0)",TPC.Get_GXe_maxZ(),TPC.Get_GXe_minZ(),TPC.Get_LXe_minRR(),TPC.Get_LXe_maxRR());
	file_input_tree->Draw("xp_pri/10. : yp_pri/10. >> xy_det_top", draw_selection, "goff");
	h_xy_det_top->Draw("colz");
	if (file_outplot) c_xy_det_top->Write();
	
	/*=================================================================*/
	// X vs. Y detected events (BOTTOM PMTs)
	/*=================================================================*/
	style_2D->cd();
	TCanvas *c_xy_det_bottom = new TCanvas("xy_det_bottom","xy_det_bottom",canvas_x,canvas_x);
	TH2F* h_xy_det_bottom = new TH2F("xy_det_bottom", "X vs. Y detected events (BOTTOM PMTs)", TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR(), TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR());
	h_xy_det_bottom->SetXTitle("X [cm]");
	h_xy_det_bottom->GetXaxis()->CenterTitle();
	h_xy_det_bottom->SetYTitle("Y [cm]");
	h_xy_det_bottom->GetYaxis()->CenterTitle();
	h_xy_det_bottom->SetZTitle("BOTTOM detected [#]");
	h_xy_det_bottom->GetZaxis()->CenterTitle();
	sprintf(draw_selection,"zp_pri/10<=%f && zp_pri/10>=%f && rrp_pri>=%f && rrp_pri<=%f && (nbpmthits > 0)",TPC.Get_GXe_maxZ(),TPC.Get_GXe_minZ(),TPC.Get_LXe_minRR(),TPC.Get_LXe_maxRR());
	file_input_tree->Draw("xp_pri/10. : yp_pri/10. >> xy_det_bottom", draw_selection, "goff");
	h_xy_det_bottom->Draw("colz");
	if (file_outplot) c_xy_det_bottom->Write();
	
	/*=================================================================*/
	if (!batch) {gROOT->SetBatch(kFALSE);}
	/*=================================================================*/
	/*=================================================================*/
	// LCE of X vs. Y
	/*=================================================================*/	
	style_2D->cd();
	TCanvas *c_LCE_xy = new TCanvas("LCE_xy","LCE_xy",canvas_x,canvas_x);
	TH2F* h_LCE_xy = new TH2F("LCE_xy", "LCE S2 of X vs. Y (ALL PMTs)", TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR(), TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR());
	h_LCE_xy->SetXTitle("X [cm]");
	h_LCE_xy->GetXaxis()->CenterTitle();
	h_LCE_xy->SetYTitle("Y [cm]");
	h_LCE_xy->GetYaxis()->CenterTitle();
	h_LCE_xy->SetZTitle("LCE [%]");
	h_LCE_xy->GetZaxis()->CenterTitle();
	h_LCE_xy->Sumw2();
	h_LCE_xy->Divide(h_xy_det, h_xy, 1.,1., "b");
	h_LCE_xy->Scale(100.);
	h_LCE_xy->Draw("colz");
	if (file_outplot) c_LCE_xy->Write();
	sprintf(canvasfile,"%s/%s_S2_xy_LCE.%s", output_dir.c_str(), rawdatafilename.c_str(), export_format.c_str());
	if (!(export_format=="")) c_LCE_xy->SaveAs(canvasfile);
	
	/*=================================================================*/
	// LCE of X vs. Y (TOP PMTs)
	/*=================================================================*/
	style_2D->cd();
	TCanvas *c_LCE_xy_top = new TCanvas("LCE_xy_top","LCE_xy_top",canvas_x,canvas_x);
	TH2F* h_LCE_xy_top = new TH2F("LCE_xy_top", "LCE S2 of X vs. Y (TOP PMTs)", TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR(), TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR());
	h_LCE_xy_top->SetXTitle("X [cm]");
	h_LCE_xy_top->GetXaxis()->CenterTitle();
	h_LCE_xy_top->SetYTitle("Y [cm]");
	h_LCE_xy_top->GetYaxis()->CenterTitle();
	h_LCE_xy_top->SetZTitle("TOP LCE [%]");
	h_LCE_xy_top->GetZaxis()->CenterTitle();
	h_LCE_xy_top->Sumw2();
	h_LCE_xy_top->Divide(h_xy_det_top, h_xy, 1.,1., "b");
	h_LCE_xy_top->Scale(100.);
	h_LCE_xy_top->Draw("colz");
	if (file_outplot) c_LCE_xy_top->Write();
	sprintf(canvasfile,"%s/%s_S2_xy_LCE_top.%s", output_dir.c_str(), rawdatafilename.c_str(), export_format.c_str());
	if (!(export_format=="")) c_LCE_xy_top->SaveAs(canvasfile);
	
	/*=================================================================*/
	gROOT->SetBatch(kTRUE);
	/*=================================================================*/
	/*=================================================================*/
	// LCE of X vs. Y (BOTTOM PMTs)
	/*=================================================================*/
	style_2D->cd();
	TCanvas *c_LCE_xy_bottom = new TCanvas("LCE_xy_bottom","LCE_xy_bottom",canvas_x,canvas_x);
	TH2F* h_LCE_xy_bottom = new TH2F("LCE_xy_bottom", "LCE S2 of X vs. Y (BOTTOM PMTs)", TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR(), TPC.Get_nbinsR(), -TPC.Get_LXe_maxR(), TPC.Get_LXe_maxR());
	h_LCE_xy_bottom->SetXTitle("X [cm]");
	h_LCE_xy_bottom->GetXaxis()->CenterTitle();
	h_LCE_xy_bottom->SetYTitle("Y [cm]");
	h_LCE_xy_bottom->GetYaxis()->CenterTitle();
	h_LCE_xy_bottom->SetZTitle("BOTTOM LCE [%]");
	h_LCE_xy_bottom->GetZaxis()->CenterTitle();
	h_LCE_xy_bottom->Sumw2();
	h_LCE_xy_bottom->Divide(h_xy_det_bottom, h_xy, 1.,1., "b");
	h_LCE_xy_bottom->Scale(100.);
	h_LCE_xy_bottom->Draw("colz");
	if (file_outplot) c_LCE_xy_bottom->Write();
	sprintf(canvasfile,"%s/%s_S2_xy_LCE_bottom.%s", output_dir.c_str(), rawdatafilename.c_str(), export_format.c_str());
	if (!(export_format=="")) c_LCE_xy_bottom->SaveAs(canvasfile);
	
	file_outstat << "============================================================" << "\n";
	file_outstat.close();	
	if (batch) {file_outplot->Close();}  
}

/*=================================================================*/