//
// Test code to illustrate problem when drawing histograms 
// 
//

#include <iostream>
#include "TTree.h"
#include "TFile.h"
#include <string>
#include <vector>
#include <stdio.h>
#include <stdlib.h>

//using namespace std;

int nHistogramBins = 100;
float upperMassBin = 0.005;
float mottLowBin = 0;

TH1F* genericMottHist = new TH1F("genericMottHist","Generic Mott Distribution",nHistogramBins,mottLowBin,upperMassBin);
TH1F* gsMottHist = new TH1F("gsMottHist","GS Mott Distribution",nHistogramBins,mottLowBin,upperMassBin);

float calcGSDistribution(float m, float mAvg, float caseMass, float j);

void testProblem(){

    gStyle->SetOptStat(0000);
    float caseMass = 271.7/1000; //kg
    float rogersAvgMass = 0.154443175/1000; //kg
	
    for(int j = 1; j < gsMottHist->GetNbinsX()+1; j++){
	float bin1 = gsMottHist->GetBinLowEdge(j);
	float bin2 = gsMottHist->GetBinLowEdge(j+1);
	float avgBinMass = (bin1 + bin2)/2.;

	if(j < 3){
	    cout << "bin # = " << j << ", value = " << calcGSDistribution(avgBinMass, rogersAvgMass, caseMass, 2) << ", lowEdge = " << bin1 << ", highEdge = " << bin2 << endl;
	}
	gsMottHist->SetBinContent(j, calcGSDistribution(avgBinMass, rogersAvgMass, caseMass, 2) );
    }

    TCanvas* c1 = new TCanvas;
    
    //option 1 = INCORRECT
    //option 2 = CORRECT
    int option = 1;

    if(option == 1){
    // OPTION 1
    // Will display leftmost bin (bin index #2, lowEdge = 5e-5, highEdge = 1e-4) as 995 <- INCORRECT
    //////////////////////////////////////

    //gsMottHist->GetXaxis()->SetRangeUser(5e-5,upperMassBin);
    genericMottHist->SetMinimum(1);
    genericMottHist->SetMaximum(10000);
    genericMottHist->GetXaxis()->SetRangeUser(5e-5,upperMassBin);
    
    }else if(option ==  2){
    // OPTION 2
    // Will display leftmost bin (bin index #2, lowEdge = 5e-5, highEdge = 1e-4) as 656 <- CORRECT
    //////////////////////////////////////

    gsMottHist->GetXaxis()->SetRangeUser(5e-5,upperMassBin);
    genericMottHist->SetMinimum(1);
    genericMottHist->SetMaximum(10000);
    genericMottHist->GetXaxis()->SetRangeUser(5e-5,upperMassBin);
    }

    genericMottHist->Draw();
    gsMottHist->Draw("same");

    gPad->SetLogx();
    gPad->SetLogy();

    if(option == 1){
	c1->SaveAs("/mnt/d/Incorrect.png","png");
    }else if(option == 2){
	c1->SaveAs("/mnt/d/Correct.png","png");
    }
}

float calcGSDistribution(float m, float mAvg, float caseMass, float j){
    float mu = mAvg/TMath::Factorial(j);
    float nm = caseMass/mAvg * exp(-pow(m/mu,1/j));
    return nm;
}