/*
 * analyzer.cxx
 *
 *  Created on: Jan 9, 2013
 *      Author: yqiang
 */

#include <time.h>
#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.h"
#include "pxiData.h"

using namespace std;

ClassImp(TObjTimeStamp);
ClassImp(MyMainFrame);

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

int DoStartupFile;
char InputFileList[256];

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);
  TGaxis::SetExponentOffset(-0.06, -0.04,"y");
  // Grid
  gStyle->SetPadGridX(kTRUE);
  gStyle->SetPadGridY(kTRUE);
  
  for (int i = 0; i < 100; i++)
    datafiles[i] = 0;
  nfile = 0;
  datatree = 0;
  h_current = 0;
  h_shunt = 0;
  h_field = 0;
  h_voltage = 0;
  h_history = 0;
  h_qd = 0;
  h_gfd = 0;
  offset_current = 0.;
  offset_shunt = 0.;
  offset_field = 0;
  offset_voltage = 0;
  offset_qd = 0;
  offset_gfd = 0;
  for (int i = 0; i < ch_vtt; i++) {
    h_vtt[i] = 0;
    offset_vtt[i] = 0.;
  }
  for (int i = 0; i < ch_vtts; i++) {
    h_vtts[i] = 0;
    offset_vtts[i] = 0.;
  }
  for (int i = 0; i < ch_vttb; i++) {
    h_vttb[i] = 0;
    offset_vttb[i] = 0.;
  }
  for (int i = 0; i < ch_pickup; i++) {
    h_pickup[i] = 0;
    offset_pickup[i] = 0.;
  }
  for (int i = 0; i < ch_lead; i++) {
    h_lead[i] = 0;
    offset_lead[i] = 0;
  }
  for (int i = 0; i < ch_acc; i++) {
    h_acc[i] = 0;
    offset_acc[i] = 0;
  }
  version = 2;
  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);
  
  // Hall probe delay, 0-10 s
  field_delay = new TGNumberEntryField(frame_cut, -1, 0.0,
				       TGNumberFormat::kNESRealTwo, TGNumberFormat::kNEANonNegative,
				       TGNumberFormat::kNELLimitMinMax, 0., 10.);
  field_delay->SetWidth(40);
  frame_cut->AddFrame(new TGLabel(frame_cut, "Hall Probe Delay"),
		      new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 2, 2, 2));
  frame_cut->AddFrame(field_delay,
		      new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 2, 5, 2, 2));
  
  // Data resample, 1-10000 Hz
  field_resample = new TGNumberEntryField(frame_cut, -1, 100,
					  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));
  
  // dI/dt window, 1-100s
  field_window = new TGNumberEntryField(frame_cut, -1, 10,
					TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative,
					TGNumberFormat::kNELLimitMinMax, 0.1, 100);
  field_window->SetWidth(40);
  frame_cut->AddFrame(new TGLabel(frame_cut, "dI/dt Window"),
		      new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 2, 2, 2));
  frame_cut->AddFrame(field_window,
		      new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 2, 10, 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 to use ZFCT current instead of shunt current
  check_ZFCT = new TGCheckButton(frame_cut, "ZFCT");
  frame_cut->AddFrame(check_ZFCT,
		      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 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));
  
  /*
   * Frame to generate plots
   */
  TGHorizontalFrame *frame_plot = new TGHorizontalFrame(this, 1000, 40);
  
  // button to do plot field
  TGTextButton *button_plotfield = new TGTextButton(frame_plot,
						    "Field/Cur.");
  button_plotfield->Connect("Clicked()", "MyMainFrame", this,
			    "DoPlotField()");
  frame_plot->AddFrame(button_plotfield,
		       new TGLayoutHints(kLHintsLeft | kLHintsExpandX | kLHintsCenterY, 5,
					 5, 2, 2));
  
  // button to do plot field to current field from pickup coils
  TGTextButton *button_plotratio = new TGTextButton(frame_plot,
						    "Pickups");
  button_plotratio->Connect("Clicked()", "MyMainFrame", this,
			    "DoPlotRatio()");
  frame_plot->AddFrame(button_plotratio,
		       new TGLayoutHints(kLHintsLeft | kLHintsExpandX | kLHintsCenterY, 5,
					 5, 2, 2));
  
  // button to do plot voltage
  TGTextButton *button_plotvoltage = new TGTextButton(frame_plot, "Voltage");
  button_plotvoltage->Connect("Clicked()", "MyMainFrame", this,
			      "DoPlotVoltage()");
  frame_plot->AddFrame(button_plotvoltage,
		       new TGLayoutHints(kLHintsLeft | kLHintsExpandX | kLHintsCenterY, 5,
					 5, 2, 2));
  
  // button to do plot voltage Coils 1 and 2
  TGTextButton *button_plotvoltageC12 = new TGTextButton(frame_plot, "V C1/2");
  button_plotvoltageC12->Connect("Clicked()", "MyMainFrame", this,
				 "DoPlotVoltageC12()");
  frame_plot->AddFrame(button_plotvoltageC12,
		       new TGLayoutHints(kLHintsLeft | kLHintsExpandX | kLHintsCenterY, 5,
					 5, 2, 2));
  
  // button to do plot voltage Coils 3 and 4
  TGTextButton *button_plotvoltageC34 = new TGTextButton(frame_plot, "V C3/4");
  button_plotvoltageC34->Connect("Clicked()", "MyMainFrame", this,
				 "DoPlotVoltageC34()");
  frame_plot->AddFrame(button_plotvoltageC34,
		       new TGLayoutHints(kLHintsLeft | kLHintsExpandX | kLHintsCenterY, 5,
					 5, 2, 2));
  
  // button to do plot voltage Coils 3 and 4
  TGTextButton *button_plotvoltageLeads = new TGTextButton(frame_plot, "Leads");
  button_plotvoltageLeads->Connect("Clicked()", "MyMainFrame", this,
				   "DoPlotVoltageLeads()");
  frame_plot->AddFrame(button_plotvoltageLeads,
		       new TGLayoutHints(kLHintsLeft | kLHintsExpandX | kLHintsCenterY, 5,
					 5, 2, 2));
  
  // button to plot quench detector and lead voltages
  TGTextButton *button_plotQD = new TGTextButton(frame_plot,
						 "Quench Det.");
  button_plotQD->Connect("Clicked()", "MyMainFrame", this, "DoPlotQD()");
  frame_plot->AddFrame(button_plotQD,
		       new TGLayoutHints(kLHintsLeft | kLHintsExpandX | kLHintsCenterY, 5,
					 10, 2, 2));
  
  // button to calculate Inductance
  TGTextButton *button_plotInductance = new TGTextButton(frame_plot,
							 "Inductance");
  button_plotInductance->Connect("Clicked()", "MyMainFrame", this,
				 "DoPlotInductance()");
  frame_plot->AddFrame(button_plotInductance,
		       new TGLayoutHints(kLHintsLeft | kLHintsExpandX | kLHintsCenterY, 10,
					 5, 2, 2));
  
  // button to calculate resistance
  TGTextButton *button_plotResistance = new TGTextButton(frame_plot,
							 "Resistance");
  button_plotResistance->Connect("Clicked()", "MyMainFrame", this,
				 "DoPlotResistance()");
  frame_plot->AddFrame(button_plotResistance,
		       new TGLayoutHints(kLHintsLeft | kLHintsExpandX | kLHintsCenterY, 5,
					 5, 2, 2));
  
  // button to calculate temperature
  TGTextButton *button_plotAccelerometers = new TGTextButton(frame_plot,
							     "Accel.");
  button_plotAccelerometers->Connect("Clicked()", "MyMainFrame", this,
				     "DoPlotAccelerometers()");
  frame_plot->AddFrame(button_plotAccelerometers,
		       new TGLayoutHints(kLHintsLeft | kLHintsExpandX | kLHintsCenterY, 5,
					 10, 2, 2));
  
  // button to calculate temperature
  TGTextButton *button_plotTemperature = new TGTextButton(frame_plot,
							  "Temp.");
  button_plotTemperature->Connect("Clicked()", "MyMainFrame", this,
				  "DoPlotTemperature()");
  frame_plot->AddFrame(button_plotTemperature,
		       new TGLayoutHints(kLHintsLeft | kLHintsExpandX | kLHintsCenterY, 5,
					 10, 2, 2));
  
  // button to print
  TGTextButton *button_print = new TGTextButton(frame_plot, "&Print");
  button_print->Connect("Clicked()", "MyMainFrame", this, "DoPrint()");
  frame_plot->AddFrame(button_print,
		       new TGLayoutHints(kLHintsRight | kLHintsExpandX | kLHintsCenterY,
					 20, 5, 2, 2));
  
  AddFrame(frame_plot, new TGLayoutHints(kLHintsExpandX));
  
  /*
   * Canvas for misc plots: zero and ratio
   */
  canvas_misc = new TRootEmbeddedCanvas("Canvas_MISC", this, 800, 800);
  AddFrame(canvas_misc,
	   new TGLayoutHints(kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
  
  canvas_misc->GetCanvas()->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)","MyMainFrame",this, 
				    "EventInfo(Int_t,Int_t,Int_t,TObject*)");
  
  // status bar
  Int_t parts[] = {20, 15, 10, 55};
  fStatusBar = new TGStatusBar(this, 50, 10, kVerticalFrame);
  fStatusBar->SetParts(parts, 4);
  fStatusBar->Draw3DCorner(kFALSE);
  AddFrame(fStatusBar, new TGLayoutHints(kLHintsExpandX, 0, 0, 10, 0));
  
  SetWindowName("Solenoid Analysis");
  MapSubwindows();
  Resize(GetDefaultSize());
  MapWindow();
  if (DoStartupFile){
    std::cout<<"Read file "<<InputFileList<<std::endl;
    DoOpen();
  }

  return;
}

