#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <numeric>
#include <chrono>

void fill_hists(string input_file_name, string output_file_name="output.root"){
	// Input file is a root tree
	// Call with: root 'fill_hists.cpp("output_061.root")'
	
	TFile * input_file = new TFile(input_file_name.c_str());
	TTree * clover_array = (TTree*) input_file->Get("clover_array");
	
	TFile * output_file = new TFile(output_file_name.c_str(),"RECREATE");

	int n_events = clover_array->GetEntries();
	cout << n_events << endl;
	
	double amplitude_beam_monitor[16];
	double amplitude_clover_135[16];
	double amplitude_clover_90[16];
	double integration_long[16];
	double time_clover_135[16];
	double time_clover_90[16];
	double channel_time[16];

	double time_p1=0.098; //in ns
	
	//CeBr3 calibration parameters 
	vector <double> cebr_p0({-13.1403447516191,3.69678252687933,-9.1702188841228,-3.47442299494794,-1.30462351564509,1.67482470826693,9.927802495021,-8.08442860771037,5.98822953506204,1.60528564159611,-7.69057784536478,-1.25458432139152});
	vector <double> cebr_p1({2.94004530168439,2.87080299062921,3.14403301258813,2.97104155682001,2.70746623846035,2.80349351708509,2.89812550447883,3.20443283370779,2.99227742725396,2.64164778119365,3.11605043416524,2.92331319375623});

	
	// Calibration params from a file
	ifstream calfile("encal_quad_56Fe.txt");

    string line, val;
    vector <double> row; 

    vector <double> clover_90_step;
    vector <double> clover_135_step;
	vector <vector<double>> clover_90_p0;
	vector <vector<double>> clover_90_p1;
	vector <vector<double>> clover_90_p2;
	vector <vector<double>> clover_135_p0;
	vector <vector<double>> clover_135_p1;
	vector <vector<double>> clover_135_p2;
    
    for (int i=0; i<2; i++){
        clover_90_p0.push_back(vector<double>());
        clover_90_p1.push_back(vector<double>());
        clover_90_p2.push_back(vector<double>());
        clover_135_p0.push_back(vector<double>());
        clover_135_p1.push_back(vector<double>());
        clover_135_p2.push_back(vector<double>());
    }
    
    string clover_90 = "clover_90";
    string clover_135 = "clover_135";
	string bad_crystal = "clover_135_5";

    int i = 0;
	while(getline(calfile, line)) {
        row.clear();
        // Check module by peeking at beginning of line
        if (strncmp(line.c_str(), clover_90.c_str(), clover_90.size()) ==0){

            stringstream s(line);
    
            // read every column data of a row and
            // store it in a string variable
            while (getline(s, val, '\t')) {
                // add all the column data
                // of a row to a vector
                row.push_back(strtold(val.c_str(), NULL));
            }
            // Write parameters from each row to their respective vectors
			clover_90_step.push_back(row[1]);
            clover_90_p0[0].push_back(row[2]);
            clover_90_p1[0].push_back(row[3]);
            clover_90_p2[0].push_back(row[4]);

            if (row.size()>7){
                clover_90_p0[1].push_back(row[8]);
                clover_90_p1[1].push_back(row[9]);
                clover_90_p2[1].push_back(row[10]);
            }
            else{
                // Not a step function -> fill with zeros 
                clover_90_p0[1].push_back(0);
                clover_90_p1[1].push_back(0);
                clover_90_p2[1].push_back(0);
            }
        }
        // Check module by peeking at beginning of line
        else if (strncmp(line.c_str(), clover_135.c_str(), clover_135.size()) ==0){
            
            stringstream s(line);
            // read every column data of a row and
            // store it in a string variable
            while (getline(s, val, '\t')) {
                // add all the column data
                // of a row to a vector
                row.push_back(strtold(val.c_str(), NULL));
            }
            // Write parameters from each row to their respective vectors
			// Write zeros for bad crystal on B4
			if (strncmp(line.c_str(), bad_crystal.c_str(), bad_crystal.size()) ==0){
				clover_135_step.push_back(70000);
				clover_135_p0[0].push_back(0);
				clover_135_p1[0].push_back(0);
				clover_135_p2[0].push_back(0);
				clover_135_p0[1].push_back(0);
				clover_135_p1[1].push_back(0);
				clover_135_p2[1].push_back(0);
			}
			else{
				clover_135_step.push_back(row[1]);
				clover_135_p0[0].push_back(row[2]);
				clover_135_p1[0].push_back(row[3]);
				clover_135_p2[0].push_back(row[4]);

				if (row.size()>7){
					clover_135_p0[1].push_back(row[8]);
					clover_135_p1[1].push_back(row[9]);
					clover_135_p2[1].push_back(row[10]);
				}
				else{
					// Not a step function -> fill with zeros 
					clover_135_p0[1].push_back(0);
					clover_135_p1[1].push_back(0);
					clover_135_p2[1].push_back(0);
				}
			}
			            
        }
        else {
        }
        i++;
	}
    
    calfile.close();

	clover_array-> SetBranchAddress("amplitude_beam_monitor",amplitude_beam_monitor);
	clover_array-> SetBranchAddress("amplitude_clover_135",amplitude_clover_135);
	clover_array-> SetBranchAddress("amplitude_clover_90",amplitude_clover_90);
	clover_array-> SetBranchAddress("integration_long",integration_long);
	clover_array-> SetBranchAddress("time_clover_135",time_clover_135);
	clover_array-> SetBranchAddress("time_clover_90",time_clover_90);
	clover_array-> SetBranchAddress("channel_time",channel_time);

	TH1D* zero_degree = new TH1D("zero_degree","zero_degree",65536,-0.5,65535.5);
	TH1D* zero_degree_cal = new TH1D("zero_degree_cal","zero_degree_cal",65536,-0.125,16383.875);

	vector <TH1D*> histograms_clover_135;
	vector <TH1D*> histograms_clover_135_cal;
	vector <TH1D*> addback_clover_135;
	vector <TH1D*> sum_clover_135;

	//vector < vector <TH1D*> > time_diff_135;
	vector <TH1D*> histograms_clover_90;
	vector <TH1D*> histograms_clover_90_cal;
	vector <TH1D*> addback_clover_90;
	vector <TH1D*> sum_clover_90;

	//vector < vector <TH1D*> > time_diff_90;

	vector <TH1D*> histograms_cebr;
	vector <TH1D*> histograms_cebr_cal;

	double percent = 0;

	double energy_addback = 0;
	double highest_energy = 0;
	int highest_energy_index = -1;

	// Modules, excluding beam monitors
	vector <string> modules({"clover_90","clover_135","cebr"});
	vector <double> energy_clover_135(16);
	vector <string> clover_135_detnames({"B1", "B4", "B5", "B6"});

	vector <double> energy_clover_90(16);
	vector <string> clover_90_detnames({"1", "3", "5", "7"});

	vector <double> energy_cebr(12);
	vector <string> cebr_detnames({"B", "C", "E", "F", "I", "M", "BB", "BC", "BD", "BK", "BL", "BM"});

	vector <string> detname_lst;
	detname_lst.insert(detname_lst.end(),clover_90_detnames.begin(),clover_90_detnames.end());
	detname_lst.insert(detname_lst.end(),clover_135_detnames.begin(), clover_135_detnames.end());
	detname_lst.insert(detname_lst.end(),cebr_detnames.begin(), cebr_detnames.end());

	vector <double> energy_addback_135(4);
	vector <double> energy_addback_90(4);
	
	vector<int> highest_energy_index_90(4);
	vector<int> highest_energy_index_135(4);

	vector<double> final_energies;
	vector<int> highest_energy_indices;
	vector<double> times;


	//135 degree clover module
	for(int i=0; i<16; i++){
		histograms_clover_135.push_back(new TH1D(("amplitude_clover_135_"+to_string(i)).c_str(),("amplitude_clover_135_"+to_string(i)).c_str(),65536,-0.5,65535.5));
		histograms_clover_135_cal.push_back(new TH1D(("amplitude_clover_135_cal_"+to_string(i)).c_str(),("amplitude_clover_135_cal_"+to_string(i)).c_str(),65536,-0.125,16383.875));
		//time_diff_135.push_back(vector<TH1D*>());
		//for(int j=i+1; j<16; j++){
			//time_diff_135[i].push_back(new TH1D(("time_diff_135_"+to_string(i)+"_"+to_string(j)).c_str(), ("time_diff_135_"+to_string(i)+"_"+to_string(j)).c_str(), 60000, -3000, 3000));
		//}
	}

	for (auto detname : clover_135_detnames){
		addback_clover_135.push_back(new TH1D(("addback_"+detname).c_str(),("addback_"+ detname).c_str(),65536,-0.125,16383.875));
		sum_clover_135.push_back(new TH1D(("sum_"+detname).c_str(),("sum_"+ detname).c_str(),65536,-0.125,16383.875));

	}

	//90 degree clover module
	for(int i=0; i<16; i++){
		histograms_clover_90.push_back(new TH1D(("amplitude_clover_90_"+to_string(i)).c_str(),("amplitude_clover_90_"+to_string(i)).c_str(),65536,-0.5,65535.5));
		histograms_clover_90_cal.push_back(new TH1D(("amplitude_clover_90_cal_"+to_string(i)).c_str(),("amplitude_clover_90_cal_"+to_string(i)).c_str(),65536,-0.125,16383.875));
		// time_diff_90.push_back(vector<TH1D*>());
		// for(int j=i+1; j<16; j++){
		// 	time_diff_90[i].push_back(new TH1D(("time_diff_90_"+to_string(i)+"_"+to_string(j)).c_str(), ("time_diff__90"+to_string(i)+"_"+to_string(j)).c_str(), 60000, -3000, 3000));
		// }
	}

	for (auto detname : clover_90_detnames){
		addback_clover_90.push_back(new TH1D(("addback_"+detname).c_str(),("addback_"+ detname).c_str(),65536,-0.125,16383.875));
		sum_clover_90.push_back(new TH1D(("sum_"+detname).c_str(),("sum_"+ detname).c_str(),65536,-0.125,16383.875));
	}

	//cebr module
	for(int i=0; i<12; i++){
		histograms_cebr.push_back(new TH1D(("integration_long_"+to_string(i)).c_str(),("integration_long_"+to_string(i)).c_str(),4096,-0.5,4095.5));
		histograms_cebr_cal.push_back(new TH1D(("integration_long_cal_"+to_string(i)).c_str(),("integration_long_cal_"+to_string(i)).c_str(),4096,-0.5,4095.5));
	}

	cout << "Finished creating histograms" << endl;
	
	cout << "Filling histograms for " << n_events << " events." << endl;

	auto start = std::chrono::steady_clock::now();

	//for (int i=32089; i<32090; ++i){
	//for(int i=1; i<=10000000; ++i){
	for(int i=1; i<=clover_array->GetEntries(); ++i){
		clover_array->GetEntry(i);

		final_energies.clear();
		highest_energy_indices.clear();
		times.clear();

		//zero degree detector
		zero_degree->Fill(amplitude_beam_monitor[0]);
		zero_degree_cal->Fill(2.792486 + 2.409417e-01*amplitude_beam_monitor[0]);

		//90 degree clover module
		for(int j=0; j<16; j++){
			if (isnan(amplitude_clover_90[j])){
				energy_clover_90[j]=0;
			}
			else if (amplitude_clover_90[j]<=clover_90_step[j]){
				// Quadratic calibration
				// First part of step function, or entire range if no step function
				energy_clover_90.at(j)=clover_90_p0[0][j] + clover_90_p1[0][j]*amplitude_clover_90[j] + clover_90_p2[0][j]*amplitude_clover_90[j]*amplitude_clover_90[j];
			}
			
			else{
				// Quadratic calibration
				// Second part of step function
				energy_clover_90[j]=clover_90_p0[1][j] + clover_90_p1[1][j]*amplitude_clover_90[j] + clover_90_p2[1][j]*amplitude_clover_90[j]*amplitude_clover_90[j];
			}
		}
		for(int det=0; det<4; det++){
			highest_energy = 0;
			highest_energy_index= -1;
			for(int j=4*(det+1)-4; j<4*(det+1); j++){
				if (energy_clover_90[j]>150){  //energy cut in keV
					histograms_clover_90[j]->Fill(amplitude_clover_90[j]);
					histograms_clover_90_cal[j]->Fill(energy_clover_90[j]);
					sum_clover_90[det]->Fill(energy_clover_90[j]);
					if (energy_clover_90[j]>highest_energy){
						highest_energy = energy_clover_90[j];
						highest_energy_index = j;
					}
					// for (int k=j+1; k<4*(det+1); k++){
					// 	//if (! isnan(amplitude_clover_90[k])){
					// 	if (energy_clover_90[k]>15){ //energy cut in keV
					// 		time_diff_90[j][k-j-1]->Fill(time_p1*(time_clover_90[j]-time_clover_90[k]));
					// 	}
					// }	
				}		
			}
			highest_energy_index_90[det]=highest_energy_index;
			energy_addback = 0;
			if (highest_energy_index > -1){
				for(int j=4*(det+1)-4; j<4*(det+1); j++){
					if(fabs(time_p1*(time_clover_90[highest_energy_index]-time_clover_90[j]))<150){ //time cut in ns
						if (energy_clover_90[j]>150){		//energy cut in keV
							energy_addback += energy_clover_90[j];
						}
					}
				}
				addback_clover_90[det]->Fill(energy_addback);
				
			}	
			//energy_addback_90[det]=energy_addback;
			final_energies.push_back(energy_addback);
			highest_energy_indices.push_back(highest_energy_index);
			times.push_back(time_p1*(time_clover_90[highest_energy_index]));
		}
		//final_energies.push_back(energy_addback_90);
		//highest_energy_indices.push_back(highest_energy_index_90);
		//135 degree clover module
		//loop through each channel (in sets of four)
		for(int j=0; j<16; j++){
			if (isnan(amplitude_clover_135[j])){
				energy_clover_135[j]=0;
			}
			else if (amplitude_clover_135[j]<=clover_135_step[j]){
				// Quadratic calibration
				// First part of step function, or entire range if no step function
				energy_clover_135[j]=clover_135_p0[0][j] + clover_135_p1[0][j]*amplitude_clover_135[j] + clover_135_p2[0][j]*amplitude_clover_135[j]*amplitude_clover_135[j];
			}
			
			else{
				// Quadratic calibration
				// Second part of step function
				energy_clover_135[j]=clover_135_p0[1][j] + clover_135_p1[1][j]*amplitude_clover_135[j] + clover_135_p2[1][j]*amplitude_clover_135[j]*amplitude_clover_135[j];
			}
		}
		for(int det=0; det<4; ++det){
			highest_energy = 0;
			highest_energy_index= -1;
			for(int j=4*det; j<4*(det+1); ++j){
				if (energy_clover_135[j]>150){  //energy cut in keV
					histograms_clover_135[j]->Fill(amplitude_clover_135[j]);
					histograms_clover_135_cal[j]->Fill(energy_clover_135[j]);
					sum_clover_135[det]->Fill(energy_clover_135[j]);
					if (energy_clover_135[j]>highest_energy){
						highest_energy = energy_clover_135[j];
						highest_energy_index = j;
					}
					// for (int k=j+1; k<4*(det+1); k++){
					// 	//if (! isnan(amplitude_clover_135[k])){
					// 	if (energy_clover_135[k]>15){ //energy cut in keV
					// 		time_diff_135[j][k-j-1]->Fill(time_p1*(time_clover_135[j]-time_clover_135[k]));
					// 	}
					// }	
				}		
			}
			highest_energy_index_135[det]=highest_energy_index;
			energy_addback = 0;
			if (highest_energy_index > -1){
				for(int j=4*det; j<4*(det+1); ++j){
					if(fabs(time_p1*(time_clover_135[highest_energy_index]-time_clover_135[j]))<150){ //time cut in ns
						if (energy_clover_135[j]>150){		//energy cut in keV
							energy_addback += energy_clover_135[j];
						}
					}
				}
				addback_clover_135[det]->Fill(energy_addback);
			}
			energy_addback_135[det]=energy_addback;
			//cout << "Final addback energy is \t" << energy_addback_135[det] << endl;
			final_energies.push_back(energy_addback);
			highest_energy_indices.push_back(highest_energy_index);
			times.push_back(time_p1*(time_clover_135[highest_energy_index]));

		}

		//final_energies.push_back(energy_addback_135);
		//highest_energy_indices.push_back(highest_energy_index_135);	

		//cebr scintillator module
		for(int j=0; j<12; j++){
			if (isnan(integration_long[j])){
				//index_cebr[j]=-1;
				//energy_cebr[j]=0;
			}
			else{
				//index_cebr[j]=j;
				histograms_cebr[j]->Fill(integration_long[j]);
				histograms_cebr_cal[j]->Fill(cebr_p0[j] + cebr_p1[j]*(5+integration_long[j])); //plus 5 bins for hdtv offset
				//energy_cebr[j]=cebr_p0[j] + cebr_p1[j]*(5+integration_long[j]); //plus 5 bins for hdtv offset
			}
			final_energies.push_back(cebr_p0[j] + cebr_p1[j]*(5+integration_long[j]));
			highest_energy_indices.push_back(j);
			times.push_back(time_p1*(channel_time[j]));

		}
		
		if (i%100000 == 0){
			auto end = std::chrono::steady_clock::now();
			std::chrono::duration<double> diff = end - start;
			percent = float(i)/n_events;
			cout << "\r" << 100*percent << "%           \t";
			cout <<diff.count()*(n_events/float(i) - 1) << "s remaining.           " << flush;
		}
	}
	
	cout << "Finished filling histograms" << endl;

	for(auto hist : addback_clover_90){
		hist->Write();
	}
	for(auto hist : addback_clover_135){
		hist->Write();
	}
	for(auto hist : sum_clover_90){
		hist->Write();
	}
	for(auto hist : sum_clover_135){
		hist->Write();
	}
	for(auto hist : histograms_clover_90_cal){
		hist->Write();
	}
	for(auto hist : histograms_clover_135_cal){
		hist->Write();
	}
	for(auto hist : histograms_clover_90){
		hist->Write();
	}
	for(auto hist : histograms_clover_135){
		hist->Write();
	}
	for(auto hist : histograms_cebr_cal){
		hist->Write();
	}
	for(auto hist : histograms_cebr){
		hist->Write();
	}
	zero_degree->Write();
	zero_degree_cal->Write();
	
	// for(int i=0; i<16; i++){
	// 	for(int j=i+1; j<16; j++){
	// 		time_diff_90[i][j-i-1]->Write();
	// 	}
	// }

	cout << "Finished writing histograms" << endl;

	input_file->Close();
	output_file->Close();
}