/* Code for checking cross-talk, to be compiled into an executable. Compilation:
   INC=`root-config --incdir`
   LIBS=`root-config --libs`
   g++ -std=c++11 -I${INC} ${LIBS} -o measure_xtalk.exe xtalk_compiled.cpp
- see usage with arguments in main() below
*/

#include "TF1.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}};

bool analyse_tree(const char* fname) {
  ifstream in(fname);
  if (in.fail()) {
    cerr << "Cannot open file named " << fname << "." << endl;
    return false;
  }
  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
    }
  }
  return true;
}

/*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, "a");
  fprintf(outfile, "Peak #");
  for (int i = parity; i < nchan; i += 2)
    fprintf(outfile, " %5s %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, " %8.1e", peaks[i][j] / peak_avg[i]);
    fprintf(outfile, "\n");
  }
}

int main(int argc, char** argv) {
  if (argc < 4) {
    printf("Not enough arguments: argc = %d\nSyntax should be:\n./measure_xtalk.exe chip_string config adcfilename\n", argc);
    printf("Example:\n./measure_xtalk.exe 100000 20mV_calib-odd /home/tpc/robot/20180725/123510_sampa100000/123510_sampa100000_trorc00_link00_20mV_calib-odd_adc.txt\n", argc);
    return -1;
  }

  char* infile = argv[3];
  // find dirname (last /)
  string str(infile);
  size_t found;
  found=str.find_last_of("/");
  char* outfile = Form("%s/xtalk_sampa%s_%s.log",  str.substr(0,found).c_str(), argv[1], argv[2]);

  if (!analyse_tree(infile))
    return 1;

  bool parity;
  if (strstr(argv[2], Form("odd"))) {
    cout << " odd selected " << endl; 
    parity = false;
  }
  else if (strstr(argv[2], Form("even"))) {
    cout << " even selected " << endl; 
    parity = true;
  }
  else {
    cerr << "Parity must be either 'even' or 'odd'" << endl;
    return 1;
  }
  print_cross_talk(outfile, parity);

  return 0;
}