void MyMainFrame::SetStatusText(const char *txt, Int_t pi)
{
  // Set text in status bar.
  fStatusBar->SetText(txt,pi);
  
  return;
}

void MyMainFrame::EventInfo(Int_t event, Int_t px, Int_t py, TObject *selected)
{
  //  Writes the event status in the status bar parts
  
  const char *text0, *text1, *text3;
  char text2[50];
  text0 = selected->GetTitle();
  SetStatusText(text0,0);
  text1 = selected->GetName();
  SetStatusText(text1,1);
  if (event == kKeyPress)
    sprintf(text2, "%c", (char) px);
  else
    sprintf(text2, "%d,%d", px, py);
  SetStatusText(text2,2);
  text3 = selected->GetObjectInfo(px,py);
  SetStatusText(text3,3);
  
  return;
}


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

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

  time_t t = time(0); // get current time 
  struct tm *now = localtime(&t);
  
  std::cout << std::endl;
  // get list of files from GUI
  TList myli;
  TGFileInfo fi;
  TObjString str1;
  fi.fFileTypes = filetypes;
  fi.fIniDir = StrDup(workingdir);
  fi.fMultipleSelection = true;
  TList *filelist;
  if (!DoStartupFile){
    new TGFileDialog(gClient->GetRoot(), this, kFDOpen, &fi);
    if (!fi.fFileNamesList) {
      std::cout << "No file selected!" << std::endl;
      return;
    }
    filelist = fi.fFileNamesList;
  } else {
    DoStartupFile = 0;
    str1.SetString(InputFileList);
    myli.Add(&str1);
    filelist = &myli;
  }
  
  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) {
      h_shunt = ReadTree(datatree, name_shunt, 1. / scale_shunt.coeff);
      if (!h_shunt)
	h_shunt = ReadTree(datatree, name_shunt_old,
			   1. / scale_shunt.coeff);
      h_voltage = ReadTree(datatree, name_voltage,
			   1. / scale_voltage.coeff);
      h_field = ReadTree(datatree, name_field, 1. / scale_field.coeff);
      
      // read voltage tabbs data
      for (int i = 0; i < ch_vtt; i++)
	h_vtt[i] = ReadTree(datatree, name_vtt[i],
			    1. / scale_vtt[i].coeff);


      for (int i = 0; i < ch_vtts; i++)
	h_vtts[i] = ReadTree(datatree, name_vtts[i],
			    1. / scale_vtts[i].coeff);


      for (int i = 0; i < ch_vttb; i++)
	h_vttb[i] = ReadTree(datatree, name_vttb[i],
			    1. / scale_vttb[i].coeff);

      // read pick up data
      for (int i = 0; i < ch_pickup; i++)
	h_pickup[i] = ReadTree(datatree, name_pickup[i],
			       1. / scale_pickup[i].coeff);
      // read accelerometer data 
      for (int i = 0; i < ch_acc; i++){
	h_acc[i] = ReadTree(datatree, name_acc[i], 1.);
      }
      
      h_qd = ReadTree(datatree, name_qd, 1. / scale_qd.coeff);
      if (!h_qd)
	h_qd = ReadTree(datatree, name_qd_old, 1. / scale_qd.coeff);
      h_gfd = ReadTree(datatree, name_gfd, 1. / scale_gfd.coeff);
      
      // version id, 0 w/o QD/GFD/VTT3.., 1 w/o GFD, 2 all
      if (!h_qd)
	version = 0;
      else {
	for (int i = 0; i < ch_lead; i++)
	  h_lead[i] = ReadTree(datatree, name_lead[i],
			       1. / scale_lead[i].coeff);
	if (!h_gfd)
	  version = 1;
	else
	  version = 2;
      }
      
      // 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("pxi", "");
	
	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 *current_data = new pxiData(filelist,
					name_prefix + name_current);
    if (current_data->entries == 0) {
      name_prefix = name_epics;
      current_data = new pxiData(filelist, name_prefix + name_current);
    }
    starttime = current_data->starttime;
    Double_t freq = 0;
    Int_t length = current_data->GetLength(freq);
    h_current = current_data->GetHistogram(starttime, length, freq);
    if (h_current) {
      pxiData *shunt_data = new pxiData(filelist,
					name_prefix + name_shunt);
      if (shunt_data->entries == 0)
	shunt_data = new pxiData(filelist,
				 name_prefix + name_shunt_old);
      h_shunt = shunt_data->GetHistogram(starttime, length, freq);
      pxiData *field_data = new pxiData(filelist,
					name_prefix + name_field);
      h_field = field_data->GetHistogram(starttime, length, freq);
      pxiData *voltage_data = new pxiData(filelist,
					  name_prefix + name_voltage);
      h_voltage = voltage_data->GetHistogram(starttime, length, freq);
      
      pxiData *vtt_data[ch_vtt], *pickup_data[ch_pickup], *acc_data[ch_acc],
	*lead_data[ch_lead], *vtts_data[ch_vtts], *vttb_data[ch_vttb];
      // read vtt voltage tab data
      for (int i = 0; i < ch_vtt; i++) {
	vtt_data[i] = new pxiData(filelist, name_prefix + name_vtt[i]);
	h_vtt[i] = vtt_data[i]->GetHistogram(starttime, length, freq);
      }
      // read vtts voltage tab data
      for (int i = 0; i < ch_vtts; i++) {
	vtts_data[i] = new pxiData(filelist, name_prefix + name_vtts[i]);
	h_vtts[i] = vtts_data[i]->GetHistogram(starttime, length, freq);
      }
      // read vtt voltage tab data
      for (int i = 0; i < ch_vttb; i++) {
	vttb_data[i] = new pxiData(filelist, name_prefix + name_vttb[i]);
	h_vttb[i] = vttb_data[i]->GetHistogram(starttime, length, freq);
      }
      // read pick-up coil data
      for (int i = 0; i < ch_pickup; i++) {
	pickup_data[i] = new pxiData(filelist,
				     name_prefix + name_pickup[i]);
	h_pickup[i] = pickup_data[i]->GetHistogram(starttime, length,
						   freq);
      }
      // read accelerometer data
      for (int i = 0; i < ch_acc; i++) {
	acc_data[i] = new pxiData(filelist,
				  name_prefix + name_acc[i]);
	h_acc[i] = acc_data[i]->GetHistogram(starttime, length,
					     freq);
      }
      
      pxiData *qd_data = new pxiData(filelist, name_prefix + name_qd);
      if (qd_data->entries == 0)
	qd_data = new pxiData(filelist, name_prefix + name_qd_old);
      pxiData *gfd_data = new pxiData(filelist, name_prefix + name_gfd);
      
      // version id
      if (qd_data->entries == 0) {
	version = 0;
	for (int i = 0; i < ch_pickup; i++) {
	  h_pickup[i]->Scale(pickup_pol_corr[i]);
	}
      } else {
	h_qd = qd_data->GetHistogram(starttime, length, freq);
	for (int i = 0; i < ch_lead; i++) {
	  lead_data[i] = new pxiData(filelist,
				     name_prefix + name_lead[i]);
	  h_lead[i] = lead_data[i]->GetHistogram(starttime, length,
						 freq);
	}
	if (gfd_data->entries == 0)
	  version = 1;
	else {
	  version = 2;
	  h_gfd = gfd_data->GetHistogram(starttime, length, freq);
	}
      }
      int offsetT = 5;
      std::cout<<"Daylight Savings time is: "<<now->tm_isdst<<std::endl;
      if (now->tm_isdst){
	offsetT = 4;
      }
      starttime.SetSec(starttime.GetSec() - offsetT * 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();
    
    // plot history
    if (h_shunt) {
      Double_t xmin = h_shunt->GetXaxis()->GetXmin();
      Double_t xmax = h_shunt->GetXaxis()->GetXmax();
      h_history = Resample(h_shunt, 100, 0, 0, kFALSE, kTRUE, xmin, xmax,
			   kBlack, scale_shunt);
      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");
	//			ComboPlot(fCanvas->GetSelectedPad(), &h_history, 1, "History", true);
	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, 3);

	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 first for all ! after thatt apply offsets and plot
	 */
// field and current
	TH1F *hr[30];
	Int_t nhr = 0;
	hr[nhr] = Resample(h_current, freq, 0, 0, kFALSE, kTRUE, xmin, xmax,
			col_current);
	Double_t dummy_error;
	offset_current = -Average(hr[nhr], dummy_error);
	nhr++;
	hr[nhr] = Resample(h_shunt, freq, 0, 0, kFALSE, kTRUE, xmin, xmax,
			col_shunt);
	offset_shunt = -Average(hr[nhr], dummy_error);
	nhr++;
	hr[nhr] = Resample(h_field, freq, 0, field_delay->GetNumber(), kFALSE,
			kTRUE, xmin, xmax, col_field);
	offset_field = -Average(hr[nhr], dummy_error);
	nhr++;
	ComboPlot(fCanvas->GetPad(1), hr, nhr, "Field and Current",
			check_legend->IsOn());

// voltage
	nhr = 0;
	hr[nhr] = Resample(h_voltage, freq, 0, 0, kFALSE, kTRUE, xmin, xmax,
			col_voltage);
	offset_voltage = -Average(hr[nhr], dummy_error);
	nhr++;

	nhr++;
	for (int i = 0; i < ch_vtt; i++) {
		hr[nhr] = Resample(h_vtt[i], freq, 0, field_delay->GetNumber(), kFALSE,
				kTRUE, xmin, xmax, col_vtt[i]);
		offset_vtt[i] = -Average(hr[nhr], dummy_error);
		nhr++;
	}
	if (version > 0) {
		hr[nhr] = Resample(h_qd, freq, 0, 0, kFALSE, kTRUE, xmin, xmax, col_qd);
		offset_qd = -Average(hr[nhr], dummy_error);
		nhr++;
		for (int i = 0; i < ch_lead; i++) {
			hr[nhr] = Resample(h_lead[i], freq, 0, field_delay->GetNumber(),
					kFALSE, kTRUE, xmin, xmax, col_lead[i]);
			offset_lead[i] = -Average(hr[nhr], dummy_error);
			nhr++;
		}
	} else {
		offset_qd = 0.;
		for (int i = 0; i < ch_lead; i++) {
			offset_lead[i] = 0.;
		}
	}
	if (version > 1) {
		hr[nhr] = Resample(h_gfd, freq, 0, 0, kFALSE, kTRUE, xmin, xmax,
				col_gfd);
		offset_gfd = -Average(hr[nhr], dummy_error);
		nhr++;
	} else
		offset_gfd = 0.;

        ComboPlot(fCanvas->GetPad(3), hr, nhr, "Voltage Taps",
                  check_legend->IsOn());

// voltage vtts
	for (int i = 0; i < ch_vtts; i++) {
		hr[nhr] = Resample(h_vtts[i], freq, 0, field_delay->GetNumber(), kFALSE,
				kTRUE, xmin, xmax, col_vtts[i]);
		offset_vtts[i] = -Average(hr[nhr], dummy_error);
		nhr++;
	}

// voltage vttb
	for (int i = 0; i < ch_vttb; i++) {
		hr[nhr] = Resample(h_vttb[i], freq, 0, field_delay->GetNumber(), kFALSE,
				kTRUE, xmin, xmax, col_vttb[i]);
		offset_vttb[i] = -Average(hr[nhr], dummy_error);
		nhr++;
	}

	//ComboPlot(fCanvas->GetPad(3), hr, nhr, "Voltage Taps",
        //        check_legend->IsOn());
	

// pickup coils
	nhr = 0;
	for (int i = 0; i < ch_pickup; i++) {
		hr[nhr] = Resample(h_pickup[i], freq, 0, field_delay->GetNumber(),
				kFALSE, kTRUE, xmin, xmax, col_pickup[i]);
		offset_pickup[i] = -Average(hr[nhr], dummy_error);
		nhr++;
	}
	ComboPlot(fCanvas->GetPad(5), hr, nhr, "Pickup Coils",
			check_legend->IsOn());

// Accelerometers:
	
	nhr = 0;
	for (int i = 0; i < ch_acc; i++) {
	  hr[nhr] = Resample(h_acc[i], freq, 0, field_delay->GetNumber(),
			     kFALSE, kTRUE, xmin, xmax, col_acc[i]);
	  offset_acc[i] = -Average(hr[nhr], dummy_error);
	  nhr++;
	}
	// do not plot accelerometers here and later.

	/*
	 * Now Plot signal after zeroing
	 */

// field and current
	nhr = 0;
	hr[nhr] = Resample(h_current, freq, offset_current, 0, kFALSE, kTRUE, xmin,
			xmax, col_current);
	nhr++;
	hr[nhr] = Resample(h_shunt, freq, offset_shunt, 0, kFALSE, kTRUE, xmin,
			xmax, col_shunt);
	nhr++;
	hr[nhr] = Resample(h_field, freq, offset_field, field_delay->GetNumber(),
			kFALSE, kTRUE, xmin, xmax, col_field);
	nhr++;
	ComboPlot(fCanvas->GetPad(2), hr, nhr);

// voltage tabs
	nhr = 0;
	hr[nhr] = Resample(h_voltage, freq, offset_voltage, 0, kFALSE, kTRUE, xmin, xmax,
			col_voltage);
	nhr++;

	for (int i = 0; i < ch_vtt; i++) {
	  hr[nhr] = Resample(h_vtt[i], freq, offset_vtt[i], field_delay->GetNumber(), kFALSE,
			     kTRUE, xmin, xmax, col_vtt[i]);
	  nhr++;
	}
	for (int i = 0; i < ch_vtts; i++) {
	  hr[nhr] = Resample(h_vtts[i], freq, offset_vtts[i], field_delay->GetNumber(), kFALSE,
			     kTRUE, xmin, xmax, col_vtts[i]);
	  nhr++;
	}
	for (int i = 0; i < ch_vttb; i++) {
	  hr[nhr] = Resample(h_vttb[i], freq, offset_vttb[i], field_delay->GetNumber(), kFALSE,
			     kTRUE, xmin, xmax, col_vttb[i]);
	  nhr++;
	}
	ComboPlot(fCanvas->GetPad(4), hr, nhr);

// pickup coils
	nhr = 0;
	for (int i = 0; i < ch_pickup; i++) {
		hr[nhr] = Resample(h_pickup[i], freq, offset_pickup[i], 0, kFALSE,
				kTRUE, xmin, xmax, col_pickup[i]);
		nhr++;
	}
	ComboPlot(fCanvas->GetPad(6), hr, nhr);

// accelerometers
	nhr = 0;
	for (int i = 0; i < ch_acc; i++) {
		hr[nhr] = Resample(h_acc[i], freq, offset_acc[i], 0, kFALSE,
				kTRUE, xmin, xmax, col_acc[i]);
		nhr++;
	}

	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()) {
      offsetfile << offset_current << "\n";
      offsetfile << offset_shunt << "\n";
      offsetfile << offset_field << "\n";
      offsetfile << offset_voltage << "\n";
      for (int i = 0; i < ch_vtt; i++)
	offsetfile << offset_vtt[i] << "\n";
      for (int i = 0; i < ch_pickup; i++)
	offsetfile << offset_pickup[i] << "\n";
      offsetfile << offset_qd << "\n";
      for (int i = 0; i < ch_lead; i++)
	offsetfile << offset_lead[i] << "\n";
      offsetfile << offset_gfd << "\n";
      for (int i = 0; i < ch_acc; i++)
	offsetfile << offset_acc[i] << "\n";
      for (int i = 0; i < ch_vtts; i++)
	offsetfile << offset_vtts[i] << "\n";
      for (int i = 0; i < ch_vttb; i++)
	offsetfile << offset_vttb[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()) {
			offsetfile >> offset_current;
			offsetfile >> offset_shunt;
			offsetfile >> offset_field;
			offsetfile >> offset_voltage;
			for (int i = 0; i < ch_vtt; i++)
				offsetfile >> offset_vtt[i];
			for (int i = 0; i < ch_pickup; i++)
				offsetfile >> offset_pickup[i];
			offsetfile >> offset_qd;
			for (int i = 0; i < ch_lead; i++)
				offsetfile >> offset_lead[i];
			offsetfile >> offset_gfd;
			for (int i = 0; i < ch_acc; i++)
				offsetfile >> offset_acc[i];
			for (int i = 0; i < ch_vtts; i++)
				offsetfile >> offset_vtts[i];
			for (int i = 0; i < ch_vttb; i++)
				offsetfile >> offset_vttb[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() {
	offset_current = 0.;
	offset_shunt = 0.;
	offset_field = 0;
	offset_voltage = 0;
	for (int i = 0; i < ch_vtt; i++) {
		offset_vtt[i] = 0.;
	}
	for (int i = 0; i < ch_pickup; i++) {
		offset_pickup[i] = 0.;
	}
	offset_qd = 0;
	for (int i = 0; i < ch_lead; i++) {
		offset_lead[i] = 0;
	}
	offset_gfd = 0;
	for (int i = 0; i < ch_acc; i++) {
		offset_acc[i] = 0;
	}
	for (int i = 0; i < ch_vtts; i++) {
		offset_vtts[i] = 0;
	}
	for (int i = 0; i < ch_vttb; i++) {
		offset_vttb[i] = 0;
	}
	if (isTree) {
		label_offset->ChangeText("Not Needed");
	} else {
		label_offset->ChangeText("Empty");
	}
	this->Layout();
}
/*
 * Action to make field plot and ratio
 */
void MyMainFrame::DoPlotAccelerometers() {
	if (!h_history)
		return;

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

	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 Accelerometer voltages
	TH1F *vhr[30];
	for (int i = 0; i < ch_acc; i++) {
		vhr[i] = Resample(h_acc[i], field_resample->GetNumber(), offset_acc[i],
				0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_acc[i],
				scale_acc[i]);
	}
	
	ComboPlot(fCanvas->GetPad(1), vhr, 2, "Accelerometers Coil 2",
		  check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
		  kFALSE, kFALSE, -10, 10);
	
	ComboPlot(fCanvas->GetPad(2), vhr+2, 2, "Accelerometers Coil 1",
		  check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
		  kFALSE, kFALSE, -10, 10);
	
	ComboPlot(fCanvas->GetPad(3), vhr+4, 2, "Accelerometers Coil 3",
		  check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
		  kFALSE, kFALSE, -10, 10);
	
	ComboPlot(fCanvas->GetPad(4), vhr+6, 2, "Accelerometers Coil 4",
		  check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
		  kFALSE, kFALSE, -10, 10);

	fCanvas->Update();
}


/*
 * Action to make field plot and ratio
 */
void MyMainFrame::DoPlotField() {

  if (!h_history)
    return;
  
  TCanvas *fCanvas = canvas_misc->GetCanvas();
  fCanvas->Clear();
  fCanvas->SetTitle("Field");
  fCanvas->Divide(1, 3);
  
  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 Current
  TH1F *ihr[30];
  Int_t nhr = 0;
  ihr[nhr] = Resample(h_current, field_resample->GetNumber(), offset_current,
		      0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_current,
		      scale_current);
  nhr++;
  ihr[nhr] = Resample(h_shunt, field_resample->GetNumber(), offset_shunt, 0,
		      check_denoise->IsOn(), kTRUE, xmin, xmax, col_shunt, scale_shunt);
  nhr++;
  ComboPlot(fCanvas->GetPad(1), ihr, nhr,
	    starttime.AsString("s") + TString(" Current"), check_legend->IsOn(),
	    check_error->IsOn(), check_mean->IsOn(), kFALSE, kFALSE, -10, 10);
  
  // Plot Field
  TH1F *fhr[30];
  nhr = 0;
  fhr[nhr] = Resample(h_field, field_resample->GetNumber(), offset_field,
		      field_delay->GetNumber(), check_denoise->IsOn(), kTRUE, xmin, xmax,
		      col_field, scale_field);
  nhr++;
  ComboPlot(fCanvas->GetPad(2), fhr, nhr, "Field", check_legend->IsOn(),
	    check_error->IsOn(), check_mean->IsOn(), kFALSE, kFALSE, -10, 10);
  
  // Plot Ratio
  TH1F *rhr[30];
  nhr = 0;
  //	rhr[nhr] = Ratio(fhr[nhr], ihr[0], scale_ratio);
  rhr[nhr] = Ratio(fhr[nhr], ihr[1], scale_ratio);
  nhr++;
  ComboPlot(fCanvas->GetPad(3), rhr, nhr, "Ratio of Field to Current",
	    check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
	    kFALSE, kFALSE, -10., 10.);
  
  fCanvas->Update();

  return;
}

/*
 * Action to make field plot and ratio
 */
void MyMainFrame::DoPlotRatio() {

  
  //std::cout<<"History: "<<h_history<<std::endl;

  if (!h_history)
    return;
  
  TCanvas *fCanvas = canvas_misc->GetCanvas();
  fCanvas->Clear();
  fCanvas->SetTitle("Ratio");
  fCanvas->Divide(1, 2);
  
  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_pickup; i++) {
    phr[nhr] = Resample(h_pickup[i], field_resample->GetNumber(),
			offset_pickup[i], 0, check_denoise->IsOn(), kTRUE, xmin, xmax,
			col_pickup[i], scale_pickup[i]);
    nhr++;
  }
  ComboPlot(fCanvas->GetPad(1), phr, nhr,
	    starttime.AsString("s") + TString(" Pickup Signal"),
	    check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
	    kFALSE, kFALSE, -10, 10);
  
  // Plot Pickup Integral
  TH1F *fhr[30];
  nhr = 0;
  for (int i = 0; i < ch_pickup; i++) {
    fhr[nhr] = IntHist(phr[i], scale_pickupfield);
    nhr++;
  }
  ComboPlot(fCanvas->GetPad(2), fhr, nhr, "Pickup Signal Integral",
	    check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
	    kFALSE, kFALSE, -10, 10);
  
  // Plot Ratio
  //	TH1F *ihr = Resample(h_current, field_resample->GetNumber(), offset_current,
  //			0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_current);
  TH1F *ihr = Resample(h_shunt, field_resample->GetNumber(), offset_shunt, 0,
		       check_denoise->IsOn(), kTRUE, xmin, xmax, col_shunt, scale_shunt);
  /*  
  TH1F *rhr[30];
  nhr = 0;
  for (int i = 0; i < ch_pickup; i++) {
    rhr[nhr] = Ratio(fhr[i], ihr, scale_ratio);
    nhr++;
  }
  ComboPlot(fCanvas->GetPad(3), rhr, nhr,
	    "Ratio of Pickup Integral to Current", check_legend->IsOn(),
	    check_error->IsOn(), check_mean->IsOn(), kFALSE, kFALSE, -10., 10.);
  */

  fCanvas->Update();

  return;
}

/*
 * 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, 3);
  
  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 VTT voltages
  TH1F *vhr[30];
  for (int i = 0; i < ch_vtt; i++) {
    vhr[i] = Resample(h_vtt[i], field_resample->GetNumber(), offset_vtt[i],
		      0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vtt[i],
		      scale_vtt[i]);
  }
  
  // Plot Total Voltage
  TH1F *vthr[30];
  Int_t nhr = 0;
  vthr[nhr] = Resample(h_voltage, field_resample->GetNumber(), offset_voltage,
		       0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_voltage,
		       scale_voltage);
  nhr++;
  vthr[nhr] = AddHist(vhr, ch_vtt, "hv_total", "Total Voltage", col_softvol);
  nhr++;
  
  ComboPlot(fCanvas->GetPad(1), vthr, nhr,
	    starttime.AsString("s") + TString(" Total Voltage"),
	    check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
	    kFALSE, kFALSE, -10, 10);
  
  ComboPlot(fCanvas->GetPad(2), vhr, 9, "Voltages in Coil 1 and 2",
	    check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
	    kFALSE, kFALSE, -10, 10);
  
  ComboPlot(fCanvas->GetPad(3), vhr + 9, ch_vtt - 9,
	    "Voltages in Coil 3 and 4", check_legend->IsOn(),
	    check_error->IsOn(), check_mean->IsOn(), kFALSE, kFALSE, -10, 10);
  
  fCanvas->Update();



}

/*
 * Action to make voltage plot for Coil 1 and 2
 */
void MyMainFrame::DoPlotVoltageC12() {

  if (!h_history)
    return;
  
  TCanvas *fCanvas = canvas_misc->GetCanvas();
  fCanvas->Clear();
  fCanvas->SetTitle("Voltage Coil 1_2");
  fCanvas->Divide(1, 2);
  
  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 VTT voltages
  TH1F *vhrC1[30];
  int C1 = 0;
  for (int i = 0; i < 4; i++) {
    vhrC1[i] = Resample(h_vtt[i], field_resample->GetNumber(), offset_vtt[i],
			0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vtt[i],
			scale_vtt[i]);
    C1++;
  }
  
  for (int i = 0; i < 2; i++) {
    vhrC1[C1] = Resample(h_vttb[i], field_resample->GetNumber(), offset_vttb[i],
			 0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vttb[i],
			 scale_vttb[i]);
    C1++;
  }
  
  
  TH1F *vhrC2[30];
  int C2 = 0;
  for (int i = 4; i < 9; i++) {
    vhrC2[C2++] = Resample(h_vtt[i], field_resample->GetNumber(), offset_vtt[i],
			   0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vtt[i],
			   scale_vtt[i]);
  }
  
  for (int i = 2; i < 4; i++) {
    vhrC2[C2++] = Resample(h_vttb[i], field_resample->GetNumber(), offset_vttb[i],
			   0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vttb[i],
			   scale_vttb[i]);
  }
  // Plot Total Voltage
  ComboPlot(fCanvas->GetPad(1), vhrC1, C1, "Voltages in Coil 2",
	    check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
	    kFALSE, kFALSE, -10, 10);
  
  ComboPlot(fCanvas->GetPad(2), vhrC2, C2,
	    "Voltages in Coil 1", check_legend->IsOn(),
	    check_error->IsOn(), check_mean->IsOn(), kFALSE, kFALSE, -10, 10);
  
  fCanvas->Update();


  return;
}


/*
 * Action to make voltage plot for Coil 3 and 4
 */
void MyMainFrame::DoPlotVoltageC34() {
	if (!h_history)
		return;

	TCanvas *fCanvas = canvas_misc->GetCanvas();
	fCanvas->Clear();
	fCanvas->SetTitle("Voltage Coil 3_4");
	fCanvas->Divide(1, 2);

	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 VTT voltages
	TH1F *vhrC3[30];
	int C3 = 0;
	for (int i = 9; i < 13; i++) {
		vhrC3[C3++] = Resample(h_vtt[i], field_resample->GetNumber(), offset_vtt[i],
				0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vtt[i],
				scale_vtt[i]);
	}

	for (int i = 4; i < 5; i++) {
		vhrC3[C3++] = Resample(h_vttb[i], field_resample->GetNumber(), offset_vttb[i],
				0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vttb[i],
				scale_vttb[i]);
	}

	TH1F *vhrC4[30];
	int C4 = 0;
	for (int i = 13; i < 16; i++) {
		vhrC4[C4++] = Resample(h_vtt[i], field_resample->GetNumber(), offset_vtt[i],
				0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vtt[i],
				scale_vtt[i]);
	}

	for (int i = 5; i < 7; i++) {
		vhrC4[C4++] = Resample(h_vttb[i], field_resample->GetNumber(), offset_vttb[i],
				0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vttb[i],
				scale_vttb[i]);
	}

// Plot Total Voltage
	ComboPlot(fCanvas->GetPad(1), vhrC3, C3, "Voltages in Coil 3",
			check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
			kFALSE, kFALSE, -10, 10);

	ComboPlot(fCanvas->GetPad(2), vhrC4, C4,
			"Voltages in Coil 4", check_legend->IsOn(),
			check_error->IsOn(), check_mean->IsOn(), kFALSE, kFALSE, -10, 10);

	fCanvas->Update();
}

/*
 * Action to make voltage plot for splices and leads
 */
void MyMainFrame::DoPlotVoltageLeads() {
	if (!h_history)
		return;

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

	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 VTT voltages
	TH1F *vhrC1[30];
	int C1=0;
	for (int i = 0; i < ch_vtts; i++) {
		vhrC1[C1++] = Resample(h_vtts[i], field_resample->GetNumber(), offset_vtts[i],
				0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vtts[i],
				scale_vtts[i]);
	}

	TH1F *vhrC2[30];
	int C2 = 0;
	for (int i = 0; i < ch_lead; i++) {
		vhrC2[C2++] = Resample(h_lead[i], field_resample->GetNumber(), offset_lead[i],
				0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_lead[i],
				scale_lead[i]);
	}

// Plot Total Voltage
	ComboPlot(fCanvas->GetPad(1), vhrC1, C1, "Voltages with Splices",
			check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
			kFALSE, kFALSE, -10, 10);

	ComboPlot(fCanvas->GetPad(2), vhrC2, C2,
			"Voltages in Leads", check_legend->IsOn(),
			check_error->IsOn(), check_mean->IsOn(), kFALSE, kFALSE, -10, 10);

	fCanvas->Update();
}

/*
 * Action to plot inductance
 */
void MyMainFrame::DoPlotInductance() {
	if (!h_history)
		return;

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

	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);

// Resampled VTT voltages
	TH1F *vhr[30];
	for (int i = 0; i < ch_vtt; i++) {
		vhr[i] = Resample(h_vtt[i], field_resample->GetNumber(), offset_vtt[i],
				0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vtt[i],
				scale_vtt[i]);
	}

// Total voltage and current
	TH1F *vtotal = AddHist(vhr, ch_vtt, "hv_total", "Total Voltage");
	TH1F *ihr;
	if (check_ZFCT->IsOn())
		ihr = Resample(h_current, field_resample->GetNumber(), offset_current,
				0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_current,
				scale_current);
	else
		ihr = Resample(h_shunt, field_resample->GetNumber(), offset_shunt, 0,
				check_denoise->IsOn(), kTRUE, xmin, xmax, col_shunt,
				scale_shunt);
	TH1F *ihd = Derivative(ihr, field_window->GetNumber());
	ihd->SetTitle("Current Derivative");
	ihd->GetYaxis()->SetTitle("Current Change Rate (A/s)");

// VTT voltages to total voltage, and ratio vs current
	TH1F *rhr[30];
	TH1F *rvi[30];
	for (int i = 0; i < ch_vtt; i++) {
		rhr[i] = Ratio(vhr[i], ihd, CreateScale(1, "Inductance (H)"));
		rvi[i] = CorrHist(rhr[i], ihr);
	}

	ComboPlot(fCanvas->GetPad(1), &ihd, 1,
			starttime.AsString("s") + TString(" Current Change Rate"),
			check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn());

	ComboPlot(fCanvas->GetPad(2), rvi, 9, "Inductance in Coil 1 and 2",
			check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
			check_fit->IsOn());

	ComboPlot(fCanvas->GetPad(3), rvi + 9, ch_vtt - 9,
			"Inductance in Coil 3 and 4", check_legend->IsOn(),
			check_error->IsOn(), check_mean->IsOn(), check_fit->IsOn());

	fCanvas->Update();
}

/*
 * Plot Quench Detector Error signal and Voltages across leads
 */
void MyMainFrame::DoPlotQD() {
	if (!h_history || version == 0)
		return;

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

	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 Quench Detector and Ground Fault Detector Error
	TH1F *dhr[2];
	int ndhr = 1;
	dhr[0] = Resample(h_qd, field_resample->GetNumber(), offset_qd, 0,
			check_denoise->IsOn(), kTRUE, xmin, xmax, col_qd, scale_qd);
	if (version > 1) {
		dhr[ndhr] = Resample(h_gfd, field_resample->GetNumber(), offset_gfd, 0,
				check_denoise->IsOn(), kTRUE, xmin, xmax, col_gfd, scale_gfd);
		ndhr++;
	}

// Plot Lead Voltages
	TH1F *lhr[ch_lead];
	for (int i = 0; i < ch_lead; i++) {
		lhr[i] = Resample(h_lead[i], field_resample->GetNumber(),
				offset_lead[i], 0, check_denoise->IsOn(), kTRUE, xmin, xmax,
				col_lead[i], scale_lead[i]);
	}

	ComboPlot(fCanvas->GetPad(1), dhr, ndhr,
			"Quench Detector and Ground Fault Detector Error Signals",
			check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
			kFALSE, kFALSE, -10, 10);

	ComboPlot(fCanvas->GetPad(2), lhr, ch_lead, "Voltages across Leads",
			check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
			kFALSE, kFALSE, -10, 10);

	fCanvas->Update();
}

/*
 * Plot resistance calculated from current, inductance and voltage
 */
void MyMainFrame::DoPlotResistance() {
	if (!h_history)
		return;
// find range
	TCanvas *fCanvas = canvas_misc->GetCanvas();
	fCanvas->Clear();
	fCanvas->SetTitle("Resistance");
	fCanvas->Divide(1, 3);
	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);

	TH1F *ihr = Resample(h_shunt, field_resample->GetNumber(), offset_shunt, 0,
			check_denoise->IsOn(), kTRUE, xmin, xmax, col_shunt, scale_shunt);
	TH1F *ihd[3];
	ihd[0] = Derivative(ihr, field_window->GetNumber());
	ihd[0]->SetTitle("Current Derivative");
	ihd[0]->GetYaxis()->SetTitle("Current Change Rate (A/s)");
	ihd[1] = Resample(h_pickup[0], field_resample->GetNumber(),
			offset_pickup[0], 0, check_denoise->IsOn(), kTRUE, xmin, xmax,
			col_pickup[0], scale_pickup[0]);
	ihd[1]->Scale(13.8 / 1.71);
	ihd[2] = Resample(h_pickup[1], field_resample->GetNumber(),
			offset_pickup[1], 0, check_denoise->IsOn(), kTRUE, xmin, xmax,
			col_pickup[1], scale_pickup[1]);
	ihd[2]->Scale(13.8 / 1.93);

// Resampled VTT voltages
	TH1F *vhr[30];
	for (int i = 0; i < ch_vtt; i++) {
		vhr[i] = Resample(h_vtt[i], field_resample->GetNumber(), offset_vtt[i],
				0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vtt[i],
				scale_vtt[i]);
	}

	TH1F *rhr[30];
	for (int i = 0; i < ch_vtt; i++) {
		TString hname = TString("resistance_") + name_vtt[i];
		rhr[i] = (TH1F*) gROOT->FindObject(hname);
		if (rhr[i])
			delete rhr[i];
		rhr[i] = (TH1F*) vhr[i]->Clone(hname);
		rhr[i]->GetYaxis()->SetTitle("Resistance (#Omega)");
		Int_t nbin = rhr[i]->GetNbinsX();
		for (int j = 1; j <= nbin; j++) {
			Double_t tmpy = 0;
			if (ihr->GetBinContent(j) != 0)
				tmpy = (vhr[i]->GetBinContent(j)
						- Inductance[i] * ihd[0]->GetBinContent(j))
						/ ihr->GetBinContent(j);
			rhr[i]->SetBinContent(j, tmpy);
		}
	}

	ComboPlot(fCanvas->GetPad(1), ihd, 3,
			starttime.AsString("s") + TString(" Current Change Rate"),
			check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn());

	ComboPlot(fCanvas->GetPad(2), rhr, 9, "Resistance in Coil 1 and 2",
			check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
			kFALSE, kFALSE, -10, 10);

	ComboPlot(fCanvas->GetPad(3), rhr + 9, ch_vtt - 9,
			"Resistance in Coil 3 and 4", check_legend->IsOn(),
			check_error->IsOn(), check_mean->IsOn(), kFALSE, kFALSE, -10, 10);

	fCanvas->Update();
}

/*
 * Plot temperature calculated from resistance
 */
void MyMainFrame::DoPlotTemperature() {
	if (!h_history)
		return;
// find range
	TCanvas *fCanvas = canvas_misc->GetCanvas();
	fCanvas->Clear();
	fCanvas->SetTitle("Temperature");
	fCanvas->Divide(1, 3);
	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);

	TH1F *ihr = Resample(h_shunt, field_resample->GetNumber(), offset_shunt, 0,
			check_denoise->IsOn(), kTRUE, xmin, xmax, col_shunt, scale_shunt);
	TH1F *ihd[3];
	ihd[0] = Derivative(ihr, field_window->GetNumber());
	ihd[0]->SetTitle("Current Derivative");
	ihd[0]->GetYaxis()->SetTitle("Current Change Rate (A/s)");
	ihd[1] = Resample(h_pickup[0], field_resample->GetNumber(),
			offset_pickup[0], 0, check_denoise->IsOn(), kTRUE, xmin, xmax,
			col_pickup[0], scale_pickup[0]);
	ihd[1]->Scale(13.8 / 1.71);
	ihd[2] = Resample(h_pickup[1], field_resample->GetNumber(),
			offset_pickup[1], 0, check_denoise->IsOn(), kTRUE, xmin, xmax,
			col_pickup[1], scale_pickup[1]);
	ihd[2]->Scale(13.8 / 1.93);

// Resampled VTT voltages
	TH1F *vhr[30];
	for (int i = 0; i < ch_vtt; i++) {
		vhr[i] = Resample(h_vtt[i], field_resample->GetNumber(), offset_vtt[i],
				0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vtt[i],
				scale_vtt[i]);
	}

	TH1F *rhr[30];
	for (int i = 0; i < ch_vtt; i++) {
		TString hname = TString("temperature_") + name_vtt[i];
		rhr[i] = (TH1F*) gROOT->FindObject(hname);
		if (rhr[i])
			delete rhr[i];
		rhr[i] = (TH1F*) vhr[i]->Clone(hname);
		rhr[i]->GetYaxis()->SetTitle("Temperature (K)");
		Int_t nbin = rhr[i]->GetNbinsX();
		for (int j = 1; j <= nbin; j++) {
			Double_t tmpy = (vhr[i]->GetBinContent(j)
					- Inductance[i] * ihd[0]->GetBinContent(j))
					/ ihr->GetBinContent(j);
			tmpy = Temperature(tmpy / Resistance[i]);
			rhr[i]->SetBinContent(j, tmpy);
		}
	}

	ComboPlot(fCanvas->GetPad(1), &ihr, 1,
			starttime.AsString("s") + TString(" Magnet Current"),
			check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn());

	ComboPlot(fCanvas->GetPad(2), rhr, 9, "Temperature in Coil 1 and 2",
			check_legend->IsOn(), check_error->IsOn(), check_mean->IsOn(),
			kFALSE, kFALSE, -10, 10);

	ComboPlot(fCanvas->GetPad(3), rhr + 9, ch_vtt - 9,
			"Temperature in Coil 3 and 4", check_legend->IsOn(),
			check_error->IsOn(), check_mean->IsOn(), kFALSE, kFALSE, -10, 10);

	fCanvas->Update();
	return;
}

/*
 * 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[60];
  Float_t value[60];
  Double_t time;
  
  tree->Branch("time", &time, "time/D");
  
  Int_t nhr = 0;
  std::cout << "start data resampling ..." << std::endl;
  // Resampled Current, V_MPS, Field
  hr[nhr] = Resample(h_current, field_resample->GetNumber(), offset_current,
		     0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_current,
		     scale_current);
  tree->Branch(name_current, value + nhr, name_current + TString("/F"));
  nhr++;
  
  hr[nhr] = Resample(h_shunt, field_resample->GetNumber(), offset_shunt, 0,
		     check_denoise->IsOn(), kTRUE, xmin, xmax, col_shunt, scale_shunt);
  tree->Branch(name_shunt, value + nhr, name_shunt + TString("/F"));
  nhr++;
  
  hr[nhr] = Resample(h_voltage, field_resample->GetNumber(), offset_voltage,
		     0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_voltage,
		     scale_voltage);
  tree->Branch(name_voltage, value + nhr, name_voltage + TString("/F"));
  nhr++;
  
  hr[nhr] = Resample(h_field, field_resample->GetNumber(), offset_field, 0,
		     check_denoise->IsOn(), kTRUE, xmin, xmax, col_field, scale_field);
  tree->Branch(name_field, value + nhr, name_field + TString("/F"));
  nhr++;
  
  // Resampled VTT voltages
  for (int i = 0; i < ch_vtt; i++) {
    hr[nhr] = Resample(h_vtt[i], field_resample->GetNumber(), offset_vtt[i],
		       0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vtt[i],
		       scale_vtt[i]);
    tree->Branch(name_vtt[i], value + nhr, name_vtt[i] + TString("/F"));
    nhr++;
  }
  
  // Resampled VTTS and VTTB voltage
  for (int i = 0; i < ch_vtts; i++) {
    hr[nhr] = Resample(h_vtts[i], field_resample->GetNumber(), offset_vtts[i],
		       0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vtts[i],
		       scale_vtts[i]);
    tree->Branch(name_vtts[i], value + nhr, name_vtts[i] + TString("/F"));
    nhr++;
  }

  for (int i = 0; i < ch_vttb; i++) {
    hr[nhr] = Resample(h_vttb[i], field_resample->GetNumber(), offset_vttb[i],
		       0, check_denoise->IsOn(), kTRUE, xmin, xmax, col_vttb[i],
		       scale_vttb[i]);
    tree->Branch(name_vttb[i], value + nhr, name_vttb[i] + TString("/F"));
    nhr++;
  }

  // Resampled Pickup coil signals
  for (int i = 0; i < ch_pickup; i++) {
    hr[nhr] = Resample(h_pickup[i], field_resample->GetNumber(),
		       offset_pickup[i], 0, check_denoise->IsOn(), kTRUE, xmin, xmax,
		       col_pickup[i], scale_pickup[i]);
    tree->Branch(name_pickup[i], value + nhr,
		 name_pickup[i] + TString("/F"));
    nhr++;
  }
  
  // Resampled Accelerometer signals
  for (int i = 0; i < ch_acc; i++) {
    hr[nhr] = Resample(h_acc[i], field_resample->GetNumber(),
		       offset_acc[i], 0, check_denoise->IsOn(), kTRUE, xmin, xmax,
		       col_acc[i], scale_acc[i]);
    tree->Branch(name_acc[i], value + nhr,
		 name_pickup[i] + TString("/F"));
    nhr++;
  }
  
  if (version > 0) {
    hr[nhr] = Resample(h_qd, field_resample->GetNumber(), offset_qd, 0,
		       check_denoise->IsOn(), kTRUE, xmin, xmax, col_qd, scale_qd);
    tree->Branch(name_qd, value + nhr, name_qd + TString("/F"));
    nhr++;
    // Resampled Leads Voltages
    for (int i = 0; i < ch_lead; i++) {
      hr[nhr] = Resample(h_lead[i], field_resample->GetNumber(),
			 offset_lead[i], 0, check_denoise->IsOn(), kTRUE, xmin, xmax,
			 col_lead[i], scale_lead[i]);
      tree->Branch(name_lead[i], value + nhr,
		   name_lead[i] + TString("/F"));
      nhr++;
    }
  }
  if (version > 1) {
    hr[nhr] = Resample(h_gfd, field_resample->GetNumber(), offset_gfd, 0,
		       check_denoise->IsOn(), kTRUE, xmin, xmax, col_gfd, scale_gfd);
    tree->Branch(name_gfd, value + nhr, name_gfd + TString("/F"));
    nhr++;
  }
  
  std::cout << "finished resampling data and now saving to a ROOT Tree ..."<< nhr<<"  "
	    << 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");
  /*	TObjTimeStamp *ts2 = (TObjTimeStamp*) gROOT->FindObject("StartTime");
	ts2->GetTime().Print();*/
  rootfile->Close();
  std::cout << "rootfile " << filename << " has been created." << std::endl;

  return;
}

