/*  General function which generates a file of cut values for a given run and sample data type.
 *  Started 03/09/18 by John Russell. Update on 02/27/18 to account for proper
 *  wrapping of deviations in phi. Edited 05/20/18 to include WAIS VPol template.
 */


#include "FFTtools.h"
//#include "FFTtoolsRev.h"
#include "AnitaDataset.h"
//#include "UCMapStatTools.h"
#include "Hical2.h"
#include "qualityCutLib.h"


/*  Structure corresponding to selected cut data for signal waveforms. Two structures for coherently summed and deconvolved waveforms.
 *  Each entry is 2-dimensional, with the first entry corresponding to horizontal polarization, and the second vertical.
 */
struct {

	double maxWaveform[2], maxWaveformTime[2];
	double minWaveform[2], minWaveformTime[2];
	double waveformRMS[2], waveformZScore[2];
	double waveformPeak[2], waveformPeakTime[2], waveformPeakZScore[2];
	double analyticPeak[2], analyticPeakTime[2], analyticPeakZScore[2];
	double snr[2];
	double linearPolFrac[2], linearPolAngle[2];
	double totalPolFrac[2];
	double instantaneousLinearPolFrac[2], instantaneousLinearPolAngle[2];
	double instantaneousTotalPolFrac[2];
	double impulsivity[2];
	double autocorr5nsMean[2];
	double autocorrWaveform5nsRMS[2], autocorrWaveform5nsValZScore[2];
	double autocorrAnalytic5nsRMS[2], autocorrAnalytic5nsValZScore[2];
	double autocorrWaveformFullRMS[2], autocorrWaveformFullValZScore[2];
	double autocorrAnalyticFullRMS[2], autocorrAnalyticFullValZScore[2];
	double TIMdB[2];
	double maxCorr[2], maxCorrTime[2];
	double minCorr[2], minCorrTime[2];
	double corrRMS[2], corrValZScore[2];
	double corrPeak[2], corrPeakTime[2], corrPeakZScore[2];
	double analyticCorrPeak[2], analyticCorrPeakTime[2], analyticCorrPeakZScore[2];
} coherent, deconvolved;

struct {

	double mapPeak[2];
	double rectMapRMS[2], sphMapRMS[2];
	double rectMapSNR[2], sphMapSNR[2];
	double rectMapPeakZScore[2], sphMapPeakZScore[2];
	double mapPeakNegTheta[2];
} mapStruct;


