#include "TH1F.h"
#include "TProfile.h"
#include "TMath.h"
#include "TF1.h"
#include "TLegend.h"
#include "TCanvas.h"
 //
#include "TROOT.h"
//#include "TStyle.h"
#include "TTree.h"
#include "TFile.h"
#include "TGraphErrors.h"
#include "TStyle.h"
#include <string>
#include <cmath>
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
//

using namespace std;

//const int debug = 0;
//const bool regeneratePlots = false;
const int debug = 3; //DS test
const bool regeneratePlots = true; //DS test

const double kScaleFactorADC = 2.15; // for amp, from ADC to mV

const int kNSteps = 13;
const double kmVSteps[kNSteps] = { 179.9316, 170.0439, 157.5928, 150.0244, 140.0146,
				   130.0049, 119.9951, 109.9854, 90.08789, 69.94629,
				   50.04883, 30.0293, 9.887695 };

// for 4 mV = MuCh settings
const int kNStepsMuCh = 12;
/* 
   const int AMuCh[kNStepsMuCh] = { 2063, 1950, 
                                    1836, 1719, 1606, 1491, 1376, 
				    1261, 1034, 802, 574, 346 };
   kmVStepsMuCh = 5000 * (AMuCh / 16384.0) ; 
*/
const double kmVStepsMuCh[kNStepsMuCh] = { 629.5776, 595.0928,
					   560.3027, 524.5972, 490.1123, 455.0171, 419.9219,
					   384.8267, 315.5518, 244.7510, 175.1709, 105.5908 };


/* // 10 pF settings:
// const double kScaleFactor = 2.35/(2.35+47.0) * 10.0; // theory, single
//const double kScaleFactor = 2.42/(49.5) * 10.0; // measured
const double kScaleFactor = 2.05/(2.05+47.0) * 10.0; // odd/even pulser
*/
// 1 pF settings:
//const double kPulseGenVoltScale = 1.5; // for ALICE relative to settings for 10 pF above
const double kPulseGenVoltScale = 1.25; // for sPHENIX : 1.5*2.5/3.0
const double kScaleFactor = kPulseGenVoltScale * 11.0/(11.0 + 39.0) * 1.172; // 1.172 factor from calibration in lab (before odd/even pulser) 
const double kScaleFactorMuCh = 11.0/(11.0 + 39.0) * 1.172; // 1.172 factor from calibration in lab (before odd/even pulser) 

const int kMaxPeakVal = 1000; // values should not be larger than this
const int kLimit = 5;
const double kPrecisionfC = 0.1; // there is some uncertainty on the fC scale, but mostly on the absolute level

const int kNChan = 32;
//const int kMaxLineLength = kNChan * 10; // max 4 digits + 1 space per channel + few extra for margin
const int kMaxLineLength = 100000;

const int kFirstPed = 3; // first sample for baseline calculation
const int kLastPed = 12; // last sample for baseline calculation
const int kSpacing = 14; // distance between (later) peaks
const int kSpacingWindow = 2; // +-2 for spacing tolerance 
const int kMaxSample = 300 ; // one full staircase should be about 231 samples, i.e. significantly less than 300 samples
const int kMinZeroSamples = 12; // should be at least 12 zeroes in a row (old 25)
const int kMinZerosLine = 8; // half the channels on a line should be zeros (odd or even) for a zeroline

