/*
 * analyzer.cxx
 *
 *  Created on: Jan 9, 2013
 *      Author: yqiang
 *
 *      Modified: Oct 14, 2014 by yqiang for Active Collimator
 *
 */

#include <TApplication.h>
#include <TStyle.h>
#include <TGaxis.h>
#include <TCanvas.h>
#include <TGFrame.h>
#include <TGButton.h>
#include <TGLabel.h>
#include <TGFileDialog.h>
#include <TText.h>
#include <TLegend.h>
#include <TSystem.h>
#include <TH2.h>
#include <TF1.h>
#include <TProfile.h>
#include <TObject.h>
#include <TVirtualFFT.h>
#include <TMath.h>
#include <TList.h>
#include <stdio.h>
#include <iostream>
#include <fstream>
#include <math.h>
#include <TGenericClassInfo.h>
#include "analyzer.hh"
#include "pxiData.h"

ClassImp(TObjTimeStamp);
ClassImp(MyMainFrame);

const char *filetypes[] = { "ROOT files", "*.root", 0, 0 };
const char *offsetfiletypes[] = { "Offset files", "*.txt", 0, 0 };

using namespace std;

MyMainFrame::MyMainFrame(const TGWindow *p, UInt_t w, UInt_t h) :
		TGMainFrame(p, w, h) {

	// Fill color
	gStyle->SetStatColor(0);
	gStyle->SetTitleFillColor(0);
	gStyle->SetCanvasColor(0);
	gStyle->SetPadColor(0);
	gStyle->SetFrameFillColor(0);
	// Border mode
	gStyle->SetCanvasBorderMode(0);
	gStyle->SetPadBorderMode(0);
	gStyle->SetFrameBorderMode(0);
	// Margin
	gStyle->SetPadLeftMargin(0.12);
	gStyle->SetPadRightMargin(0.04);
	gStyle->SetPadTopMargin(0.06);
	gStyle->SetPadBottomMargin(0.12);
	// Font
	gStyle->SetTextFont(132);
	gStyle->SetLabelFont(132, "XYZ");
	gStyle->SetTitleFont(132, "XYZ");
	gStyle->SetStatFont(132);
	gStyle->SetLegendFont(132);
	// Fontsize
	gStyle->SetLabelSize(0.05, "XYZ");
	gStyle->SetTitleSize(0.06, "XYZ");
	gStyle->SetTitleOffset(0.9, "XY");
	// Opt
	gStyle->SetOptTitle(0);
	gStyle->SetOptFit(1);
	gStyle->SetOptStat(0);
	gStyle->SetStatX(0.97);
	gStyle->SetStatY(0.98);
	// Axis
	TGaxis::SetMaxDigits(3);
	// Grid
	gStyle->SetPadGridX(kTRUE);
	gStyle->SetPadGridY(kTRUE);

	for (int i = 0; i < 100; i++)
		datafiles[i] = 0;
	nfile = 0;
	datatree = 0;
	h_history = 0;
	for (int i = 0; i < ch_voltage; i++) {
		h_voltage[i] = 0;
		offset_voltage[i] = 0.;
	}
	version = 0;
	isTree = false;
	basename = TString("");
	workingdir = TString(gSystem->WorkingDirectory());

	/*
	 * Horizontal frame for file opening
	 */
	TGHorizontalFrame *frame_dir = new TGHorizontalFrame(this, 1000, 40);

	// button to open file
	TGTextButton *button_open = new TGTextButton(frame_dir, "&Open");
	button_open->Connect("Clicked()", "MyMainFrame", this, "DoOpen()");
	frame_dir->AddFrame(button_open,
			new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 5, 2, 2));
	label_dir = new TGLabel(frame_dir, "[input data file] ");
	label_dir->SetTextJustify(kTextLeft | kLHintsExpandX | kTextCenterY);
	label_dir->SetWrapLength(-1);
	frame_dir->AddFrame(label_dir,
			new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 5, 2, 2));

	// button to exit
	TGTextButton *button_exit = new TGTextButton(frame_dir, "&Exit",
			"gApplication->Terminate(0)");
	frame_dir->AddFrame(button_exit,
			new TGLayoutHints(kLHintsRight | kLHintsCenterY, 10, 5, 2, 2));

	// button to reset offset to 0
	TGTextButton *button_reset = new TGTextButton(frame_dir, "Reset");
	button_reset->Connect("Clicked()", "MyMainFrame", this, "DoResetOffset()");
	frame_dir->AddFrame(button_reset,
			new TGLayoutHints(kLHintsRight | kLHintsCenterY, 5, 10, 2, 2));

	// button to load offset
	TGTextButton *button_load = new TGTextButton(frame_dir, "Load");
	button_load->Connect("Clicked()", "MyMainFrame", this, "DoLoadOffset()");
	frame_dir->AddFrame(button_load,
			new TGLayoutHints(kLHintsRight | kLHintsCenterY, 5, 5, 2, 2));

	// button to save offset
	TGTextButton *button_save = new TGTextButton(frame_dir, "Save");
	button_save->Connect("Clicked()", "MyMainFrame", this, "DoSaveOffset()");
	frame_dir->AddFrame(button_save,
			new TGLayoutHints(kLHintsRight | kLHintsCenterY, 5, 5, 2, 2));

	// button to zero signal
	TGTextButton *button_zero = new TGTextButton(frame_dir, "Calculate");
	button_zero->Connect("Clicked()", "MyMainFrame", this, "DoZeroOffset()");
	frame_dir->AddFrame(button_zero,
			new TGLayoutHints(kLHintsRight | kLHintsCenterY, 5, 5, 2, 2));

	// offset file
	label_offset = new TGLabel(frame_dir, "Empty");
	label_offset->SetTextJustify(kTextLeft | kLHintsExpandX | kTextCenterY);
	label_offset->SetWrapLength(-1);
	frame_dir->AddFrame(label_offset,
			new TGLayoutHints(kLHintsRight | kLHintsCenterY, 5, 5, 2, 2));

	// offset label
	frame_dir->AddFrame(new TGLabel(frame_dir, "Offset File:"),
			new TGLayoutHints(kLHintsRight | kLHintsCenterY, 10, 5, 2, 2));

	AddFrame(frame_dir, new TGLayoutHints(kLHintsExpandX, 2, 2, 2, 2));

	/*
	 * Canvas for history plot
	 */
	canvas_history = new TRootEmbeddedCanvas("Canvas_History", this, 1000, 250);
	AddFrame(canvas_history, new TGLayoutHints(kLHintsExpandX, 2, 2, 2, 2));

	/*
	 * Horizontal frame for plot parameters
	 */
	TGHorizontalFrame *frame_cut = new TGHorizontalFrame(this, 1000, 40);

	// Data resample, 1-10000 Hz
	field_resample = new TGNumberEntryField(frame_cut, -1, 128,
			TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative,
			TGNumberFormat::kNELLimitMinMax, 1, 10000);
	field_resample->SetWidth(40);
	frame_cut->AddFrame(new TGLabel(frame_cut, "Resample Freq."),
			new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 2, 2, 2));
	frame_cut->AddFrame(field_resample,
			new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 2, 5, 2, 2));

	// check box for legend
	check_legend = new TGCheckButton(frame_cut, "Legend");
	frame_cut->AddFrame(check_legend,
			new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 10, 5, 2, 2));

	// check box for error bars
	check_error = new TGCheckButton(frame_cut, "Error");
	frame_cut->AddFrame(check_error,
			new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 5, 2, 2));

	// check box for mean value
	check_mean = new TGCheckButton(frame_cut, "Mean");
	frame_cut->AddFrame(check_mean,
			new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 5, 2, 2));

	// check box to do pol2 fit
	check_fit = new TGCheckButton(frame_cut, "Fit");
	frame_cut->AddFrame(check_fit,
			new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 5, 2, 2));

	// check box for denoise
	check_denoise = new TGCheckButton(frame_cut, "DeNoise");
	frame_cut->AddFrame(check_denoise,
			new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 20, 2, 2));

	// button to do plot voltage
	TGTextButton *button_plotvoltage = new TGTextButton(frame_cut, "Voltage");
	button_plotvoltage->Connect("Clicked()", "MyMainFrame", this,
			"DoPlotVoltage()");
	frame_cut->AddFrame(button_plotvoltage,
			new TGLayoutHints(kLHintsLeft | kLHintsExpandX | kLHintsCenterY, 5,
					5, 2, 2));

	// button to create root file
	TGTextButton *button_root = new TGTextButton(frame_cut, "Save ROOT File");
	button_root->Connect("Clicked()", "MyMainFrame", this, "DoRoot()");
	frame_cut->AddFrame(button_root,
			new TGLayoutHints(kLHintsRight | kLHintsCenterY, 20, 5, 2, 2));

	AddFrame(frame_cut, new TGLayoutHints(kLHintsExpandX));

	// button to print
	TGTextButton *button_print = new TGTextButton(frame_cut, "&Print");
	button_print->Connect("Clicked()", "MyMainFrame", this, "DoPrint()");
	frame_cut->AddFrame(button_print,
			new TGLayoutHints(kLHintsRight | kLHintsExpandX | kLHintsCenterY,
					20, 5, 2, 2));
	/*
	 * Frame to generate plots
	 */
