{
//This program produces histograms displaying the IBD xsec and IBD energy+baseline in KamLAND. Also, an example of simulating 500 IBD events is given. 
//Program works with root version 5 or version 6
//run with the following command
//root ibd.C

double RATE_NORMALIZATION=29756.78;//nuebar-electron elastic events with the following assumptions
//total flux:    1.147e+23 
//nyrs(with DF):  4.00 
//protons/yr:     1.97e+24
//isotopes/proton:   0.01456 
//isotope_br:  1.00
//beta-decay endpoint =   15.000 MeV
//mass:    897.0 tons 
//radius:   6.50 m
//shield:    3.0 m
//distance:  16.00 m

//nuebar-e xsec. Note that changes here are not automatically implemented in the rate calculation.
double sinsqthw = 0.238;
double eminus = -0.5-sinsqthw;
double eplus = -sinsqthw;
//sig0 = 2Gmu^2*me*enu/(pi) in cm^2
double sig0 = 1.72328e-44;//*enu with enu in MeV
//xsec_nuebar_e = sig0*enu*(eplus**2+eminus**2/3.0)*8.0 // 8 electrons per CH2
TF1 *xsec_elastic = new TF1("fa","[0]*x*8.*(pow([1],2)+(pow([2],2)/3.0))",0.,16.);
xsec_elastic->SetParameter(0,sig0);
xsec_elastic->SetParameter(1,eplus);
xsec_elastic->SetParameter(2,eminus);
TH1D *hist_elasticxsec=new TH1D("elasticxsec","elasticxsec",1000,0,16);

//input file and histograms
TFile *file=new TFile("elastic_inputs.root");
TH2D *energies= (TH2D*)file->Get("Graph2D");

//normalize input histograms
energies->Scale(RATE_NORMALIZATION/energies->Integral());

//initialize histograms to put simulated events in
TH2D *energy_sim=new TH2D("","",100,0,16,100,0,16);//100 bins from 0-16 MeV
TH1D *neutrino_energy_sim=new TH1D("","",100,0,16);//100 bins from 0-16 MeV
TH1D *electron_energy_sim=new TH1D("","",45,3,12);//100 bins from 0-16 MeV

//Example of simulating 5000 nuebar-electron elastic events and plotting neutrino energy and electron energy
int TRIALS=50000000;
double neutrino_energy,electron_energy=0.;
for(int i=0;i<TRIALS;i++)
{
energies->GetRandom2(neutrino_energy,electron_energy);
energy_sim->Fill(neutrino_energy,electron_energy);
neutrino_energy_sim->Fill(neutrino_energy);
electron_energy_sim->Fill(electron_energy);
}

//Produce a bunch of plots
new TCanvas;
energies->SetTitle("nuebar-electron elastic input");
energies->SetXTitle("neutrino energy (MeV)");
energies->SetYTitle("electron energy (MeV)");
energies->Draw("COLZ");

new TCanvas;
energy_sim->SetTitle("Simulated nuebar-electron energies");
energy_sim->SetXTitle("neutrino energy (MeV)");
energy_sim->SetYTitle("electron energy (MeV)");
energy_sim->Draw("COLZ");

new TCanvas;
neutrino_energy_sim->SetTitle("Simulated nuebar-electron neutrino energy");
neutrino_energy_sim->SetXTitle("neutrino energy (MeV)");
neutrino_energy_sim->SetYTitle("events/bin");
neutrino_energy_sim->Draw();

new TCanvas;
electron_energy_sim->SetTitle("Simulated nuebar-electron electron energy");
electron_energy_sim->SetXTitle("electron energy (MeV)");
electron_energy_sim->SetYTitle("events/bin");
electron_energy_sim->Draw();


//put nuebar-elastic xsec into easier-to-manipulate histogram
for(int i=1;i<=hist_elasticxsec->GetNbinsX();i++)
{
if(xsec_elastic->Eval(hist_elasticxsec->GetBinCenter(i))>0.)
hist_elasticxsec->SetBinContent(hist_elasticxsec->GetXaxis()->FindBin(hist_elasticxsec->GetBinCenter(i)),xsec_elastic->Eval(hist_elasticxsec->GetBinCenter(i)));
}

new TCanvas;
hist_elasticxsec->SetTitle("nuebar-electron elastic xsec");
hist_elasticxsec->SetXTitle("neutrino energy (MeV)");
hist_elasticxsec->SetYTitle("#sigma (cm^{2})");
hist_elasticxsec->Draw();


TFile *out=new TFile("mark_input_2.root", "RECREATE"); 
electron_energy_sim->Write();
out->Close();

}