int analyse_file(const int iChip, const char *config,
		 const char *infile, const char *outfile)
{
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0111);
  gStyle->SetStatX(0.5); 

  FILE* fout = fopen(outfile, "w");
  fprintf(fout, "%9s %13s\n", "Channel #", " Calib (mV/fC)");

  char outdir[50];
  char command[100];

  if (regeneratePlots) {
    sprintf(outdir, "png/calib/%d/%s", iChip, config);
    sprintf(command, ".! mkdir -p %s", outdir);
    gROOT->ProcessLine(command);
    cout << " executing command " << command << endl;
  }

  const char * detString[] = { "TPC", "MuCh" };
  int idetConfig = 0; // default TPC
  int kMaxBin = 40; // location of first peak from block of zeros end (old: 27)
  int kMinSignalSamples = 225; // should be 231 non-zero values in a row for 20 or 30 mV, 205 for 4 mV
  int kMaxWindow = 5; // +/- window around expected maxbin
  int kNPeaks = kNSteps;
  int kNPeaksSkip = 0;
  int kMinFit = 1; // min fC value
  int kMaxFit = 99; // max fC value

  if ( strstr(config, "4mV") ) {
    idetConfig = 1; // 4mV is for MuCh config
    kMaxBin = 25;
    kMaxWindow = 10;
    kMinSignalSamples = 200;
    kNPeaks = kNStepsMuCh; 
    kNPeaksSkip = 2; // skip first 2 peaks for MuCh setting
    kMinFit = 1; // min fC value
    kMaxFit = 150; // max fC value
  }
  cout << " idetConfig = " << idetConfig << " " << detString[idetConfig] << endl;

  ifstream inADC;
  inADC.open(infile);
  cout << " reading file " << infile << endl;
  cout << " writing file " << outfile << endl;
    
  // signal stuff
  int ichan = 0;

  // local var in loop
  double amp = 0;
  double base = 0;
  int a = 0;
  int ipeak = 0;
  
  char hname[100];
  char hdesc[100];
  TProfile* hPeaks[32];
  TProfile* hMax[32];
  for (ichan=0; ichan<32; ichan++) {
    sprintf(hname, "hPeaks_%d", ichan);
    sprintf(hdesc, "Peaks - Chan %d", ichan);
    hPeaks[ichan] = new TProfile(hname, hdesc, 
				 kNPeaks+1, -0.5, kNPeaks + 0.5);

    sprintf(hname, "hMax_%d", ichan);
    sprintf(hdesc, "Max - Chan %d", ichan);
    hMax[ichan] = new TProfile(hname, hdesc, 
			       kNPeaks+1, -0.5, kNPeaks + 0.5);
  }
  TProfile* hBase = new TProfile("hBase", "hBase", 32, -0.5, 31.5);
  
  int nlines  = 0; // over whole file

  // sample loop, inside block of interesting values, between zeros
  int nsamples = 0;
  char line[kMaxLineLength];
  int samples[kNChan];
  int samplesMem[kNChan][kMaxSample];
  
  int nZeros = 0; // no zero values seen yet

  while ( inADC.good() ) {
    // cout << "nlines " << nlines << endl;
    inADC.getline(line, kMaxLineLength);
    if (! inADC.good()) break;
    nlines++;      

    // split line into the expected 32 pieces 
    int ns = 0;
    int nsZeros = 0; // this line
    char * pch;
    pch = strtok (line," ");
    while ( (pch != NULL) && (ns<kNChan) ) {
      // printf ("ns %d pch %s\n",ns, pch);
      samples[ns] = atoi(pch);
      //      if (samples[ns] == 0) nsZeros++; // actually zero
      if (samples[ns] < 40) nsZeros++; // effectively below any ped values
      pch = strtok (NULL, " ");
      ns++;
    }
    // cout << " ns " << ns << endl;

    // check for sufficient number of zeros for this sample line (should be at least close to half; pulsing odd or even) 
    if ( nsZeros > kMinZerosLine ) {    
      if (nsamples > kMinSignalSamples)  { // some previous interesting data; analyze this sample block
	
	for (ichan=0; ichan<32; ichan++) {

	  //--get baseline
	  base = 0;                          
	  int nbase = 0;
	  if (idetConfig == 0) {
	    // Take baseline among first samples before pulse (not quite the first))                                
	    for(Int_t a = kFirstPed; a <= kLastPed; a++) {
	      base += samplesMem[ichan][a];
	      nbase++;
	    }
	  }
	  else {
	    // Take baseline among last samples (not quite the last))                                
	    for(Int_t a = kFirstPed; a <= kLastPed; a++) {
	      int ab = nsamples - a;
	      base += samplesMem[ichan][ab];
	      nbase++;
	    }
	  }
	  base /= nbase;
	  if (debug>2) printf("ichan %d baseline estimate %3.1f\n", ichan, base);
	  hBase->Fill(ichan, base);

	  int iMaxBin = kMaxBin;
	  // allow for some wiggle room around the expected maxbin position (+- 2 positions)
	  int iMaxVal = samplesMem[ichan][iMaxBin];
	  int ibin = - kMaxWindow;
	  while (ibin <= kMaxWindow) {
	    if (samplesMem[ichan][kMaxBin + ibin] > iMaxVal) {
	      iMaxBin = kMaxBin + ibin;
	      iMaxVal = samplesMem[ichan][iMaxBin];
	    }
	    ibin++;
	  }

	  // get peak info
	  for (ipeak=0; ipeak<kNPeaks; ipeak++) {
	    int iPeakMaxBin = iMaxBin + ipeak*kSpacing;
	    iMaxVal = samplesMem[ichan][iPeakMaxBin];

	    // do a max search around suggested ibin
	    int locPeakMaxBin = iPeakMaxBin;
	    ibin = - kSpacingWindow;
	    while (ibin <= kSpacingWindow) {
	      if (samplesMem[ichan][iPeakMaxBin + ibin] > iMaxVal) {
		locPeakMaxBin = iPeakMaxBin + ibin;
		iMaxVal = samplesMem[ichan][locPeakMaxBin];
	      }
	      ibin++;
	    }

	    ibin = locPeakMaxBin; 	    // for easy indexing below 
	    amp = samplesMem[ichan][ibin] - base;
	    if (ipeak>= kNPeaksSkip) {
	      hPeaks[ichan]->Fill(ipeak, amp);
	      hMax[ichan]->Fill(ipeak, samplesMem[ichan][ibin] );
	      if (debug>2) printf("ipeak %d ibin %d maxVal %d amp %3.0f\n", ipeak, ibin, samplesMem[ichan][ibin], amp);
	    } 
	  }
	    
	  int baseBin = iMaxBin + kNPeaks*kSpacing + kSpacing/2;
	  if (baseBin<nsamples) {
	    amp = samplesMem[ichan][baseBin] - base;
	    hPeaks[ichan]->Fill(kNPeaks, amp);
	    hMax[ichan]->Fill(kNPeaks, samplesMem[ichan][baseBin] );
	    if (debug>2) {
	      printf("baseBin %d amp %3.0f\n",
		     baseBin, amp);
	    }
	  }	    
	} // ichan
	
	nsamples = 0; // analysis done for this block
      } // samples > kMinSignalSamples 
      nZeros++;      
    }
    else { // non-zero values
      if (nZeros >= kMinZeroSamples) { // aha, looks like start of new block
	nsamples = 0;
      }
      if (nsamples >= kMaxSample) { // if we have crappy data, we may not find the zero-samples - then just reset when memory is full
	nsamples = 0;
      }
      nZeros = 0;
      // store samples from now
      for (ichan=0; ichan<32; ichan++) {
	samplesMem[ichan][nsamples] = samples[ichan];
      }
      nsamples++;
    }
    if (debug>2) printf("line %d sample0 adc %d nsamples %d nZerosInARow %d\n", nlines, samples[0], nsamples, nZeros);
  } // inADC

    // go through all the maxVal and baseline info and fill graphs
  TGraphErrors *g1[kNChan];
  int nOK[kNChan] = {0};
  double rms = 0;

  for (ichan=0; ichan<32; ichan++) { 

    // make a graph for each channel
    g1[ichan] = new TGraphErrors(kNPeaks + 2); 
    g1[ichan]->SetMarkerStyle(20);
    //   g1[ichan]->SetMarkerColor(ichan + 1);
    g1[ichan]->SetTitle("Max - Baseline");
    
    if (debug>0) { 
      cout << " ichan " << ichan << " : " 
	   << hPeaks[ichan]->GetEntries() << " values found" << endl;
    }
    
    if ( hPeaks[ichan]->GetEntries() > 1 ) {

      base = hBase->GetBinContent(ichan+1);

      /* -- skip adding zero point = baseline for now..
      // add baseline info
      amp = hPeaks[ichan]->GetBinContent(kNPeaks+1);
      rms = hPeaks[ichan]->GetBinError(kNPeaks+1);
      g1[ichan]->SetPoint(nOK[ichan], 0, amp);
      g1[ichan]->SetPointError(nOK[ichan], kPrecisionfC, rms);
      if (debug>2) {
      printf("ichan %d base amp %4.1f rms %4.1f\n",
      ichan, amp, rms);
      }
      nOK[ichan]++;
      */
      g1[ichan]->SetPoint(nOK[ichan], 0, 0);
      g1[ichan]->SetPointError(nOK[ichan], kPrecisionfC, 0.1);
      nOK[ichan]++;

      // figure out amp and rms for each step a.k.a. peak
      for (ipeak=kNPeaks-1; ipeak>=0; ipeak--) {	
	//	if ( (hPeaks[ichan]->GetBinContent(ipeak+1) + base) < kMaxPeakVal ) {
	if ( hMax[ichan]->GetBinContent(ipeak+1) < kMaxPeakVal ) {
	  // add points to graph
	  double fC_estimate = kmVSteps[ipeak] * kScaleFactor;
	  if (idetConfig == 1) { // special MuCh 
	    fC_estimate = kmVStepsMuCh[ipeak] * kScaleFactorMuCh;
	  }
	  amp = kScaleFactorADC * hPeaks[ichan]->GetBinContent(ipeak+1);
	  rms = kScaleFactorADC * (hPeaks[ichan]->GetBinError(ipeak+1) + 0.1); // smear a little 
	  if (debug>1) {
	    printf("ichan %d ipeak %d fC_estimate %4.1f amp %4.1f rms %4.1f\n",
		   ichan, ipeak, fC_estimate, amp, rms);
	  }
	  // add points to graph
	  g1[ichan]->SetPoint(nOK[ichan], fC_estimate, amp); 
	  g1[ichan]->SetPointError(nOK[ichan], kPrecisionfC, rms);
	  nOK[ichan]++;
	}
      }
    }      // some fillings for this channel
  } // ichan loop

    // graphs are filled - now we do the fitting...
  
    // fit channel info, and store fit parameters (with errors) and chi2 values
  char name[50], description[100];
  TCanvas* c1 = new TCanvas("c1", "Amplitude vs Sig", 900, 600);
  TCanvas* c2 = new TCanvas("c2", "Amplitude vs Sig - Chip", 900, 600);
  
  TGraphErrors *gp0 = new TGraphErrors(kNChan);
  gp0->SetMarkerStyle(20);
  gp0->SetTitle("Fit Par 0");
  TGraphErrors *gp1 = new TGraphErrors(kNChan);
  gp1->SetMarkerStyle(20);
  gp1->SetTitle("Fit Par 1");
  
  gStyle->SetOptFit(0111);
  
  for (ichan=0; ichan<kNChan; ichan++) {
    if (debug>1) cout << " ichan " << ichan << " nOK " << nOK[ichan] << endl; 
    if (nOK[ichan] < 2) continue; // don't try to fit an empty plot
    TF1 *func1 = new TF1("func1","[0] + [1]*x", kMinFit, kMaxFit); // skip 0 point
    func1->SetParameter(0, 0);
    func1->SetParameter(1, 20);

    c1->cd();
    g1[ichan]->Fit("func1","RQ");

    fprintf(fout, "%2d%7s %5.2f +/- %3.2f\n", ichan, "", func1->GetParameter(1), func1->GetParError(1) );
    
    if (regeneratePlots) {
      g1[ichan]->Draw("A*");
      sprintf(description, "Amplitude vs Sig for chip %d chan %d", iChip, ichan);
      g1[ichan]->GetXaxis()->SetTitle("signal (fC)");
      if (kScaleFactorADC>2) {
	g1[ichan]->GetYaxis()->SetTitle("amplitude (mV)");
      }
      else {
	g1[ichan]->GetYaxis()->SetTitle("amplitude (ADC)");
      }
      g1[ichan]->SetTitle(description);
      c1->Modified();
      sprintf(name, "%s/amp_sig_chip_%d_chan_%02d.png", outdir, iChip, ichan);
      c1->SaveAs(name);
    }
    gp0->SetPoint(ichan, ichan, func1->GetParameter(0));
    gp0->SetPointError(ichan, 0, func1->GetParError(0));
    gp1->SetPoint(ichan, ichan, func1->GetParameter(1));
    gp1->SetPointError(ichan, 0, func1->GetParError(1));
    
    if (func1) delete func1;
  } // ichan

  gStyle->SetOptFit(0);
  
  if (regeneratePlots) {
    // plot fit info for all channels in chip
    c2->cd();
    sprintf(description, "Fit Par 0 : Chip %d", iChip);
    gp0->SetTitle(description);
    //    gp0->Fit("pol0","Q");
    gp0->Draw("A*");
    gp0->GetXaxis()->SetTitle("Chan #");
    gp0->GetYaxis()->SetTitle("Fit Par0");
    c2->Modified();
    sprintf(name, "%s/amp_sig_chip_%d_fitpar0.png", outdir, iChip);
    c2->SaveAs(name);
    
    c2->cd();
    sprintf(description, "Fit Par 1 : Chip %d", iChip);
    gp1->SetTitle(description);
    //    gp1->Fit("pol0","");
    gp1->Draw("A*");
    gp1->GetXaxis()->SetTitle("Chan #");
    gp1->GetYaxis()->SetTitle("Fit Par1");
    c2->Modified();
    sprintf(name, "%s/amp_sig_chip_%d_fitpar1.png", outdir, iChip);
    c2->SaveAs(name);
  }

  // clean up
  for (ichan=0; ichan<kNChan; ichan++) {
    if (g1[ichan]) delete g1[ichan];
  }
  if (gp0) delete gp0;
  if (gp1) delete gp1;
  
  fclose(fout);
  return 0;
}

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

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

  int iret = analyse_file(iChip, argv[2], infile, outfile);
  return iret;
}