//	TGHorizontalFrame *frame_plot = new TGHorizontalFrame(this, 1000, 40);
//	AddFrame(frame_plot, new TGLayoutHints(kLHintsExpandX));
	/*
	 * Canvas for misc plots
	 */
	canvas_misc = new TRootEmbeddedCanvas("Canvas_MISC", this, 800, 400);
	AddFrame(canvas_misc,
			new TGLayoutHints(kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));

	SetWindowName("AC ROOT Analysis");
	MapSubwindows();
	Resize(GetDefaultSize());
	MapWindow();
}

MyMainFrame::~MyMainFrame() {
	Cleanup();
	gApplication->Terminate(0);
}

/*
 * Action to open a rootfile, fill histograms and plot history chart
 */
void MyMainFrame::DoOpen() {

	std::cout << std::endl;

	// get list of files from GUI
	TGFileInfo fi;
	fi.fFileTypes = filetypes;
	fi.fIniDir = StrDup(workingdir);
	fi.fMultipleSelection = true;
	new TGFileDialog(gClient->GetRoot(), this, kFDOpen, &fi);
	if (!fi.fFileNamesList) {
		std::cout << "No file selected!" << std::endl;
		return;
	}
	TList *filelist = fi.fFileNamesList;
//	filelist->Print();

// reset existing root tree
	if (datatree)
		delete datatree;

	// reset data files
	if (nfile > 0) {
		for (int i = 0; i < nfile; i++) {
			if (datafiles[i]->IsOpen()) {
				datafiles[i]->Close();
			}
			if (datafiles[i]) {
				delete datafiles[i];
				datafiles[i] = 0;
			}
		}
		nfile = 0;
	}

	// reset history plot
	h_history = (TH1F*) gROOT->FindObject("h_history");
	if (h_history)
		delete h_history;

	// get link to first object in the list
	TObjLink *lnk = filelist->FirstLink();
	// update GUI with first filename
	label_dir->ChangeText(lnk->GetObject()->GetName());
	this->Layout();
	// open all files, max 100 files
	while (lnk && nfile <= 100) {
		datafiles[nfile] = new TFile(lnk->GetObject()->GetName(), "r");
		std::cout << nfile + 1 << " " << datafiles[nfile]->GetName() << " added"
				<< std::endl;
		lnk = lnk->Next();
		nfile++;
	}

	// if the data are saved as processed histograms, only the first file will be checked
	datafiles[0]->cd();
	datatree = (TTree *) gROOT->FindObject("T");
	if (datatree) {
		isTree = true;
		h_current = ReadTree(datatree, name_current, 1. / scale_current.coeff);
		if (h_current) {
			for (int i = 0; i < ch_voltage; i++)
				h_voltage[i] = ReadTree(datatree, name_voltage[i],
						1. / scale_voltage[i].coeff);

			// read TimeStamp object if available, otherwise get time info from filename
			TObjTimeStamp *ts = (TObjTimeStamp*) gROOT->FindObject("StartTime");
			if (ts) {
				starttime = ts->GetTime();
			} else {
				TString timestring = label_dir->GetText()->GetString();
				timestring.Replace(0, timestring.Last('/') + 1, 0, 0);
				timestring.ReplaceAll(".root", "");
				timestring.ReplaceAll(":", "");
				timestring.ReplaceAll("-", "");
				timestring.ReplaceAll("_", "");
				timestring.ReplaceAll("ac", "");

				Int_t year, month, day, hour, min, sec;
				sscanf(timestring.Data(), "%04d%02d%02d%02d%02d%02d", &year,
						&month, &day, &hour, &min, &sec);
				starttime.Set(year, month, day, hour, min, sec, 0, kTRUE, 0);
			}
			//		std::cout << "time string is " << starttime.AsString("s") << std::endl;
		}
	}

	// Data are in raw format
	else {
		isTree = false;
		TString name_prefix = "";
		pxiData *voltage_data[ch_voltage];
		voltage_data[0] = new pxiData(filelist, name_prefix + name_voltage[0]);
		starttime = voltage_data[0]->starttime;
		Double_t freq = 0;
		Int_t length = voltage_data[0]->GetLength(freq);
		h_voltage[0] = voltage_data[0]->GetHistogram(starttime, length, freq);
		if (h_voltage[0]) {
			h_current = voltage_data[0]->GetHistogram(starttime, length, freq,
					true);
			for (int i = 1; i < ch_voltage; i++) {
				voltage_data[i] = new pxiData(filelist,
						name_prefix + name_voltage[i]);
				h_voltage[i] = voltage_data[i]->GetHistogram(starttime, length,
						freq);
			}
			starttime.SetSec(starttime.GetSec() - 4 * 60 * 60);
		}
	}
//	delete file;

// plot current history
	if (h_current) {
		UInt_t tmpyear;
		UInt_t tmpmonth;
		UInt_t tmpday;
		UInt_t tmphour;
		UInt_t tmpmin;
		UInt_t tmpsec;
		starttime.GetDate(true, 0, &tmpyear, &tmpmonth, &tmpday);
		starttime.GetTime(true, 0, &tmphour, &tmpmin, &tmpsec);
		starttime.Print();

		// add time offset to translate all X axis into absolute time
		TDatime da(tmpyear, tmpmonth, tmpday, tmphour, tmpmin, tmpsec);
		Double_t timeoffset = Double_t(da.Convert())
				+ Double_t(starttime.GetNanoSec()) * 1.e-9;
		gStyle->SetTimeOffset(timeoffset);

		std::cout << "The data file is in version " << version << std::endl;

		if (isTree) {
			label_offset->ChangeText("Not Needed");
		} else {
			label_offset->ChangeText("Empty");
		}
		this->Layout();

		Double_t xmin = h_voltage[0]->GetXaxis()->GetXmin();
		Double_t xmax = h_voltage[0]->GetXaxis()->GetXmax();
		h_history = Resample(h_current, 128, 0, 0, kFALSE, kTRUE, xmin, xmax,
				kBlack, scale_current);
		if (h_history) {
			h_history->SetName("h_history");
			TCanvas *fCanvas = canvas_history->GetCanvas();
			fCanvas->cd();
			h_history->GetXaxis()->SetTimeDisplay(1);
			h_history->GetXaxis()->SetTimeFormat("%H:%M:%S");
			h_history->Draw("hist");
			fCanvas->Update();
			std::cout << "All data loaded" << std::endl;
			DoResetOffset();
		}

		basename = label_dir->GetText()->GetString();
		basename.Replace(0, basename.Last('/') + 1, 0, 0);
		basename.ReplaceAll(".root", "");
		basename.ReplaceAll(":", "");
		basename.ReplaceAll("-", "");
	} else
		std::cout << "File is corrupted or not compatible with this analyzer!"
				<< std::endl;
}

