#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 kNPadSweepX = 1;
const int kNPadSweepY = 4;
const int kNPlotSweep = kNPadSweepX * kNPadSweepY;

// for normal running:
const int debug = 0; 
const bool regeneratePlots = false;

// for debug running:
/*
const int debug = 1; //DS test
const bool regeneratePlots = true; //DS test
*/

const char *plotInfo = "last"; const int kSelectedChanSweep[kNPlotSweep] = { 28, 29, 30, 31}; // last, most sensitive
//const char *plotInfo = "first"; const int kSelectedChanSweep[kNPlotSweep] = { 0, 1, 2, 3}; // first, less sensitive?
//const char *plotInfo = "mid"; const int kSelectedChanSweep[kNPlotSweep] = { 14, 15, 16, 17}; // first, less sensitive?

//const char *plotFormat = "pdf";
const char *plotFormat = "png";

// constants; reading data
const int kNChan = 32;
const int kMaxLineLength = kNChan * 6; // max 4 digits + 1 space per channel + 1 extra for margin
const int kNSkip = 0; // skip some first samples

const int kADCWindow = 10; // +- window size to use for fit around mean ADC value

// for sweep plots
const int kNBinsPerSweep = 1000; // 0.2 us per bin  
//const int kNSweep = 36;   
const int kNSweep = 6;   

