#include <TFile.h>
#include <TH2.h>

#include <BCAL/DBCALGeometry.h>
#include <units.h>

void bcalInit();
int getDarkHits();

extern double BCAL_mevPerPE;


int main(int narg, char *argv[])
{
	bcalInit();
	
	// To avoid distracting binning effects when plotting energy, we determine
	// a good bin size in MeV here.
	double binSize = (1.0 * BCAL_mevPerPE * k_MeV)*1000.0;

	TFile *f = new TFile("darkhits.root", "RECREATE");

	TH1D *Ndarkhits = new TH1D("Ndarkhits", "Number of dark hits per event in single SiPM", 100, -0.5, 99.5);
	Ndarkhits->SetXTitle("Number of darkhits (+cross-talk)");

	TH1D *Edarkhits = new TH1D("Edarkhits", "Energy(attenuated) equivalent dark hits per event in single SiPM", 60, 0.0, 60.0*binSize);
	Edarkhits->SetXTitle("Energy (MeV)");
	Edarkhits->SetFillColor(kRed);
	Edarkhits->SetFillStyle(3000);

	TH1D *NdarkhitsSum3 = new TH1D("NdarkhitsSum3", "Number of dark hits per event for sum of 3 SiPMs", 100, -0.5, 99.5);
	NdarkhitsSum3->SetXTitle("Number of darkhits (+cross-talk) for sum of 3 SiPMs");

	TH1D *EdarkhitsSum3 = new TH1D("EdarkhitsSum3", "Energy(attenuated) equivalent dark hits per event for sum of 3 SiPMs", 60, 0.0, 60.0*binSize);
	EdarkhitsSum3->SetXTitle("Energy (MeV)");
	EdarkhitsSum3->SetFillColor(kBlue);
	EdarkhitsSum3->SetFillStyle(3000);

	for(int i=0; i<10000; i++){
		int Nhits = getDarkHits();
		Ndarkhits->Fill((double)Nhits);
		
		double Edark = (double)Nhits * BCAL_mevPerPE * k_MeV;
		Edarkhits->Fill(Edark*1000.0);
		
		Nhits = 0;
		for(int i=0; i<3; i++)Nhits += getDarkHits();
		NdarkhitsSum3->Fill((double)Nhits);
		EdarkhitsSum3->Fill((double)Nhits * BCAL_mevPerPE * k_MeV*1000.0);
	}	

	f->Write();
	delete f;

	return 0;
}