/*
 *  makeBkgEstGraphs.C
 *
 *  Created on: Aug 21, 2021
 *      Author: John Russell
 *
 *	With the inclusion of iceMC-only signal responses now included, the work to produce
 *	produce background estimate graphs in Excel-like software now exceeds the work to
 *	produce such graphs through a ROOT macro. This is that macro.
 */


int numPassCutIceMCThermal(const char * responseName, TMatrixD * recallNormConfMat, TChain * classChain, double probLeaked) {

	TMatrixDColumn recallsIceMC(* recallNormConfMat, 0);
	TMatrixDColumn recallsAboveHorizontalThermal(* recallNormConfMat, 1);

	struct {

		double iceMC;
		double aboveHorizontalThermal;

	} responseStruct;

	classChain -> SetBranchAddress(TString::Format("%sCDF", responseName), & responseStruct);

	int numPass = 0;

	for (int i = 0; i < classChain -> GetEntries(); ++i) {

		classChain -> GetEntry(i);

		double arrAct[2] = {responseStruct.iceMC, responseStruct.aboveHorizontalThermal};
		int idxAct = TMath::LocMax(2, arrAct);

		if (idxAct > 0) continue;

		double probPredIceMC = recallsIceMC(0) * responseStruct.iceMC + recallsIceMC(1) * responseStruct.aboveHorizontalThermal;
		double probPredAboveHorizontalThermal = recallsAboveHorizontalThermal(0) * responseStruct.iceMC + recallsAboveHorizontalThermal(1) * responseStruct.aboveHorizontalThermal;

		double arrPred[2] = {probPredIceMC, probPredAboveHorizontalThermal};
		int idxPred = TMath::LocMax(2, arrPred);

		if (idxPred != idxAct) continue;

		//  Referring to the normalized confusion matrix, we want to minimize signal false positives.
		//  What is the probability of actually being any background classification assuming predicted maximum signal classification?
		TMatrixDColumn recallMaxSig(* recallNormConfMat, idxPred);
		double probPredMaxSig = arrPred[idxPred];
		double probActAboveHorizontalThermalWhenPredMaxSig = recallMaxSig(1) * responseStruct.aboveHorizontalThermal / probPredMaxSig;

		if (probActAboveHorizontalThermalWhenPredMaxSig > probLeaked) continue;

		++numPass;
	}

	classChain -> ResetBranchAddresses();

	return numPass;
}


