#include "hitreader.h"
#include <string>
#include <stdio.h>
#include <iostream>
#include <numeric>
#include <vector>
#include <utility>
#include <TH2.h>
#include <TF1.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"
#include "TMatrixD.h"
#include "TApplication.h"
#include "TRint.h"
#include "TSystem.h"
#include "TMultiGraph.h"
#include "TObject.h"

//function declarations
int analyse(int argc, char **argv);
int opendatafiles(int argc, char ** argv);
void histograms(int argc, char ** argv);
int closedatafiles();
int set_background_v2(int start_frame, int frames);


//////global variables

//opendatafiles
long int  fileframesize;
char * filename;
ifstream file;
TFile *rootFile; //output root file
TTree *rootTree; 
char * rootfilename;
char * offsetfilename;
ifstream offsetfile;

struct bpm_frame_v2 {
  double channel_amp[320];
  double channel_position[320];
  double avg_position;
  double avg_width;
  double integratedsignalamp;
  int maxchannel;
  double maxchannel_amp;
};


bpm_frame_v2 board_b[4];
bpm_frame_v2 board_b_bkg[4];


bpm_frame_v2 readboard(Fullframe frame, int boardnumber);
int framestart = 0;

struct beamRecon  {

  double Position;
  double Focus;
  double Peak;
  double Rsqr;
  double Skew;
  double Kurtosis;
  double Sum;
  int n_channels;
};


beamRecon beamreconstruction(bpm_frame_v2 frametoanalyse, double threshold);
beamRecon BPMbeamrecon_0,BPMbeamrecon_1,BPMbeamrecon_2,BPMbeamrecon_3;

beamRecon BPMbeamrecon_Zeroed;

    

TGraph *  board_graph = new TGraph(320);
TF1 * beamrecon_function1 = new TF1("beamrecon_function1", "gaus(0)",0.,320.);


class Channel{
 public:
  double amplitude;;
  double position;
  int chnumber;
  int last_neighbour;
  int next_neighbour;
};


bool CompareChannels(Channel A, Channel B);

double avg ( vector<Channel> v );
double variance ( vector<Channel> v , double mean ); 

bool visualize = false;
int eventID;

////histograms
TH2D * TH2D_b0_signal_vs_channel;
TH2D * TH2D_b1_signal_vs_channel;
TH2D * TH2D_b2_signal_vs_channel;
TH2D * TH2D_b3_signal_vs_channel;