/*
 * Read branch from tree
 */
TH1F *ReadTree(TTree *tree, TString branchname, Double_t scale) {

	if (!tree->GetBranch(branchname))
		return 0;
	if (!tree->GetBranch("time"))
		return 0;

	TString hname = branchname;
	TString htitle = TString("History plot of ") + hname
			+ TString(";Time (s);Value (V)");
	hname.Prepend("hist_");

	TH1F *h1 = (TH1F*) gROOT->FindObject(hname);
	if (h1)
		delete h1;

	tree->SetBranchStatus("*", 0);

	Double_t tmpt;
	tree->SetBranchStatus("time", 1);
	tree->SetBranchAddress("time", &tmpt);
	Int_t entries = tree->GetEntries();
	if (entries < 2)
		return 0;

	tree->GetEntry(0);
	Double_t t_start = tmpt;
	tree->GetEntry(entries - 1);
	Double_t t_end = tmpt;
	Double_t t_hw = (t_end - t_start) / Double_t(entries - 1) / 2.;
	t_start -= t_hw;
	t_end += t_hw;

	Float_t tmpy;
	tree->SetBranchStatus(branchname, 1);
	tree->SetBranchAddress(branchname, &tmpy);
	h1 = new TH1F(hname, htitle, entries, t_start, t_end);
	for (int i = 0; i < entries; i++) {
		tree->GetEntry(i);
		h1->SetBinContent(i + 1, Double_t(tmpy) * scale);
	}
	std::cout << branchname << " was read" << std::endl;
	return h1;
}

