#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 compare(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         fRmsAna20mV1[32];
  Float_t         fPedAna20mV1[32];
  Float_t         fGain20mV1[32];
  Float_t         fGainRT20mV1[32];
  Float_t         fRT20mV1[32];
  Int_t           iSync1[kNSync];

  // Set branch addresses.
  tree1->SetBranchAddress("iChip",&iChip1);
  tree1->SetBranchAddress("iSta",&iSta1);
  tree1->SetBranchAddress("fPedAna20mV",fPedAna20mV1);
  tree1->SetBranchAddress("fRmsAna20mV",fRmsAna20mV1);
  tree1->SetBranchAddress("fGain20mV",fGain20mV1);
  tree1->SetBranchAddress("fGainRT20mV",fGainRT20mV1);
  tree1->SetBranchAddress("fRT20mV",fRT20mV1);
  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         fRmsAna20mV2[32];
  Float_t         fPedAna20mV2[32];
  Float_t         fGain20mV2[32];
  Float_t         fGainRT20mV2[32];
  Float_t         fRT20mV2[32];
  Int_t           iSync2[kNSync];

  // Set branch addresses.
  tree2->SetBranchAddress("iChip",&iChip2);
  tree2->SetBranchAddress("iSta",&iSta2);
  tree2->SetBranchAddress("fPedAna20mV",fPedAna20mV2);
  tree2->SetBranchAddress("fRmsAna20mV",fRmsAna20mV2);
  tree2->SetBranchAddress("fGain20mV",fGain20mV2);
  tree2->SetBranchAddress("fGainRT20mV",fGainRT20mV2);
  tree2->SetBranchAddress("fRT20mV",fRT20mV2);
  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 1,2 or one in each (case 0)
  TH1F *h1DiffPedAna[kNSta+1];
  TH1F *h1DiffPedAnaEven[kNSta+1];
  TH1F *h1DiffPedAnaOdd[kNSta+1];
  TH2F *h2DiffPedAna[kNSta+1];
  TH2F *h2DiffPedAnaChip[kNSta+1];
  const float kMinDiffPedAna = -5;
  const float kMaxDiffPedAna = 5;

  TH1F *h1DiffRmsAna[kNSta+1];
  TH1F *h1DiffRmsAnaEven[kNSta+1];
  TH1F *h1DiffRmsAnaOdd[kNSta+1];
  TH2F *h2DiffRmsAna[kNSta+1];
  TH2F *h2DiffRmsAnaChip[kNSta+1];
  const float kMinDiffRmsAna = -2;
  const float kMaxDiffRmsAna = 2;

  TH1F *h1DiffGain[kNSta+1];
  TH1F *h1DiffGainEven[kNSta+1];
  TH1F *h1DiffGainOdd[kNSta+1];
  TH2F *h2DiffGain[kNSta+1];
  TH2F *h2DiffGainChip[kNSta+1];
  const float kMinDiffGain = -2;
  const float kMaxDiffGain = 2;

  TH1F *h1DiffGainRT[kNSta+1];
  TH1F *h1DiffGainRTEven[kNSta+1];
  TH1F *h1DiffGainRTOdd[kNSta+1];
  TH2F *h2DiffGainRT[kNSta+1];
  TH2F *h2DiffGainRTChip[kNSta+1];
  const float kMinDiffGainRT = -2;
  const float kMaxDiffGainRT = 2;

  TH1F *h1DiffRT[kNSta+1];
  TH1F *h1DiffRTEven[kNSta+1];
  TH1F *h1DiffRTOdd[kNSta+1];
  TH2F *h2DiffRT[kNSta+1];
  TH2F *h2DiffRTChip[kNSta+1];
  const float kMinDiffRT = -5;
  const float kMaxDiffRT = 5;

  int igain = 0;
  char title[100];
  for (ista=0; ista<=kNSta; ista++) {
    sprintf(title,"DiffPedAna : Sta %d", ista);
    h1DiffPedAna[ista] = new TH1F(Form("h1DiffPedAna%d", ista), title,
				  100, kMinDiffPedAna, kMaxDiffPedAna );
    sprintf(title,"DiffPedAnaEven : Sta %d", ista);
    h1DiffPedAnaEven[ista] = new TH1F(Form("h1DiffPedAnaEven%d", ista), title,
				      100, kMinDiffPedAna, kMaxDiffPedAna );
    sprintf(title,"DiffPedAnaOdd : Sta %d", ista);
    h1DiffPedAnaOdd[ista] = new TH1F(Form("h1DiffPedAnaOdd%d", ista), title,
				     100, kMinDiffPedAna, kMaxDiffPedAna );
    sprintf(title,"DiffPedAna vs Chan : Sta %d", ista);
    h2DiffPedAna[ista] = new TH2F( Form("h2DiffPedAna%d", ista), title,
				   32, -0.5, 31.5,
				   100, kMinDiffPedAna, kMaxDiffPedAna );        
    sprintf(title,"DiffPedAna vs Chip : Sta %d", ista);
    h2DiffPedAnaChip[ista] = new TH2F(Form("h2DiffPedAnaChip%d", ista), title,
				      chipcount+1, -0.5, chipcount+0.5, 				
				      100, kMinDiffPedAna, kMaxDiffPedAna );

    sprintf(title,"DiffRmsAna : Sta %d", ista);
    h1DiffRmsAna[ista] = new TH1F(Form("h1DiffRmsAna%d", ista), title,
				  100, kMinDiffRmsAna, kMaxDiffRmsAna );
    sprintf(title,"DiffRmsAnaEven : Sta %d", ista);
    h1DiffRmsAnaEven[ista] = new TH1F(Form("h1DiffRmsAnaEven%d", ista), title,
				      100, kMinDiffRmsAna, kMaxDiffRmsAna );
    sprintf(title,"DiffRmsAnaOdd : Sta %d", ista);
    h1DiffRmsAnaOdd[ista] = new TH1F(Form("h1DiffRmsAnaOdd%d", ista), title,
				     100, kMinDiffRmsAna, kMaxDiffRmsAna );
    sprintf(title,"DiffRmsAna vs Chan : Sta %d", ista);
    h2DiffRmsAna[ista] = new TH2F( Form("h2DiffRmsAna%d", ista), title,
				   32, -0.5, 31.5,
				   100, kMinDiffRmsAna, kMaxDiffRmsAna );        
    sprintf(title,"DiffRmsAna vs Chip : Sta %d", ista);
    h2DiffRmsAnaChip[ista] = new TH2F(Form("h2DiffRmsAnaChip%d", ista), title,
				      chipcount+1, -0.5, chipcount+0.5, 				
				      100, kMinDiffRmsAna, kMaxDiffRmsAna );
    
    sprintf(title,"DiffGain : Sta %d", ista);
    h1DiffGain[ista] = new TH1F(Form("h1DiffGain%d", ista), title,
				100, kMinDiffGain, kMaxDiffGain );
    sprintf(title,"DiffGainEven : Sta %d", ista);
    h1DiffGainEven[ista] = new TH1F(Form("h1DiffGainEven%d", ista), title,
				    100, kMinDiffGain, kMaxDiffGain );
    sprintf(title,"DiffGainOdd : Sta %d", ista);
    h1DiffGainOdd[ista] = new TH1F(Form("h1DiffGainOdd%d", ista), title,
				   100, kMinDiffGain, kMaxDiffGain );
    sprintf(title,"DiffGain vs Chan : Sta %d", ista);
    h2DiffGain[ista] = new TH2F( Form("h2DiffGain%d", ista), title,
				 32, -0.5, 31.5,
				 100, kMinDiffGain, kMaxDiffGain );    
    sprintf(title,"DiffGain vs Chip : Sta %d", ista);
    h2DiffGainChip[ista] = new TH2F(Form("h2DiffGainChip%d", ista), title,
				    chipcount+1, -0.5, chipcount+0.5, 				
				    100, kMinDiffGain, kMaxDiffGain );

    sprintf(title,"DiffGainRT : Sta %d", ista);
    h1DiffGainRT[ista] = new TH1F(Form("h1DiffGainRT%d", ista), title,
				  100, kMinDiffGainRT, kMaxDiffGainRT );
    sprintf(title,"DiffGainRTEven : Sta %d", ista);
    h1DiffGainRTEven[ista] = new TH1F(Form("h1DiffGainRTEven%d", ista), title,
				      100, kMinDiffGainRT, kMaxDiffGainRT );
    sprintf(title,"DiffGainRTOdd : Sta %d", ista);
    h1DiffGainRTOdd[ista] = new TH1F(Form("h1DiffGainRTOdd%d", ista), title,
				     100, kMinDiffGainRT, kMaxDiffGainRT );
    sprintf(title,"DiffGainRT vs Chan : Sta %d", ista);
    h2DiffGainRT[ista] = new TH2F( Form("h2DiffGainRT%d", ista), title,
				   32, -0.5, 31.5,
				   100, kMinDiffGainRT, kMaxDiffGainRT );    
    sprintf(title,"DiffGainRT vs Chip : Sta %d", ista);
    h2DiffGainRTChip[ista] = new TH2F(Form("h2DiffGainRTChip%d", ista), title,
				      chipcount+1, -0.5, chipcount+0.5, 				
				      100, kMinDiffGainRT, kMaxDiffGainRT );
    
    sprintf(title,"DiffRT : Sta %d", ista);
    h1DiffRT[ista] = new TH1F(Form("h1DiffRT%d", ista), title,
			      100, kMinDiffRT, kMaxDiffRT );
    sprintf(title,"DiffRTEven : Sta %d", ista);
    h1DiffRTEven[ista] = new TH1F(Form("h1DiffRTEven%d", ista), title,
				  100, kMinDiffRT, kMaxDiffRT );
    sprintf(title,"DiffRTOdd : Sta %d", ista);
    h1DiffRTOdd[ista] = new TH1F(Form("h1DiffRTOdd%d", ista), title,
				 100, kMinDiffRT, kMaxDiffRT );
    sprintf(title,"DiffRT vs Chan : Sta %d", ista);
    h2DiffRT[ista] = new TH2F( Form("h2DiffRT%d", ista), title,
			       32, -0.5, 31.5,
			       100, kMinDiffRT, kMaxDiffRT );        
    sprintf(title,"DiffRT vs Chip : Sta %d", ista);
    h2DiffRTChip[ista] = new TH2F(Form("h2DiffRTChip%d", ista), title,
				  chipcount+1, -0.5, chipcount+0.5, 				
				  100, kMinDiffRT, kMaxDiffRT );
  }
  
  int nOK[kNSta+1][kNChips] = {0};
  int nFail[kNSta+1][kNChips] = {0};
  float meanGain[kNSta+1][kNChips] = {0};
  float rmsGain[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; // meaning, 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;
	
	h1DiffGain[ista]->Fill(dG);
	if ( (ichan%2) == 0 ) h1DiffGainEven[ista]->Fill(dG); 
	if ( (ichan%2) == 1 ) h1DiffGainOdd[ista]->Fill(dG); 
	h2DiffGain[ista]->Fill(ichan, dG);
	
	sum += dG;
	sum2 += dG * dG;
	nSum ++;
	
	if ( (kMinDiffGain<dG) && (dG<kMaxDiffGain) ) {
	  h2DiffGainChip[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;
      }

      // same for gain from RiseTime measurement
      for (ichan = 0; ichan<32; ichan++) {
	float dG = fGainRT20mV1[ichan] - fGainRT20mV2[ichan];
	
	h1DiffGainRT[ista]->Fill(dG);
	if ( (ichan%2) == 0 ) h1DiffGainRTEven[ista]->Fill(dG); 
	if ( (ichan%2) == 1 ) h1DiffGainRTOdd[ista]->Fill(dG); 
	h2DiffGainRT[ista]->Fill(ichan, dG);
	if ( (kMinDiffGainRT<dG) && (dG<kMaxDiffGainRT) ) {
	  h2DiffGainRTChip[ista]->Fill(chipcount, dG);
	}
      }
      

      // also fill other plots
      for (ichan = 0; ichan<32; ichan++) {
	float dPed = fPedAna20mV1[ichan] - fPedAna20mV2[ichan];
	float dRms = fRmsAna20mV1[ichan] - fRmsAna20mV2[ichan];
	float dRT = fRT20mV1[ichan] - fRT20mV2[ichan];
	
	h1DiffPedAna[ista]->Fill(dPed);
	if ( (ichan%2) == 0 ) h1DiffPedAnaEven[ista]->Fill(dPed); 
	if ( (ichan%2) == 1 ) h1DiffPedAnaOdd[ista]->Fill(dPed); 
	h2DiffPedAna[ista]->Fill(ichan, dPed);
	if ( (kMinDiffPedAna<dPed) && (dPed<kMaxDiffPedAna) ) {
	  h2DiffPedAnaChip[ista]->Fill(chipcount, dPed);
	}

	h1DiffRmsAna[ista]->Fill(dRms);
	if ( (ichan%2) == 0 ) h1DiffRmsAnaEven[ista]->Fill(dRms); 
	if ( (ichan%2) == 1 ) h1DiffRmsAnaOdd[ista]->Fill(dRms); 
	h2DiffRmsAna[ista]->Fill(ichan, dRms);
	if ( (kMinDiffRmsAna<dRms) && (dRms<kMaxDiffRmsAna) ) {
	  h2DiffRmsAnaChip[ista]->Fill(chipcount, dRms);
	}

	h1DiffRT[ista]->Fill(dRT);
	if ( (ichan%2) == 0 ) h1DiffRTEven[ista]->Fill(dRT); 
	if ( (ichan%2) == 1 ) h1DiffRTOdd[ista]->Fill(dRT); 
	h2DiffRT[ista]->Fill(ichan, dRT);
	if ( (kMinDiffRT<dRT) && (dRT<kMaxDiffRT) ) {
	  h2DiffRTChip[ista]->Fill(chipcount, dRT);
	}

      }
      
    } // entries for this chip
    chipcount++;
  }
  fout.close();


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

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

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

      h1DiffPedAna[ista]->SetMarkerStyle(20);
      h1DiffPedAna[ista]->GetXaxis()->SetTitle("PedAna Diff");
      h1DiffPedAna[ista]->GetYaxis()->SetTitle("Yield");
      //      h1DiffPedAna[ista]->SetTitle( Form("PedAna Diff : Sta %d", kSta[ista]) );
      h1DiffPedAna[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_peddiff.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_peddiff.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);

      h1DiffPedAnaEven[ista]->SetMarkerStyle(20);
      h1DiffPedAnaEven[ista]->GetXaxis()->SetTitle("PedAna-Diff-Even");
      h1DiffPedAnaEven[ista]->GetYaxis()->SetTitle("Yield");
      //      h1DiffPedAnaEven[ista]->SetTitle( Form("PedAna Diff Even : Sta %d", kSta[ista]) );
      h1DiffPedAnaEven[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_peddiff_even.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_peddiff_even.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);

      h1DiffPedAnaOdd[ista]->SetMarkerStyle(20);
      h1DiffPedAnaOdd[ista]->GetXaxis()->SetTitle("PedAna-Diff-Odd");
      h1DiffPedAnaOdd[ista]->GetYaxis()->SetTitle("Yield");
      //      h1DiffPedAnaOdd[ista]->SetTitle( Form("PedAna Diff Odd : Sta %d", kSta[ista]) );
      h1DiffPedAnaOdd[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_peddiff_odd.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_peddiff_odd.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);
      
      gStyle->SetOptStat(0);      
      c1->SetLogy(0);
      c1->cd();
      
      c1->cd();
      h2DiffPedAna[ista]->GetYaxis()->SetTitle("PedAna Diff");
      h2DiffPedAna[ista]->GetXaxis()->SetTitle("Chan #");
      //      h2DiffPedAna[ista]->SetTitle( Form("PedAna Diff vs Chan : Sta %d", kSta[ista]) );
      h2DiffPedAna[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_peddiff_vs_chan.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_peddiff_vs_chan.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);
      
      c1->cd();
      h2DiffPedAnaChip[ista]->GetYaxis()->SetTitle("PedAna Diff");
      h2DiffPedAnaChip[ista]->GetXaxis()->SetTitle("Chip #");
      //      h2DiffPedAnaChip[ista]->SetTitle( Form("PedAna Diff vs Chip : Sta %d", kSta[ista]) );
      h2DiffPedAnaChip[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_peddiff_vs_chip.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_peddiff_vs_chip.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);


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

      h1DiffRmsAna[ista]->SetMarkerStyle(20);
      h1DiffRmsAna[ista]->GetXaxis()->SetTitle("RmsAna Diff");
      h1DiffRmsAna[ista]->GetYaxis()->SetTitle("Yield");
      //      h1DiffRmsAna[ista]->SetTitle( Form("RmsAna Diff : Sta %d", kSta[ista]) );
      h1DiffRmsAna[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_rmsdiff.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_rmsdiff.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);

      h1DiffRmsAnaEven[ista]->SetMarkerStyle(20);
      h1DiffRmsAnaEven[ista]->GetXaxis()->SetTitle("RmsAna-Diff-Even");
      h1DiffRmsAnaEven[ista]->GetYaxis()->SetTitle("Yield");
      //      h1DiffRmsAnaEven[ista]->SetTitle( Form("RmsAna Diff Even : Sta %d", kSta[ista]) );
      h1DiffRmsAnaEven[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_rmsdiff_even.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_rmsdiff_even.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);

      h1DiffRmsAnaOdd[ista]->SetMarkerStyle(20);
      h1DiffRmsAnaOdd[ista]->GetXaxis()->SetTitle("RmsAna-Diff-Odd ");
      h1DiffRmsAnaOdd[ista]->GetYaxis()->SetTitle("Yield");
      //      h1DiffRmsAnaOdd[ista]->SetTitle( Form("RmsAna Diff Odd : Sta %d", kSta[ista]) );
      h1DiffRmsAnaOdd[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_rmsdiff_odd.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_rmsdiff_odd.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);
      
      gStyle->SetOptStat(0);      
      c1->SetLogy(0);
      c1->cd();
      
      c1->cd();
      h2DiffRmsAna[ista]->GetYaxis()->SetTitle("RmsAna Diff");
      h2DiffRmsAna[ista]->GetXaxis()->SetTitle("Chan #");
      //      h2DiffRmsAna[ista]->SetTitle( Form("RmsAna Diff vs Chan : Sta %d", kSta[ista]) );
      h2DiffRmsAna[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_rmsdiff_vs_chan.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_rmsdiff_vs_chan.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);
      
      c1->cd();
      h2DiffRmsAnaChip[ista]->GetYaxis()->SetTitle("RmsAna Diff");
      h2DiffRmsAnaChip[ista]->GetXaxis()->SetTitle("Chip #");
      //      h2DiffRmsAnaChip[ista]->SetTitle( Form("RmsAna Diff vs Chip : Sta %d", kSta[ista]) );
      h2DiffRmsAnaChip[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_rmsdiff_vs_chip.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_rmsdiff_vs_chip.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);

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

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

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

      h1DiffGainOdd[ista]->SetMarkerStyle(20);
      h1DiffGainOdd[ista]->GetXaxis()->SetTitle("Gain-Diff-Odd (mV/fC)");
      h1DiffGainOdd[ista]->GetYaxis()->SetTitle("Yield");
      //      h1DiffGainOdd[ista]->SetTitle( Form("Gain Diff Odd : Sta %d", kSta[ista]) );
      h1DiffGainOdd[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_gaindiff_odd.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gaindiff_odd.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);
      
      gStyle->SetOptStat(0);      
      c1->SetLogy(0);
      c1->cd();
      
      c1->cd();
      h2DiffGain[ista]->GetYaxis()->SetTitle("Gain Diff (mV/fC)");
      h2DiffGain[ista]->GetXaxis()->SetTitle("Chan #");
      //      h2DiffGain[ista]->SetTitle( Form("Gain Diff vs Chan : Sta %d", kSta[ista]) );
      h2DiffGain[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_gaindiff_vs_chan.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gaindiff_vs_chan.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);
      
      c1->cd();
      h2DiffGainChip[ista]->GetYaxis()->SetTitle("Gain Diff (mV/fC)");
      h2DiffGainChip[ista]->GetXaxis()->SetTitle("Chip #");
      //      h2DiffGainChip[ista]->SetTitle( Form("Gain Diff vs Chip : Sta %d", kSta[ista]) );
      h2DiffGainChip[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_gaindiff_vs_chip.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gaindiff_vs_chip.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);

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

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

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

      h1DiffGainRTOdd[ista]->SetMarkerStyle(20);
      h1DiffGainRTOdd[ista]->GetXaxis()->SetTitle("GainRT-Diff-Odd (mV/fC)");
      h1DiffGainRTOdd[ista]->GetYaxis()->SetTitle("Yield");
      //      h1DiffGainRTOdd[ista]->SetTitle( Form("GainRT Diff Odd : Sta %d", kSta[ista]) );
      h1DiffGainRTOdd[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_gainrtdiff_odd.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gainrtdiff_odd.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);
      
      gStyle->SetOptStat(0);      
      c1->SetLogy(0);
      c1->cd();
      
      c1->cd();
      h2DiffGainRT[ista]->GetYaxis()->SetTitle("GainRT Diff (mV/fC)");
      h2DiffGainRT[ista]->GetXaxis()->SetTitle("Chan #");
      //      h2DiffGainRT[ista]->SetTitle( Form("GainRT Diff vs Chan : Sta %d", kSta[ista]) );
      h2DiffGainRT[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_gainrtdiff_vs_chan.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gainrtdiff_vs_chan.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);
      
      c1->cd();
      h2DiffGainRTChip[ista]->GetYaxis()->SetTitle("GainRT Diff (mV/fC)");
      h2DiffGainRTChip[ista]->GetXaxis()->SetTitle("Chip #");
      //      h2DiffGainRTChip[ista]->SetTitle( Form("GainRT Diff vs Chip : Sta %d", kSta[ista]) );
      h2DiffGainRTChip[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_gainrtdiff_vs_chip.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_gainrtdiff_vs_chip.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);

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

      h1DiffRT[ista]->SetMarkerStyle(20);
      h1DiffRT[ista]->GetXaxis()->SetTitle("RT Diff");
      h1DiffRT[ista]->GetYaxis()->SetTitle("Yield");
      //      h1DiffRT[ista]->SetTitle( Form("RT Diff : Sta %d", kSta[ista]) );
      h1DiffRT[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_risetimediff.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_risetimediff.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);

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

      h1DiffRTOdd[ista]->SetMarkerStyle(20);
      h1DiffRTOdd[ista]->GetXaxis()->SetTitle("RT-Diff-Odd (mV/fC)");
      h1DiffRTOdd[ista]->GetYaxis()->SetTitle("Yield");
      //      h1DiffRTOdd[ista]->SetTitle( Form("RT Diff Odd : Sta %d", kSta[ista]) );
      h1DiffRTOdd[ista]->Draw();
      c1->Modified();
      sprintf(name, "%s/sta%d_risetimediff_odd.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_risetimediff_odd.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);
      
      gStyle->SetOptStat(0);      
      c1->SetLogy(0);
      c1->cd();
      
      c1->cd();
      h2DiffRT[ista]->GetYaxis()->SetTitle("RT Diff (mV/fC)");
      h2DiffRT[ista]->GetXaxis()->SetTitle("Chan #");
      //      h2DiffRT[ista]->SetTitle( Form("RT Diff vs Chan : Sta %d", kSta[ista]) );
      h2DiffRT[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_risetimediff_vs_chan.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_risetimediff_vs_chan.%s", outdir, ista, suffix2); 
      c1->SaveAs(name);
      
      c1->cd();
      h2DiffRTChip[ista]->GetYaxis()->SetTitle("RT Diff (mV/fC)");
      h2DiffRTChip[ista]->GetXaxis()->SetTitle("Chip #");
      //      h2DiffRTChip[ista]->SetTitle( Form("RT Diff vs Chip : Sta %d", kSta[ista]) );
      h2DiffGainChip[ista]->Draw("colz");
      c1->Modified();
      sprintf(name, "%s/sta%d_risetimediff_vs_chip.%s", outdir, ista, suffix1); 
      c1->SaveAs(name);
      sprintf(name, "%s/sta%d_risetimediff_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
}