void makeBkgEstIceMCThermalGraph() {

	//  Open file containing relevant confusion matrices, then access them to be passed to other functions.
	TFile confMatFile("../filterEvents/confusionMatricesIceMCThermal.root");

	TMatrixD * PengRecallNormConfMat = (TMatrixD *) confMatFile.Get("Peng/recallNormConfMat");
	TMatrixD * AndrewRecallNormConfMat = (TMatrixD *) confMatFile.Get("Andrew/recallNormConfMat");
	TMatrixD * StokesHybridPlusRecallNormConfMat = (TMatrixD *) confMatFile.Get("StokesHybridPlus/recallNormConfMat");

	//  Create TChains of response values per classification to be passed to other functions.
	TChain survivingChain("responseTree"), RCPChain("responseTree"), LCPChain("responseTree"), payloadBlastChain("responseTree"), sunChain("responseTree"), aboveHorizontalThermalChain("responseTree");

	survivingChain.Add("~/jwrussWork/responseIceMCThermalValues/belowHorizontal/surviving/*.root");
	RCPChain.Add("~/jwrussWork/responseIceMCThermalValues/other/RCP/*.root");
	LCPChain.Add("~/jwrussWork/responseIceMCThermalValues/other/LCP/*.root");
	payloadBlastChain.Add("~/jwrussWork/responseIceMCThermalValues/other/payloadBlast/*.root");
	sunChain.Add("~/jwrussWork/responseIceMCThermalValues/other/sun/*.root");
	aboveHorizontalThermalChain.Add("~/jwrussWork/responseIceMCThermalValues/aboveHorizontal/thermal/*.root");

	vector<TChain *> classChains(6);
	classChains[0] = & survivingChain;
	classChains[1] = & RCPChain;
	classChains[2] = & LCPChain;
	classChains[3] = & payloadBlastChain;
	classChains[4] = & sunChain;
	classChains[5] = & aboveHorizontalThermalChain;

	//  Create arrays in which to store background estimates.
	double sigmaVals[13];
	double numPeng[6][13], numAndrew[6][13], numStokesHybridPlus[6][13];
	double numPengErr[13], numAndrewErr[13], numStokesHybridPlusErr[13];

	for (int i = 0; i < 13; ++i) {

		sigmaVals[i] = 0.5 * i;

		double probLeaked = ROOT::Math::normal_cdf_c(sigmaVals[i]);

		int numPengSurviving = numPassCutIceMCThermal("Peng", PengRecallNormConfMat, classChains[0], probLeaked);
		int numAndrewSurviving = numPassCutIceMCThermal("Andrew", AndrewRecallNormConfMat, classChains[0], probLeaked);
		int numStokesHybridPlusSurviving = numPassCutIceMCThermal("StokesHybridPlus", StokesHybridPlusRecallNormConfMat, classChains[0], probLeaked);

		numPeng[0][i] = numPengSurviving * probLeaked;
		numAndrew[0][i] = numAndrewSurviving * probLeaked;
		numStokesHybridPlus[0][i] = numStokesHybridPlusSurviving * probLeaked;

		numPengErr[i] = sqrt(numPengSurviving) * probLeaked;
		numAndrewErr[i] = sqrt(numAndrewSurviving) * probLeaked;
		numStokesHybridPlusErr[i] = sqrt(numStokesHybridPlusSurviving) * probLeaked;

		for (int j = 1; j < 6; ++j) {

			numPeng[j][i] = numPassCutIceMCThermal("Peng", PengRecallNormConfMat, classChains[j], probLeaked);
			numAndrew[j][i] = numPassCutIceMCThermal("Andrew", AndrewRecallNormConfMat, classChains[j], probLeaked);
			numStokesHybridPlus[j][i] = numPassCutIceMCThermal("StokesHybridPlus", StokesHybridPlusRecallNormConfMat, classChains[j], probLeaked);
		}
	}

	//  Create TGraphs.
	TGraphErrors PengSurvivingGraph(13, sigmaVals, numPeng[0], 0, numPengErr);
	TGraphErrors AndrewSurvivingGraph(13, sigmaVals, numAndrew[0], 0, numAndrewErr);
	TGraphErrors StokesHybridPlusSurvivingGraph(13, sigmaVals, numStokesHybridPlus[0], 0, numStokesHybridPlusErr);

	PengSurvivingGraph.SetTitle("surviving");
	AndrewSurvivingGraph.SetTitle("surviving");
	StokesHybridPlusSurvivingGraph.SetTitle("surviving");

	vector<TGraph> PengBkgGraphs(5), AndrewBkgGraphs(5), StokesHybridPlusBkgGraphs(5);

	vector<TString> bkgNames = {"RCP", "LCP", "payload blast", "sun", "above horizontal thermal"};

	for (int i = 0; i < 5; ++i) {

		PengBkgGraphs[i] = TGraph(13, sigmaVals, numPeng[i + 1]);
		AndrewBkgGraphs[i] = TGraph(13, sigmaVals, numAndrew[i + 1]);
		StokesHybridPlusBkgGraphs[i] = TGraph(13, sigmaVals, numStokesHybridPlus[i + 1]);

		PengBkgGraphs[i].SetTitle(bkgNames[i]);
		AndrewBkgGraphs[i].SetTitle(bkgNames[i]);
		StokesHybridPlusBkgGraphs[i].SetTitle(bkgNames[i]);
	}

	TMultiGraph PengGraphs("PengGraphs", "Peng; Rejection #sigma; Estimated background size");
	TMultiGraph AndrewGraphs("AndrewGraphs", "Andrew; Rejection #sigma; Estimated background size");
	TMultiGraph StokesHybridPlusGraphs("StokesHybridPlusGraphs", "Stokes Hybrid Plus; Rejection #sigma; Estimated background size");

	PengGraphs.Add(& PengSurvivingGraph);
	AndrewGraphs.Add(& AndrewSurvivingGraph);
	StokesHybridPlusGraphs.Add(& StokesHybridPlusSurvivingGraph);

	for (int i = 0; i < 5; ++i) {

		PengGraphs.Add(& PengBkgGraphs[i]);
		AndrewGraphs.Add(& AndrewBkgGraphs[i]);
		StokesHybridPlusGraphs.Add(& StokesHybridPlusBkgGraphs[i]);
	}

	//  Save images of TMultiGraph objects.
	TCanvas PengCanvas("PengCanvas", "", 1);
	PengCanvas.cd();
	PengGraphs.Draw("A L PLC goff");
	gPad -> SetLogy();
	PengGraphs.GetXaxis() -> SetRangeUser(0, 6);
	PengGraphs.GetYaxis() -> SetRangeUser(1e-6, 1e6);
	PengCanvas.BuildLegend();
	PengCanvas.SaveAs("PengIceMCThermal.png");

	TCanvas AndrewCanvas("AndrewCanvas", "", 1);
	AndrewCanvas.cd();
	AndrewGraphs.Draw("A L PLC goff");
	gPad -> SetLogy();
	AndrewGraphs.GetXaxis() -> SetRangeUser(0, 6);
	AndrewGraphs.GetYaxis() -> SetRangeUser(1e-6, 1e6);
	AndrewCanvas.BuildLegend();
	AndrewCanvas.SaveAs("AndrewIceMCThermal.png");

	TCanvas StokesHybridPlusCanvas("StokesHybridPlusCanvas", "", 1);
	StokesHybridPlusCanvas.cd();
	StokesHybridPlusGraphs.Draw("A L PLC goff");
	gPad -> SetLogy();
	StokesHybridPlusGraphs.GetXaxis() -> SetRangeUser(0, 6);
	StokesHybridPlusGraphs.GetYaxis() -> SetRangeUser(1e-6, 1e6);
	StokesHybridPlusCanvas.BuildLegend();
	StokesHybridPlusCanvas.SaveAs("StokesHybridPlusIceMCThermal.png");

	//  Close the confusion matrix file.
	confMatFile.Close();
}