/*
 * Action to zero signals
 */
void MyMainFrame::DoZeroOffset() {

	if (!h_history)
		return;

	TCanvas *fCanvas = canvas_misc->GetCanvas();
	fCanvas->Clear();
	fCanvas->Divide(2, 1);

	Int_t bin_min = h_history->GetXaxis()->GetFirst();
	Int_t bin_max = h_history->GetXaxis()->GetLast();
	Double_t xmin = h_history->GetXaxis()->GetBinLowEdge(bin_min);
	Double_t xmax = h_history->GetXaxis()->GetBinUpEdge(bin_max);
	std::cout << "Offsets are calculated using data between " << xmin << " and "
			<< xmax << std::endl;

	Double_t freq = 10;

	/*
	 * Calculate offset
	 */
	// voltage coils
	int nhr = 0;
	TH1F *hr[10];
	Double_t dummy_error;
	for (int i = 0; i < ch_voltage; i++) {
		hr[nhr] = Resample(h_voltage[i], freq, 0, 0, kFALSE, kTRUE, xmin, xmax,
				col_voltage[i]);
		offset_voltage[i] = -Average(hr[nhr], dummy_error);
		nhr++;
	}
	ComboPlot(fCanvas->GetPad(1), hr, nhr, "Before", check_legend->IsOn());

	/*
	 * Plot signal after zeroing
	 */

	// voltage coils
	nhr = 0;
	for (int i = 0; i < ch_voltage; i++) {
		hr[nhr] = Resample(h_voltage[i], freq, offset_voltage[i], 0, kFALSE,
				kTRUE, xmin, xmax, col_voltage[i]);
		nhr++;
	}
	ComboPlot(fCanvas->GetPad(2), hr, nhr, "After", check_legend->IsOn());

	fCanvas->Update();

	label_offset->ChangeText("In Memory");
	this->Layout();
}

/*
 * Action to save offset
 */
void MyMainFrame::DoSaveOffset() {
	if (!h_history)
		return;
	TString offsetdir = workingdir + TString("/offsets/");
	gSystem->MakeDirectory(offsetdir);
	TGFileInfo fi;
	fi.fFileTypes = offsetfiletypes;
	fi.fIniDir = StrDup(offsetdir);
	new TGFileDialog(gClient->GetRoot(), this, kFDSave, &fi);
	if (!fi.fFilename)
		return;
	else {
		ofstream offsetfile;
		offsetfile.open(fi.fFilename);
		if (offsetfile.is_open()) {
			for (int i = 0; i < ch_voltage; i++)
				offsetfile << offset_voltage[i] << "\n";
			offsetfile.close();
			std::cout << "offset file " << fi.fFilename << " was saved."
					<< std::endl;
		}
	}
	TString filename = fi.fFilename;
	filename.Replace(0, filename.Last('/') + 1, 0, 0);
	label_offset->ChangeText(filename);
	this->Layout();
}

/*
 * Action to load offset
 */
void MyMainFrame::DoLoadOffset() {
	TString offsetdir = workingdir + TString("/offsets/");
	TGFileInfo fi;
	fi.fFileTypes = offsetfiletypes;
	fi.fIniDir = StrDup(offsetdir);
	new TGFileDialog(gClient->GetRoot(), this, kFDOpen, &fi);
	if (!fi.fFilename)
		return;
	else {
		ifstream offsetfile;
		offsetfile.open(fi.fFilename);
		if (offsetfile.is_open()) {
			for (int i = 0; i < ch_voltage; i++)
				offsetfile >> offset_voltage[i];
			offsetfile.close();
			std::cout << "offsets were read from file " << fi.fFilename << "."
					<< std::endl;
		}
	}
	TString filename = fi.fFilename;
	filename.Replace(0, filename.Last('/') + 1, 0, 0);
	label_offset->ChangeText(filename);
	this->Layout();
}

/*
 * Action to reset offset
 */
void MyMainFrame::DoResetOffset() {
	for (int i = 0; i < ch_voltage; i++) {
		offset_voltage[i] = 0.;
	}
	this->Layout();
}

/*
 * Action to make voltage plot
 */
void MyMainFrame::DoPlotVoltage() {
	if (!h_history)
		return;

	TCanvas *fCanvas = canvas_misc->GetCanvas();
	fCanvas->Clear();
	fCanvas->SetTitle("Voltage");
//	fCanvas->Divide(1, 1);

	Int_t bin_min = h_history->GetXaxis()->GetFirst();
	Int_t bin_max = h_history->GetXaxis()->GetLast();
	Double_t xmin = h_history->GetXaxis()->GetBinLowEdge(bin_min);
	Double_t xmax = h_history->GetXaxis()->GetBinUpEdge(bin_max);
//	std::cout << "X axis range is chosen from " << xmin << " to " << xmax
//	<< std::endl;

// Plot Pickup Signal
	TH1F *phr[30];
	Int_t nhr = 0;
	for (int i = 0; i < ch_voltage; i++) {
		phr[nhr] = Resample(h_voltage[i], field_resample->GetNumber(),
				offset_voltage[i], 0, check_denoise->IsOn(), kTRUE, xmin, xmax,
				col_voltage[i], scale_voltage[i]);
		nhr++;
	}
	ComboPlot(fCanvas->GetPad(0), phr, nhr,
			starttime.AsString("s") + TString(" Voltage Readings"),
			check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
			kFALSE, kFALSE, -10, 10);

	fCanvas->Update();
}

