/// to do:
/// 1. correct for gap between arrays

#define convert2_cxx
#include "hitutils.h"
#include <string>
#include <stdio.h>
#include <iostream>
#include <vector>
#include <utility>
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TTree.h>
#include <TSystemDirectory.h>
#include <gsl/gsl_statistics.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <gsl/gsl_sort.h>
#include <TVector.h>
#include "TStopwatch.h"
#include "Math/MinimizerOptions.h"
#include "TVirtualFitter.h"

using namespace std;

bool smooth_on = false;
bool doFit = false;


int eventID = 0;
int framestart=0;
int frameend = 0;
double board_b0_ch[128];
double board_b1_ch[128];
double board_b2_ch[128];
double board_b3_ch[128];
double board_b0_ch_bkg[128];
double board_b1_ch_bkg[128];
double board_b2_ch_bkg[128];
double board_b3_ch_bkg[128];
double signal_b0 = 0.;
double signal_b1 = 0.;
double signal_b2 = 0.;
double signal_b3 = 0.;

double signal_b0_left = 0.;
double signal_b1_left = 0.;
double signal_b2_left = 0.;
double signal_b3_left = 0.;

double signal_b0_right = 0.;
double signal_b1_right = 0.;
double signal_b2_right = 0.;
double signal_b3_right = 0.;

double signal_gsl_b0[128];
double signal_gsl_b1[128];
double signal_gsl_b2[128];
double signal_gsl_b3[128];


double channellist_gsl_b0[128];
double channellist_gsl_b1[128];
double channellist_gsl_b2[128];
double channellist_gsl_b3[128];
double channellist[128];
double pos[128];

double maxchannelamp_b0 =0.;
double maxchannelamp_b1 =0.;
double maxchannelamp_b2 =0.;
double maxchannelamp_b3 =0.;
int maxchannel_b0 =0;
int maxchannel_b1 =0;
int maxchannel_b2 =0;
int maxchannel_b3 =0;
int numtocalc_b0 =0;
int numtocalc_b1 =0;
int numtocalc_b2 =0;
int numtocalc_b3 =0;
int nfwhm_b0 = 0;
int nfwhm_b1 = 0;
int nfwhm_b2 = 0;
int nfwhm_b3 = 0;

double	beamPosX_b0,beamPosX_b1,beamPosX_b2,beamPosX_b3;
double	beamFocusX_b0,beamFocusX_b1,beamFocusX_b2,beamFocusX_b3;
double	beamPosX_fit_b0,beamPosX_fit_b1,beamPosX_fit_b2,beamPosX_fit_b3;
double	beamFocusX_fit_b0,beamFocusX_fit_b1,beamFocusX_fit_b2,beamFocusX_fit_b3; 
double	beamSkewX_b0,beamSkewX_b1,beamSkewX_b2,beamSkewX_b3;
double	beamKurtX_b0,beamKurtX_b1,beamKurtX_b2,beamKurtX_b3; 
double	beamNumX_b0,beamNumX_b1,beamNumX_b2,beamNumX_b3;

double beamSidebandNoise_b0, beamSidebandNoise_b1, beamSidebandNoise_b2, beamSidebandNoise_b3 ;
int sidenumtocalc_b0,sidenumtocalc_b1,sidenumtocalc_b2,sidenumtocalc_b3;
size_t size = 5;
vector<double> boxcar;
vector<double> boxcarsort;
vector<double> boxcarweight{1,3,5,3,1};
vector<double> data(128);
TVector  sumvector_b0(128);
TVector  sumvector_b1(128);
TVector  sumvector_b2(128);
TVector  sumvector_b3(128);

TVector * sumvector_b0_ptr = &sumvector_b0; 
TVector * sumvector_b1_ptr =  &sumvector_b1; 
TVector * sumvector_b2_ptr =  &sumvector_b2;
TVector * sumvector_b3_ptr = &sumvector_b3; 

const int length = 100; //length of the rolling average
double array_b0[length][128] = {{0.}};
double array_b1[length][128] = {{0.}};
double array_b2[length][128] = {{0.}};
double array_b3[length][128] = {{0.}};

double arrayavg_b0[128] = {0.};
double arrayavg_b1[128] = {0.};
double arrayavg_b2[128] = {0.};
double arrayavg_b3[128] = {0.};



double board_b0_smooth_ch[128];
double board_b1_smooth_ch[128];
double board_b2_smooth_ch[128];
double board_b3_smooth_ch[128];	    

bool graphsaved_b0 = false;
bool graphsaved_b1 = false;
bool graphsaved_b2 = false;
bool graphsaved_b3 = false;

TVector beamontime(0,3,0.,0.,0.,0.,"END");
TVector * beamontime_ptr = &beamontime;

double calibration_b0[128] = {0.};
double calibration_b1[128] = {0.};
double calibration_b2[128] = {0.};
double calibration_b3[128] = {0.};

TF1 * gausfunc_b0 = new TF1("gausfunc_b0","gaus(0)+[3]");
TF1 * gausfunc_b1 = new TF1("gausfunc_b1","gaus(0)+[3]");



TGraphErrors * gausgraph_b0;
TGraphErrors * gausgraph_b1;
Int_t lastfit_b0 = 0;
Int_t lastfit_b1 = 0;
double errorx[128];
double errory[128];
      


/// compile with:
//// $ make clean; make
/// run with:
//// $ ./convert <path to data> <run>
//// $ ./convert /work/leverington/beamprofilemonitor/hitdata/HIT_17_12_2017/ run2



