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

void staircaseWithTree(const char *listname="staircaseWithTreelist.txt", 
		       const char *outname="staircaseWithTree.root")
{
  gStyle->SetOptStat(0);
  //  gStyle->SetOptFit(0);
  gStyle->SetOptFit(0111);
  gStyle->SetStatX(0.5); 

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

  const int zoomWindow=5; // plot zoom in on peak
  const double kScaleFactorADC = 2.15; // for amp, from ADC to mV
  //const double kScaleFactorADC = 1; // we do this conversion in fitStudy.C instead

  const int kNSteps = 13;
  /*
  // 1 pF (effective cap seems more like 1.172; tuned with reference chip), with 11 and 39 Ohm; fewer steps for V3 
  const int kmVSteps[kNSteps] = { 360, 340, 320, 300, 280, 260, 240, 220, 180, 140, 100, 60, 20 }; 
  const double kScaleFactor = 11.0/(11.0 + 39.0) * 1.172; // multiply pulsegenerator mV with 
  */
  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 };
  const double kScaleFactor = 2.35/(2.35+47.0) * 10.0; // theory
  //  const double kScaleFactor = 2.42/(49.5) * 10.0; // measured

  const int kMaxPeaks = 20; // should really only be kNSteps intervals but use a lot of extras to avoid crashes 
  const int kLimit = 5;

  const double kPrecisionfC = 1; // there is some uncertainty on the fC scale

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

  const int kFirstPed = 3; // first sample for baseline calculation
  const int kLastPed = 12; // last sample for baseline calculation
  const int kMaxBin = 27; // location of first peak
  const int kSpacing = 14; // distance between (later) peaks

  const int kMaxSample = 300 ; // one full staircase should be about 231 samples, i.e. significantly less than 300 samples
  const int kMinSignalSamples = 225; // should be 231 non-zero values in a row
  const int kMinZeroSamples = 20; // should be 25 zeroes in a row

  char infostr[100];
 
  // variables common for all channels
  int iChip = 0;

  // fit variables
  int iChan = 0;
  float fP0 = 0;
  float fP0Err = 0;
  float fP1 = 0;
  float fP1Err = 0;
  float fChi2 = 0;
  int iNDF = 0;

  TFile *fout = new TFile(outname, "RECREATE");
  TTree* tree = new TTree("tree", "fit info");
 
  tree->Branch("iChip",&iChip,"iChip/I");
  tree->Branch("iChan",&iChan,"iChan/I");
  tree->Branch("fP0",&fP0,"fP0/F");
  tree->Branch("fP0Err",&fP0Err,"fP0Err/F");
  tree->Branch("fP1",&fP1,"fP1/F");
  tree->Branch("fP1Err",&fP1Err,"fP1Err/F");
  tree->Branch("fChi2",&fChi2,"fChi2/F");
  tree->Branch("iNDF",&iNDF,"iNDF/I");

  ifstream input;
  input.open(listname);
  char filename[200];
  char outdir[50];
  char command[100];

  int nfiles = 0;
  while (1) {
    input >> iChip >> filename;
    if (! input.good()) break;
    nfiles++;
        
    cout << " iChip " << iChip 
	 << " filename " << filename << endl;

    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, "png/calib/%s/%d", chipVer, iChip);

      sprintf(command, ".! mkdir -p %s", outdir);
      gROOT->ProcessLine(command);
      cout << " executing command " << command << endl;
    }

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

    // local var in loop
    int maxBin = 0;
    int maxVal = 0;
    double amp = 0;
    double base = 0;
    int a = 0;
    int ipeak = 0;

    char hname[100];
    char hdesc[100];
    TProfile* hPeaks[32];
    for (ichan=0; ichan<32; ichan++) {
      sprintf(hname, "hPeaks_%d", ichan);
      sprintf(hdesc, "Peaks - Chan %d", ichan);
      hPeaks[ichan] = new TProfile(hname, hdesc, 
				   kNSteps+1, -0.5, kNSteps + 0.5);
    }

    int nlines  = 0; // over whole file

    // sample loop, inside block of interesting values, between zeros
    int i; // sample counter
    int nsamples = 0;
    char line[kMaxLineLength];
    int samples[kNChan];
    int samplesMem[kNChan][kMaxSample];

    int nZeros = 0; // no zero values seen yet

    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++;
      }
      // cout << " ns " << ns << endl;

      // check for zeroes for first and last channels
      if ( (samples[0] == 0) && (samples[31] == 0) ) {    
	if (nsamples > kMinSignalSamples)  { // some previous interesting data; analyze this sample block

	  for (ichan=0; ichan<32; ichan++) {
	    maxVal = 0;
	    maxBin = 0;
	    for(Int_t l = 0; l < nsamples; l++) {
	      if (maxVal<samplesMem[ichan][l]) {
		maxVal = samplesMem[ichan][l];
		maxBin = l;
	      }
	    }

	    if (debug>1) {
	      cout << " maxBin " << maxBin << endl;
	    }

	    //--get baseline
	    base = 0;                          
	    int nbase = 0;
	    // Take baseline among first samples before pulse (not quite the first))                                
	    if(maxBin > 20){ 
	      for(Int_t a = kFirstPed; a <= kLastPed; a++) {
		base += samplesMem[ichan][a];
		nbase++;
	      }
	    }
	    base /= nbase;
	    if (debug>1) printf("ichan %d baseline estimate %3.1f\n", ichan, base);
	
	    // get peak info
	    int maxBinPeak[kMaxPeaks] = { 0 };
	    int maxValPeak[kMaxPeaks] = { 0 };

	    int nPeaks = 0;
	    int ibin = maxBin - 3; // first pulse in staircase should be the maxBin 
	    while (ibin < nsamples) {
	      // a maximum should be bigger than the earlier samples and later samples          
	      // and the neighbor samples should also be bigger than their neighbors            
	      if ( ((samplesMem[ichan][ibin-1]+kLimit) < samplesMem[ichan][ibin]) &&
		   ((samplesMem[ichan][ibin+1]+kLimit) < samplesMem[ichan][ibin]) &&
		   ((samplesMem[ichan][ibin+2]+kLimit) < samplesMem[ichan][ibin+1]) ) {

		bool spacingOK = false;
		if (nPeaks == 0) {
		  if (ibin == kMaxBin) spacingOK = true;
		}
		else if (nPeaks > 0) {
		  int spacing = ibin - maxBinPeak[nPeaks-1];
		  if (spacing == kSpacing) {
		    spacingOK = true;
		  }
		}

		if (spacingOK) {
		  maxBinPeak[nPeaks] = ibin;
		  maxValPeak[nPeaks] = samplesMem[ichan][ibin];		
		  nPeaks++;
		}
	      }
	      ibin++;
	    } // ibin loop

	    // if expected amount of peaks, fill TProfile
	    if (debug>1) printf("ichan %d nPeaks %d\n", ichan, nPeaks);
	    if (nPeaks == kNSteps) {
	      for (ipeak=0; ipeak<nPeaks; ipeak++) {
		amp = maxValPeak[ipeak] - base;
		hPeaks[ichan]->Fill(ipeak, amp);
		if (debug>1) {
		  printf("ipeak %d maxBin %d amp %3.0f\n",
			 ipeak, maxBinPeak[ipeak], amp);
		}
	      } // peak loop

	      int baseBin = maxBinPeak[nPeaks-1] + kSpacing/2;
	      if (baseBin<nsamples) {
		amp = samplesMem[ichan][baseBin] - base;
		hPeaks[ichan]->Fill(kNSteps, amp);
		if (debug>2) {
		  printf("baseBin %d amp %3.0f\n",
			 baseBin, amp);
		}
	      }

	    } // all peaks found
	  } // ichan
 
	} // samples > kMinSignalSamples 
	nZeros++;      
      }
      else { // non-zero values
	if (nZeros > kMinZeroSamples) { // aha, looks like start of new block
	  nsamples = 0;
	}
	nZeros = 0;
	// store samples from now
	for (ichan=0; ichan<32; ichan++) {
	  samplesMem[ichan][nsamples] = samples[ichan];
	}
	nsamples++;
      }
      if (debug>1) printf("line %d sample0 adc %d nsamples %d nZerosInARow %d\n", nlines, samples[0], nsamples, nZeros);
    }

    // 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(kNSteps + 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() > 0 ) {

	/*
	// add baseline info
	amp = hPeaks[ichan]->GetBinContent(kNSteps+1);
	rms = hPeaks[ichan]->GetBinError(kNSteps+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=kNSteps-1; ipeak>=0; ipeak--) {	
	// add points to graph
	  double fC_estimate = kmVSteps[ipeak] * kScaleFactor;
	  amp = kScaleFactorADC * hPeaks[ichan]->GetBinContent(ipeak+1);
	  rms = kScaleFactorADC * (hPeaks[ichan]->GetBinError(ipeak+1) + 0.1); // smear a little 
	  if (debug>2) {
	    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", 1, 99); // skip 0 point
      func1->SetParameter(0, 0);
      func1->SetParameter(1, 9.5);

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

      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);
	
	cout << " ichip " << iChip << " ichan " << ichan 
	     << " fit par1 " << func1->GetParameter(1)
	     << " error " << func1->GetParError(1) << endl;
      }

      iChan = ichan;
      fP0 = func1->GetParameter(0);
      fP0Err = func1->GetParError(0);
      fP1 = func1->GetParameter(1);
      fP1Err = func1->GetParError(1);
      fChi2 = func1->GetChisquare();
      iNDF = func1->GetNDF();
      tree->Fill();

      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;

  }
  cout << " nfiles " << nfiles << endl;

  fout->Write();
}