//void makeBkgEstIceMCGraph() {
//
//	//  Open file containing relevant confusion matrices, then access them to be passed to other functions.
//	TFile confMatFile("../filterEvents/confusionMatricesIceMC.root");
//
//	TMatrixD * PengRecallNormConfMat = (TMatrixD *) confMatFile.Get("Peng/recallNormConfMat");
//	TMatrixD * AndrewRecallNormConfMat = (TMatrixD *) confMatFile.Get("Andrew/recallNormConfMat");
//	TMatrixD * StokesHybridPlusRecallNormConfMat = (TMatrixD *) confMatFile.Get("StokesHybridPlus/recallNormConfMat");
//
//	//  Create corresponding vector of these objects.
//	vector<TMatrixD *> recallNormConfMats(3);
//	recallNormConfMats[0] = PengRecallNormConfMat;
//	recallNormConfMats[1] = AndrewRecallNormConfMat;
//	recallNormConfMats[2] = StokesHybridPlusRecallNormConfMat;
//
//	//  Create TChains of response values per classification to be passed to other functions.
//	TChain RCPChain("responseTree"), LCPChain("responseTree"), payloadBlastChain("responseTree"), sunChain("responseTree"), aboveHorizontalThermalChain("responseTree"), survivingChain("responseTree");
//
//	RCPChain.Add("~/jwrussWork/responseIceMCValues/other/RCP/*.root");
//	LCPChain.Add("~/jwrussWork/responseIceMCValues/other/LCP/*.root");
//	payloadBlastChain.Add("~/jwrussWork/responseIceMCValues/other/payloadBlast/*.root");
//	sunChain.Add("~/jwrussWork/responseIceMCValues/other/sun/*.root");
//	aboveHorizontalThermalChain.Add("~/jwrussWork/responseIceMCValues/aboveHorizontal/thermal/*.root");
//	survivingChain.Add("~/jwrussWork/responseIceMCValues/belowHorizontal/surviving/*.root");
//
//	vector<TChain *> classChains(6);
//	classChains[0] = & RCPChain;
//	classChains[1] = & LCPChain;
//	classChains[2] = & payloadBlastChain;
//	classChains[3] = & sunChain;
//	classChains[4] = & aboveHorizontalThermalChain;
//	classChains[5] = & survivingChain;
//
//	//  Point to relevant responses from TChains.
//	struct responsesStruct {
//
//		double iceMC;
//		double RCP;
//		double LCP;
//		double payloadBlast
//		double sun;
//		double aboveHorizontalThermal;
//	};
//
//	vector<responsesStruct *> responsesPeng(6), responsesAndrew(6), responsesStokesHybridPlus(6);
//
//	//  Point to relevant responses from TChains.
//	for (int i = 0; i < 6; ++i) {
//
//		classChains[i] -> SetBranchAddress("PengCDF", responsesPeng[i]);
//		classChains[i] -> SetBranchAddress("AndrewCDF", responsesAndrew[i]);
//		classChains[i] -> SetBranchAddress("StokesHybridPlusCDF", responsesStokesHybridPlus[i]);
//	}
//
//
//
//	//  Close the confusion matrix file.
//	confMatFile.Close();
//}
//
//
//void makeBkgEstThermalGraph() {
//
//	//  Open file containing relevant confusion matrices, then access them to be passed to other functions.
//	TFile confMatFile("../filterEvents/confusionMatricesThermal.root");
//
//	TMatrixD * PengRecallNormConfMat = (TMatrixD *) confMatFile.Get("Peng/recallNormConfMat");
//	TMatrixD * AndrewRecallNormConfMat = (TMatrixD *) confMatFile.Get("Andrew/recallNormConfMat");
//	TMatrixD * StokesHybridPlusRecallNormConfMat = (TMatrixD *) confMatFile.Get("StokesHybridPlus/recallNormConfMat");
//
//	//  Create corresponding vector of these objects.
//	vector<TMatrixD *> recallNormConfMats(3);
//	recallNormConfMats[0] = PengRecallNormConfMat;
//	recallNormConfMats[1] = AndrewRecallNormConfMat;
//	recallNormConfMats[2] = StokesHybridPlusRecallNormConfMat;
//
//	//  Create TChains of response values per classification to be passed to other functions.
//	TChain RCPChain("responseTree"), LCPChain("responseTree"), payloadBlastChain("responseTree"), sunChain("responseTree"), aboveHorizontalThermalChain("responseTree"), survivingChain("responseTree");
//
//	RCPChain.Add("~/jwrussWork/responseThermalValues/other/RCP/*.root");
//	LCPChain.Add("~/jwrussWork/responseThermalValues/other/LCP/*.root");
//	payloadBlastChain.Add("~/jwrussWork/responseThermalValues/other/payloadBlast/*.root");
//	sunChain.Add("~/jwrussWork/responseThermalValues/other/sun/*.root");
//	aboveHorizontalThermalChain.Add("~/jwrussWork/responseThermalValues/aboveHorizontal/thermal/*.root");
//	survivingChain.Add("~/jwrussWork/responseThermalValues/belowHorizontal/surviving/*.root");
//
//	vector<TChain *> classChains(6);
//	classChains[0] = & RCPChain;
//	classChains[1] = & LCPChain;
//	classChains[2] = & payloadBlastChain;
//	classChains[3] = & sunChain;
//	classChains[4] = & aboveHorizontalThermalChain;
//	classChains[5] = & survivingChain;
//
//	//  Point to relevant responses from TChains.
//	struct responsesStruct {
//
//		double WAISHPol;
//		double WAISVPol;
//		double HiCal2A;
//		double HiCal2B;
//		double iceMC;
//		double aboveHorizontalThermal;
//	};
//
//	vector<responsesStruct *> responsesPeng(6), responsesAndrew(6), responsesStokesHybridPlus(6);
//
//	//  Point to relevant responses from TChains.
//	for (int i = 0; i < 6; ++i) {
//
//		classChains[i] -> SetBranchAddress("PengCDF", responsesPeng[i]);
//		classChains[i] -> SetBranchAddress("AndrewCDF", responsesAndrew[i]);
//		classChains[i] -> SetBranchAddress("StokesHybridPlusCDF", responsesStokesHybridPlus[i]);
//	}
//
//
//
//	//  Close the confusion matrix file.
//	confMatFile.Close();
//}
//
//
//void makeBkgEstGraph() {
//
//	//  Open file containing relevant confusion matrices, then access them to be passed to other functions.
//	TFile confMatFile("../filterEvents/confusionMatrices.root");
//
//	TMatrixD * PengRecallNormConfMat = (TMatrixD *) confMatFile.Get("Peng/recallNormConfMat");
//	TMatrixD * AndrewRecallNormConfMat = (TMatrixD *) confMatFile.Get("Andrew/recallNormConfMat");
//	TMatrixD * StokesHybridPlusRecallNormConfMat = (TMatrixD *) confMatFile.Get("StokesHybridPlus/recallNormConfMat");
//
//	//  Create corresponding vector of these objects.
//	vector<TMatrixD *> recallNormConfMats(3);
//	recallNormConfMats[0] = PengRecallNormConfMat;
//	recallNormConfMats[1] = AndrewRecallNormConfMat;
//	recallNormConfMats[2] = StokesHybridPlusRecallNormConfMat;
//
//	//  Create TChains of response values per classification to be passed to other functions.
//	TChain RCPChain("responseTree"), LCPChain("responseTree"), payloadBlastChain("responseTree"), sunChain("responseTree"), aboveHorizontalThermalChain("responseTree"), survivingChain("responseTree");
//
//	RCPChain.Add("~/jwrussWork/responseValues/other/RCP/*.root");
//	LCPChain.Add("~/jwrussWork/responseValues/other/LCP/*.root");
//	payloadBlastChain.Add("~/jwrussWork/responseValues/other/payloadBlast/*.root");
//	sunChain.Add("~/jwrussWork/responseValues/other/sun/*.root");
//	aboveHorizontalThermalChain.Add("~/jwrussWork/responseValues/aboveHorizontal/thermal/*.root");
//	survivingChain.Add("~/jwrussWork/responseValues/belowHorizontal/surviving/*.root");
//
//	vector<TChain *> classChains(6);
//	classChains[0] = & RCPChain;
//	classChains[1] = & LCPChain;
//	classChains[2] = & payloadBlastChain;
//	classChains[3] = & sunChain;
//	classChains[4] = & aboveHorizontalThermalChain;
//	classChains[5] = & survivingChain;
//
//	//  Point to relevant responses from TChains.
//	struct responsesStruct {
//
//		double WAISHPol;
//		double WAISVPol;
//		double HiCal2A;
//		double HiCal2B;
//		double iceMC;
//		double RCP;
//		double LCP;
//		double payloadBlast;
//		double sun;
//		double aboveHorizontalThermal;
//		double signalLike;
//	};
//
//	vector<responsesStruct *> responsesPeng[6], responsesAndrew[6], responsesStokesHybridPlus[6];
//
//	//  Point to relevant responses from TChains.
//	for (int i = 0; i < 6; ++i) {
//
//		classChains[i] -> SetBranchAddress("PengCDF", responsesPeng[i]);
//		classChains[i] -> SetBranchAddress("AndrewCDF", responsesAndrew[i]);
//		classChains[i] -> SetBranchAddress("StokesHybridPlusCDF", responsesStokesHybridPlus[i]);
//	}
//
//
//
//	//  Close the confusion matrix file.
//	confMatFile.Close();
//}


void makeBkgEstGraphs() {

	makeBkgEstIceMCThermalGraph();
//	makeBkgEstIceMCGraph();
//	makeBkgEstThermalGraph();
//	makeBkgEstGraph();
}