void doCutFileA4(int runNum = 119, int sampleType = 0) {
	
	//  To point towards the correct event files.
	TString ANITA4_ROOT_DATA = gSystem -> Getenv("ANITA4_ROOT_DATA");
	gSystem -> Setenv("ANITA_ROOT_DATA", ANITA4_ROOT_DATA);
	
	// Butterworth low-pass filter in case interpolation introduces unphysical high frequency terms.
	double deltaT = 0.1;  //  Step size in resampling.
	FFTtools::ButterworthFilter but(FFTtools::LOWPASS, 4, 2 * deltaT * 0.65);  //  650 MHz 3dB point for first Butterworth low-pass filter in terms of Nyquist frequency multiplier.

	//  Get pointer addresses to WAIS template TGraph objects.
	TFile templateFile("polCutFiles/waisPolCutFilesA4Update/interpWaisTemplatesA4Update.root");
	TGraph * cohTemplateGraphH = (TGraph *) templateFile.Get("waisHPolTemplateA4_1");
	TGraph * cohTemplateGraphV = (TGraph *) templateFile.Get("waisVPolTemplateA4_1");
	TGraph * deconvTemplateGraphH = (TGraph *) templateFile.Get("deconvolvedWaisHPolTemplateA4_1");
	TGraph * deconvTemplateGraphV = (TGraph *) templateFile.Get("deconvolvedWaisVPolTemplateA4_1");

	//  Open the headFile corresponding to run so that we can iterate over the corresponding run numbers.
	TString sampleName;
	if (sampleType == 0) sampleName = "wais";
	else if (sampleType == 1) sampleName = "decimatedThermal";
	else if (sampleType == 2) sampleName = "decimatedSignal";
	else if (sampleType == 3) sampleName = "minBias";
	else return;

	TString sampleFileName;
	if (sampleType != 2) sampleFileName = "runFiles/" + sampleName + "/" + TString::Itoa(runNum, 10) + "_sinsub_10_3_ad_2.root";  //  Change this later, both signal and thermal use the same decimated data.
	else sampleFileName = "runFiles/decimatedThermal/" + TString::Itoa(runNum, 10) + "_sinsub_10_3_ad_2.root";

	TFile sampleFile(sampleFileName);
	TTree * sampleTree;
	if (sampleType != 2) sampleTree = (TTree *) sampleFile.Get(sampleName);  //  Change this later, both signal and thermal use the same decimated data.
	else sampleTree = (TTree *) sampleFile.Get("decimatedThermal");
	AnitaEventSummary * sampleSummary = 0;
	sampleTree -> SetBranchAddress("summary", & sampleSummary);
	int & run = sampleSummary -> run;
	unsigned int & eventNumber = sampleSummary -> eventNumber;
	
	//  Create TFile of polarization cuts.
	TFile cutFile("cutFiles/" + sampleName + "CutFiles/" + sampleName + "CutFileRun" + TString::Itoa(runNum, 10) + ".root", "recreate");

	//  Tree storing all the cut results.
	TTree cutTree("cutTree", "TTree storing all cut results.");

	cutTree.Branch("run", & run);
	cutTree.Branch("eventNumber", & eventNumber);
	
	TString waveformCutStruct = "maxWaveform[2]/D:maxWaveformTime[2]";
	waveformCutStruct += ":minWaveform[2]:minWaveformTime[2]";
	waveformCutStruct += ":waveformRMS[2]:waveformZScore[2]";
	waveformCutStruct += ":waveformPeak[2]:waveformPeakTime[2]:waveformPeakZScore[2]";
	waveformCutStruct += ":analyticPeak[2]:analyticPeakTime[2]:analyticPeakZScore[2]";
	waveformCutStruct += ":snr[2]";
	waveformCutStruct += ":linearPolFrac[2]:linearPolAngle[2]";
	waveformCutStruct += ":totalPolFrac[2]";
	waveformCutStruct += ":instantaneousLinearPolFrac[2]:instantaneousLinearPolAngle[2]";
	waveformCutStruct += ":instantaneousTotalPolFrac[2]";
	waveformCutStruct += ":impulsivity[2]";
	waveformCutStruct += ":autocorr5nsMean[2]";
	waveformCutStruct += ":autocorrWaveform5nsRMS[2]:autocorrWaveform5nsValZScore[2]";
	waveformCutStruct += ":autocorrAnalytic5nsRMS[2]:autocorrAnalytic5nsValZScore[2]";
	waveformCutStruct += ":autocorrWaveformFullRMS[2]:autocorrWaveformFullValZScore[2]";
	waveformCutStruct += ":autocorrAnalyticFullRMS[2]:autocorrAnalyticFullValZScore[2]";
	waveformCutStruct += ":TIMdB[2]";
	waveformCutStruct += ":maxCorr[2]:maxCorrTime[2]";
	waveformCutStruct += ":minCorr[2]:minCorrTime[2]";
	waveformCutStruct += ":corrRMS[2]:corrValZScore[2]";
	waveformCutStruct += ":corrPeak[2]:corrPeakTime[2]:corrPeakZScore[2]";
	waveformCutStruct += ":analyticCorrPeak[2]:analyticCorrPeakTime[2]:analyticCorrPeakZScore[2]";

	cutTree.Branch("cohWaveformCuts", & coherent, waveformCutStruct);
	cutTree.Branch("deconvWaveformCuts", & deconvolved, waveformCutStruct);

	TString mapCutStruct = "mapPeak[2]/D";
	mapCutStruct += ":rectMapRMS[2]:sphMapRMS[2]";
	mapCutStruct += ":rectMapSNR[2]:sphMapSNR[2]";
	mapCutStruct += ":rectMapPeakZScore[2]:sphMapPeakZScore[2]";
	mapCutStruct += ":mapPeakNegTheta[2]";

	cutTree.Branch("mapCuts", & mapStruct, mapCutStruct);

	//  Set up AnitaDataset object for given runNum.
	AnitaDataset d(runNum);
	
	//  UCorrelator stuff. Set up Analyzer objects.
	UCorrelator::AnalysisConfig cfg;
	cfg.deconvolution_method = new AnitaResponse::AllPassDeconvolution;  //  Rather than use the band-limited default, use all-pass to account for lack of elephant trunk corrections.
	cfg.response_option = UCorrelator::AnalysisConfig::ResponseA4;  //  The new A4 response.
	UCorrelator::Analyzer analyzer(& cfg, true);  //  "true" is set so "analyzer" doesn't reset each time.
	FilterStrategy * fStrat = UCorrelator::getStrategyWithKey("sinsub_10_3_ad_2");  //  Basic sine subtraction filter strategy.

	//  Use the sine subtract cache to speed up processing!
//	UCorrelator::SineSubtractFilter::setUseCache(true);

	//  Fill polarityTestsTree with data.
//	int nEntries = 100;
	int nEntries = sampleTree -> GetEntries();
	for (int entryNum = 0; entryNum < nEntries; ++entryNum) {
		
		sampleTree -> GetEntry(entryNum);

		//  Make quality cuts. General.
		if (sampleType != 1 && !sampleSummary -> flags.isGood) continue;  //  Already assumed for WAIS. Only exception are minBias events.
		if (sampleSummary -> flags.isPayloadBlast) continue;
		if (sampleSummary -> flags.isStepFunction) continue;
		if (sampleSummary -> flags.hasGlitch) continue;
		if (Hical2::isHical(sampleSummary)) continue;  //  Should only apply to thermal, because this is hard coded for known runs where HiCal was present.
		//  Specific to WAIS.
		if (sampleType == 0 && (run >= 129 && run <= 132)) continue;
		int triggeredSum = 0;
		for (int i = 0; i < 2; ++i) {

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

				if (sampleSummary -> peak[i][j].triggered) ++triggeredSum;
			}
		}
		if ((sampleType == 0 || sampleType == 2) && !triggeredSum) continue;
		double waisPhiDispH = FFTtools::wrap(sampleSummary -> peak[0][0].phi - sampleSummary -> wais.phi, 360, 0);
		double waisPhiDispV = FFTtools::wrap(sampleSummary -> peak[1][0].phi - sampleSummary -> wais.phi, 360, 0);
		double waisThetaDispH = sampleSummary -> peak[0][0].theta - sampleSummary -> wais.theta;
		double waisThetaDispV = sampleSummary -> peak[1][0].theta - sampleSummary -> wais.theta;
		if (sampleType == 0 && pow(waisPhiDispH / 3, 2) + pow(waisThetaDispH, 2) > 1 && pow(waisPhiDispV / 3, 2) + pow(waisThetaDispV, 2) > 1) continue;  //  Check pointing on largest peak in HPol.
		//  Specific to thermal.
		if (sampleType == 1 && sampleSummary -> peak[0][0].theta > 0) continue;
		if (sampleType == 1 && sampleSummary -> peak[1][0].theta > 0) continue;
		if ((sampleType == 1 || sampleType == 2 || sampleType == 3) && (sampleSummary -> flags.maxBottomToTopRatio[0] <= 1 || sampleSummary -> flags.maxBottomToTopRatio[0] >= 2.6)) continue;  //  Check HPol ratio.
		if ((sampleType == 1 || sampleType == 2 || sampleType == 3) && (sampleSummary -> flags.maxBottomToTopRatio[1] <= 1 || sampleSummary -> flags.maxBottomToTopRatio[1] >= 2.6)) continue;  //  Check VPol ratio.
//		int maskedSum = 0, masked_xpolSum = 0;
//		for (int i = 0; i < 2; ++i) {
//
//			for (int j = 0; j < 5; ++j) {  //  For no masked events.
//
//				if (sampleSummary -> peak[i][j].masked) ++maskedSum;
//				if (sampleSummary -> peak[i][j].masked_xpol) ++masked_xpolSum;
//			}
//		}
//		if ((sampleType == 1 || sampleType == 2) && (maskedSum || masked_xpolSum)) continue;
		//  Specific to signal.
		if (sampleType == 2 && sampleSummary -> peak[0][0].theta <= 0) continue;
		if (sampleType == 2 && sampleSummary -> peak[1][0].theta <= 0) continue;

		//  Use the analyzer declared above to acquire coherently summed waveform for given run and event number.
		d.getEvent(eventNumber);  //  In the run, pointing to the correct eventNumber.
		FilteredAnitaEvent filtEvent(d.useful(), fStrat, d.gps(), d.header());

		analyzer.analyze(& filtEvent, sampleSummary);

		//  Objects with which to calculate quantities to fill TTrees with.
		TGraph * cohWaveformH = FFTtools::getInterpolatedGraph(analyzer.getCoherent(AnitaPol::kHorizontal, 0) -> even(), deltaT);
		but.filterGraph(cohWaveformH);

		TGraph * cohWaveformV = FFTtools::getInterpolatedGraph(analyzer.getCoherent(AnitaPol::kVertical, 0) -> even(), deltaT);
		but.filterGraph(cohWaveformV);

		TGraph * deconvWaveformH = FFTtools::getInterpolatedGraph(analyzer.getDeconvolved(AnitaPol::kHorizontal, 0) -> even(), deltaT);
		but.filterGraph(deconvWaveformH);

		TGraph * deconvWaveformV = FFTtools::getInterpolatedGraph(analyzer.getDeconvolved(AnitaPol::kVertical, 0) -> even(), deltaT);
		but.filterGraph(deconvWaveformV);

		//  Create correlation graphs, checking that they're created correctly and regenerating them if they aren't.
		TGraph * cohAutocorrH = FFTtools::getCorrGraph(cohWaveformH, cohWaveformH);
		while (!isfinite(cohAutocorrH -> GetHistogram() -> GetMaximum())) cohAutocorrH = FFTtools::getCorrGraph(cohWaveformH, cohWaveformH);

		TGraph * cohAutocorrV = FFTtools::getCorrGraph(cohWaveformV, cohWaveformV);
		while (!isfinite(cohAutocorrV -> GetHistogram() -> GetMaximum())) cohAutocorrV = FFTtools::getCorrGraph(cohWaveformV, cohWaveformV);

		TGraph * cohTemplateCorrH = FFTtools::getCorrGraph(cohTemplateGraphH, cohWaveformH);
		while (!isfinite(max(abs(cohTemplateCorrH -> GetHistogram() -> GetMaximum()), abs(cohTemplateCorrH -> GetHistogram() -> GetMinimum())))) cohTemplateCorrH = FFTtools::getCorrGraph(cohTemplateGraphH, cohWaveformH);

		TGraph * cohTemplateCorrV = FFTtools::getCorrGraph(cohTemplateGraphV, cohWaveformV);
		while (!isfinite(max(abs(cohTemplateCorrV -> GetHistogram() -> GetMaximum()), abs(cohTemplateCorrV -> GetHistogram() -> GetMinimum())))) cohTemplateCorrV = FFTtools::getCorrGraph(cohTemplateGraphV, cohWaveformV);

		TGraph * cohTemplateAnalyticCorrAbsH = getAnalyticCorrAbs(cohTemplateGraphH, cohWaveformH);
		while (!isfinite(max(abs(cohTemplateAnalyticCorrAbsH -> GetHistogram() -> GetMaximum()), abs(cohTemplateAnalyticCorrAbsH -> GetHistogram() -> GetMinimum())))) cohTemplateAnalyticCorrAbsH = getAnalyticCorrAbs(cohTemplateGraphH, cohWaveformH);

		TGraph * cohTemplateAnalyticCorrAbsV = getAnalyticCorrAbs(cohTemplateGraphV, cohWaveformV);
		while (!isfinite(max(abs(cohTemplateAnalyticCorrAbsV -> GetHistogram() -> GetMaximum()), abs(cohTemplateAnalyticCorrAbsV -> GetHistogram() -> GetMinimum())))) cohTemplateAnalyticCorrAbsV = getAnalyticCorrAbs(cohTemplateGraphV, cohWaveformV);

		TGraph * deconvAutocorrH = FFTtools::getCorrGraph(deconvWaveformH, deconvWaveformH);
		while (!isfinite(deconvAutocorrH -> GetHistogram() -> GetMaximum())) deconvAutocorrH = FFTtools::getCorrGraph(deconvWaveformH, deconvWaveformH);

		TGraph * deconvAutocorrV = FFTtools::getCorrGraph(deconvWaveformV, deconvWaveformV);
		while (!isfinite(deconvAutocorrV -> GetHistogram() -> GetMaximum())) deconvAutocorrV = FFTtools::getCorrGraph(deconvWaveformV, deconvWaveformV);

		TGraph * deconvTemplateCorrH = FFTtools::getCorrGraph(deconvTemplateGraphH, deconvWaveformH);
		while (!isfinite(max(abs(deconvTemplateCorrH -> GetHistogram() -> GetMaximum()), abs(deconvTemplateCorrH -> GetHistogram() -> GetMinimum())))) deconvTemplateCorrH = getAnalyticCorrAbs(deconvTemplateGraphH, deconvWaveformH);

		TGraph * deconvTemplateCorrV = FFTtools::getCorrGraph(deconvTemplateGraphV, deconvWaveformV);
		while (!isfinite(max(abs(deconvTemplateCorrV -> GetHistogram() -> GetMaximum()), abs(deconvTemplateCorrV -> GetHistogram() -> GetMinimum())))) deconvTemplateCorrV = getAnalyticCorrAbs(deconvTemplateGraphV, deconvWaveformV);

		TGraph * deconvTemplateAnalyticCorrAbsH = getAnalyticCorrAbs(deconvTemplateGraphH, deconvWaveformH);
		while (!isfinite(max(abs(deconvTemplateAnalyticCorrAbsH -> GetHistogram() -> GetMaximum()), abs(deconvTemplateAnalyticCorrAbsH -> GetHistogram() -> GetMinimum())))) deconvTemplateAnalyticCorrAbsH = getAnalyticCorrAbs(deconvTemplateGraphH, deconvWaveformH);

		TGraph * deconvTemplateAnalyticCorrAbsV = getAnalyticCorrAbs(deconvTemplateGraphV, deconvWaveformV);
		while (!isfinite(max(abs(deconvTemplateAnalyticCorrAbsV -> GetHistogram() -> GetMaximum()), abs(deconvTemplateAnalyticCorrAbsV -> GetHistogram() -> GetMinimum())))) deconvTemplateAnalyticCorrAbsV = getAnalyticCorrAbs(deconvTemplateGraphV, deconvWaveformV);

		TH2D * corrMapH = new TH2D(* analyzer.getCorrelationMap(AnitaPol::kHorizontal));
		TH2D * corrMapV = new TH2D(* analyzer.getCorrelationMap(AnitaPol::kVertical));

		//  Fill branches corresponding to waveform values.
		int NCohH = cohWaveformH -> GetN();
		double * XCohH = cohWaveformH -> GetX();
		double * YCohH = cohWaveformH -> GetY();
		int maxIndCohH = TMath::LocMax(NCohH, YCohH);
		int minIndCohH = TMath::LocMin(NCohH, YCohH);

		int NCohV = cohWaveformV -> GetN();
		double * XCohV = cohWaveformV -> GetX();
		double * YCohV = cohWaveformV -> GetY();
		int maxIndCohV = TMath::LocMax(NCohV, YCohV);
		int minIndCohV = TMath::LocMin(NCohV, YCohV);

		int NCohCorrH = cohTemplateCorrH -> GetN();
		double * XCohCorrH = cohTemplateCorrH -> GetX();
		double * YCohCorrH = cohTemplateCorrH -> GetY();
		int maxIndCohCorrH = TMath::LocMax(NCohCorrH, YCohCorrH);
		int minIndCohCorrH = TMath::LocMin(NCohCorrH, YCohCorrH);

		int NCohCorrV = cohTemplateCorrV -> GetN();
		double * XCohCorrV = cohTemplateCorrV -> GetX();
		double * YCohCorrV = cohTemplateCorrV -> GetY();
		int maxIndCohCorrV = TMath::LocMax(NCohCorrV, YCohCorrV);
		int minIndCohCorrV = TMath::LocMin(NCohCorrV, YCohCorrV);

		int NCohAnalyticCorrAbsH = cohTemplateAnalyticCorrAbsH -> GetN();
		double * XCohAnalyticCorrAbsH = cohTemplateAnalyticCorrAbsH -> GetX();
		double * YCohAnalyticCorrAbsH = cohTemplateAnalyticCorrAbsH -> GetY();
		int maxIndCohAnalyticCorrAbsH = TMath::LocMax(NCohAnalyticCorrAbsH, YCohAnalyticCorrAbsH);

		int NCohAnalyticCorrAbsV = cohTemplateAnalyticCorrAbsV -> GetN();
		double * XCohAnalyticCorrAbsV = cohTemplateAnalyticCorrAbsV -> GetX();
		double * YCohAnalyticCorrAbsV = cohTemplateAnalyticCorrAbsV -> GetY();
		int maxIndCohAnalyticCorrAbsV = TMath::LocMax(NCohAnalyticCorrAbsV, YCohAnalyticCorrAbsV);

		coherent.maxWaveform[0] = YCohH[maxIndCohH];
		coherent.maxWaveform[1] = YCohV[maxIndCohV];
		coherent.maxWaveformTime[0] = XCohH[maxIndCohH];
		coherent.maxWaveformTime[1] = XCohV[maxIndCohV];
		coherent.minWaveform[0] = YCohH[minIndCohH];
		coherent.minWaveform[1] = YCohV[minIndCohV];
		coherent.minWaveformTime[0] = XCohH[minIndCohH];
		coherent.minWaveformTime[1] = XCohV[minIndCohV];
		coherent.waveformRMS[0] = cohWaveformH -> GetRMS(2);
		coherent.waveformRMS[1] = cohWaveformV -> GetRMS(2);
		coherent.waveformZScore[0] = abs(coherent.maxWaveform[0]) >= abs(coherent.minWaveform[0]) ? coherent.maxWaveform[0] : coherent.minWaveform[0];
		coherent.waveformZScore[0] /= coherent.waveformRMS[0];
		coherent.waveformZScore[1] = abs(coherent.maxWaveform[1]) >= abs(coherent.minWaveform[1]) ? coherent.maxWaveform[1] : coherent.minWaveform[1];
		coherent.waveformZScore[1] /= coherent.waveformRMS[1];
		coherent.waveformPeak[0] = max(abs(coherent.maxWaveform[0]), abs(coherent.minWaveform[0]));
		coherent.waveformPeak[1] = max(abs(coherent.maxWaveform[1]), abs(coherent.minWaveform[1]));
		coherent.waveformPeakTime[0] = abs(coherent.maxWaveform[0]) >= abs(coherent.minWaveform[0]) ? coherent.maxWaveformTime[0] : coherent.minWaveformTime[0];
		coherent.waveformPeakTime[1] = abs(coherent.maxWaveform[1]) >= abs(coherent.minWaveform[1]) ? coherent.maxWaveformTime[1] : coherent.minWaveformTime[1];
		coherent.waveformPeakZScore[0] = abs(coherent.waveformZScore[0]);
		coherent.waveformPeakZScore[1] = abs(coherent.waveformZScore[1]);
		coherent.analyticPeak[0] = sampleSummary -> coherent[0][0].peakHilbert;
		coherent.analyticPeak[1] = sampleSummary -> coherent[1][0].peakHilbert;
		coherent.analyticPeakTime[0] = sampleSummary -> coherent[0][0].peakTime;
		coherent.analyticPeakTime[1] = sampleSummary -> coherent[1][0].peakTime;
		coherent.analyticPeakZScore[0] = coherent.analyticPeak[0] / (sqrt(2) * coherent.waveformRMS[0]);  //  Scaling of RMS comes from analytic waveform definition and Parseval's theorem.
		coherent.analyticPeakZScore[1] = coherent.analyticPeak[1] / (sqrt(2) * coherent.waveformRMS[1]);
		coherent.snr[0] = sampleSummary -> coherent[0][0].snr;
		coherent.snr[1] = sampleSummary -> coherent[1][0].snr;
		coherent.linearPolFrac[0] = sampleSummary -> coherent[0][0].linearPolFrac();
		coherent.linearPolFrac[1] = sampleSummary -> coherent[1][0].linearPolFrac();
		coherent.linearPolAngle[0] = sampleSummary -> coherent[0][0].linearPolAngle();
		coherent.linearPolAngle[1] = sampleSummary -> coherent[1][0].linearPolAngle();
		coherent.totalPolFrac[0] = sampleSummary -> coherent[0][0].totalPolFrac();
		coherent.totalPolFrac[1] = sampleSummary -> coherent[1][0].totalPolFrac();
		coherent.instantaneousLinearPolFrac[0] = sampleSummary -> coherent[0][0].instantaneousLinearPolFrac();
		coherent.instantaneousLinearPolFrac[1] = sampleSummary -> coherent[1][0].instantaneousLinearPolFrac();
		coherent.instantaneousLinearPolAngle[0] = sampleSummary -> coherent[0][0].instantaneousLinearPolAngle();
		coherent.instantaneousLinearPolAngle[1] = sampleSummary -> coherent[1][0].instantaneousLinearPolAngle();
		coherent.instantaneousTotalPolFrac[0] = sampleSummary -> coherent[0][0].instantaneousTotalPolFrac();
		coherent.instantaneousTotalPolFrac[1] = sampleSummary -> coherent[1][0].instantaneousTotalPolFrac();
		coherent.impulsivity[0] = sampleSummary -> coherent[0][0].impulsivityMeasure;
		coherent.impulsivity[1] = sampleSummary -> coherent[1][0].impulsivityMeasure;
		coherent.autocorr5nsMean[0] = get5nsMean(cohAutocorrH);
		coherent.autocorr5nsMean[1] = get5nsMean(cohAutocorrV);
		coherent.autocorrWaveform5nsRMS[0] = get5nsRMS(cohAutocorrH);
		coherent.autocorrWaveform5nsRMS[1] = get5nsRMS(cohAutocorrV);
		coherent.autocorrWaveform5nsValZScore[0] = (1 - coherent.autocorr5nsMean[0]) / coherent.autocorrWaveform5nsRMS[0];
		coherent.autocorrWaveform5nsValZScore[1] = (1 - coherent.autocorr5nsMean[1]) / coherent.autocorrWaveform5nsRMS[1];
		coherent.autocorrAnalytic5nsRMS[0] =  getAutocorrAnalytic5nsRMS(cohWaveformH);
		coherent.autocorrAnalytic5nsRMS[1] =  getAutocorrAnalytic5nsRMS(cohWaveformV);
		coherent.autocorrAnalytic5nsValZScore[0] = (1 - coherent.autocorr5nsMean[0]) / coherent.autocorrAnalytic5nsRMS[0];
		coherent.autocorrAnalytic5nsValZScore[1] = (1 - coherent.autocorr5nsMean[1]) / coherent.autocorrAnalytic5nsRMS[1];
		coherent.autocorrWaveformFullRMS[0] = cohAutocorrH -> GetRMS(2);
		coherent.autocorrWaveformFullRMS[1] = cohAutocorrV -> GetRMS(2);
		coherent.autocorrWaveformFullValZScore[0] = 1 / coherent.autocorrWaveformFullRMS[0];
		coherent.autocorrWaveformFullValZScore[1] = 1 / coherent.autocorrWaveformFullRMS[1];
		coherent.autocorrAnalyticFullRMS[0] =  sqrt(2) * coherent.autocorrWaveformFullRMS[0];
		coherent.autocorrAnalyticFullRMS[1] =  sqrt(2) * coherent.autocorrWaveformFullRMS[1];
		coherent.autocorrAnalyticFullValZScore[0] = 1 / coherent.autocorrAnalyticFullRMS[0];
		coherent.autocorrAnalyticFullValZScore[1] = 1 / coherent.autocorrAnalyticFullRMS[1];
		coherent.TIMdB[0] = getTIMdB(cohWaveformH);
		coherent.TIMdB[1] = getTIMdB(cohWaveformV);
		coherent.maxCorr[0] = YCohCorrH[maxIndCohCorrH];
		coherent.maxCorr[1] = YCohCorrV[maxIndCohCorrV];
		coherent.maxCorrTime[0] = XCohCorrH[maxIndCohCorrH];
		coherent.maxCorrTime[1] = XCohCorrV[maxIndCohCorrV];
		coherent.minCorr[0] = YCohCorrH[minIndCohCorrH];
		coherent.minCorr[1] = YCohCorrV[minIndCohCorrV];
		coherent.minCorrTime[0] = XCohCorrH[minIndCohCorrH];
		coherent.minCorrTime[1] = XCohCorrV[minIndCohCorrV];
		coherent.corrRMS[0] = cohTemplateCorrH -> GetRMS(2);
		coherent.corrRMS[1] = cohTemplateCorrV -> GetRMS(2);
		coherent.corrValZScore[0] = abs(coherent.maxCorr[0]) >= abs(coherent.minCorr[0]) ? coherent.maxCorr[0] : coherent.minCorr[0];
		coherent.corrValZScore[0] /= coherent.corrRMS[0];
		coherent.corrValZScore[1] = abs(coherent.maxCorr[1]) >= abs(coherent.minCorr[1]) ? coherent.maxCorr[1] : coherent.minCorr[1];
		coherent.corrValZScore[1] /= coherent.corrRMS[1];
		coherent.corrPeak[0] = max(abs(coherent.maxCorr[0]), abs(coherent.minCorr[0]));
		coherent.corrPeak[1] = max(abs(coherent.maxCorr[1]), abs(coherent.minCorr[1]));
		coherent.corrPeakTime[0] = abs(coherent.maxCorr[0]) >= abs(coherent.minCorr[0]) ? coherent.maxCorrTime[0] : coherent.minCorrTime[0];
		coherent.corrPeakTime[1] = abs(coherent.maxCorr[1]) >= abs(coherent.minCorr[1]) ? coherent.maxCorrTime[1] : coherent.minCorrTime[1];
		coherent.corrPeakZScore[0] = abs(coherent.corrValZScore[0]);
		coherent.corrPeakZScore[1] = abs(coherent.corrValZScore[1]);
		coherent.analyticCorrPeak[0] = YCohAnalyticCorrAbsH[maxIndCohAnalyticCorrAbsH];
		coherent.analyticCorrPeak[1] = YCohAnalyticCorrAbsV[maxIndCohAnalyticCorrAbsV];
		coherent.analyticCorrPeakTime[0] = XCohAnalyticCorrAbsH[maxIndCohAnalyticCorrAbsH];
		coherent.analyticCorrPeakTime[1] = XCohAnalyticCorrAbsV[maxIndCohAnalyticCorrAbsV];
		coherent.analyticCorrPeakZScore[0] = coherent.analyticCorrPeak[0] / (sqrt(2) * coherent.corrRMS[0]);
		coherent.analyticCorrPeakZScore[1] = coherent.analyticCorrPeak[1] / (sqrt(2) * coherent.corrRMS[1]);

		int NDeconvH = deconvWaveformH -> GetN();
		double * XDeconvH = deconvWaveformH -> GetX();
		double * YDeconvH = deconvWaveformH -> GetY();
		int maxIndDeconvH = TMath::LocMax(NDeconvH, YDeconvH);
		int minIndDeconvH = TMath::LocMin(NDeconvH, YDeconvH);

		int NDeconvV = deconvWaveformV -> GetN();
		double * XDeconvV = deconvWaveformV -> GetX();
		double * YDeconvV = deconvWaveformV -> GetY();
		int maxIndDeconvV = TMath::LocMax(NCohV, YCohV);
		int minIndDeconvV = TMath::LocMin(NCohV, YCohV);

		int NDeconvCorrH = deconvTemplateCorrH -> GetN();
		double * XDeconvCorrH = deconvTemplateCorrH -> GetX();
		double * YDeconvCorrH = deconvTemplateCorrH -> GetY();
		int maxIndDeconvCorrH = TMath::LocMax(NDeconvCorrH, YDeconvCorrH);
		int minIndDeconvCorrH = TMath::LocMin(NDeconvCorrH, YDeconvCorrH);

		int NDeconvCorrV = deconvTemplateCorrV -> GetN();
		double * XDeconvCorrV = deconvTemplateCorrV -> GetX();
		double * YDeconvCorrV = deconvTemplateCorrV -> GetY();
		int maxIndDeconvCorrV = TMath::LocMax(NDeconvCorrV, YDeconvCorrV);
		int minIndDeconvCorrV = TMath::LocMin(NDeconvCorrV, YDeconvCorrV);

		int NDeconvAnalyticCorrAbsH = deconvTemplateAnalyticCorrAbsH -> GetN();
		double * XDeconvAnalyticCorrAbsH = deconvTemplateAnalyticCorrAbsH -> GetX();
		double * YDeconvAnalyticCorrAbsH = deconvTemplateAnalyticCorrAbsH -> GetY();
		int maxIndDeconvAnalyticCorrAbsH = TMath::LocMax(NDeconvAnalyticCorrAbsH, YDeconvAnalyticCorrAbsH);

		int NDeconvAnalyticCorrAbsV = deconvTemplateAnalyticCorrAbsV -> GetN();
		double * XDeconvAnalyticCorrAbsV = deconvTemplateAnalyticCorrAbsV -> GetX();
		double * YDeconvAnalyticCorrAbsV = deconvTemplateAnalyticCorrAbsV -> GetY();
		int maxIndDeconvAnalyticCorrAbsV = TMath::LocMax(NDeconvAnalyticCorrAbsV, YDeconvAnalyticCorrAbsV);

		deconvolved.maxWaveform[0] = YDeconvH[maxIndDeconvH];
		deconvolved.maxWaveform[1] = YDeconvV[maxIndDeconvV];
		deconvolved.maxWaveformTime[0] = XDeconvH[maxIndDeconvH];
		deconvolved.maxWaveformTime[1] = XDeconvV[maxIndDeconvV];
		deconvolved.minWaveform[0] = YDeconvH[minIndDeconvH];
		deconvolved.minWaveform[1] = YDeconvV[minIndDeconvV];
		deconvolved.minWaveformTime[0] = XDeconvH[minIndDeconvH];
		deconvolved.minWaveformTime[1] = XDeconvV[minIndDeconvV];
		deconvolved.waveformRMS[0] = deconvWaveformH -> GetRMS(2);
		deconvolved.waveformRMS[1] = deconvWaveformV -> GetRMS(2);
		deconvolved.waveformZScore[0] = abs(deconvolved.maxWaveform[0]) >= abs(deconvolved.minWaveform[0]) ? deconvolved.maxWaveform[0] : deconvolved.minWaveform[0];
		deconvolved.waveformZScore[0] /= deconvolved.waveformRMS[0];
		deconvolved.waveformZScore[1] = abs(deconvolved.maxWaveform[1]) >= abs(deconvolved.minWaveform[1]) ? deconvolved.maxWaveform[1] : deconvolved.minWaveform[1];
		deconvolved.waveformZScore[1] /= deconvolved.waveformRMS[1];
		deconvolved.waveformPeak[0] = max(abs(deconvolved.maxWaveform[0]), abs(deconvolved.minWaveform[0]));
		deconvolved.waveformPeak[1] = max(abs(deconvolved.maxWaveform[1]), abs(deconvolved.minWaveform[1]));
		deconvolved.waveformPeakTime[0] = abs(deconvolved.maxWaveform[0]) >= abs(deconvolved.minWaveform[0]) ? deconvolved.maxWaveformTime[0] : deconvolved.minWaveformTime[0];
		deconvolved.waveformPeakTime[1] = abs(deconvolved.maxWaveform[1]) >= abs(deconvolved.minWaveform[1]) ? deconvolved.maxWaveformTime[1] : deconvolved.minWaveformTime[1];
		deconvolved.waveformPeakZScore[0] = abs(deconvolved.waveformZScore[0]);
		deconvolved.waveformPeakZScore[1] = abs(deconvolved.waveformZScore[1]);
		deconvolved.analyticPeak[0] = sampleSummary -> deconvolved[0][0].peakHilbert;
		deconvolved.analyticPeak[1] = sampleSummary -> deconvolved[1][0].peakHilbert;
		deconvolved.analyticPeakTime[0] = sampleSummary -> deconvolved[0][0].peakTime;
		deconvolved.analyticPeakTime[1] = sampleSummary -> deconvolved[1][0].peakTime;
		deconvolved.analyticPeakZScore[0] = deconvolved.analyticPeak[0] / (sqrt(2) * deconvolved.waveformRMS[0]);
		deconvolved.analyticPeakZScore[1] = deconvolved.analyticPeak[1] / (sqrt(2) * deconvolved.waveformRMS[1]);
		deconvolved.snr[0] = sampleSummary -> deconvolved[0][0].snr;
		deconvolved.snr[1] = sampleSummary -> deconvolved[1][0].snr;
		deconvolved.linearPolFrac[0] = sampleSummary -> deconvolved[0][0].linearPolFrac();
		deconvolved.linearPolFrac[1] = sampleSummary -> deconvolved[1][0].linearPolFrac();
		deconvolved.linearPolAngle[0] = sampleSummary -> deconvolved[0][0].linearPolAngle();
		deconvolved.linearPolAngle[1] = sampleSummary -> deconvolved[1][0].linearPolAngle();
		deconvolved.totalPolFrac[0] = sampleSummary -> deconvolved[0][0].totalPolFrac();
		deconvolved.totalPolFrac[1] = sampleSummary -> deconvolved[1][0].totalPolFrac();
		deconvolved.instantaneousLinearPolFrac[0] = sampleSummary -> deconvolved[0][0].instantaneousLinearPolFrac();
		deconvolved.instantaneousLinearPolFrac[1] = sampleSummary -> deconvolved[1][0].instantaneousLinearPolFrac();
		deconvolved.instantaneousLinearPolAngle[0] = sampleSummary -> deconvolved[0][0].instantaneousLinearPolAngle();
		deconvolved.instantaneousLinearPolAngle[1] = sampleSummary -> deconvolved[1][0].instantaneousLinearPolAngle();
		deconvolved.instantaneousTotalPolFrac[0] = sampleSummary -> deconvolved[0][0].instantaneousTotalPolFrac();
		deconvolved.instantaneousTotalPolFrac[1] = sampleSummary -> deconvolved[1][0].instantaneousTotalPolFrac();
		deconvolved.impulsivity[0] = sampleSummary -> deconvolved[0][0].impulsivityMeasure;
		deconvolved.impulsivity[1] = sampleSummary -> deconvolved[1][0].impulsivityMeasure;
		deconvolved.autocorr5nsMean[0] = get5nsMean(deconvAutocorrH);
		deconvolved.autocorr5nsMean[1] = get5nsMean(deconvAutocorrV);
		deconvolved.autocorrWaveform5nsRMS[0] = get5nsRMS(deconvAutocorrH);
		deconvolved.autocorrWaveform5nsRMS[1] = get5nsRMS(deconvAutocorrV);
		deconvolved.autocorrWaveform5nsValZScore[0] = (1 - deconvolved.autocorr5nsMean[0]) / deconvolved.autocorrWaveform5nsRMS[0];
		deconvolved.autocorrWaveform5nsValZScore[1] = (1 - deconvolved.autocorr5nsMean[1]) / deconvolved.autocorrWaveform5nsRMS[1];
		deconvolved.autocorrAnalytic5nsRMS[0] = getAutocorrAnalytic5nsRMS(deconvWaveformH);
		deconvolved.autocorrAnalytic5nsRMS[1] = getAutocorrAnalytic5nsRMS(deconvWaveformV);
		deconvolved.autocorrAnalytic5nsValZScore[0] = (1 - deconvolved.autocorr5nsMean[0]) / deconvolved.autocorrAnalytic5nsRMS[0];
		deconvolved.autocorrAnalytic5nsValZScore[1] = (1 - deconvolved.autocorr5nsMean[1]) / deconvolved.autocorrAnalytic5nsRMS[1];
		deconvolved.autocorrWaveformFullRMS[0] = deconvAutocorrH -> GetRMS(2);
		deconvolved.autocorrWaveformFullRMS[1] = deconvAutocorrV -> GetRMS(2);
		deconvolved.autocorrWaveformFullValZScore[0] = 1 / deconvolved.autocorrWaveformFullRMS[0];
		deconvolved.autocorrWaveformFullValZScore[1] = 1 / deconvolved.autocorrWaveformFullRMS[1];
		deconvolved.autocorrAnalyticFullRMS[0] = sqrt(2) * deconvolved.autocorrWaveformFullRMS[0];
		deconvolved.autocorrAnalyticFullRMS[1] = sqrt(2) * deconvolved.autocorrWaveformFullRMS[1];
		deconvolved.autocorrAnalyticFullValZScore[0] = 1 / deconvolved.autocorrAnalyticFullRMS[0];
		deconvolved.autocorrAnalyticFullValZScore[1] = 1 / deconvolved.autocorrAnalyticFullRMS[1];
		deconvolved.TIMdB[0] = getTIMdB(deconvWaveformH);
		deconvolved.TIMdB[1] = getTIMdB(deconvWaveformV);
		deconvolved.maxCorr[0] = YDeconvCorrH[maxIndDeconvCorrH];
		deconvolved.maxCorr[1] = YDeconvCorrV[maxIndDeconvCorrV];
		deconvolved.maxCorrTime[0] = XDeconvCorrH[maxIndDeconvCorrH];
		deconvolved.maxCorrTime[1] = XDeconvCorrV[maxIndDeconvCorrV];
		deconvolved.minCorr[0] = YDeconvCorrH[minIndDeconvCorrH];
		deconvolved.minCorr[1] = YDeconvCorrV[minIndDeconvCorrV];
		deconvolved.minCorrTime[0] = XDeconvCorrH[minIndDeconvCorrH];
		deconvolved.minCorrTime[1] = XDeconvCorrV[minIndDeconvCorrV];
		deconvolved.corrRMS[0] = deconvTemplateCorrH -> GetRMS(2);
		deconvolved.corrRMS[1] = deconvTemplateCorrV -> GetRMS(2);
		deconvolved.corrValZScore[0] = abs(deconvolved.maxCorr[0]) >= abs(deconvolved.minCorr[0]) ? deconvolved.maxCorr[0] : deconvolved.minCorr[0];
		deconvolved.corrValZScore[0] /= deconvolved.corrRMS[0];
		deconvolved.corrValZScore[1] = abs(deconvolved.maxCorr[1]) >= abs(deconvolved.minCorr[1]) ? deconvolved.maxCorr[1] : deconvolved.minCorr[1];
		deconvolved.corrValZScore[1] /= deconvolved.corrRMS[1];
		deconvolved.corrPeak[0] = max(abs(deconvolved.maxCorr[0]), abs(deconvolved.minCorr[0]));
		deconvolved.corrPeak[1] = max(abs(deconvolved.maxCorr[1]), abs(deconvolved.minCorr[1]));
		deconvolved.corrPeakTime[0] = abs(deconvolved.maxCorr[0]) >= abs(deconvolved.minCorr[0]) ? deconvolved.maxCorrTime[0] : deconvolved.minCorrTime[0];
		deconvolved.corrPeakTime[1] = abs(deconvolved.maxCorr[1]) >= abs(deconvolved.minCorr[1]) ? deconvolved.maxCorrTime[1] : deconvolved.minCorrTime[1];
		deconvolved.corrPeakZScore[0] = abs(deconvolved.corrValZScore[0]);
		deconvolved.corrPeakZScore[1] = abs(deconvolved.corrValZScore[1]);
		deconvolved.analyticCorrPeak[0] = YDeconvAnalyticCorrAbsH[maxIndDeconvAnalyticCorrAbsH];
		deconvolved.analyticCorrPeak[1] = YDeconvAnalyticCorrAbsV[maxIndDeconvAnalyticCorrAbsV];
		deconvolved.analyticCorrPeakTime[0] = XDeconvAnalyticCorrAbsH[maxIndDeconvAnalyticCorrAbsH];
		deconvolved.analyticCorrPeakTime[1] = XDeconvAnalyticCorrAbsV[maxIndDeconvAnalyticCorrAbsV];
		deconvolved.analyticCorrPeakZScore[0] = deconvolved.analyticCorrPeak[0] / (sqrt(2) * deconvolved.corrRMS[0]);
		deconvolved.analyticCorrPeakZScore[1] = deconvolved.analyticCorrPeak[1] / (sqrt(2) * deconvolved.corrRMS[1]);

		mapStruct.mapPeak[0] = sampleSummary -> peak[0][0].value;
		mapStruct.mapPeak[1] = sampleSummary -> peak[1][0].value;
		mapStruct.rectMapRMS[0] = sampleSummary -> peak[0][0].mapRMS;
		mapStruct.rectMapRMS[1] = sampleSummary -> peak[1][0].mapRMS;
		mapStruct.sphMapRMS[0] = UCorrelator::getMapRMS(corrMapH);
		mapStruct.sphMapRMS[1] = UCorrelator::getMapRMS(corrMapV);
		mapStruct.rectMapSNR[0] = sampleSummary -> peak[0][0].snr;
		mapStruct.rectMapSNR[1] = sampleSummary -> peak[1][0].snr;
		mapStruct.sphMapSNR[0] = UCorrelator::getMapSNR(corrMapH);
		mapStruct.sphMapSNR[1] = UCorrelator::getMapSNR(corrMapV);
		mapStruct.rectMapPeakZScore[0] = UCorrelator::getMapPeakZScore(corrMapH, false);
		mapStruct.rectMapPeakZScore[1] = UCorrelator::getMapPeakZScore(corrMapV, false);
		mapStruct.sphMapPeakZScore[0] = UCorrelator::getMapPeakZScore(corrMapH, true);
		mapStruct.sphMapPeakZScore[1] = UCorrelator::getMapPeakZScore(corrMapV, true);
		mapStruct.mapPeakNegTheta[0] = -sampleSummary -> peak[0][0].theta;
		mapStruct.mapPeakNegTheta[1] = -sampleSummary -> peak[1][0].theta;

		//  Filling the trees.
		cutTree.Fill();

		//  Further annoyance with the necessity of memory clearance.
		analyzer.clearInteractiveMemory();
		delete cohWaveformH, delete cohWaveformV;
		delete deconvWaveformH, delete deconvWaveformV;
		delete cohAutocorrH, delete cohAutocorrV;
		delete deconvAutocorrH, delete deconvAutocorrV;
		delete cohTemplateCorrH, delete cohTemplateCorrV;
		delete deconvTemplateCorrH, delete deconvTemplateCorrV;
		delete corrMapH, delete corrMapV;
	}
	
	//  Writing trees to cutFile.
	cutFile.cd();
	cutTree.Write();
	
	//  Closing files.
	cutFile.Close();
	sampleFile.Close();
	templateFile.Close();
}