#include "TCanvas.h"
#include "TGraph.h"
#include "TStyle.h"
#include "TAxis.h"
#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <fstream>

using namespace std;

const int nchan = 32, npeaks = 15;
double pedestals[nchan] = {0};
double peaks[npeaks][nchan] = {{0}};

void analyse_tree(const char* fname) {
  ifstream in(fname);
  gStyle->SetOptStat(0);
  const int spacing = 14, len = 218, ped_length = 17;
  int ped_start = len - (ped_length + 3);
  int n = 0; // number of pulse trains
  bool first_zero = false, gap_reached = false;
  for (int i = 0; in.good(); ++i) {
    int sample[nchan];
    int nzeroes = 0;
    for (int j = 0; j < nchan; ++j) {
      in >> sample[j];
      if (sample[j] < 30)
	++nzeroes;
    }
    if (!gap_reached && nzeroes > 10) { // dip in most channels
      if (!first_zero) {
	for (int k = 0; k < 3; ++k) {
	  nzeroes = 0;
	  for (int j = 0; j < nchan; ++j) {
	    in >> sample[j];
	    if (sample[j] < 30)
	      ++nzeroes;
	  }
	  if (nzeroes <= 10)
	    break;
	}
	if (nzeroes > 10) { // quadruple row of zeroes
	  first_zero = true;
	  i = 0;
	}
      }
      if (i == 0)
	i = -1; // outside pulse train
      else if (first_zero && i < 30) {
	gap_reached = true;
	i = 0; // start of pulse train
      }
    }
    if (!gap_reached || !first_zero)
      continue;
    ++n;
    if (i == len) { // end of pulse train
      gap_reached = false;
      first_zero = false;
      i = -1;
      continue;
    }
    if (i >= ped_start && i < ped_start + ped_length) {
      for (int j = 0; j < nchan; ++j)
	pedestals[j] += sample[j];
      continue;
    }
    if (i % spacing == 0)
      for (int j = 0; j < nchan; ++j)
	peaks[i / spacing][j] += sample[j];
  }
  in.close();
  int nrem = n % len; // remainder of n
  n /= len; // integer part
  int ped_n = n * ped_length;
  if (nrem >= ped_start) {
    if (nrem - ped_start < ped_length)
      ped_n += nrem - ped_start + 1;
    else
      ped_n += ped_length;
  }
  for (int j = 0; j < nchan; ++j)
    pedestals[j] /= ped_n;
  for (int i = 0; i < npeaks; ++i) {
    for (int j = 0; j < nchan; ++j) {
      if (i * spacing < nrem)
	peaks[i][j] /= (n + 1);
      else
	peaks[i][j] /= n;
      peaks[i][j] -= pedestals[j]; // baseline subtraction
    }
  }
}

// displays peak amplitude as a function of peak number for the given channel
void display_peak_amp(int channel) {
  TGraph* graph = new TGraph(npeaks);
  for (int i = 0; i < npeaks; ++i)
    graph->SetPoint(i, i, peaks[i][channel]);
  TCanvas* c1 = new TCanvas(Form("amp_canvas_ch%d", channel), 
			    Form("Amplitude of pulse / cross-talk in channel %d"
				 , channel), 900, 600);
  graph->SetMarkerStyle(20);
  graph->SetMarkerSize(2);
  graph->SetMarkerColor(1);
  graph->Draw("AP");
  graph->GetXaxis()->SetTitle("peak #");
  graph->GetYaxis()->SetTitle("amplitude (ADC)");
  graph->SetTitle(Form("Amplitude as a function of peak # in channel %d", 
		       channel));
  c1->Modified();
}

/*Prints the cross-talk in all non-pulsed channels in fraction of the peak 
  amplitude to the specified file.
  parity == true -> even pulsing, parity == false -> odd pulsing*/
void print_cross_talk(char* fname, bool parity) {
  FILE* outfile = fopen(fname, "w");
  fprintf(outfile, "Peak #");
  for (int i = parity; i < nchan; i += 2)
    fprintf(outfile, " %6s %2d", "ch", i);
  fprintf(outfile, "\n");
  double peak_avg[npeaks] = {0};
  for (int i = 0; i < npeaks; ++i) {
    for (int j = !parity; j < nchan; j += 2)
      peak_avg[i] += peaks[i][j];
    peak_avg[i] /= (nchan / 2);
  }
  for (int i = 0; i < npeaks; ++i) {
    fprintf(outfile, "%6d", i);
    for (int j = parity; j < nchan; j += 2)
      fprintf(outfile, " %9.2e", peaks[i][j] / peak_avg[i]);
    fprintf(outfile, "\n");
  }
}