#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>
#include <stdio.h>

void gainCompare(const int date1=20180822, const int date2=20180821)
{
  char file1[100];
  sprintf(file1,"log/%d/result.root", date1);
  char file2[100];
  sprintf(file2, "log/%d/result.root", date2);
  const int kNFiles = 2;

  const int debug = 1;

  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, "compare/%d-%d", date1, date2);
  char command[100];
  sprintf(command, ".! mkdir -p %s", outdir);
  gROOT->ProcessLine(command);
  cout << " executing command " << command << endl;

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

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

  // file1 info
  TFile *f1 = new TFile(file1);
  TTree *tree1 = (TTree*) f1->Get("tree");
  //Declaration of leaves types
  Int_t           iChip1;
  Int_t           iSta1;
  Float_t         fGain20mV1[32];
  Int_t           iSync1[kNSync];

  // Set branch addresses.
  tree1->SetBranchAddress("iChip",&iChip1);
  tree1->SetBranchAddress("iSta",&iSta1);
  tree1->SetBranchAddress("fGain20mV",fGain20mV1);
  tree1->SetBranchAddress("iSync",iSync1);

  int nentries1 = tree1->GetEntries();

  // file2 info
  TFile *f2 = new TFile(file2);
  TTree *tree2 = (TTree*) f2->Get("tree");
  //Declaration of leaves types
  Int_t           iChip2;
  Int_t           iSta2;
  Float_t         fGain20mV2[32];
  Int_t           iSync2[kNSync];

  // Set branch addresses.
  tree2->SetBranchAddress("iChip",&iChip2);
  tree2->SetBranchAddress("iSta",&iSta2);
  tree2->SetBranchAddress("fGain20mV",fGain20mV2);
  tree2->SetBranchAddress("iSync",iSync2);

  int nentries2 = tree2->GetEntries();

  cout << " nentries1 " << nentries1 << endl;
  cout << " nentries2 " << nentries2 << endl;

  int ifile = 0;
  int ista = 0;
  Long64_t nbytes = 0;

  int minChip = kNChips;
  int maxChip = 0;
  int iEntry1[kNChips] = {0};
  int iEntry2[kNChips] = {0};
  int locChip = 0;  
  int ichip = 0;
  int ichan = 0;
  int k = 0;

  ifile = 0;
  for (int i=0; i<nentries1;i++) {
    nbytes += tree1->GetEntry(i);

    int nSyncFail = 0;
    int nSyncOK = 0;
    for (k=0; k<kNSync; k++) {
      if ( (200<iSync1[k]) && (iSync1[k]<300) ) { nSyncOK++; }
      else { nSyncFail++; }
    }
    if (debug>1) cout << " nSyncOK " << nSyncOK << endl;
    if ( nSyncOK == kNSync ) {
      
      locChip = iChip1 - kChipOffset;  
      if ( (locChip>=0) && (locChip < kNChips) ) { // right range of chips
	
	if (minChip>locChip) minChip = locChip;
	if (maxChip<locChip) maxChip = locChip;    
	
	iEntry1[locChip] = i; // keep track of latest OK entry for each chip	
	if (debug>1) cout << " ok entry1 " << i << " chip1 " << iChip1 << endl;
      }
    }
  }

  // we now know our range of chip numbers, for the first file
  int nChips = maxChip - minChip + 1;
  int chipcount = nentries1; // can't have more chip matches than that

  // 2nd files
  ifile = 1;
  for (int i=0; i<nentries2;i++) {
    nbytes += tree2->GetEntry(i);

    int nSyncFail = 0;
    int nSyncOK = 0;
    for (k=0; k<kNSync; k++) {
      if ( (200<iSync2[k]) && (iSync2[k]<300) ) { nSyncOK++; }
      else { nSyncFail++; }
    }
    if (debug>1) cout << " nSyncOK " << nSyncOK << endl;
    if ( nSyncOK == kNSync ) {
      
      locChip = iChip2 - kChipOffset;  
      if ( (locChip>=0) && (locChip < kNChips) ) { // right range of chips	
	iEntry2[locChip] = i; // keep track of latest entry for each chip
	if (debug>1) cout << " ok entry2 " << i << " chip2 " << iChip2 << endl;
      }
    }
  }


  ofstream fout;
  char fname[100];
  sprintf(fname, "%s/aveGain.txt", outdir);
  fout.open(fname);
  char fstring[100];

  // let's separate the difference plots for when the tests were in sta 0,1 or one in each (case 2)
  TH1F *h1Diff[kNSta+1];
  TH1F *h1DiffEven[kNSta+1];
  TH1F *h1DiffOdd[kNSta+1];
  TH2F *h2Diff[kNSta+1];
  TH2F *h2DiffChip[kNSta+1];

  const float kMinDiff = -2;
  const float kMaxDiff = 2;

  int igain = 0;
  for (ista=0; ista<=kNSta; ista++) {
    h1Diff[ista] = new TH1F(Form("h1Diff%d", ista),
			    Form("GainDiff %dmV : Sta%d", kGain[igain], kSta[ista]),
			    100, kMinDiff, kMaxDiff );
    h1DiffEven[ista] = new TH1F(Form("h1DiffEven%d", ista),
				Form("GainDiffEven %dmV : Sta%d", kGain[igain], kSta[ista]),
				100, kMinDiff, kMaxDiff );
    h1DiffOdd[ista] = new TH1F(Form("h1DiffOdd%d", ista),
			       Form("GainDiffOdd %dmV : Sta%d", kGain[igain], kSta[ista]),
			       100, kMinDiff, kMaxDiff );
    h2Diff[ista] = new TH2F( Form("h2Diff%d", ista),
			     Form("Gain (Diff) %dmV : Sta%d vs Chan #", kGain[igain], kSta[ista]),
			     32, -0.5, 31.5,
			     100, kMinDiff, kMaxDiff );    
    
    h2DiffChip[ista] = new TH2F(Form("h2DiffChip%d", ista),
				Form("Gain (Diff) %dmV : Sta%d vs Chip #", kGain[igain], kSta[ista]),			     
				chipcount+1, -0.5, chipcount+0.5, 				100, kMinDiff, kMaxDiff );
  }
  
  int nOK[kNSta+1][kNChips] = {0};
  float meanGain[kNSta+1][kNChips] = {0};
  float rmsGain[kNSta+1][kNChips] = {0};
  int nFail[kNSta+1][kNChips] = {0};
  
  chipcount = 0;
  // OK, here we do the Diff loop  						 
  for (ichip = minChip; ichip<=maxChip; ichip++) {
    if ( (iEntry1[ichip] > 0) && (iEntry2[ichip] > 0) ) {
      nbytes += tree1->GetEntry(iEntry1[ichip]);
      nbytes += tree2->GetEntry(iEntry2[ichip]);
      locChip = iChip1 - kChipOffset;  
      if (iSta1 == iSta2) ista = iSta1;
      else ista = 0; // not matching station      
      
      int nSum = 0;      
      float sum = 0;
      float sum2 = 0;
      
      for (ichan = 0; ichan<32; ichan++) {
	float dG = fGain20mV1[ichan] - fGain20mV2[ichan];
	//	cout << " line " << __LINE__ << ": ichan " << ichan << " gain " << fGain20mV[ichan] << " avecorr " << avecorr[ista][ichan] << endl;
	
	h1Diff[ista]->Fill(dG);
	if ( (ichan%2) == 0 ) h1DiffEven[ista]->Fill(dG); 
	if ( (ichan%2) == 1 ) h1DiffOdd[ista]->Fill(dG); 
	h2Diff[ista]->Fill(ichan, dG);

	sum += dG;
	sum2 += dG * dG;
	nSum ++;
	
	if ( (kMinDiff<dG) && (dG<kMaxDiff) ) {
	  h2DiffChip[ista]->Fill(chipcount, 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]));
	sprintf(fstring, "%d %02d %02d %5.2f %5.2f", iChip1, nOK[ista][locChip], nFail[ista][locChip], 
		meanGain[ista][locChip], meanGain[ista][locChip] );
	fout << fstring << endl;
      }
    } // entries for this chip
    chipcount++;
  }
  fout.close();


  cout << " histos filled " << endl;
  TCanvas* c1 = new TCanvas("c1", "c1", 900, 600);

  for (ista=0; ista<=kNSta; ista++) {	
    if (h1Diff[ista]->GetEntries() > 0) {

      c1->SetLogy(1);
      c1->cd();
      
      gStyle->SetOptStat(1);

      h1Diff[ista]->SetMarkerStyle(20);
      h1Diff[ista]->GetXaxis()->SetTitle("Gain Diff (mV/fC)");
      h1Diff[ista]->GetYaxis()->SetTitle("Yield");
      h1Diff[ista]->SetTitle( Form ("Sta %d Gain Diff (mV/fC)", ista) );
      h1Diff[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_diff.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_diff.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);

      h1DiffEven[ista]->SetMarkerStyle(20);
      h1DiffEven[ista]->GetXaxis()->SetTitle("Gain-Diff-Even (mV/fC)");
      h1DiffEven[ista]->GetYaxis()->SetTitle("Yield");
      h1DiffEven[ista]->SetTitle( Form ("Sta %d Gain-Diff-Even (mV/fC)", ista) );
      h1DiffEven[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_diff_even.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_diff_even.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);

      h1DiffOdd[ista]->SetMarkerStyle(20);
      h1DiffOdd[ista]->GetXaxis()->SetTitle("Gain-Diff-Odd (mV/fC)");
      h1DiffOdd[ista]->GetYaxis()->SetTitle("Yield");
      h1DiffOdd[ista]->SetTitle(Form ("Sta %d Gain-Diff-Even (mV/fC)", ista) );
      h1DiffOdd[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_diff_odd.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_diff_odd.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);
      
      gStyle->SetOptStat(0);
      
      c1->SetLogy(0);
      c1->cd();
      
      c1->cd();
      h2Diff[ista]->GetYaxis()->SetTitle("Gain Diff (mV/fC)");
      h2Diff[ista]->GetXaxis()->SetTitle("Chan #");
      h2Diff[ista]->SetTitle( Form("Sta %d Gain Diff (mV/fC) vs Chan #", ista) );
      h2Diff[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_diff_vs_chan.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_diff_vs_chan.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);
      
      c1->cd();
      h2DiffChip[ista]->GetYaxis()->SetTitle("Gain Diff (mV/fC)");
      h2DiffChip[ista]->GetXaxis()->SetTitle("Chip #");
      h2DiffChip[ista]->SetTitle( Form("Sta %d Gain Diff (mV/fC) vs Chip #", ista) );
      h2DiffChip[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_diff_vs_chip.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_diff_vs_chip.%s", outdir, ista, 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) {
	  iChip1 = locChip + kChipOffset;  
	  printf("iSta %d iChip %d meanGain %4.3f rmsGain %4.3f\n", 
		 kSta[ista], iChip1, 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) {
	  iChip1 = locChip + kChipOffset;  
	  printf("iSta %d iChip %d nOK %02d nFail %02d meanGain %4.3f rmsGain %4.3f\n", 
		 kSta[ista], iChip1, 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) {
	  iChip1 = locChip + kChipOffset;  
	  printf("iSta %d iChip %d nOK %02d nFail %02d meanGain %4.3f rmsGain %4.3f\n", 
		 kSta[ista], iChip1, 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) {
	  iChip1 = locChip + kChipOffset;  
	  printf("iSta %d iChip %d nOK %02d nFail %02d meanGain %4.3f rmsGain %4.3f\n", 
		 kSta[ista], iChip1, nOK[ista][locChip], nFail[ista][locChip], meanGain[ista][locChip], rmsGain[ista][locChip]);
	  nproblem++;
	}
      }
      cout << " - number of more problematic chips : " << nproblem << endl;
      
    } // some entries
  } // ista

}