int analyse(int argc, char **argv)
{

  TVirtualFitter::SetDefaultFitter("Minuit2");
  TVirtualFitter::SetPrecision(0.01);

  std::cout << "Default Fitter:" << TVirtualFitter::GetDefaultFitter() << std::endl;

  //initialise some values;
  for (int i = 0;i<128;i++){
    board_b0_ch_bkg[i] = 0.;
    board_b1_ch_bkg[i] = 0.;
    board_b2_ch_bkg[i] = 0.;
    board_b3_ch_bkg[i] = 0.;
    channellist_gsl_b0[i] = 0;
    channellist_gsl_b1[i] = 0;
    channellist_gsl_b2[i] = 0;
    channellist_gsl_b3[i] = 0;
    sumvector_b0[i] = 0;
    sumvector_b1[i] = 0;
    sumvector_b2[i] = 0;
    sumvector_b3[i] = 0;
    calibration_b0[i] = 1.0;
    calibration_b1[i] = 1.0;
    calibration_b2[i] = 1.0;
    calibration_b3[i] = 1.0;
    errorx[i]=17.;
    errory[i]=17.;

    if (i<64) {pos[i] = double(i);}
    else {pos[i] = double(i) + 0.2; }
  }
  if (argc > 1){
    //open the file names specified in the command line
    string ethercatfile = argv[3];
    //  argv[1]= "/work/leverington/beamprofilemonitor/hitdata/HIT_17_12_2017/";
    //  argv[2]= "Run11"
    string filename = Form("%s%s.dat",argv[1],argv[2]);
    string timestampfilename =  Form("%s%s_timestamp.csv",argv[1],argv[2]); 
    string rootfilename = Form("%s/root/%s.root",argv[1],argv[2]); // Giulia needs to change the path %s/root/ of the written root file. 
    ifstream file(filename, ifstream::in | ifstream::binary);
    ifstream timestampfile(timestampfilename, ifstream::in);
    if (!file.is_open()){
      printf(".dat did not open.\n");
      return -1;	//file could not be opened
    }
        if (!timestampfile.is_open()){
      printf("timestamp.csv did not open.\n");
      return -1;	//file could not be opened
    }
    
    
   
    ///some variables to get the data
    // int number_of_records_to_read = 4*10;
    //   BufferData* buffer = new BufferData[number_of_records_to_read];
    //if(argc > 4)  {
    //  framestart = atoi(argv[3]);
    //  frameend = atoi(argv[4]);
	// }


    int board_b = -1;
    long int  fileframesize =   getFileSize(filename.c_str()) / ( 4*sizeof(BufferData) );
    BufferData* dataptr = new BufferData(); 
    
    if (fileframesize>0){
      std::cout << "Number of frames:" << fileframesize << std::endl;
      //open ROOT file for saving histos and TTree.
      TFile *rootFile = new TFile(rootfilename.c_str(),"RECREATE");
      if ( rootFile->IsOpen() ) printf("ROOT file opened successfully\n");
      TH2D * th2_signal_vs_channel_b0 = new TH2D("th2_signal_vs_channel_b0","th2_signal_vs_channel_b0",128,0,128,2200,-1000,10000);
      TH2D * th2_signal_vs_channel_b1 = new TH2D("th2_signal_vs_channel_b1","th2_signal_vs_channel_b1",128,0,128,2200,-1000,10000);
      TH2D * th2_signal_vs_channel_b2 = new TH2D("th2_signal_vs_channel_b2","th2_signal_vs_channel_b2",128,0,128,2200,-1000,10000);
      TH2D * th2_signal_vs_channel_b3 = new TH2D("th2_signal_vs_channel_b3","th2_signal_vs_channel_b3",128,0,128,2200,-1000,10000);
      TH2D * th2_signal_vs_channel_sub_b0 = new TH2D("th2_signal_vs_channel_sub_b0","th2_signal_vs_channel_sub_b0",128,0,128,2200,-1000,10000);
      TH2D * th2_signal_vs_channel_sub_b1 = new TH2D("th2_signal_vs_channel_sub_b1","th2_signal_vs_channel_sub_b1",128,0,128,2200,-1000,10000);
      TH2D * th2_signal_vs_channel_sub_b2 = new TH2D("th2_signal_vs_channel_sub_b2","th2_signal_vs_channel_sub_b2",128,0,128,2200,-1000,10000);
      TH2D * th2_signal_vs_channel_sub_b3 = new TH2D("th2_signal_vs_channel_sub_b3","th2_signal_vs_channel_sub_b3",128,0,128,2200,-1000,10000);
      TH1D * th1_signal_b0 = new TH1D("th1_signal_b0","th1_signal_b0",2200,-10000,50000);
      TH1D * th1_signal_b1 = new TH1D("th1_signal_b1","th1_signal_b1",2200,-10000,50000);
      TH1D * th1_signal_b2 = new TH1D("th1_signal_b2","th1_signal_b2",2200,-10000,50000);
      TH1D * th1_signal_b3 = new TH1D("th1_signal_b3","th1_signal_b3",2200,-10000,50000);
      TGraph * graph_bkg_b0 = new TGraph();
      TGraph * graph_bkg_b1 = new TGraph();
      TGraph * graph_bkg_b2 = new TGraph();
      TGraph * graph_bkg_b3 = new TGraph();
      TGraph * gr_b0;
      TGraph * gr_sm_b0;
      TGraph * gr_b1;
      TGraph * gr_sm_b1;
      TGraph * gr_b2;
      TGraph * gr_sm_b2;
      TGraph * gr_b3;
      TGraph * gr_sm_b3;

      TH2D * th2_signal_vs_channel_bkg_b0 = new TH2D("th2_signal_vs_channel_bkg_b0","th2_signal_vs_channel_bkg_b0",128,0,128,1000,-500,500);
      TH2D * th2_signal_vs_channel_bkg_b1 = new TH2D("th2_signal_vs_channel_bkg_b1","th2_signal_vs_channel_bkg_b1",128,0,128,1000,-500,500);
      TH2D * th2_signal_vs_channel_bkg_b2 = new TH2D("th2_signal_vs_channel_bkg_b2","th2_signal_vs_channel_bkg_b2",128,0,128,1000,-500,500);
      TH2D * th2_signal_vs_channel_bkg_b3 = new TH2D("th2_signal_vs_channel_bkg_b3","th2_signal_vs_channel_bkg_b3",128,0,128,1000,-500,500);

      TTree *rootTree = new TTree("t","HIT Data Root Tree");
      rootTree->Branch("beamPosX_b0",&beamPosX_b0,"beamPosX_b0/D");
      rootTree->Branch("beamPosX_b1",&beamPosX_b1,"beamPosX_b1/D");
      rootTree->Branch("beamPosX_b2",&beamPosX_b2,"beamPosX_b2/D");
      rootTree->Branch("beamPosX_b3",&beamPosX_b3,"beamPosX_b3/D");

      rootTree->Branch("beamPosX_fit_b0",&beamPosX_fit_b0,"beamPosX_fit_b0/D");
      rootTree->Branch("beamPosX_fit_b1",&beamPosX_fit_b1,"beamPosX_fit_b1/D");
      rootTree->Branch("beamPosX_fit_b2",&beamPosX_fit_b2,"beamPosX_fit_b2/D");
      rootTree->Branch("beamPosX_fit_b3",&beamPosX_fit_b3,"beamPosX_fit_b3/D");
      
      rootTree->Branch("beamFocusX_b0",&beamFocusX_b0,"beamFocusX_b0/D");
      rootTree->Branch("beamFocusX_b1",&beamFocusX_b1,"beamFocusX_b1/D");
      rootTree->Branch("beamFocusX_b2",&beamFocusX_b2,"beamFocusX_b2/D");
      rootTree->Branch("beamFocusX_b3",&beamFocusX_b3,"beamFocusX_b3/D");

      rootTree->Branch("beamFocusX_fit_b0",&beamFocusX_fit_b0,"beamFocusX_fit_b0/D");
      rootTree->Branch("beamFocusX_fit_b1",&beamFocusX_fit_b1,"beamFocusX_fit_b1/D");
      rootTree->Branch("beamFocusX_fit_b2",&beamFocusX_fit_b2,"beamFocusX_fit_b2/D");
      rootTree->Branch("beamFocusX_fit_b3",&beamFocusX_fit_b3,"beamFocusX_fit_b3/D");


      
      rootTree->Branch("beamSkewX_b0",&beamSkewX_b0,"beamSkewX_b0/D");
      rootTree->Branch("beamSkewX_b1",&beamSkewX_b1,"beamSkewX_b1/D");
      rootTree->Branch("beamSkewX_b2",&beamSkewX_b2,"beamSkewX_b2/D");
      rootTree->Branch("beamSkewX_b3",&beamSkewX_b3,"beamSkewX_b3/D");
      rootTree->Branch("beamKurtX_b0",&beamKurtX_b0,"beamKurtX_b0/D");
      rootTree->Branch("beamKurtX_b1",&beamKurtX_b1,"beamKurtX_b1/D");
      rootTree->Branch("beamKurtX_b2",&beamKurtX_b2,"beamKurtX_b2/D");
      rootTree->Branch("beamKurtX_b3",&beamKurtX_b3,"beamKurtX_b3/D");
      rootTree->Branch("beamNumX_b0",&beamNumX_b0,"beamNumX_b0/D");
      rootTree->Branch("beamNumX_b1",&beamNumX_b1,"beamNumX_b1/D");
      rootTree->Branch("beamNumX_b2",&beamNumX_b2,"beamNumX_b2/D");
      rootTree->Branch("beamNumX_b3",&beamNumX_b3,"beamNumX_b3/D");
      rootTree->Branch("beamSignal_b0",&signal_b0,"beamSignal_b0/D");
      rootTree->Branch("beamSignal_b1",&signal_b1,"beamSignal_b1/D");
      rootTree->Branch("beamSignal_b2",&signal_b2,"beamSignal_b2/D");
      rootTree->Branch("beamSignal_b3",&signal_b3,"beamSignal_b3/D");

      rootTree->Branch("beamSignal_b0_left",&signal_b0_left,"beamSignal_b0_left/D");
      rootTree->Branch("beamSignal_b1_left",&signal_b1_left,"beamSignal_b1_left/D");
      rootTree->Branch("beamSignal_b2_left",&signal_b2_left,"beamSignal_b2_left/D");
      rootTree->Branch("beamSignal_b3_left",&signal_b3_left,"beamSignal_b3_left/D");
      rootTree->Branch("beamSignal_b0_right",&signal_b0_right,"beamSignal_b0_right/D");
      rootTree->Branch("beamSignal_b1_right",&signal_b1_right,"beamSignal_b1_right/D");
      rootTree->Branch("beamSignal_b2_right",&signal_b2_right,"beamSignal_b2_right/D");
      rootTree->Branch("beamSignal_b3_right",&signal_b3_right,"beamSignal_b3_right/D");

      rootTree->Branch("eventID",&eventID,"eventID/I");
      //rootTree->Branch("board_b0_ch",&board_b0_ch,"board_b0_ch[128]/D");
      //rootTree->Branch("board_b1_ch",&board_b1_ch,"board_b1_ch[128]/D");
      //rootTree->Branch("board_b2_ch",&board_b2_ch,"board_b2_ch[128]/D");
      //rootTree->Branch("board_b3_ch",&board_b3_ch,"board_b3_ch[128]/D");

      rootTree->Branch("beamSidebandNoise_b0",&beamSidebandNoise_b0,"beamSidebandNoise_b0/D");
      rootTree->Branch("beamSidebandNoise_b1",&beamSidebandNoise_b1,"beamSidebandNoise_b1/D");
      rootTree->Branch("beamSidebandNoise_b2",&beamSidebandNoise_b2,"beamSidebandNoise_b2/D");
      rootTree->Branch("beamSidebandNoise_b3",&beamSidebandNoise_b3,"beamSidebandNoise_b3/D");
      //rootTree->Branch("arrayavg_b0",&arrayavg_b0,"arrayavg_b0[128]/D");
      //rootTree->Branch("arrayavg_b1",&arrayavg_b1,"arrayavg_b1[128]/D");
      // rootTree->Branch("arrayavg_b2",&arrayavg_b2,"arrayavg_b2[128]/D");
      //rootTree->Branch("arrayavg_b3",&arrayavg_b3,"arrayavg_b3[128]/D");
      double time;
      
      // import and align matching ethercat data
      TTree *tree2 = new TTree("t2", "t2");
      std::cout << " Loading Ethercat data." << std::endl;
      tree2->ReadFile(ethercatfile.c_str(), "RELTIME2/D:IC1/D:MW1_POSX/D:MW1_POSY/D:ANALOG_IN1/D:ENERGY/D:INTENSITY/D:ION-SORT/D:TIME2/D", '\t');
      std::cout << "Ethercat data loaded." << std::endl;
      tree2->Print();
      int k = 0;
      // int count = 0;
      double ic1, ic2, mw1_focusx, mw1_focusy, mw2_focusx, mw2_focusy, mw1_posx, mw1_posy, mw2_posx, mw2_posy;
      double ic1_avg, ic2_avg, mw1_focusx_avg, mw1_focusy_avg, mw2_focusx_avg, mw2_focusy_avg, mw1_posx_avg, mw1_posy_avg, mw2_posx_avg, mw2_posy_avg;
      
      double analog_in1;
      double energy;
      double intensity;
      double ionsort;

      double rel_time2,time2;

      double timeoffset;
      int mwoffset;
      double timewindow;
      double timeoffset2;
      double timewindow2;


      rootTree->Branch("ic1_avg", &ic1_avg, "ic1_avg/D");
      rootTree->Branch("mw1_posx", &mw1_posx_avg, "mw1_posx/D");
      rootTree->Branch("mw1_posy", &mw1_posy_avg, "mw1_posy/D");
      rootTree->Branch("energy", &energy, "energy/D");
      rootTree->Branch("intensity", &intensity, "intensity/D");
      rootTree->Branch("ionsort", &ionsort, "ionsort/D");


      timeoffset = 0.10; //offset between ic and bpm readouts
      mwoffset = 2; // offset for timestamped event. MW chamber position correlation seems to be better in other windows
      timewindow = -0.0999; //should be a negative number. window size in time to average over.
      timeoffset2 = -0.00; //offset between ic and bpm readouts
      timewindow2 = -0.05; //should be a negative number. window size in time to average over.
      //////////////////////////////////////////////////////////////////
      ///// ETHERCAT DATA
      //////////////////////////////////////////////////////////////////
      //		tree->SetBranchAddress("time", &time);

      tree2->SetBranchAddress("RELTIME2", &rel_time2);
      tree2->SetBranchAddress("TIME2", &time2);
      tree2->SetBranchAddress("IC1", &ic1);
      //    tree2->SetBranchAddress("IC2", &ic2);
      //tree2->SetBranchAddress("MW1_FOCUSX", &mw1_focusx);
      //tree2->SetBranchAddress("MW1_FOCUSY", &mw1_focusy);
      ///tree2->SetBranchAddress("MW2_FOCUSX", &mw2_focusx);
      ///tree2->SetBranchAddress("MW2_FOCUSY", &mw2_focusy);
      tree2->SetBranchAddress("MW1_POSX", &mw1_posx);
      tree2->SetBranchAddress("MW1_POSY", &mw1_posy);
      //tree2->SetBranchAddress("MW2_POSX", &mw2_posx);
      //tree2->SetBranchAddress("MW2_POSY", &mw2_posy);
      tree2->SetBranchAddress("ENERGY", &energy);
      tree2->SetBranchAddress("INTENSITY", &intensity);
      tree2->SetBranchAddress("ION-SORT", &ionsort);
      tree2->SetBranchAddress("ANALOG_IN1", &analog_in1);





      int currentEntryTree2 = 1;

      //	int nevents = tree->GetEntries();
      int nevents2 = tree2->GetEntries();

      int icCounter = 0;
      
      int count = 0;
      int count2 = 0;
      
      //loop through recorded data frames
      for (int frame = 0; frame<fileframesize; frame++){
	//	if (icCounter>10000) break;
	eventID = frame;
	if (timestampfile) timestampfile >> time ;  
	//printf("%i %f\n", eventID, time); 
	
	count= 0;    
	count2 = 0;
	ic1_avg = 0.;
	mw1_posx_avg = 0.;
	mw1_posy_avg = 0.; 
	tree2->GetEntry(currentEntryTree2);
	
        
	/*  if (frame % 100 == 0)
	  {

	  printf("merging event %d ,", frame);
	  printf("Time %f \n", time);
	  printf("Entry hit %d ,", currentEntryTree2);
	  printf("Time hit %f \n", time2);
	  }
	*/
	
	while (time2 <  time + timeoffset && currentEntryTree2 < nevents2 )
	  {
	    if (time2 - time - timeoffset >  timewindow)
	      {
		tree2->GetEntry(currentEntryTree2);
			    
		if (ic1>0.0) ic1_avg += ic1;
			    
		if (ic1>0.0) count++;
		//	        if (frame % 200 == 0) printf("%i %2.3f %2.3f  %2.3f  %2.3f %2.3f %i \n", count, ic1, time2, time, timeoffset, time2 - time - timeoffset, mwoffset);
	      }
	  	
	    tree2->GetEntry(currentEntryTree2);
			
	    if ( time2 - time - timeoffset2 > timewindow2)
	      {
		tree2->GetEntry(currentEntryTree2 + mwoffset); //why currentEtryTree2-4?
		mw1_posx_avg += mw1_posx;
		mw1_posy_avg += mw1_posy;
		count2++;
		//  if (i % 2000 == 0) printf("%i  %2.3f  %2.3f  %2.3f  %2.3f\n", count2, time, time2, ic1, mw1_posx);
		//if (i % 2001 == 0) printf("%i  %2.3f  %2.3f  %2.3f  %2.3f\n", count2, time, time2, ic1, mw1_posx);
		// if (i % 2002 == 0) printf("%i  %2.3f  %2.3f  %2.3f  %2.3f\n", count2, time, time2, ic1, mw1_posx);
		// printf("%i  %2.3f  %2.3f  %2.3f  %2.3f\n", count2, time, time2, ic1, mw1_posx);
	      }
	    // currentEntryTree2++;
	    tree2->GetEntry(currentEntryTree2);
	    currentEntryTree2++;
			
	  
		
		
	    if (count2>0){
	      mw1_posx_avg /= double(count2); //the positions are weighted by the charge
	      mw1_posy_avg /= double(count2);
	    }
	    if(count>0){
	      ic1_avg /= double(count);
	      if (ic1_avg>1.) icCounter++;
	      //   if (frame % 2000 == 0) printf("%i %f.2 %i \n", count, ic1_avg, icCounter);  
	    }
	    //	    std::cout << ic1_avg << " " << icCounter << std::endl;
	   

	  }
	/////////end of ethercat matching
	if(frame >=2000&& ic1_avg<0.5) continue;
	

	///////////////////////////////////
	//// HIT DATA
	///////////////////////////////////


	//board_b 0
	board_b=0;
	signal_b0 = 0.;
	maxchannelamp_b0 = 0.;

	file.seekg(framestart*sizeof(BufferData)+board_b*sizeof(BufferData)+4*frame*sizeof(BufferData));
	file.read ((char*)dataptr ,sizeof(BufferData));
	if (dataptr->sync_frame.device_nr==board_b){
	  //	  std::cout << frame << " " << dataptr->sync_frame.device_nr <<  " " << dataptr->sync_frame.local_ctr << " " << dataptr->sync_frame.global_ctr << " " << dataptr->sync_frame.sma_state << std::endl;
	  if (frame%1000==0) std::cout << "Frame: " <<  frame << " (" <<double(frame)/double(fileframesize)*100.0 << "%)" << std::endl;
	  for (int i = 0;i<128;i++){
	    board_b0_ch[i] = dataptr->sensor_data[i];
	    if (frame<1000){ 
	      board_b0_ch_bkg[i] +=  board_b0_ch[i]/1000.; //find baseline from the average of first 1000 events
	    
	    }
	    else if (frame>=1000&& frame<2000){ 
	      board_b0_ch[i] -=  board_b0_ch_bkg[i]; //histogram the subtracted  baseline in the next 1000 events
	      th2_signal_vs_channel_bkg_b0->Fill(i,board_b0_ch[i]);
	      if (i==1) graph_bkg_b0->SetPoint(graph_bkg_b0->GetN(), frame, board_b0_ch[i]);
	    }
	    else if (frame>=2000) {
	      board_b0_ch[i]-=board_b0_ch_bkg[i]; // the rest background subtracted
	      board_b0_ch[i]*=calibration_b0[i]; //calibration factor
	      th2_signal_vs_channel_b0->Fill(i,board_b0_ch[i]);
	      //  signal_b0 +=board_b0_ch[i] ;
	      if (board_b0_ch[i]> maxchannelamp_b0) {
		maxchannel_b0 = i;
		maxchannelamp_b0 = board_b0_ch[i];
	      }
	      //calculate a rolling average of the signal
	      arrayavg_b0[i] = 0;
	      for (int j = 1; j<length;j++){
		array_b0[j-1][i] = array_b0[j][i];
		arrayavg_b0[i] += array_b0[j-1][i]/double(length);

	      }
	      if( board_b0_ch[i]>-1000){ 
		array_b0[length-1][i] =  board_b0_ch[i];
		arrayavg_b0[i] +=  array_b0[length-1][i]/double(length);
	      }
	      ////////////////////////////////////////////
	    }
	  }
	  //	th1_signal_b0->Fill(signal_b0);
	}
	else {
	  std::cout << "Error." << std::endl;
	}
    
	//board_b 1
	board_b=1;
	signal_b1=0.;
	maxchannelamp_b1 = 0.;

	file.seekg(framestart*sizeof(BufferData)+board_b*sizeof(BufferData)+4*frame*sizeof(BufferData));
	file.read ((char*)dataptr ,sizeof(BufferData));
	if (dataptr->sync_frame.device_nr==board_b){
	  //std::cout << frame << " " << dataptr->sync_frame.device_nr << std::endl;
	  for (int i = 0;i<128;i++){
	    board_b1_ch[i] = dataptr->sensor_data[i];
	    if (frame<1000){ 
	      board_b1_ch_bkg[i] +=  board_b1_ch[i]/1000.; //find baseline from the average of first 1000 events
	    }
	    else if (frame>=1000&& frame<2000){ 
	      board_b1_ch[i] -=  board_b1_ch_bkg[i]; //histogram the subtracted  baseline in the next 1000 events
	      th2_signal_vs_channel_bkg_b1->Fill(i,board_b1_ch[i]);
	      if (i==1) graph_bkg_b1->SetPoint(graph_bkg_b1->GetN(), frame, board_b1_ch[i]);

	    }
	    else if (frame>=2000) {
	      board_b1_ch[i]-=board_b1_ch_bkg[i]; // the rest are background subtracted
	      board_b1_ch[i]*=calibration_b1[i]; //calibration factor
	      
	      th2_signal_vs_channel_b1->Fill(i,board_b1_ch[i]);
	      // signal_b1 +=board_b1_ch[i] ;
	      if (board_b1_ch[i]> maxchannelamp_b1) {
		maxchannel_b1 = i;
		maxchannelamp_b1 = board_b1_ch[i];
	      }
	      //calculate a rolling average of the signal
	      arrayavg_b1[i] = 0;
	      for (int j = 1; j<length;j++){
		array_b1[j-1][i] = array_b1[j][i];
		arrayavg_b1[i] += array_b1[j-1][i]/double(length);

	      }
	      if( board_b1_ch[i]>-1000){ 
		array_b1[length-1][i] =  board_b1_ch[i];
		arrayavg_b1[i] +=  array_b1[length-1][i]/double(length);
	      }
	      ////////////////////////////////////////////
	    }
	  }
	  //	th1_signal_b1->Fill(signal_b1);
	}
	else {
	  std::cout << "Error." << std::endl;
	}
	//board_b 2
	board_b=2;
	signal_b2=0.;
	maxchannelamp_b2 = 0.;
	file.seekg(framestart*sizeof(BufferData)+board_b*sizeof(BufferData)+4*frame*sizeof(BufferData));
	file.read ((char*)dataptr ,sizeof(BufferData));
	if (dataptr->sync_frame.device_nr==board_b){
	 //std::cout << frame << " " << dataptr->sync_frame.device_nr << std::endl;
	  for (int i = 0;i<128;i++){
	    board_b2_ch[i] = dataptr->sensor_data[i];
	    if (frame<1000){ 
	      board_b2_ch_bkg[i] +=  board_b2_ch[i]/1000.; //find baseline from the average of first 1000 events
	    }
	    else if (frame>=1000&& frame<2000){ 
	      board_b2_ch[i] -=  board_b2_ch_bkg[i]; //histogram the subtracted  baseline in the next 1000 events
	      th2_signal_vs_channel_bkg_b2->Fill(i,board_b2_ch[i]);
	      if (i==1) graph_bkg_b2->SetPoint(graph_bkg_b2->GetN(), frame, board_b2_ch[i]);

	    }
	    else if (frame>=2000) {
	      board_b2_ch[i]-=board_b2_ch_bkg[i]; // the rest background subtracted
	      board_b2_ch[i]*=calibration_b2[i]; //calibration factor
	      
	      th2_signal_vs_channel_b2->Fill(i,board_b2_ch[i]);
	      //  signal_b2 +=board_b2_ch[i] ;
	      if (board_b2_ch[i]> maxchannelamp_b2) {
		maxchannel_b2 = i;
		maxchannelamp_b2 = board_b2_ch[i];
	      }
	      //calculate a rolling average of the signal
	      arrayavg_b2[i] = 0;
	      for (int j = 1; j<length;j++){
		array_b2[j-1][i] = array_b2[j][i];
		arrayavg_b2[i] += array_b2[j-1][i]/double(length);

	      }
	      if( board_b2_ch[i]>-1000){ 
		array_b2[length-1][i] =  board_b2_ch[i];
		arrayavg_b2[i] +=  array_b2[length-1][i]/double(length);
	      }
	      ////////////////////////////////////////////
	    }
	  }
	  //	th1_signal_b2->Fill(signal_b2);
	}
	else {
	  std::cout << "Error." << std::endl;
	}

	//board_b 3
	board_b=3;
	signal_b3=0.;
	maxchannelamp_b3 = 0.;
	file.seekg(framestart*sizeof(BufferData)+board_b*sizeof(BufferData)+4*frame*sizeof(BufferData));
	file.read ((char*)dataptr ,sizeof(BufferData));
	if (dataptr->sync_frame.device_nr==board_b){
	  //std::cout << frame << " " << dataptr->sync_frame.device_nr << std::endl;
	  for (int i = 0;i<128;i++){
	    board_b3_ch[i] = dataptr->sensor_data[i];
	    if (frame<1000){ 
	      board_b3_ch_bkg[i] +=  board_b3_ch[i]/1000.; //find baseline from the average of first 1000 events
	    }
	    else if (frame>=1000&& frame<2000){ 
	      board_b3_ch[i] -=  board_b3_ch_bkg[i]; //histogram the subtracted  baseline in the next 1000 events
	      th2_signal_vs_channel_bkg_b3->Fill(i,board_b3_ch[i]);
	      if (i==1) graph_bkg_b3->SetPoint(graph_bkg_b3->GetN(), frame, board_b3_ch[i]);

	    }
	    else if (frame>=2000) {
	      board_b3_ch[i]-=board_b3_ch_bkg[i]; // the rest of the events are background subtracted
	      board_b3_ch[i]*=calibration_b3[i]; //with a calibration factor
	      
	      th2_signal_vs_channel_b3->Fill(i,board_b3_ch[i]);
	      //  signal_b3 +=board_b3_ch[i] ;
	      if (board_b3_ch[i]> maxchannelamp_b3) {
		maxchannel_b3 = i;
		maxchannelamp_b3 = board_b3_ch[i];
	      }
	      //calculate a rolling average of the signal
	      arrayavg_b3[i] = 0;
	      for (int j = 1; j<length;j++){
		array_b3[j-1][i] = array_b3[j][i];
		arrayavg_b3[i] += array_b3[j-1][i]/double(length);

	      }
	      if( board_b3_ch[i]>-1000){ 
		array_b3[length-1][i] =  board_b3_ch[i];
		arrayavg_b3[i] +=  array_b3[length-1][i]/double(length);
	      }
	      ////////////////////////////////////////////
	    }
	  }
	  //	th1_signal_b3->Fill(signal_b3);
	}
	else {
	  std::cout << "Error." << std::endl;
	}




	////////////////////////////////////////////////////////////
	////////////////////////////////////////////////////////////
	//start the signal analysis
	////////////////////////////////////////////////////////////
	////////////////////////////////////////////////////////////
	
	if (frame>=2000){




	  ////////////////////////////////////////////////////////////
	  //boxcar smoothing filter
	  ////////////////////////////////////////////////////////////
	  if (smooth_on){
	    maxchannelamp_b0 = 0.0;
	    maxchannelamp_b1 = 0.0;
	    maxchannelamp_b2 = 0.0;
	    maxchannelamp_b3 = 0.0;
	  }
	  boxcar.clear();
	  boxcar.push_back(board_b0_ch[0]);
	  boxcar.push_back(board_b0_ch[1]);
	  //     std::cout << frame << std::endl;
	  for (int i = 0;i<128;i++){
	    channellist[i] = i;
	    if (i<126) {
	      boxcar.push_back(board_b0_ch[i+2]);
	    }
	    if ( (i<126&&boxcar.size()>5) || (i>=126&&boxcar.size()>3) ) boxcar.erase(boxcar.cbegin());
	    board_b0_smooth_ch[i] = gsl_stats_mean(boxcar.data(), 1, 5 );
	    if (smooth_on && board_b0_smooth_ch[i]> maxchannelamp_b0) {
		maxchannel_b0 = i;
		maxchannelamp_b0 = board_b0_smooth_ch[i];
	      }
	  }
	  if (eventID>2001&&maxchannelamp_b0>100&&!graphsaved_b0) {
	    gr_b0 = new TGraph(128, channellist, board_b0_ch);
	    gr_b0->SetTitle("RawData_b0");
	    gr_b0->SetName("RawData_b0");
	    gr_b0->Write();
	    gr_sm_b0 = new TGraph(128, channellist, board_b0_smooth_ch);
	    gr_sm_b0->SetTitle("SmoothedData_b0");
	    gr_sm_b0->SetName("SmoothedData_b0");
	    gr_sm_b0->Write();
	    graphsaved_b0=true;
	  }

	  boxcar.clear();
	  boxcar.push_back(board_b1_ch[0]);
	  boxcar.push_back(board_b1_ch[1]);
	  //     std::cout << frame << std::endl;
	  for (int i = 0;i<128;i++){
	    channellist[i] = i;
	    if (i<126) {
	      boxcar.push_back(board_b1_ch[i+2]);
	    }
	    if ( (i<126&&boxcar.size()>5) || (i>=126&&boxcar.size()>3) ) boxcar.erase(boxcar.cbegin());
	    board_b1_smooth_ch[i] = gsl_stats_mean(boxcar.data(), 1, 5 );
	    if (smooth_on && board_b1_smooth_ch[i]> maxchannelamp_b1) {
		maxchannel_b1 = i;
		maxchannelamp_b1 = board_b1_smooth_ch[i];
	      }
	  }
	  if (eventID>2001&&maxchannelamp_b1>100&&!graphsaved_b1) {
	    gr_b1 = new TGraph(128, channellist, board_b1_ch);
	    gr_b1->SetTitle("RawData_b1");
	    gr_b1->SetName("RawData_b1");
	    gr_b1->Write();
	    gr_sm_b1 = new TGraph(128, channellist, board_b1_smooth_ch);
	    gr_sm_b1->SetTitle("SmoothedData_b1");
	    gr_sm_b1->SetName("SmoothedData_b1");
	    gr_sm_b1->Write();
	    graphsaved_b1=true;
	  }

	  boxcar.clear();
	  boxcar.push_back(board_b2_ch[0]);
	  boxcar.push_back(board_b2_ch[1]);
	  //     std::cout << frame << std::endl;
	  for (int i = 0;i<128;i++){
	    channellist[i] = i;
	    if (i<126) {
	      boxcar.push_back(board_b2_ch[i+2]);
	    }
	    if ( (i<126&&boxcar.size()>5) || (i>=126&&boxcar.size()>3) ) boxcar.erase(boxcar.cbegin());
	    board_b2_smooth_ch[i] = gsl_stats_mean(boxcar.data(), 1, 5 );
	    if (smooth_on && board_b2_smooth_ch[i]> maxchannelamp_b2) {
		maxchannel_b2 = i;
		maxchannelamp_b2 = board_b2_smooth_ch[i];
	      }
	  }
	  if (eventID>2001&&maxchannelamp_b2>100&&!graphsaved_b2) {
	    gr_b2 = new TGraph(128, channellist, board_b2_ch);
	    gr_b2->SetTitle("RawData_b2");
	    gr_b2->SetName("RawData_b2");
	    gr_b2->Write();
	    gr_sm_b2 = new TGraph(128, channellist, board_b2_smooth_ch);
	    gr_sm_b2->SetTitle("SmoothedData_b2");
	    gr_sm_b2->SetName("SmoothedData_b2");
	    gr_sm_b2->Write();
	    graphsaved_b2=true;
	  }


	  boxcar.clear();
	  boxcar.push_back(board_b3_ch[0]);
	  boxcar.push_back(board_b3_ch[1]);
	  //     std::cout << frame << std::endl;
	  for (int i = 0;i<128;i++){
	    channellist[i] = i;
	    if (i<126) {
	      boxcar.push_back(board_b3_ch[i+2]);
	    }
	    if ( (i<126&&boxcar.size()>5) || (i>=126&&boxcar.size()>3) ) boxcar.erase(boxcar.cbegin());
	    board_b3_smooth_ch[i] = gsl_stats_mean(boxcar.data(), 1, 5 );
	    if (smooth_on && board_b3_smooth_ch[i]> maxchannelamp_b3) {
		maxchannel_b3 = i;
		maxchannelamp_b3 = board_b3_smooth_ch[i];
	      }
	  }
	  if (eventID>2001&&maxchannelamp_b3>100&&!graphsaved_b3) {
	    gr_b3 = new TGraph(128, channellist, board_b3_ch);
	    gr_b3->SetTitle("RawData_b3");
	    gr_b3->SetName("RawData_b3");
	    gr_b3->Write();
	    gr_sm_b3 = new TGraph(128, channellist, board_b3_smooth_ch);
	    gr_sm_b3->SetTitle("SmoothedData_b3");
	    gr_sm_b3->SetName("SmoothedData_b3");
	    gr_sm_b3->Write();
	    graphsaved_b3=true;
	  }



	////////////////////////////////////////////////////////////
	//find the approx FWHM
	////////////////////////////////////////////////////////////
	
	  nfwhm_b0 =0;
	  nfwhm_b1 =0;
	  nfwhm_b2 =0;
	  nfwhm_b3 =0;
	  if (smooth_on){
	    for (int i = 0;i<128;i++){
	      if (board_b0_smooth_ch[i] > maxchannelamp_b0/2.) nfwhm_b0++;
	      if (board_b1_smooth_ch[i] > maxchannelamp_b1/2.) nfwhm_b1++;
	      if (board_b2_smooth_ch[i] > maxchannelamp_b2/2.) nfwhm_b2++;
	      if (board_b3_smooth_ch[i] > maxchannelamp_b3/2.) nfwhm_b3++;
	      signal_gsl_b0[i] = 0.;
	      signal_gsl_b1[i] = 0.;
	      signal_gsl_b2[i] = 0.;
	      signal_gsl_b3[i] = 0.;

	    }

	  }
	  else {
	    for (int i = 0;i<128;i++){
	      if (board_b0_ch[i] > maxchannelamp_b0/2.) nfwhm_b0++;
	      if (board_b1_ch[i] > maxchannelamp_b1/2.) nfwhm_b1++;
	      if (board_b2_ch[i] > maxchannelamp_b2/2.) nfwhm_b2++;
	      if (board_b3_ch[i] > maxchannelamp_b3/2.) nfwhm_b3++;
	      signal_gsl_b0[i] = 0.;
	      signal_gsl_b1[i] = 0.;
	      signal_gsl_b2[i] = 0.;
	      signal_gsl_b3[i] = 0.;

	    }
	  }
	  

	  ////////////////////////////////////////////////////////////
	  //integrate the sidebands first to check for baseline shift
	  ////////////////////////////////////////////////////////////
	 
	  beamSidebandNoise_b0 = 0.;
	  sidenumtocalc_b0 = 0;
	  for (int i =0; i<maxchannel_b0-nfwhm_b0;  i++){
	    if (i>=0 && i<=127){
	      beamSidebandNoise_b0 +=  board_b0_ch[i]; //integrate the noise outside the peak
	      sidenumtocalc_b0++;
	    }
	  }
	    for (int i = maxchannel_b0+nfwhm_b0; i <  128;  i++){
	      if (i>=0 && i<=127){
		beamSidebandNoise_b0 +=  board_b0_ch[i]; //integrate the noise outside the peak
		sidenumtocalc_b0++;
	      }
	    }
	    if (sidenumtocalc_b0>0) beamSidebandNoise_b0 = beamSidebandNoise_b0 /double(sidenumtocalc_b0);  // channel baseline shift

	    beamSidebandNoise_b1 = 0.;
	    sidenumtocalc_b1 = 0;
	    for (int i =0; i<maxchannel_b1-nfwhm_b1;  i++){
	      if (i>=0 && i<=127){
		beamSidebandNoise_b1 +=  board_b1_ch[i]; //integrate the noise outside the peak
		sidenumtocalc_b1++;
	      }
	    }
	      for (int i = maxchannel_b1+nfwhm_b1; i <  128;  i++){
		if (i>=0 && i<=127){
		  beamSidebandNoise_b1 +=  board_b1_ch[i]; //integrate the noise outside the peak
		  sidenumtocalc_b1++;
		}
	      }
	      if (sidenumtocalc_b1>0) beamSidebandNoise_b1 = beamSidebandNoise_b1 /double(sidenumtocalc_b1);  // channel baseline shift

	      beamSidebandNoise_b2 = 0.;
	      sidenumtocalc_b2 = 0;
	      for (int i =0; i<maxchannel_b2-nfwhm_b2;  i++){
		if (i>=0 && i<=127){
		  beamSidebandNoise_b2 +=  board_b2_ch[i]; //integrate the noise outside the peak
		  sidenumtocalc_b2++;
		}
	      }
		for (int i = maxchannel_b2+nfwhm_b2; i <  128;  i++){
		  if (i>=0 && i<=127){
		    beamSidebandNoise_b2 +=  board_b2_ch[i]; //integrate the noise outside the peak
		    sidenumtocalc_b2++;
		  }
		}
		if (sidenumtocalc_b2>0) beamSidebandNoise_b2 = beamSidebandNoise_b2 /double(sidenumtocalc_b2);  // channel baseline shift

		beamSidebandNoise_b3 = 0.;
		sidenumtocalc_b3 = 0;
		for (int i =0; i<maxchannel_b3-nfwhm_b3;  i++){
		  if (i>=0 && i<=127){
		    beamSidebandNoise_b3 +=  board_b3_ch[i]; //integrate the noise outside the peak
		    sidenumtocalc_b3++;
		  }
		}
		  for (int i = maxchannel_b3+nfwhm_b3; i <  128;  i++){
		    if (i>=0 && i<=127){
		      beamSidebandNoise_b3 +=  board_b3_ch[i]; //integrate the noise outside the peak
		      sidenumtocalc_b3++;
		    }
		  }
		  if (sidenumtocalc_b3>0) beamSidebandNoise_b3 = beamSidebandNoise_b3 /double(sidenumtocalc_b3);  // channel baseline shift


	  


	  ////////////////////////////////////////////////////////////
	  //integrate under the approximate peak
	  //build the channel list for statistical analysis
	  ////////////////////////////////////////////////////////////

	  numtocalc_b0 = 0;
	  numtocalc_b1 = 0;
	  numtocalc_b2 = 0;
	  numtocalc_b3 = 0;
	  beamSidebandNoise_b0=0.;
	  beamSidebandNoise_b1=0.;
	  beamSidebandNoise_b2=0.;
	  beamSidebandNoise_b3=0.;

	  signal_b0_left = 0.;
	  signal_b0_right = 0.;
	  for (int i = maxchannel_b0-nfwhm_b0 ; i <= maxchannel_b0 + nfwhm_b0; i++){
	    if (i>=0 && i<=127 && board_b0_ch[i]>50){
	      signal_b0 +=board_b0_ch[i]-beamSidebandNoise_b0  ;
	      if (i<64) {signal_b0_left +=board_b0_ch[i]-beamSidebandNoise_b0  ;}
	      else {signal_b0_right +=board_b0_ch[i]-beamSidebandNoise_b0  ;}
	      signal_gsl_b0[numtocalc_b0]=board_b0_ch[i]-beamSidebandNoise_b0 ;
	      // channellist_gsl_b0[numtocalc_b0] = i;
	      channellist_gsl_b0[numtocalc_b0] = pos[i];
	      numtocalc_b0++;
	    }
	  }
	  th1_signal_b0->Fill(signal_b0);

	  signal_b1_left = 0.;
	  signal_b1_right = 0.;
	  for (int i = maxchannel_b1-nfwhm_b1 ; i <= maxchannel_b1 + nfwhm_b1; i++){
	    if (i>=0 && i<=127&&board_b1_ch[i]>50){
	      signal_b1 +=board_b1_ch[i]-beamSidebandNoise_b1  ;
	      if (i<64) {signal_b1_left +=board_b1_ch[i]-beamSidebandNoise_b1  ;}
	      else {signal_b1_right +=board_b1_ch[i]-beamSidebandNoise_b1  ;}
	      signal_gsl_b1[numtocalc_b1]=board_b1_ch[i]-beamSidebandNoise_b1 ;
	      // channellist_gsl_b1[numtocalc_b1] = i;
	      channellist_gsl_b1[numtocalc_b1] = pos[i];
	      numtocalc_b1++;
	    }
	  }
	  th1_signal_b0->Fill(signal_b1);

	  signal_b2_left = 0.;
	  signal_b2_right = 0.;
	  for (int i = maxchannel_b2-nfwhm_b2 ; i <= maxchannel_b2 + nfwhm_b2; i++){
	    if (i>=0 && i<=127&&board_b2_ch[i]>50){
	      signal_b2 +=board_b2_ch[i]-beamSidebandNoise_b2  ;
	      if (i<64) {signal_b2_left +=board_b2_ch[i]-beamSidebandNoise_b2  ;}
	      else {signal_b2_right +=board_b2_ch[i]-beamSidebandNoise_b2  ;}
	      signal_gsl_b2[numtocalc_b2]=board_b2_ch[i]-beamSidebandNoise_b2 ;
	      // channellist_gsl_b2[numtocalc_b2] = i;
	      channellist_gsl_b2[numtocalc_b2] = pos[i];
	      numtocalc_b2++;
	    }
	  }
	  th1_signal_b0->Fill(signal_b2);

	  signal_b3_left = 0.;
	  signal_b3_right = 0.;
	  for (int i = maxchannel_b3-nfwhm_b3 ; i <= maxchannel_b3 + nfwhm_b3; i++){
	    if (i>=0 && i<=127&&board_b3_ch[i]>50){
	      signal_b3 +=board_b3_ch[i]-beamSidebandNoise_b3 ;
	      if (i<64) {signal_b3_left +=board_b3_ch[i]-beamSidebandNoise_b3  ;}
	      else {signal_b3_right +=board_b3_ch[i]-beamSidebandNoise_b3  ;}
	      signal_gsl_b3[numtocalc_b3]=board_b3_ch[i]-beamSidebandNoise_b3 ;
	      //channellist_gsl_b3[numtocalc_b3] = i;
	      channellist_gsl_b3[numtocalc_b3] = pos[i];
	      numtocalc_b3++;
	    }
	  }
	  th1_signal_b0->Fill(signal_b3);


	  for (int i=0;i<128;i++){
	    if(signal_b0>700) {
	      th2_signal_vs_channel_sub_b0->Fill(i,board_b0_ch[i]-beamSidebandNoise_b0);
	      sumvector_b0[i] += board_b0_ch[i];//-beamSidebandNoise_b0;
     
	    }
	    if(signal_b1>700) {
	      th2_signal_vs_channel_sub_b1->Fill(i,board_b1_ch[i]-beamSidebandNoise_b1);
	      sumvector_b1[i] += board_b1_ch[i];//-beamSidebandNoise_b1;
	    }
	    if(signal_b2>700){
	      th2_signal_vs_channel_sub_b2->Fill(i,board_b2_ch[i]-beamSidebandNoise_b2);
	      sumvector_b2[i] += board_b2_ch[i];//-beamSidebandNoise_b2;
	    }
	    if(signal_b3>700) {
	      th2_signal_vs_channel_sub_b3->Fill(i,board_b3_ch[i]-beamSidebandNoise_b3);
	      sumvector_b3[i] += board_b3_ch[i];//-beamSidebandNoise_b3;
	    }
	  }
	  if(signal_b0>100) beamontime[0]+=1.0;
	  if(signal_b1>100) beamontime[1]+=1.0;
	  if(signal_b2>100) beamontime[2]+=1.0;
	  if(signal_b3>100) beamontime[3]+=1.0;



	  ///add gsl stuff here.
	  /*	std::cout << maxchannel_b0 << " " << maxchannel_b1 << " "<< maxchannel_b2 << " "<< maxchannel_b3 << " " << std::endl;
		std::cout << maxchannelamp_b0 << " " << maxchannelamp_b1 << " "<< maxchannelamp_b2 << " "<< maxchannelamp_b3 << " " << std::endl;
		std::cout << nfwhm_b0 << " " << nfwhm_b1 << " " << nfwhm_b2 << " " << nfwhm_b3 << " " << std::endl;
		std::cout << std::endl;*/

	  beamPosX_b0 = gsl_stats_wmean(signal_gsl_b0,1,channellist_gsl_b0,1,numtocalc_b0); //calculate the weighted mean
	  beamFocusX_b0 = gsl_stats_wsd_with_fixed_mean(signal_gsl_b0,1,channellist_gsl_b0,1,numtocalc_b0,beamPosX_b0); //Standard Deviation
	  beamSkewX_b0 = gsl_stats_wskew_m_sd(signal_gsl_b0,1,channellist_gsl_b0,1,numtocalc_b0,beamPosX_b0,beamFocusX_b0); //skewness (symmetry)
	  beamKurtX_b0 = gsl_stats_wkurtosis_m_sd(signal_gsl_b0,1,channellist_gsl_b0,1,numtocalc_b0,beamPosX_b0,beamFocusX_b0); //excess kurtosis (well behaved tails)
	  beamNumX_b0 = numtocalc_b0;
	  beamFocusX_b0 *=2.3548;//SD-->FWHM

	  beamPosX_b1 = gsl_stats_wmean(signal_gsl_b1,1,channellist_gsl_b1,1,numtocalc_b1); //calculate the weighted mean
	  beamFocusX_b1 = gsl_stats_wsd_with_fixed_mean(signal_gsl_b1,1,channellist_gsl_b1,1,numtocalc_b1,beamPosX_b1); //Standard Deviation
	  beamSkewX_b1 = gsl_stats_wskew_m_sd(signal_gsl_b1,1,channellist_gsl_b1,1,numtocalc_b1,beamPosX_b1,beamFocusX_b1); //skewness (symmetry)
	  beamKurtX_b1 = gsl_stats_wkurtosis_m_sd(signal_gsl_b1,1,channellist_gsl_b1,1,numtocalc_b1,beamPosX_b1,beamFocusX_b1); //excess kurtosis (well behaved tails)
	  beamNumX_b1 = numtocalc_b1;
	  beamFocusX_b1 *=2.3548;//SD-->FWHM

	  beamPosX_b2 = gsl_stats_wmean(signal_gsl_b2,1,channellist_gsl_b2,1,numtocalc_b2); //calculate the weighted mean
	  beamFocusX_b2 = gsl_stats_wsd_with_fixed_mean(signal_gsl_b2,1,channellist_gsl_b2,1,numtocalc_b2,beamPosX_b2); //Standard Deviation
	  beamSkewX_b2 = gsl_stats_wskew_m_sd(signal_gsl_b2,1,channellist_gsl_b2,1,numtocalc_b2,beamPosX_b2,beamFocusX_b2); //skewness (symmetry)
	  beamKurtX_b2 = gsl_stats_wkurtosis_m_sd(signal_gsl_b2,1,channellist_gsl_b2,1,numtocalc_b2,beamPosX_b2,beamFocusX_b2); //excess kurtosis (well behaved tails)
	  beamNumX_b2 = numtocalc_b2;
	  beamFocusX_b2 *=2.3548;//SD-->FWHM

	  beamPosX_b3 = gsl_stats_wmean(signal_gsl_b3,1,channellist_gsl_b3,1,numtocalc_b3); //calculate the weighted mean
	  beamFocusX_b3 = gsl_stats_wsd_with_fixed_mean(signal_gsl_b3,1,channellist_gsl_b3,1,numtocalc_b3,beamPosX_b3); //Standard Deviation
	  beamSkewX_b3 = gsl_stats_wskew_m_sd(signal_gsl_b3,1,channellist_gsl_b3,1,numtocalc_b3,beamPosX_b3,beamFocusX_b3); //skewness (symmetry)
	  beamKurtX_b3 = gsl_stats_wkurtosis_m_sd(signal_gsl_b3,1,channellist_gsl_b3,1,numtocalc_b3,beamPosX_b3,beamFocusX_b3); //excess kurtosis (well behaved tails)
	  beamNumX_b3 = numtocalc_b3;
	  beamFocusX_b3 *=2.3548;//SD-->FWHM


	  /////////////////////////////////////////////////
	  /////Add fitting algorithm here
	  /////////////////////////////////////////////////
	  /*

	    use board_b0_ch[i] to fill a TGraph. 
	    fit with a gaussian
	    subtract the difference
	    fill a THF2 with the difference, delta_amp vs channel
	   */
	
	  //	  std::cout << ic1_avg << std::endl;


	  //fit with a gaussian function;
	  if (doFit&&ic1_avg>0.1){
	    if ((lastfit_b0!=0||lastfit_b1!=0)&&numtocalc_b0>1&&numtocalc_b0<50){
	      gausfunc_b0->SetParameters(signal_b0/(sqrt(2)*beamFocusX_b0/2.3548),beamPosX_b0,4.,0.);
	      gausfunc_b1->SetParameters(signal_b1/(sqrt(2)*beamFocusX_b1/2.3548),beamPosX_b1,4.,0.);
    
	    }
	    else if (signal_b0<=0||signal_b1<=0) {
	      gausfunc_b0->SetParameters(100.,35.,10.,0.);
	      gausfunc_b1->SetParameters(100.,35.,10.,0.);
	    }

	
	    gausfunc_b0->SetParLimits(0,0.,5000.);
	    gausfunc_b0->SetParLimits(1,0.,128.);
	    gausfunc_b0->SetParLimits(2,0.5,40.);
	    gausfunc_b0->SetParLimits(3,-150.,150.);

	    gausfunc_b1->SetParLimits(0,0.,5000.);
	    gausfunc_b1->SetParLimits(1,0,128);
	    gausfunc_b1->SetParLimits(2,0.5,40.);
	    gausfunc_b1->SetParLimits(3,-150.,150.);
	    //signal_gsl_b0,1,channellist_gsl_b0,1,numtocalc_b0)
	    gausgraph_b0 = new TGraphErrors(128,channellist_gsl_b0,signal_gsl_b0,errorx,errory);
	    gausgraph_b1 = new TGraphErrors(128,channellist_gsl_b1,signal_gsl_b1,errorx,errory);
	    

	    lastfit_b0 = gausgraph_b0->Fit(gausfunc_b0,"QRN","",1,128); // single gaussian fit
	    lastfit_b1 = gausgraph_b1->Fit(gausfunc_b1,"QRN","",1,128); // single gaussian fit
	    beamPosX_fit_b0 = gausfunc_b0->GetParameter(1);
	    beamFocusX_fit_b0 =2.3548* gausfunc_b0->GetParameter(2);
	    //   beamChi2_fit_b0 = gausfunc_b0->GetChisquare()/gausfunc_b0->GetNDF();
	    //   beamPeakX_fit_b0 = gausfunc_b0->GetParameter(0);
	    
	    beamPosX_fit_b1 = gausfunc_b1->GetParameter(1);
	    beamFocusX_fit_b1 =2.3548* gausfunc_b1->GetParameter(2);
	    // beamChi2_fit_b1 = gausfunc_b1->GetChisquare()/gausfunc_b1->GetNDF();
	    // beamPeakX_fit_b1 = gausfunc_b1->GetParameter(0);
	    //	    cout << lastfit_b0 << " " << lastfit_b1 << " " << ic1_avg << " " << maxchannelamp_b0 << " " << nfwhm_b0<< " " << nfwhm_b1 <<" " << beamPosX_b0 << " " <<beamPosX_b1 << " " << beamPosX_fit_b0 << " " << beamPosX_fit_b1 << " " << signal_b0 << " " << signal_b1 << " " << energy << " " << intensity << endl;
	    
	    	   
	    /*
	      for (int ch = 4; ch < 64; ch++){
	      th2d_fitdiff_channel->Fill(ch, double(channelamp[ch]-gausfunc->Eval(ch)));
	      }
	      for (int ch = 64; ch < 128; ch++){
	      th2d_fit2diff_channel->Fill(ch, double(channelamp[ch]-gausfunc2->Eval(ch)));
	      }
	    */
	    
	  }
      

	  
	  rootTree->Fill();

	 
	 
	}
	if (frameend>0 && frame+framestart>=frameend) break;
      }//end of loop over frames
      /*    th2_signal_vs_channel_b0->Write();
      th2_signal_vs_channel_b1->Write();
      th2_signal_vs_channel_b2->Write();
      th2_signal_vs_channel_b3->Write();

      th2_signal_vs_channel_bkg_b0->Write();
      th2_signal_vs_channel_bkg_b1->Write();
      th2_signal_vs_channel_bkg_b2->Write();
      th2_signal_vs_channel_bkg_b3->Write();

      th1_signal_b0->Write();
      th1_signal_b1->Write();
      th1_signal_b2->Write();
      th1_signal_b3->Write();*/
      graph_bkg_b0->SetName("graph_bkg_b0"); // graph_bkg_b0->Write();
      graph_bkg_b1->SetName("graph_bkg_b1"); // graph_bkg_b1->Write();
      graph_bkg_b2->SetName("graph_bkg_b2"); // graph_bkg_b2->Write();
      graph_bkg_b3->SetName("graph_bkg_b3"); // graph_bkg_b3->Write();
      sumvector_b0_ptr->Write("sumvector_b0");
      sumvector_b1_ptr->Write("sumvector_b1");
      sumvector_b2_ptr->Write("sumvector_b2");
      sumvector_b3_ptr->Write("sumvector_b3");
      beamontime_ptr->Write("beamontime");
      // bic1->Fill();
      // bmw1_posx->Fill();
      //bmw1_posy->Fill();
      rootFile->Write();
      rootFile->Close();
    }
    std::cout << eventID << " frames analysed." << std::endl;
    return 0;
  }
  else {
    std::cerr << "Error 1" << std::endl;
    return 1;
  }
  
}

int main(int argc, char **argv){

  analyse(argc, argv);
  
  }