// for exclusion of periodic noise due to readout clock
const int kReadoutFreqNoise = 125; // 125 timebins = 25 us (5 MHz clock)
const int kExcludeTimebinWindow = 20; // exclude +- timebins around readout noise areas

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

  FILE* fout = fopen(outfile, "a");
  fprintf(fout, "%9s %12s %12s\n", "Channel #", " gPed     gRMS", "   hPed     hRMS");

  string data(outfile);
  string toSearch("noise");
  string replaceStr("debug");
  // Get the first occurrence
  size_t pos = data.find(toSearch);
  // Repeat till end is reached
  while( pos != std::string::npos)
    {
      // Replace this occurrence of Sub String
      data.replace(pos, toSearch.size(), replaceStr);
      // Get the next occurrence from the current position
      pos = data.find(toSearch, pos + toSearch.size());
    }

  FILE* foutDebug;
  if (debug > 0) {
    foutDebug = fopen(data.c_str(), "w");
    fprintf(foutDebug, "Chan #    Ped+/-Err          RMS+/-Err\n");
  }

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

  if (regeneratePlots) {
    char chipVer[10];
    if (iChip<2000) {
      sprintf(chipVer, "V2");
    }
    else if (iChip>=2000 && iChip<6000) {
      sprintf(chipVer, "V3");
    }
    else if (iChip>=6000) {
      sprintf(chipVer, "V4");
    }

    sprintf(outdir, "%s/calib/%s/%d", plotFormat, chipVer, iChip);
    sprintf(command, ".! mkdir -p %s", outdir);
    gROOT->ProcessLine(command);
    cout << " executing command " << command << endl;
  }

  ifstream inADC;
  inADC.open(infile);
  cout << " reading file " << infile << endl;
    
  // counter stuff
  int ichan = 0;
  int isweep = 0;
  int ibin = 0; // in sweep

  char hname[100];
  char hdesc[100];
  TH1F* hSample[32];
  TH1F* hSweep[32][kNSweep];

  for (ichan=0; ichan<32; ichan++) {
    sprintf(hname, "hSample_%d", ichan);
    sprintf(hdesc, "All ADC values - Chan %d", ichan);
    hSample[ichan] = new TH1F(hname, hdesc, 1024, -0.5, 1023.5);

    for (isweep=0; isweep<kNSweep; isweep++) {
      sprintf(hname, "hSweep_%d_%d", ichan, isweep);
      sprintf(hdesc, "TimeSweep %d (%d bins) - Chan %d", isweep, kNSweep, ichan);
      hSweep[ichan][isweep] = new TH1F(hname, hdesc, 
				       kNBinsPerSweep, -0.5, kNBinsPerSweep - 0.5);
    }
  }
  TProfile* hBase = new TProfile("hBase", "hBase", 32, -0.5, 31.5);
  
  int nlines  = 0; // over whole file
  int samples[kNChan];
  int nsamples = 0;
  char line[kMaxLineLength];

  // 1st sample loop: find min and max values for all channels
  int sampleMax[kNChan];
  int sampleMin[kNChan];
  int sampleMaxPos[kNChan];
  int sampleMinPos[kNChan];
  int modLine = 0;

  for (ichan=0; ichan<kNChan; ichan++) {
    sampleMax[ichan] = 0;
    sampleMin[ichan] = 1023;
    sampleMaxPos[ichan] = 0;
    sampleMinPos[ichan] = 0;
  }

  while ( inADC.good() && nlines<1000 ) {
    inADC.getline(line, kMaxLineLength);
    if (! inADC.good()) break;
    nlines++;      

    // split line into the expected 32 pieces 
    int ns = 0;
    char * pch;
    pch = strtok (line," ");
    while (pch != NULL) {
      // printf ("ns %d pch %s\n",ns, pch);
      samples[ns] = atoi(pch);
      pch = strtok (NULL, " ");
      ns++;
    }
    if (debug>2) { cout << "nlines " << nlines << " ns " << ns << endl; }

    for (ichan=0; ichan<kNChan; ichan++) {
      if (sampleMax[ichan] < samples[ichan]) { 
	sampleMax[ichan] = samples[ichan];
	sampleMaxPos[ichan] = nlines;
      }
      if (sampleMin[ichan] > samples[ichan]) { 
	sampleMin[ichan] = samples[ichan];
	sampleMinPos[ichan] = nlines;
      }
    }
  } // inADC
  inADC.close();

  // next: define exclusion window
  // where's the noise at? Fill histo with mean and max timebins (% kReadoutFreqNoise)
  TH1I *h1MaxMod = new TH1I("h1MaxMod", "MaxMod locations", kReadoutFreqNoise, -0.5, kReadoutFreqNoise - 0.5);
  TH1I *h1MinMod = new TH1I("h1MinMod", "MinMod locations", kReadoutFreqNoise, -0.5, kReadoutFreqNoise - 0.5);
  int modLineMin = 0;
  for (ichan=0; ichan<kNChan; ichan++) {
    modLine = sampleMaxPos[ichan] % kReadoutFreqNoise;
    modLineMin = sampleMinPos[ichan] % kReadoutFreqNoise;
    if (debug>-1) cout << " ichan " << ichan << " MaxPos " << sampleMaxPos[ichan] << " MinPos " << sampleMinPos[ichan] << " modLine " << modLine << " modLineMin " << modLineMin << endl;
    h1MaxMod->Fill( modLine );
    h1MinMod->Fill( modLineMin );
  }
  int excludeTimebin = h1MaxMod->GetMaximumBin();
  if (debug > -1) {
    cout << " exclusion info: " << endl
	 << " h1MaxMod MaximumBin " << h1MaxMod->GetMaximumBin()   
	 << " Mean " << h1MaxMod->GetMean() << endl  
	 << " h1MinMod MaximumBin " << h1MinMod->GetMaximumBin()   
	 << " Mean " << h1MinMod->GetMean() << endl;
  }

  // 2nd sample loop; exclude some samples

  inADC.open(infile);
  cout << " reading file 2nd loop " << infile << endl;

  nlines = 0; // start again
  while ( inADC.good() ) {
    inADC.getline(line, kMaxLineLength);
    if (! inADC.good()) break;
    nlines++;      

    // split line into the expected 32 pieces 
    int ns = 0;
    char * pch;
    pch = strtok (line," ");
    while (pch != NULL) {
      // printf ("ns %d pch %s\n",ns, pch);
      samples[ns] = atoi(pch);
      pch = strtok (NULL, " ");
      ns++;
    }

    bool keep = true;
    if (nlines < kNSkip) { keep = false; }
    /*
    modLine = nlines % kReadoutFreqNoise;
    int diffLine = TMath::Abs(modLine - excludeTimebin);
    if ( diffLine < kExcludeTimebinWindow ) { keep = false; } 
    */

    if (keep) {
      if (debug>2) { cout << "Fill hSample nlines " << nlines  << endl; }
      for (ichan=0; ichan<32; ichan++) {
	hSample[ichan]->Fill(samples[ichan]);
	hBase->Fill(ichan, samples[ichan]);
      }
      int iline = nlines - kNSkip - 1;
      isweep = iline/kNBinsPerSweep;
      ibin = iline % kNBinsPerSweep;
      if (debug>2) { cout << " iline " << iline << " isweep " << isweep << " ibin " << ibin  << endl; }
      if (isweep < kNSweep) {
	for (ichan=0; ichan<32; ichan++) {
	  hSweep[ichan][isweep]->SetBinContent(ibin+1, samples[ichan]);
	}
      }
    }
  } // inADC
  inADC.close();

  float mean = 0; 
  float rms = 0; 
  for (ichan=0; ichan<32; ichan++) { 
    hSample[ichan]->GetXaxis()->SetRangeUser(hBase->GetBinContent(ichan+1) - kADCWindow,
					     hBase->GetBinContent(ichan+1) + kADCWindow );
    mean = hSample[ichan]->GetMean();
    rms = hSample[ichan]->GetRMS();

    if (mean > 0){
      hSample[ichan]->Fit("gaus","Q"," ", mean-kADCWindow, mean+kADCWindow);
      TF1 *fit_gaus = hSample[ichan]->GetFunction("gaus");
      float sigma = fit_gaus->GetParameter(2);
      if (debug > 0) {
	printf("%2d%7s %5.1f     %3.2f   %5.1f     %3.2f\n", ichan, "", 
	       fit_gaus->GetParameter(1), sigma, mean, rms);
      }

      fprintf(fout, "%2d%7s %5.1f     %3.2f   %5.1f     %3.2f\n", ichan, "", 
	      fit_gaus->GetParameter(1), sigma, mean, rms);
    }
    else {
      fprintf(fout, "%2d%7s %5.1f     %3.2f\n", ichan, "", mean, rms);
    }

    if (debug > 0) {
      fprintf(foutDebug, "%2d%7s %6.3f +/- %4.3f   %4.3f +/- %4.3f\n", ichan, "", 
	      hSample[ichan]->GetMean(), hSample[ichan]->GetMeanError(), 
	      hSample[ichan]->GetRMS(), hSample[ichan]->GetRMSError() );
    }

  }

  if (regeneratePlots) {
    // make plots
    // fit channel info, and store fit parameters (with errors) and chi2 values
    char name[100], description[100];
    TCanvas* cBase = new TCanvas("cBase", "Baseline vs Chan #", 900, 600);
    TCanvas* cSample = new TCanvas("cSample", "Sample distribution", 1800, 1200);
    TCanvas* cSweep = new TCanvas("cSweep", "Swep distribution", 1800, 1200);
    
    cSample->Divide(8,4);
    cSweep->Divide(kNPadSweepX, kNPadSweepY);
    
    cBase->cd();
    hBase->Draw();
    sprintf(name, "%s/base_chip_%d_%s_%s.%s", outdir, iChip, config, plotInfo, plotFormat);
    cBase->SaveAs(name);
    
    for (ichan=0; ichan<32; ichan++) { 
      hSample[ichan]->SetMarkerStyle(20);
      cSample->cd(ichan + 1);
      hSample[ichan]->GetXaxis()->SetRangeUser(hBase->GetBinContent(ichan+1) - 10,
					       hBase->GetBinContent(ichan+1) + 10 );
      
      sprintf(description, "Sample distr for chip %d chan %d", iChip, ichan);
      hSample[ichan]->GetXaxis()->SetTitle("ADC");
      hSample[ichan]->GetYaxis()->SetTitle("Yield");
      hSample[ichan]->SetTitle(description);
      hSample[ichan]->Draw();
    } // ichan
    cSample->Modified();
    sprintf(name, "%s/sample_chip_%d_%s_%s.%s", outdir, iChip, config, plotInfo, plotFormat);
    cSample->SaveAs(name);
    
    for (isweep=0; isweep<kNSweep; isweep++) {
      int isel = 0;
      for (ichan=0; ichan<32; ichan++) { 
	if (kSelectedChanSweep[isel] == ichan) {
	  hSweep[ichan][isweep]->SetMarkerStyle(20);
	  isel++;

	  cSweep->cd(isel);    
	  sprintf(description, "Sweep %d for chip %d chan %d", isweep, iChip, ichan);
	  cout << description << " cont2 " << hSweep[ichan][isweep]->GetBinContent(2) << endl;
	  hSweep[ichan][isweep]->GetYaxis()->SetTitle("ADC");
	  hSweep[ichan][isweep]->GetXaxis()->SetTitle("Timebin");
	  hSweep[ichan][isweep]->SetTitle(description);
	  hSweep[ichan][isweep]->Draw();
	}
      } // ichan
      cout << " isweep " << isweep << " isel " << isel << endl;
      cSweep->Modified();
      sprintf(name, "%s/sweep_%d_chip_%d_%s_%s.%s", outdir, isweep, iChip, config, plotInfo, plotFormat);
      cSweep->SaveAs(name);
    } // isweep
  } // plots

  fclose(fout);
  if (debug > 0) fclose(foutDebug);
  return 0;
}

int main (int argc, char *argv[])
{
  if (argc<4) {
    printf("not enough arguments: argc = %d\nSyntax should be:\n./noise.exe chipstring config adcfilename\n",argc);
    printf("Example:\n./noise.exe 100000 30mV_noise /home/tpc/robot/20180725/123510_sampa100000/123510_sampa100000_trorc00_link00_30mV_noise_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/noiseAna_sampa%s_%s.log",  str.substr(0,found).c_str(), argv[1], argv[2]);

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