#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 <fstream>
#include <vector>
#include <algorithm>

#include <cstring>
#include <cerrno> 
#include <ctime>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

void fitStudy(const char *filename="result.root") 
{
  const int kNChan = 32;
  const int kNSync = 8;
  const int kNSta = 2;
  const int kSta[kNSta] = { 1, 2 };
  const int kNGain = 2;
  const int kGain[kNGain] = { 20, 30 };

  // acceptable limits for P1/Gain 
  const float kMinGain[kNGain] = { 18, 27 }; 
  const float kMaxGain[kNGain] = { 22, 33 };
  const float kMinGainAll = 0;
  const float kMaxGainAll = 50;

  char name[100];
  char outdir[100];
  sprintf(outdir, "png/calib");

  char suffix1[20];
  char suffix2[20];
  sprintf(suffix1, "pdf");
  sprintf(suffix2, "png");

  const bool regen_corrplots = false;
  // const bool regen_corrplots = true;

  const int kNChips = 100000; // max number of chips so far, for any station
  const int kChipOffset = 100000;

  int nOK[kNSta][kNChips] = {0};
  float meanGain[kNSta][kNChips] = {0};
  float rmsGain[kNSta][kNChips] = {0};
  int nFail[kNSta][kNChips] = {0};

  TFile *f = new TFile(filename);
  cout << " reading file " << filename << endl;
  TTree *tree = (TTree*) f->Get("tree");

  //Declaration of leaves types
  Int_t           iChip;
  Int_t           iSta;
  Float_t         fGain20mV[32];
  Float_t         fGain30mV[32];
  Int_t           iSync[kNSync];

  // Set branch addresses.
  tree->SetBranchAddress("iChip",&iChip);
  tree->SetBranchAddress("iSta",&iSta);
  tree->SetBranchAddress("fGain20mV",fGain20mV);
  tree->SetBranchAddress("fGain30mV",fGain30mV);
  tree->SetBranchAddress("iSync",iSync);

  int nentries = tree->GetEntries();

  // make 2 versions of histograms, for Gain 0 (20 mV) and Gain 1 (30 mV)
  TH1F *h1[kNSta][kNGain];
  TH2F *h2[kNSta][kNGain];
  TProfile *hp1[kNSta][kNGain];
  
  int ista = 0;
  int igain = 0;
  for (ista=0; ista<kNSta; ista++) {
    for (igain=0; igain<kNGain; igain++) {
      h1[ista][igain] = new TH1F( Form("h1%d%d", ista, igain),
				  Form("Gain %dmV : Sta%d", kGain[igain], kSta[ista]),
				  100, kMinGain[igain], kMaxGain[igain] );
      h2[ista][igain] = new TH2F( Form("h2%d%d", ista, igain),
				  Form("Gain %dmV : Sta%d vs Chan #", kGain[igain], kSta[ista]),
				  32, -0.5, 31.5,
				  100, kMinGain[igain], kMaxGain[igain] );

      hp1[ista][igain] = new TProfile(Form("hp1%d%d", ista, igain),
				      Form("Gain %dmV : Sta%d vs Chan #", kGain[igain], kSta[ista]),
				      32, -0.5, 31.5);
    }
  }

  cout << " nentries " << nentries << endl;
  Long64_t nbytes = 0;

  int minChip = kNChips;
  int maxChip = 0;
  int iEntry[kNChips] = {0};
  int locChip = 0;  
  int ichip = 0;
  int ichan = 0;

  for (int i=0; i<nentries;i++) {
    nbytes += tree->GetEntry(i);

    int nSyncFail = 0;
    int nSyncOK = 0;
    for (k=0; k<kNSync; k++) {
      if ( (200<iSync[k]) && (iSync[k]<300) ) { nSyncOK++; }
      else { nSyncFail++; }
    }
    if (debug>1) cout << " nSyncOK " << nSyncOK << endl;
    if ( nSyncOK == kNSync ) {

      locChip = iChip - kChipOffset;  
      if ( (locChip>=0) && (locChip < kNChips) ) { // right range of chips

	if (minChip>locChip) minChip = locChip;
	if (maxChip<locChip) maxChip = locChip;    

	iEntry[locChip] = i; // keep track of latest entry for each chip

	if (iSta == kSta[0]) ista = 0;
	else ista = 1;

	for (ichan=0; ichan<kNChan; ichan++) {
	  if ( (fGain20mV[ichan]>kMinGain[0]) && (fGain20mV[ichan]<kMaxGain[0]) ) {
	    h1[ista][0]->Fill(fGain20mV[ichan]);
	    h2[ista][0]->Fill(ichan, fGain20mV[ichan]);
	    hp1[ista][0]->Fill(ichan, fGain20mV[ichan]);
	  }
	  if ( (fGain30mV[ichan]>kMinGain[1]) && (fGain30mV[ichan]<kMaxGain[1]) ) {
	    h1[ista][1]->Fill(fGain30mV[ichan]);
	    h2[ista][1]->Fill(ichan, fGain30mV[ichan]);
	    hp1[ista][1]->Fill(ichan, fGain30mV[ichan]);
	  }
	}
      }
    }
  }

  // we now know our range of chip numbers
  int nChips = maxChip - minChip + 1;

  ofstream fout;
  fout.open("condcorr.txt");

  // calculate a correction factor for each channel, and station, and gain
  float dcorr[kNSta][kNGain][kNChan] = {0};
  float avecorr[kNSta][kNChan] = {0};
  int ncorr[kNSta][kNChan] = {0};

  for (ista=0; ista<kNSta; ista++) {
    for (ichan = 0; ichan<32; ichan++) {
      avecorr[ista][ichan] = 0;
      ncorr[ista][ichan] = 0;
      for (igain=0; igain<kNGain; igain++) {  
	dcorr[ista][igain][ichan] = 0;
      }
    }
  }

  for (ista=0; ista<kNSta; ista++) {
    for (igain=0; igain<kNGain; igain++) {  
      if (hp1[ista][igain]->GetEntries() > 0) {
	TF1 *func0 = new TF1("func0","[0]", 0, kNChan);
	hp1[ista][igain]->Fit("func0","RQ");
	cout << " ista " << ista << " igain " << igain << " func0 " << func0->GetParameter(0) << endl;
	
	for (ichan = 0; ichan<32; ichan++) {
	  int ibin = hp1[ista][igain]->GetBin(ichan);
	  if (hp1[ista][igain]->GetBinContent(ichan+1) > 0) {
	    dcorr[ista][igain][ichan] = func0->GetParameter(0) / hp1[ista][igain]->GetBinContent(ichan+1);
	    avecorr[ista][ichan] += dcorr[ista][igain][ichan];
	    ncorr[ista][ichan] ++;
	  }
	}
      }
    }

    for (ichan = 0; ichan<32; ichan++) {
      if (ncorr[ista][ichan] > 0) {
	avecorr[ista][ichan] /= ncorr[ista][ichan];
	cout << " kSta " << kSta[ista] << " ichan " << ichan
	     << " avecorr " << avecorr[ista][ichan] 
	     << endl;
	  
	fout << " " << kSta[ista] << " " << ichan
	     << " " << avecorr[ista][ichan] 
	     << endl;
      }
    }
  }
  fout.close();

  TH1F *h1Corr[kNSta];
  TH2F *h2Corr[kNSta];
  TProfile *hp1Corr[kNSta];
  TH2F *h2CorrChip[kNSta];
  TProfile *hp1CorrChip[kNSta];
  TH2F *h2CorrChipAll[kNSta];
  TProfile *hp1CorrChipAll[kNSta]; 
  
  igain = 0; // let's just look at one gain setting for the corrected plots... (enough already)
  for (ista=0; ista<kNSta; ista++) {
    // now, we'll do a second loop, for corrected values...
    h1Corr[ista] = new TH1F(Form("h1Corr%d", ista),
			    Form("Gain (Corr) %dmV : Sta%d", kGain[igain], kSta[ista]),
			    100, kMinGain[igain], kMaxGain[igain] );
    h2Corr[ista] = new TH2F( Form("h2Corr%d", ista),
			     Form("Gain (Corr) %dmV : Sta%d vs Chan #", kGain[igain], kSta[ista]),
			     32, -0.5, 31.5,
			     100, kMinGain[igain], kMaxGain[igain] );    
    hp1Corr[ista] = new TProfile(Form("hp1Corr%d", ista),
				 Form("Gain (Corr) %dmV : Sta%d vs Chan #", kGain[igain], kSta[ista]),
				 32, -0.5, 31.5);

    h2CorrChip[ista] = new TH2F(Form("h2CorrChip%d", ista),
				Form("Gain (Corr) %dmV : Sta%d vs Chip #", kGain[igain], kSta[ista]),			     
				nChips, minChip-0.5 + kChipOffset, maxChip+0.5 + kChipOffset, 
				100, kMinGain[igain], kMaxGain[igain] );
    hp1CorrChip[ista] = new TProfile(Form("hp1CorrChip%d", ista),
				     Form("Gain (Corr) %dmV : Sta%d vs Chip #", kGain[igain], kSta[ista]),			     
				     nChips, minChip-0.5 + kChipOffset, maxChip+0.5 + kChipOffset); 
    
    h2CorrChipAll[ista] = new TH2F( Form("h2CorrChipAll%d", ista),
				    Form("Gain (Corr, All) %dmV : Sta%d vs Chip #", kGain[igain], kSta[ista]),			     		
				    nChips, minChip-0.5 + kChipOffset, maxChip+0.5 + kChipOffset, 
				    100, kMinGainAll, kMaxGainAll );
    hp1CorrChipAll[ista] = new TProfile(Form("hp1CorrChipAll%d", ista),
					Form("Gain (Corr, All) %dmV : Sta%d vs Chip #", kGain[igain], kSta[ista]),			     
					nChips, minChip-0.5 + kChipOffset, maxChip+0.5 + kChipOffset); 
  }
  
  						 
  for (ichip = minChip; ichip<=maxChip; ichip++) {
    if (iEntry[ichip] > 0) {
      nbytes += tree->GetEntry(iEntry[ichip]);
      locChip = iChip - kChipOffset;  
      if (iSta == kSta[0]) ista = 0;
      else ista = 1;      

      int nSum = 0;      
      float sum = 0;
      float sum2 = 0;

      for (ichan = 0; ichan<32; ichan++) {
	float dG = fGain20mV[ichan] * avecorr[ista][ichan];
	//	cout << " line " << __LINE__ << ": ichan " << ichan << " gain " << fGain20mV[ichan] << " avecorr " << avecorr[ista][ichan] << endl;
	h2CorrChipAll[ista]->Fill(locChip + kChipOffset, dG);
	hp1CorrChipAll[ista]->Fill(locChip + kChipOffset, dG);

	sum += dG;
	sum2 += dG * dG;
	nSum ++;

	if ( (kMinGain[0]<dG) && (dG<kMaxGain[0]) ) {
	  h1Corr[ista]->Fill(dG);
	  h2Corr[ista]->Fill(ichan, dG);
	  hp1Corr[ista]->Fill(ichan, dG);
	  
	  h2CorrChip[ista]->Fill(locChip + kChipOffset, dG);
	  hp1CorrChip[ista]->Fill(locChip + kChipOffset, dG);

	  nOK[ista][locChip]++;
	}
	else {
	  cout << " line " << __LINE__ << " : fail ichip " << locChip << " ichan " << ichan << endl; 
	  nFail[ista][locChip]++;
	}

      }

      if (nSum>0) {
	meanGain[ista][locChip] = sum / nSum;
	rmsGain[ista][locChip] = TMath::Sqrt(TMath::Abs(sum2/nSum - meanGain[ista][locChip]*meanGain[ista][locChip]));
      }
    }
  }


  cout << " histos filled " << endl;

  TCanvas* c1 = new TCanvas("c1", "c1", 900, 600);

  for (ista=0; ista<kNSta; ista++) {
    if (hp1[ista][igain]->GetEntries() > kNChan) {
      //if ( h1Corr[ista]->GetEntries() > 0 ) {
      TF1 *func1 = new TF1("func1","[0]", 0, kNChan);
      hp1Corr[ista]->Fit("func1","RQ");
      cout << " ista " << ista << " func1 " << func1->GetParameter(0) << endl;

      c1->SetLogy(1);
      c1->cd();
      h1[ista][igain]->SetMarkerStyle(20);
      h1[ista][igain]->GetXaxis()->SetTitle("Gain (mV/fC)");
      h1[ista][igain]->GetYaxis()->SetTitle("Yield");
      h1[ista][igain]->SetTitle("Gain (mV/fC)");
      h1[ista][igain]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_gain%d.%s", outdir, kSta[ista], kGain[igain], suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gain%d.%s", outdir, kSta[ista], kGain[igain], suffix2); 
      c1->SaveAs(name);

      h1Corr[ista]->SetMarkerStyle(20);
      h1Corr[ista]->GetXaxis()->SetTitle("Gain (Corr, mV/fC)");
      h1Corr[ista]->GetYaxis()->SetTitle("Yield");
      h1Corr[ista]->SetTitle("Gain (Corr, mV/fC)");
      h1Corr[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_gain%d_corr.%s", outdir, kSta[ista], kGain[igain], suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gain%d_corr.%s", outdir, kSta[ista], kGain[igain], suffix2); 
      c1->SaveAs(name);

      gStyle->SetOptStat(0);

      c1->SetLogy(0);
      c1->cd();
      hp1[ista][igain]->SetMarkerStyle(20);
      hp1[ista][igain]->GetYaxis()->SetTitle("Gain (mV/fC)");
      hp1[ista][igain]->GetXaxis()->SetTitle("Chan #");
      hp1[ista][igain]->SetTitle("Gain (mV/fC) vs Chan #");
      hp1[ista][igain]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_gain%d_vs_chan_prof.%s", outdir, kSta[ista], kGain[igain], suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gain%d_vs_chan_prof.%s", outdir, kSta[ista], kGain[igain], suffix2); 
      c1->SaveAs(name);

      c1->cd();
      hp1Corr[ista]->SetMarkerStyle(20);
      hp1Corr[ista]->GetYaxis()->SetTitle("Gain (Corr, mV/fC)");
      hp1Corr[ista]->GetXaxis()->SetTitle("Chan #");
      hp1Corr[ista]->SetTitle("Gain (Corr, mV/fC) vs Chan #");
      hp1Corr[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_gain%d_corr_vs_chan_prof.%s", outdir, kSta[ista], kGain[igain], suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gain%d_corr_vs_chan_prof.%s", outdir, kSta[ista], kGain[igain], suffix2); 
      c1->SaveAs(name);
      
      c1->cd();
      hp1CorrChip[ista]->SetMinimum(kMinGain[igain]);
      hp1CorrChip[ista]->SetMaximum(kMaxGain[igain]);
      hp1CorrChip[ista]->SetMarkerStyle(20);
      hp1CorrChip[ista]->GetYaxis()->SetTitle("Gain (Corr, mV/fC)");
      hp1CorrChip[ista]->GetXaxis()->SetTitle("Chip #");
      hp1CorrChip[ista]->SetTitle("Gain (Corr, mV/fC) vs Chip #");
      hp1CorrChip[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_gain%d_corr_vs_chip_prof.%s", outdir, kSta[ista], kGain[igain], suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gain%d_corr_vs_chip_prof.%s", outdir, kSta[ista], kGain[igain], suffix2); 
      c1->SaveAs(name);
      
      c1->cd();
      hp1CorrChipAll[ista]->SetMarkerStyle(20);
      hp1CorrChipAll[ista]->GetYaxis()->SetTitle("Gain (Corr, All, mV/fC)");
      hp1CorrChipAll[ista]->GetXaxis()->SetTitle("Chip #");
      hp1CorrChipAll[ista]->SetTitle("Gain (Corr, All, mV/fC) vs Chip #");
      hp1CorrChipAll[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_gain%d_corr_vs_chip_prof_all.%s", outdir, kSta[ista], kGain[igain], suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gain%d_corr_vs_chip_prof_all.%s", outdir, kSta[ista], kGain[igain], suffix2); 
      c1->SaveAs(name);
      
      c1->cd();
      h2[ista][igain]->GetYaxis()->SetTitle("Gain (mV/fC)");
      h2[ista][igain]->GetXaxis()->SetTitle("Chan #");
      h2[ista][igain]->SetTitle("Gain (mV/fC) vs Chan #");
      h2[ista][igain]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_gain%d_vs_chan.%s", outdir, kSta[ista], kGain[igain], suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gain%d_vs_chan.%s", outdir, kSta[ista], kGain[igain], suffix2); 
      c1->SaveAs(name);
      
      c1->cd();
      h2Corr[ista]->GetYaxis()->SetTitle("Gain (Corr, mV/fC)");
      h2Corr[ista]->GetXaxis()->SetTitle("Chan #");
      h2Corr[ista]->SetTitle("Gain (Corr, mV/fC) vs Chan #");
      h2Corr[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_gain%d_corr_vs_chan.%s", outdir, kSta[ista], kGain[igain], suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gain%d_corr_vs_chan.%s", outdir, kSta[ista], kGain[igain], suffix2); 
      c1->SaveAs(name);
      
      c1->cd();
      h2CorrChip[ista]->GetYaxis()->SetTitle("Gain (Corr, mV/fC)");
      h2CorrChip[ista]->GetXaxis()->SetTitle("Chip #");
      h2CorrChip[ista]->SetTitle("Gain (Corr, mV/fC) vs Chip #");
      h2CorrChip[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_gain%d_corr_vs_chip.%s", outdir, kSta[ista], kGain[igain], suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gain%d_corr_vs_chip.%s", outdir, kSta[ista], kGain[igain], suffix2); 
      c1->SaveAs(name);
      
      c1->cd();
      h2CorrChipAll[ista]->GetYaxis()->SetTitle("Gain (Corr, All, mV/fC)");
      h2CorrChipAll[ista]->GetXaxis()->SetTitle("Chip #");
      h2CorrChipAll[ista]->SetTitle("Gain (Corr, All, mV/fC) vs Chip #");
      h2CorrChipAll[ista]->Draw("box");
      c1->Modified();
      sprintf(name, "%s/sta%d_gain%d_corr_vs_chip_all.%s", outdir, kSta[ista], kGain[igain], suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gain%d_corr_vs_chip_all.%s", outdir, kSta[ista], kGain[igain], suffix2); 
      c1->SaveAs(name);
      
      // Finally, the reporting...
      cout << endl << " here is a list of good chips (nOK == 32 && nFail == 0) - no failed fits: " << endl;
      int ngood = 0;
      for (locChip=minChip; locChip<=maxChip; locChip++) {
	if (nOK[ista][locChip]==32 && nFail[ista][locChip]==0) {
	  iChip = locChip + kChipOffset;  
	  printf("iSta %d iChip %d meanGain %4.3f rmsGain %4.3f\n", 
		 kSta[ista], iChip, meanGain[ista][locChip], rmsGain[ista][locChip]);
	  ngood++;
	}
      }
      cout << " - number of chips with all good fits: " << ngood << endl;
      
      cout << endl << " here is a list of grey zone chips (nOK != 32 && nFail == 0) - no failed fits but wrong number of channels: " << endl; 
      int ngrey = 0;
      for (locChip=minChip; locChip<=maxChip; locChip++) {
	if (nOK[ista][locChip]!=32 && nOK[ista][locChip]>0 && nFail[ista][locChip]==0) {
	  iChip = locChip + kChipOffset;  
	  printf("iSta %d iChip %d nOK %02d nFail %02d meanGain %4.3f rmsGain %4.3f\n", 
		 kSta[ista], iChip, nOK[ista][locChip], nFail[ista][locChip], meanGain[ista][locChip], rmsGain[ista][locChip]);
	  ngrey++;
	}
      }
      cout << " - number of chips with all good fits but not right number of channels: " << ngrey << endl;
      
      cout << endl << " next is a list of somewhat problematic chips, with 1-3 failed fits: " << endl;
      int nproblem = 0;
      for (locChip=minChip; locChip<=maxChip; locChip++) {
	if (nFail[ista][locChip]>0 && nFail[ista][locChip]<4) {
	  iChip = locChip + kChipOffset;  
	  printf("iSta %d iChip %d nOK %02d nFail %02d meanGain %4.3f rmsGain %4.3f\n", 
		 kSta[ista], iChip, nOK[ista][locChip], nFail[ista][locChip], meanGain[ista][locChip], rmsGain[ista][locChip]);
	  nproblem++;
	}
      }
      cout << " - number of somewhat problematic chips : " << nproblem << endl;
      
      
      cout << endl << " finally a list of more problematic chips, with at least 4 failed fits: " << endl;
      nproblem = 0;
      for (locChip=minChip; locChip<=maxChip; locChip++) {
	if (nFail[ista][locChip]>3) {
	  iChip = locChip + kChipOffset;  
	  printf("iSta %d iChip %d nOK %02d nFail %02d meanGain %4.3f rmsGain %4.3f\n", 
		 kSta[ista], iChip, nOK[ista][locChip], nFail[ista][locChip], meanGain[ista][locChip], rmsGain[ista][locChip]);
	  nproblem++;
	}
      }
      cout << " - number of more problematic chips : " << nproblem << endl;
      
    } // some entries
  } // ista

}