/*
 * Action to create root file
 */
void MyMainFrame::DoRoot() {

	if (!h_history)
		return;
// find range
	Int_t bin_min = h_history->GetXaxis()->GetFirst();
	Int_t bin_max = h_history->GetXaxis()->GetLast();
	Double_t xmin = h_history->GetXaxis()->GetBinLowEdge(bin_min);
	Double_t xmax = h_history->GetXaxis()->GetBinUpEdge(bin_max);

// create roottree

	TString rootdir = workingdir + TString("/rootfiles/");
	gSystem->MakeDirectory(rootdir);
	TString filename = rootdir + basename + ".root";
	TFile *rootfile = new TFile(filename, "RECREATE");
	TTree *tree = (TTree*) gROOT->FindObject("T");
	if (tree)
		delete tree;
	tree = new TTree("T", "Tree for resampled PXI data");

// create resampled histograms and branches
	TH1F *hr[50];
	Float_t value[50];
	Double_t time;

	tree->Branch("time", &time, "time/D");

	Int_t nhr = 0;
	std::cout << "start data resampling ..." << std::endl;
// Resample Voltage signals
	for (int i = 0; i < ch_voltage; i++) {
		hr[nhr] = Resample(h_voltage[i], field_resample->GetNumber(),
				offset_voltage[i], 0, check_denoise->IsOn(), kTRUE, xmin, xmax,
				col_voltage[i], scale_voltage[i]);
		tree->Branch(name_voltage[i], value + nhr,
				name_voltage[i] + TString("/F"));
		nhr++;
	}
// Resample Current reading
	hr[nhr] = Resample(h_current, field_resample->GetNumber(), 0, 0, kFALSE,
			kBlack, xmin, xmax, 0, scale_current);
	tree->Branch(name_current, value + nhr, name_current + TString("/F"));
	nhr++;
	std::cout << "finished resampling data and now saving to a ROOT Tree ..."
			<< std::endl;

// Fill ROOT Tree
	Int_t entries = hr[0]->GetNbinsX();
	for (int i = 1; i <= entries; i++) {
		time = hr[0]->GetBinCenter(i);
		for (int j = 0; j < nhr; j++) {
			value[j] = hr[j]->GetBinContent(i);
		}
		tree->Fill();
	}
	rootfile->cd();
	tree->Write("T", TObject::kOverwrite);
	TObjTimeStamp *ts = new TObjTimeStamp(starttime);
	ts->Print();
	ts->Write("StartTime");
	rootfile->Close();
	std::cout << "rootfile " << filename << " has been created." << std::endl;
}

/*
 * Action to print screen
 */
void MyMainFrame::DoPrint() {
	if (!h_history)
		return;
	TCanvas *fCanvas = canvas_misc->GetCanvas();
	if (TString(fCanvas->GetTitle()).Length() == 0)
		return;
	fCanvas->Update();
	TString pngfile = label_dir->GetText()->GetString();
	pngfile.Replace(0, pngfile.Last('/') + 1, 0, 0);
	pngfile.ReplaceAll(".root", "");
	pngfile.ReplaceAll(":", "");
	pngfile.ReplaceAll("-", "");
	TString plotdir = workingdir + TString("/plots/");
	gSystem->MakeDirectory(plotdir);
	TString pdffile = plotdir + pngfile + "_" + fCanvas->GetTitle() + ".pdf";
	pngfile = plotdir + pngfile + "_" + fCanvas->GetTitle() + ".png";
	fCanvas->Print(pngfile);
	fCanvas->Print(pdffile);
}

Scale CreateScale(Double_t coeff, TString unit) {
	Scale scale = { coeff, unit };
	return scale;
}

/*
 * Resample data to lower frequency than the base frequency (upto 10000 Hz)
 */
