#include "TSystem.h"
#include "TFile.h"
#include "TTree.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TMath.h"
#include "TH2.h"
#include <iostream>
#include <iomanip>

using namespace std;

void testScript() {

  TString infile = "/u/home/jbouffar/testfile.root";
  TFile *f = TFile::Open(infile);
  TH2F* corr = (TH2F*)f->Get("testHisto");
  gStyle->SetOptStat(0);

  cout << "   Correlation factor returns 0.08.  For other 1-bin histograms, this may return >1.0:" << endl;
  cout << "CF = " << corr->GetCorrelationFactor() << endl;

  cout << "   GetStdDev returns nonzero.  This also fluctuates with different 1 bin histograms:" << endl;
  cout << "GetStdDev(1) = " << corr->GetStdDev(1) << endl;
  cout << "GetStdDev(2) = " << corr->GetStdDev(2) << "           <--" << endl;


  cout << "   Use GetStats() function to re-calculate GetStdDev:" << endl;
  Double_t stats1[7];
  corr->GetStats(stats1);
  cout << "stats1[0] (sumw)  = " << stats1[0] << endl;
  cout << "stats1[1] (sumw2) = " << stats1[1] << endl;
  cout << "stats1[2] (sumwx) = " << stats1[2] << endl;
  cout << "stats1[3] (sumwx2)= " << stats1[3] << endl;
  cout << "stats1[4] (sumwy)  = " << stats1[4] << endl;
  cout << "stats1[5] (sumwy2)  = " << stats1[5] << endl;
  cout << "stats1[6] (sumwxy)  = " << stats1[6] << endl;

	
  Double_t mean1 = corr->GetMean(1);
  Double_t mean2 = corr->GetMean(2);
  cout << "mean_x = " << mean1 << endl;
  cout << "mean_y = " << mean2 << endl;

  cout << "   stats[4] behaves as it should, when divided by sumw, it returns mean y:" << endl;
  cout << "stats[4]/stats[0] = " << stats1[4]/stats1[0] << endl;

  //cout << "stats[3]/stats[0] - x*x = " << stats1[3]/stats1[0] - stats1[2]*stats1[2]/(stats1[0]*stats1[0]) << endl;
  cout << "   This should return zero, but doesn't!: " << endl;
  cout << "stats[5]/stats[0] - mean_y*mean_y = " << stats1[5]/stats1[0] - mean2*mean2 << endl;
  cout << "   Successfully reproduced GetStdDev(), using GetStats() and GetMean():" << endl;
  cout << "sqrt(abs(stats[5]/stats[0] - mean_y*mean_y)) = " << TMath::Sqrt(TMath::Abs(stats1[5]/stats1[0] - mean2*mean2)) << " <--" << endl;


  cout << "   Why doesn't stats[5] work? Recalculate with bin loop (overkill for single bin!)" << endl;
  //recalculate stddev
  Int_t firstBinX = corr->GetXaxis()->GetFirst();
  Int_t lastBinX = corr->GetXaxis()->GetLast();
  Int_t firstBinY = corr->GetYaxis()->GetFirst();
  Int_t lastBinY = corr->GetYaxis()->GetLast();
  Double_t myx = 0;
  Double_t myy = 0;

  cout << "   I realize that GetStats uses GetBinCenter, not GetMean, and they differ." << endl;
  Double_t mystats[7] = {0, 0, 0, 0, 0, 0, 0};
  for(Int_t biny = firstBinY; biny <= lastBinY; biny++){
    Double_t y = corr->GetYaxis()->GetBinCenter(biny);
    for(Int_t binx = firstBinX; binx < lastBinX; binx++){
      Double_t x = corr->GetXaxis()->GetBinCenter(binx);
      Int_t bin = corr->GetBin(binx, biny);
      Double_t w = corr->GetBinContent(bin);
      Double_t wx = w * x;
      Double_t wy = w * y;

      if(w != 0){
	cout << "binx = " << binx << endl;
	cout << "biny = " << biny << endl;
	cout << "x = " << setprecision(20) << x << endl;
	cout << "y = " << y << endl;
	cout << "w = " << w << endl;
	cout << "MeanX = " << corr->GetMean(1) << endl;
	cout << "MeanY = " << corr->GetMean(2) << endl;
	myy = y;
	myx = x;
      }

      mystats[0] += w;
      mystats[1] += w * w;
      mystats[2] += wx;
      mystats[3] += wx * x;
      mystats[4] += wy;
      mystats[5] += wy * y;
      mystats[6] += wx * y;
    }
  }

  cout << "   Value of mystats[5] differs from GetStats() value: " << endl;
  cout << "mystats[5] = " << mystats[5] << endl;
 
  cout << "   My version of GetStdDev returns zero: ..." << endl;
  cout << "sqrt(abs(stats[5]/stats[0] - myy*myy)) = " << TMath::Sqrt(TMath::Abs(mystats[5]/mystats[0] - myy*myy)) << endl;
  cout << "   Calculated from GetStats() and the mean y value, this matches GetStdDev()" << endl;
  cout << "sqrt(abs(stats[5]/stats[0] - meany*meany)) = " << TMath::Sqrt(TMath::Abs(stats1[5]/stats1[0] - mean2*mean2)) << endl;
	
}