/*
 * 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 fnam = fCanvas->GetTitle();
	fnam.ReplaceAll(" ","_");
	TString pdffile = plotdir + pngfile + "_" + fnam + ".pdf";
	pngfile = plotdir + pngfile + "_" + fnam + ".png";
	fCanvas->Print(pngfile);
	fCanvas->Print(pdffile);
	std::cout<<"Printing done."<<std::endl;
	return;
}

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 10 Hz
  base_frequency = floor(base_frequency / 10. + 0.5) * 10.;
  // 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");
  hh[0]->GetYaxis()->SetLabelSize(0.1);
  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.;
        // instead of using window to determine fit range use bin width
        Double_t frange = h1->GetBinWidth(1)*1.5;

	TF1 *f1 = new TF1("f1", "pol1", 0, 10);
	for (int i = 3; i <= nbin-2; 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;
}

/*
 * Temperature calculated from resistance ratio to room temperature
 */
Double_t Temperature(Double_t ratio) {
	const Int_t tbin = 23;
	Double_t t_array[tbin] = { 4, 10, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80,
			100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300 };
	Double_t r_array[tbin] = { 0.0000, 0.0032, 0.0104, 0.0146, 0.0201, 0.0267,
			0.0345, 0.0433, 0.0532, 0.0759, 0.1022, 0.1320, 0.2004, 0.2784,
			0.3637, 0.4536, 0.5455, 0.6371, 0.7257, 0.8088, 0.8839, 0.9485,
			1.0000 };
	if (ratio <= r_array[0])
		return t_array[0];
	Double_t tmpt = t_array[tbin - 1];
	for (int i = 0; i < tbin; i++) {
		if (ratio < r_array[i]) {
			tmpt = t_array[i]
					+ (ratio - r_array[i]) * (t_array[i - 1] - t_array[i])
							/ (r_array[i - 1] - r_array[i]);
			break;
		}
	}
	return tmpt;
}

/*
 * 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) {
  DoStartupFile = 0;
  for (int k=0; k<argc; k++){
    if (!strcmp(argv[k],"-F")){
      DoStartupFile = 1;
      sprintf(InputFileList,"%s",argv[++k]);
    } else if (!strcmp(argv[k],"-f")){
      DoStartupFile = 1;
      sprintf(InputFileList,"%s",argv[++k]);
    }
  }

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