TH1F *Resample(TH1F *h1, Double_t freq, Double_t offset, Double_t delay,
		Bool_t denoise, Bool_t cut, Double_t xmin, Double_t xmax, Color_t color,
		struct Scale scale) {

	if (!h1)
		return 0;

	Int_t nbin = h1->GetNbinsX();
	Int_t bin_xmin, bin_xmax;
	if (cut) {
		bin_xmin = h1->FindBin(xmin);
		bin_xmax = h1->FindBin(xmax);
	} else {
		bin_xmin = 1;
		bin_xmax = nbin;
	}

	Double_t base_frequency = nbin
			/ (h1->GetXaxis()->GetBinUpEdge(nbin)
					- h1->GetXaxis()->GetBinLowEdge(1));
// round base frequency to 1 Hz
	base_frequency = floor(base_frequency + 0.5);
// number of bins to be combined for a new bin
	Int_t nrebin = Int_t(floor(base_frequency / freq + 0.5));
	if (nrebin < 1)
		nrebin = 1;
	Int_t bin_delay = Int_t(floor(delay * base_frequency + 0.5));
// number of new bins
	Int_t nbin_new = Int_t(
			floor(Double_t(bin_xmax - bin_xmin + 1) / Double_t(nrebin)));
	Double_t xmin_new = h1->GetXaxis()->GetBinLowEdge(bin_xmin);
	Double_t xmax_new = h1->GetXaxis()->GetBinUpEdge(
			bin_xmin + nbin_new * nrebin - 1);

	TString hname = TString(h1->GetName());
	hname.Append("_res");
	TString htitle = TString(h1->GetName()) + TString(';')
			+ TString(h1->GetXaxis()->GetTitle()) + TString(';') + scale.unit;

	TH1F *h2 = (TH1F*) gROOT->FindObject(hname);
	if (h2)
		delete h2;

	h2 = new TH1F(hname, htitle, nbin_new, xmin_new, xmax_new);
	for (int i = 0; i < nbin_new; i++) {
		Double_t tmpy = 0;
		Double_t tmpy2 = 0;
		for (int j = 0; j < nrebin; j++) {
			Int_t bin_old = bin_xmin + i * nrebin + j + 1 + bin_delay;
			if (bin_old > nbin) {
				tmpy += h1->GetBinContent(nbin);
				tmpy2 += pow(h1->GetBinContent(nbin), 2.);
			} else {
				tmpy += h1->GetBinContent(
						bin_xmin + i * nrebin + j + 1 + bin_delay);
				tmpy2 += pow(
						h1->GetBinContent(
								bin_xmin + i * nrebin + j + 1 + bin_delay), 2);
			}
		}
		tmpy = tmpy / Double_t(nrebin);
		tmpy2 = sqrt(tmpy2 / Double_t(nrebin) - tmpy * tmpy);
		h2->SetBinContent(i + 1, offset + tmpy);
		h2->SetBinError(i + 1, tmpy2 / sqrt(Double_t(nrebin)));
	}
	h2->Scale(scale.coeff);
	h2->SetLineColor(color);
	if (!denoise)
		return h2;
	else {
		TH1F* h3 = DeNoise(h2);
		return h3;
	}
}

/*
 * Generate Integral Histogram
 */
TH1F * IntHist(TH1F * h1, struct Scale scale) {
	if (!h1)
		return 0;
	TString hname = TString(h1->GetName());
	hname.Append("_int");
	TString htitle = TString(h1->GetName()) + TString(';')
			+ TString(h1->GetXaxis()->GetTitle()) + TString(';') + scale.unit;
	Double_t xmin = h1->GetXaxis()->GetXmin();
	Double_t xmax = h1->GetXaxis()->GetXmax();
	Int_t nbin = h1->GetNbinsX();
	Double_t bin_width = (xmax - xmin) / Double_t(nbin);
	TH1F *h2 = (TH1F*) gROOT->FindObject(hname);
	if (h2)
		delete h2;
	h2 = new TH1F(hname, htitle, nbin, xmin, xmax);
	Double_t tmpx = 0.;
	Double_t tmpx_e2 = 0.;
	for (int i = 1; i <= nbin; i++) {
		tmpx += h1->GetBinContent(i);
		tmpx_e2 += pow(h1->GetBinError(i), 2.);
		h2->SetBinContent(i, tmpx);
		h2->SetBinError(i, sqrt(tmpx_e2));
	}
	h2->Scale(bin_width * scale.coeff);
	h2->SetLineColor(h1->GetLineColor());
	return h2;
}

/*
 * Create ratio histogram of h1/h2
 */
TH1F *Ratio(TH1F *h1, TH1F *h2, struct Scale scale) {
	if (!h1 || !h2)
		return 0;

	Double_t xmin = h1->GetXaxis()->GetXmin();
	Double_t xmin2 = h2->GetXaxis()->GetXmin();
	Double_t xmax = h1->GetXaxis()->GetXmax();
	Double_t xmax2 = h2->GetXaxis()->GetXmax();
	Int_t nbin = h1->GetNbinsX();
	Int_t nbin2 = h2->GetNbinsX();

	Double_t bin_width = (xmax - xmin) / Double_t(nbin);
	if (xmin > xmin2) {
		bin_width -= Int_t(floor((xmin - xmin2) / bin_width + 0.5));
		xmin = xmin2;
	}
	if (xmax < xmax2) {
		bin_width -= Int_t(floor((xmax2 - xmax) / bin_width + 0.5));
		xmax = xmax2;
	}
	if (bin_width <= 0)
		return 0;

	TString hname = TString(h1->GetName()) + "_to_" + TString(h2->GetName());
	hname.ReplaceAll("_res", "");
	hname.ReplaceAll("hist_", "");
	TH1F * h3 = (TH1F*) gROOT->FindObject(hname);
	if (h3)
		delete h3;
	TString htitle = TString("Ratio between ") + TString(h1->GetName())
			+ TString(" and ") + TString(h2->GetName());
	htitle.ReplaceAll("_res", 4, 0, 0);
	htitle = htitle + TString(";") + TString(h1->GetXaxis()->GetTitle())
			+ TString(";") + scale.unit;
	h3 = new TH1F(hname, htitle, nbin, xmin, xmax);
	for (int i = 1; i <= nbin; i++) {
		Double_t tmpx = h3->GetBinCenter(i);
		Double_t tmpy = h1->GetBinContent(h1->FindBin(tmpx));
		Double_t tmpy2 = h2->GetBinContent(h2->FindBin(tmpx));
		Double_t tmpy_err = h1->GetBinError(h1->FindBin(tmpx));
		Double_t tmpy2_err = h2->GetBinError(h2->FindBin(tmpx));
		if (tmpy != 0 && tmpy2 != 0)
			tmpy_err = sqrt(
					pow(tmpy_err / tmpy, 2.) + pow(tmpy2_err / tmpy2, 2.));
		if (tmpy2 == 0)
			tmpy = 0;
		else
			tmpy = tmpy / tmpy2;
		tmpy_err = tmpy * tmpy_err;
		h3->SetBinContent(i, tmpy);
		h3->SetBinError(i, tmpy_err);
	}
	h3->Scale(scale.coeff);
	h3->SetLineColor(h1->GetLineColor());
	return h3;
}

/*
 * Plot a collection of histograms on same pad
 */
void ComboPlot(TVirtualPad *pp, TH1F **hh, Int_t nh, TString title,
		Bool_t legend, Bool_t error, Bool_t mean, Bool_t fit, Bool_t manual,
		Double_t min, Double_t max) {

	pp->cd();
	Double_t ymin = 1000;
	Double_t ymax = -1000;
	if (!manual) {
		for (int i = 0; i < nh; i++) {
			Double_t lmin = hh[i]->GetBinContent(hh[i]->GetMinimumBin());
			Double_t lmax = hh[i]->GetBinContent(hh[i]->GetMaximumBin());
			if (lmin < ymin)
				ymin = lmin;
			if (lmax > ymax)
				ymax = lmax;
		}
		Double_t margin = (ymax - ymin) / 20.; // 5% margin
		ymin -= margin;
		ymax += margin;
	} else {
		ymin = min;
		ymax = max;
	}

	hh[0]->SetMinimum(ymin);
	hh[0]->SetMaximum(ymax);
	hh[0]->GetXaxis()->SetTimeDisplay(1);
	hh[0]->GetXaxis()->SetTimeFormat("%H:%M:%S");
	if (error)
		hh[0]->DrawCopy("e");
	else
		hh[0]->DrawCopy("hist");
	for (int i = 1; i < nh; i++)
		if (error)
			hh[i]->DrawCopy("e same");
		else
			hh[i]->DrawCopy("hist same");

// calculate mean value
	if (mean) {
		std::cout << "=== Average values in " << title << " ===" << std::endl;
		for (int i = 0; i < nh; i++) {
			Double_t error;
			Double_t mean = Average(hh[i], error);
			std::cout << hh[i]->GetName() << ":\t" << mean << "\t+-\t" << error
					<< std::endl;
		}
	}

// fit pol 2
	if (fit) {
		std::cout << "=== Pol2 fit results in " << title << " ===" << std::endl;
		TF1 *f1 = new TF1("f1", "pol2", 0, 1500);
		for (int i = 0; i < nh; i++) {
			hh[i]->Fit("f1", "NQ");
			Double_t par0 = f1->GetParameter(0);
			Double_t par1 = f1->GetParameter(1);
			Double_t par2 = f1->GetParameter(2);
			std::cout << hh[i]->GetName() << ":\t" << par0 << "\t" << par1
					<< "\t" << par2 << std::endl;
		}
		delete f1;
	}

// Draw Title
	TText *t1 = new TText(0.54, 0.97, title);
	t1->SetNDC(kTRUE);
	t1->SetTextAlign(22);
	t1->SetTextSize(0.06);
	t1->Draw();

// Draw legend
	if (legend) {
		Double_t loweredge = 0.96 - nh * 0.07;
		if (loweredge < 0.12)
			loweredge = 0.12;
		TLegend *l1 = new TLegend(0.77, loweredge, 0.97, 0.96, 0, "brNDC");
		for (int i = 0; i < nh; i++) {
			TString entry = hh[i]->GetName();
			entry.ReplaceAll("_res", "");
			entry.ReplaceAll("hist_", "");
			l1->AddEntry(hh[i], entry, "L");
		}
		l1->SetFillColor(kWhite);
		l1->SetTextAlign(22);
		l1->Draw();
	}

	pp->Update();
}

/*
 * Calculate average Y value
 */
Double_t Average(TH1F * hh, Double_t &error) {
	if (!hh)
		return 0.;
	Int_t nbin = hh->GetNbinsX();
	Double_t tmpy = 0;
	Double_t tmpy_error = 0;
	for (Int_t i = 1; i <= nbin; i++) {
		tmpy += hh->GetBinContent(i);
		tmpy_error += pow(hh->GetBinError(i), 2.);
	}
	tmpy = tmpy / Double_t(nbin);
	tmpy_error = sqrt(tmpy_error) / Double_t(nbin);
	error = tmpy_error;
//	std::cout << tmpy << "\t" << tmpy_error << std::endl;
	return tmpy;
}

/*
 * Create summed plot
 */
TH1F *AddHist(TH1F **hh, Int_t nhist, TString hname, TString htitle,
		Color_t color) {
	if (nhist <= 0)
		return 0;
	TH1F *h2 = (TH1F*) gROOT->FindObject(hname);
	if (h2)
		delete h2;
	h2 = (TH1F*) hh[0]->Clone(hname);
	h2->SetTitle(htitle);
	for (Int_t i = 1; i < nhist; i++) {
		h2->Add(hh[i]);
	}
	h2->SetLineColor(color);
	return h2;
}

/*
 * Populate a 1-D histogram of correlations of two plots
 */
TH1F *CorrHist(TH1F *h1, TH1F *h2, Int_t nbin, Bool_t manual, Double_t min,
		Double_t max) {

// find minimum bin number
	Int_t nx = h1->GetNbinsX();
	Int_t nx2 = h2->GetNbinsX();
	if (nx2 < nx)
		nx = nx2;
	if (nx <= 0)
		return 0;

// automatically determine range of X-axis
	if (!manual) {
		min = h2->GetBinContent(h2->GetMinimumBin());
		max = h2->GetBinContent(h2->GetMaximumBin());
	}

// determine range of Y-axis
	Double_t ymin = h1->GetBinContent(h1->GetMinimumBin());
	Double_t ymax = h1->GetBinContent(h1->GetMaximumBin());
	Int_t ybin = 500;

// create name and title, find and delete existing histograms
	TString hname = TString(h1->GetName()) + "_vs_" + TString(h2->GetName());
	hname.ReplaceAll("_res", 4, 0, 0);
	TH1F *h3 = (TH1F*) gROOT->FindObject(hname);
	if (h3)
		delete h3;
	TString htitle = TString("Correlation between ") + TString(h1->GetName())
			+ TString(" and ") + TString(h2->GetName());

// fill temporary 2-D histogram
	TH2F *dh2 = new TH2F("2dh", "", nbin, min, max, ybin, ymin, ymax);
	for (int i = 1; i <= nx; i++)
		dh2->Fill(h2->GetBinContent(i), h1->GetBinContent(i));
	dh2->ProfileX("_pfx");
	TProfile *dh2_pfx = (TProfile*) gROOT->FindObject("2dh_pfx");

// copy to 1-D histogram
	h3 = new TH1F(hname, htitle, nbin, min, max);
	for (int i = 1; i <= nbin; i++) {
		h3->SetBinContent(i, dh2_pfx->GetBinContent(i));
		h3->SetBinError(i, dh2_pfx->GetBinError(i));
	}

	delete dh2;
	delete dh2_pfx;

	h3->SetLineColor(h1->GetLineColor());
	h3->GetXaxis()->SetTitle(h2->GetYaxis()->GetTitle());
	h3->GetYaxis()->SetTitle(h1->GetYaxis()->GetTitle());
	return h3;
}

/*
 * Calculate derivative
 */
TH1F *Derivative(TH1F *h1, Double_t window) {
	if (!h1)
		return 0;

	TString htitle = TString("Derivative of") + TString(h1->GetName());
	TString hname = TString(h1->GetName()) + "_der";
	hname.ReplaceAll("_res", 4, 0, 0);
	TH1F *h2 = (TH1F*) gROOT->FindObject(hname);
	if (h2)
		delete h2;

	Double_t xmin = h1->GetXaxis()->GetXmin();
	Double_t xmax = h1->GetXaxis()->GetXmax();
	Int_t nbin = h1->GetNbinsX();
//	Double_t bin_width = (xmax - xmin) / Double_t(nbin);
	h2 = new TH1F(hname, htitle, nbin, xmin, xmax);
	Double_t frange = window / 2.;

	TF1 *f1 = new TF1("f1", "pol1", 0, 10);
	for (int i = 1; i <= nbin; i++) {
		Double_t tmpx = h1->GetBinCenter(i);
		h1->Fit(f1, "NQ", "", tmpx - frange, tmpx + frange);
		Double_t tmpy = f1->GetParameter(1); // + f1->GetParameter(2) * tmpx;
		Double_t tmpy_err = f1->GetParError(1);
		h2->SetBinContent(i, tmpy);
		h2->SetBinError(i, tmpy_err);
	}
	return h2;
}

/*
 * Remove harmonic noise
 */
TH1F *DeNoise(TH1F *h1) {

	Double_t snr_threshold = 3; // noise threhold 3 db
	Double_t filter_range = 1.02;
	Double_t noise_sup = 30; //additional noise suppresion factor in db.
	Int_t index_max = 4;
	Double_t list_noise[6] = { 60, 120, 180, 205, 240, 300 };

	Double_t xrange = h1->GetXaxis()->GetXmax() - h1->GetXaxis()->GetXmin();
	Int_t nbin = h1->GetNbinsX();
	TString hname = TString(h1->GetName());
	hname.Append("_de");
	TH1F *h2 = (TH1F*) gROOT->FindObject(hname);
	if (h2)
		delete h2;
	h2 = (TH1F*) h1->Clone(hname);

// FFT to frequency domain
	Double_t *re_full = new Double_t[nbin];
	Double_t *im_full = new Double_t[nbin];

	TVirtualFFT::SetTransform(0);
	TH1 *h1_mag = h1->FFT(0, "MAG");
	TVirtualFFT *fft = TVirtualFFT::GetCurrentTransform();
	fft->GetPointsComplex(re_full, im_full);

// Uniform filter
	Double_t *filter = new Double_t[nbin];
	for (Int_t i = 0; i < nbin; i++)
		filter[i] = 1.;
	for (Int_t index = 0; index < index_max; index++) {
		Int_t fmin = Int_t(floor(list_noise[index] * xrange / filter_range))
				- 1;
		Int_t fmax = Int_t(ceil(list_noise[index] * xrange * filter_range)) + 1;
		if (fmin < 1)
			fmin = 1;
		if (fmax >= nbin / 2)
			fmax = nbin / 2 - 1;
		for (Int_t i = fmin; i < fmax; i++) {
			filter[i] = filter[nbin - i - 1] = pow(10, -noise_sup / 10.);
		}
	}

// apply filter and transform back

	for (int i = 0; i < nbin; i++) {
		re_full[i] *= filter[i];
		im_full[i] *= filter[i];
	}

	TVirtualFFT *fft_back = TVirtualFFT::FFT(1, &nbin, "C2R M K");
	fft_back->SetPointsComplex(re_full, im_full);
	fft_back->Transform();
	TH1::TransformHisto(fft_back, (TH1*) h2, "Re");
	h2->Scale(1. / nbin);

	delete h1_mag;
	return h2;
}

int main(int argc, char **argv) {
	TApplication theApp("App", &argc, argv);
	new MyMainFrame(gClient->GetRoot(), 1000, 800);
	theApp.Run();
	return 0;
}