// ROOT macro file for plotting example B4 ntuple
//
// Can be run from ROOT session:
// root[0] .x plotNtuple.C
#include <TMath.h>
#include <TNtuple.h>
#include <TTree.h>
#include <TH2.h>
{
  gROOT->Reset();
  gROOT->SetStyle("Plain");

  // Draw histos filled by Geant4 simulation
  //

  // Open file filled by Geant4 simulation
  TFile* f = new TFile("totalN.root","READ");
  TTree *EverythingN;
  f->GetObject("EverythingN", EverythingN);
//////////////////////////reasonID == -1 /////////////////////
  Double_t w = 1000;
  Double_t h = 1000;
  TCanvas * c1 = new TCanvas("c1", "c1", w, h);
  c1->SetWindowSize(w + (w - c1->GetWw()), h + (h - c1->GetWh()));
//TCanvas* c1 = new TCanvas("c1", "", 20, 80, 800, 800);
Double_t NumberP=2000000.;
c1->Divide(2,2);
if (c1){
c1->cd(1);
TH2F *posXY = new TH2F("posXY", "The distribution of neutrons escaped from the target geometry (XY), reason ID==-1", 401, -200.5, 200.5, 401, -200.5, 200.5);
if(posXY){posXY->GetXaxis()->SetTitle("X (mm)"); posXY->GetYaxis()->SetTitle("Y (mm)"); posXY->GetZaxis()->SetLabelSize(0.025);}
EverythingN->Draw("NPosY:NPosX>>posXY", "reasonID==-1","COLZ1");
posXY->Scale(1/NumberP);
//new TCanvas();
c1->cd(2);
TH2F *posXZ = new TH2F("posXZ", "The distribution of neutrons escaped from the target geometry (XZ), reason ID==-1", 401, -200.5, 200.5, 401, -800.5, -400.5);
if(posXZ){posXZ->GetXaxis()->SetTitle("X (mm)"); posXZ->GetYaxis()->SetTitle("Z (mm)"); posXZ->GetZaxis()->SetLabelSize(0.025); }
EverythingN->Draw("NPosZ:NPosX>>posXZ", "reasonID==-1","COLZ1");
posXZ->Scale(1/NumberP);
//new TCanvas();
c1->cd(3);
TH2F *posYZ = new TH2F("posYZ", "The distribution of neutrons escaped from the target geometry (YZ), reason ID==-1", 401, -200.5, 200.5, 401, -800.5, -400.5);
if(posYZ){posYZ->GetXaxis()->SetTitle("Y (mm)"); posYZ->GetYaxis()->SetTitle("Z (mm)"); posYZ->GetZaxis()->SetLabelSize(0.025); }
EverythingN->Draw("NPosZ:NPosY>>posYZ", "reasonID==-1","COLZ1");
posYZ->Scale(1/NumberP);
c1->cd(4);
//new TCanvas();
TH2F *MonXY = new TH2F("MomentumXY", "The momentum distribution of neutrons escaped from the target geometry (XY)", 1401, -700.5, 700.5, 1401, -700.5, 700.5);
if(MonXY){MonXY->GetXaxis()->SetTitle("Px"); MonXY->GetYaxis()->SetTitle("Py"); MonXY->SetTitle("MomentumXY for reasonID==-1");}
EverythingN->Draw("mY:mX>>MonXY", "reasonID==-1","");

/*c1->cd(5);
TH1F *energy = new TH1F("Energy_full", "Neutron energy spectrum for reason ID==-1",30,0.,log10(300.));
EverythingN->Draw("log10((energyK/1000.))>>Energy_full", "reasonID==-1","H*");
energy->Scale(1/2000000.);
gPad->SetLogx();
gPad->SetLogy();

c1->cd(6);
TH1F *energy2 = new TH1F("Energy_1000kev", "Neutron energy spectrum between 0 and 1000 keV for reason ID == -1",100,0.0,log(1000.));
EverythingN->Draw("log10((energyK))>>Energy_1000kev", "reasonID==-1","H*");
energy2->Scale(1/2000000.);
gPad->SetLogx();
gPad->SetLogy();
*/
}
////////////////////////////////////////////////////////


//////////////////////////// reason ID==0 ///////////////////////////////
TCanvas* c2 = new TCanvas("c2", "c2", w, h);
  c2->SetWindowSize(w + (w - c2->GetWw()), h + (h - c2->GetWh()));
c2->Divide(2,2);

if(c2){
c2->cd(1);
TH2F *posXYr0 = new TH2F("posXYr0", "The distribution of neutrons escaped from the target geometry (XY), reasonID==0", 401, -200.5, 200.5, 401, -200.5, 200.5);
if(posXYr0){posXYr0->GetXaxis()->SetTitle("X (mm)"); posXYr0->GetYaxis()->SetTitle("Y (mm)");  posXYr0->GetZaxis()->SetLabelSize(0.025);}
EverythingN->Draw("NPosY:NPosX>>posXYr0", "reasonID==0","COLZ1");
posXYr0->Scale(1/NumberP);
//posXYr0->GetZaxis()->SetRangeUser(0.0, 0.2);
//new TCanvas();
c2->cd(2);
TH2F *posXZr0 = new TH2F("posXZr0", "The distribution of neutrons escaped from the target geometry (XZ), reasonID==0", 401, -200.5, 200.5, 401, -800.5, -400.5);
if(posXZr0){posXZr0->GetXaxis()->SetTitle("X (mm)"); posXZr0->GetYaxis()->SetTitle("Z (mm)"); posXZr0->GetZaxis()->SetLabelSize(0.025);}
EverythingN->Draw("NPosZ:NPosX>>posXZr0", "reasonID==0","COLZ1");
posXZr0->Scale(1/NumberP);
//posXZr0->GetZaxis()->SetRangeUser(0.0, 0.2);
//new TCanvas();
c2->cd(3);
TH2F *posYZr0 = new TH2F("posYZr0", "The distribution of neutrons escaped from the target geometry (YZ), reasonID==0", 401, -200.5, 200.5, 401, -800.5, -400.5);
if(posYZr0){posYZr0->GetXaxis()->SetTitle("Y (mm)"); posYZr0->GetYaxis()->SetTitle("Z (mm)"); posYZr0->GetZaxis()->SetLabelSize(0.025);}
EverythingN->Draw("NPosZ:NPosY>>posYZr0", "reasonID==0","COLZ1");
posYZr0->Scale(1/NumberP);
//posYZr0->GetZaxis()->SetRangeUser(0.0, 0.2);
c2->cd(4);
//new TCanvas();
TH2F *MonXYr0 = new TH2F("MomentumXYr0", "The momentum distribution of neutrons escaped from the target geometry (XY),reasonID==0", 1401, -700.5, 700.5, 1401, -700.5, 700.5);
if(MonXYr0){MonXYr0->GetXaxis()->SetTitle("P_{x}"); MonXYr0->GetYaxis()->SetTitle("P_{Y}");}
EverythingN->Draw("mY:mX>>MonXYr0", "reasonID==0","");


/*
c2->cd(5);
TH1F *energyr0 = new TH1F("Energy_fullr0", "Neutron energy spectrum for reason ID==0",30,0.,log10(300.));
EverythingN->Draw("log10((energyK/1000.))>>Energy_fullr0", "reasonID==0","H*");
energyr0->Scale(1/2000000.);
gPad->SetLogx();
gPad->SetLogy();

c2->cd(6);
TH1F *energy2r0 = new TH1F("Energy_100kevr0", "Neutron energy spectrum between 0 and 1000 keV for reason ID == 0",100,0.0,log10(1000.));
EverythingN->Draw("log10((energyK))>>Energy_100kevr0", "reasonID==0","H*");
energy2r0->Scale(1/2000000.);
gPad->SetLogx();
gPad->SetLogy();
*/
}
//////////////////////////////////////////reasonID==1/////////////////////////////////////////////////////
TCanvas* c3 = new TCanvas("c3", "", 50, 50, 1000, 1000);
c3->Divide(2,2);

if(c3){

c3->cd(1);
TH2F *posXYr1 = new TH2F("posXYr1", "The distribution of neutrons escaped from the target geometry (XY), reasonID==1", 401, -200.5, 200.5, 401, -200.5, 200.5);
if(posXYr1){posXYr1->GetXaxis()->SetTitle("X (mm)"); posXYr1->GetYaxis()->SetTitle("Y (mm)"); posXYr1->GetZaxis()->SetLabelSize(0.025); }
EverythingN->Draw("NPosY:NPosX>>posXYr1", "reasonID==1","COLZ1");
posXYr1->Scale(1/NumberP);
//posXYr1->GetZaxis()->SetRangeUser(0.0, 0.2);
//new TCanvas();
c3->cd(2);
TH2F *posXZr1 = new TH2F("posXZr1", "The distribution of neutrons escaped from the target geometry (XZ), reasonID==1", 401, -200.5, 200.5, 401, -800.5, -400.5);
if(posXZr1){posXZr1->GetXaxis()->SetTitle("X (mm)"); posXZr1->GetYaxis()->SetTitle("Z (mm)"); posXZr1->GetZaxis()->SetLabelSize(0.025); }
EverythingN->Draw("NPosZ:NPosX>>posXZr1", "reasonID==1","COLZ1");
posXZr1->Scale(1/NumberP);
//posXZr1->GetZaxis()->SetRangeUser(0.0, 0.2);
//new TCanvas();
c3->cd(3);
TH2F *posYZr1 = new TH2F("posYZr1", "The distribution of neutrons escaped from the target geometry (YZ), reasonID==1", 401, -200.5, 200.5, 401, -800.5, -400.5);
if(posYZr1){posYZr1->GetXaxis()->SetTitle("Y (mm)"); posYZr1->GetYaxis()->SetTitle("Z (mm)"); posYZr1->GetZaxis()->SetLabelSize(0.025);}
EverythingN->Draw("NPosZ:NPosY>>posYZr1", "reasonID==1","COLZ1");
posYZr1->Scale(1/NumberP);
//posYZr1->GetZaxis()->SetRangeUser(0.0, 0.2);
c3->cd(4);
//new TCanvas();
TH2F *MonXYr1 = new TH2F("MomentumXYr1", "The momentum distribution of neutrons escaped from the target geometry (XY),reasonID==1", 201, -100.5, 100.5, 201, -100.5, 100.5);
if(MonXYr1){MonXYr1->GetXaxis()->SetTitle("P_{x}"); MonXYr1->GetYaxis()->SetTitle("P_{y}");}
EverythingN->Draw("mY:mX>>MomentumXYr1", "reasonID==1","");
}
////////////////////////////////////////////////////////reason ID ==2 //////////////////////////////
TCanvas* c4 = new TCanvas("c4", "", 50, 50, 1000, 1000);
c4->Divide(2,2);
if(c4){

c4->cd(1);
TH2F *posXYr2 = new TH2F("posXYr2", "The distribution of neutrons escaped from the target geometry (XY), reasonID==2", 401, -200.5, 200.5, 401, -200.5, 200.5);
if(posXYr2){posXYr2->GetXaxis()->SetTitle("X (mm)"); posXYr2->GetYaxis()->SetTitle("Y (mm)"); posXYr2->GetZaxis()->SetLabelSize(0.025);}
EverythingN->Draw("NPosY:NPosX>>posXYr2", "reasonID==2","COLZ1");
posXYr2->Scale(1/NumberP);
//posXYr2->GetZaxis()->SetRangeUser(0.0, 0.2);
//new TCanvas();
c4->cd(2);
TH2F *posXZr2 = new TH2F("posXZr2", "The distribution of neutrons escaped from the target geometry (XZ), reasonID==2", 401, -200.5, 200.5, 401, -800.5, -400.5);
if(posXZr2){posXZr2->GetXaxis()->SetTitle("X (mm)"); posXZr2->GetYaxis()->SetTitle("Z (mm)"); posXZr2->GetZaxis()->SetLabelSize(0.025);}
EverythingN->Draw("NPosZ:NPosX>>posXZr2", "reasonID==2","COLZ1");
posXZr2->Scale(1/NumberP);
//posXZr2->GetZaxis()->SetRangeUser(0.0, 0.2);
//new TCanvas();
c4->cd(3);
TH2F *posYZr2 = new TH2F("posYZr2", "The distribution of neutrons escaped from the target geometry (YZ), reasonID==2", 401, -200.5, 200.5, 401, -800.5, -400.5);
if(posYZr2){posYZr2->GetXaxis()->SetTitle("Y (mm)"); posYZr2->GetYaxis()->SetTitle("Z (mm)");posYZr2->GetZaxis()->SetLabelSize(0.025);}
EverythingN->Draw("NPosZ:NPosY>>posYZr2", "reasonID==2","COLZ1");
posYZr2->Scale(1/NumberP);
//posYZr2->GetZaxis()->SetRangeUser(0.0, 0.2);
c4->cd(4);
//new TCanvas();
TH2F *MonXYr2 = new TH2F("MomentumXYr2", "The momentum distribution of neutrons escaped from the target geometry (XY), reasonID==2", 201, -100.5, 100.5, 201, -100.5, 100.5);
if(MonXYr2){MonXYr2->GetXaxis()->SetTitle("P_{x}"); MonXYr2->GetYaxis()->SetTitle("P_{y}");}
EverythingN->Draw("mY:mX>>MomentumXYr2", "reasonID==2","");
}


TCanvas* c5 = new TCanvas("c5", "", 50, 50, 1000, 1000);
c5->Divide(2,2);
if(c5){

c5->cd(1);
TH2F *posXYr3 = new TH2F("posXYr3", "The distribution of neutrons escaped from the target geometry (XY), reasonID==3", 401, -200.5, 200.5, 401, -200.5, 200.5);
if(posXYr3){posXYr3->GetXaxis()->SetTitle("X (mm)"); posXYr3->GetYaxis()->SetTitle("Y (mm)"); posXYr3->GetZaxis()->SetLabelSize(0.025);}
EverythingN->Draw("NPosY:NPosX>>posXYr3", "reasonID==3","COLZ1");
posXYr3->Scale(1/NumberP);
//posXYr3->GetZaxis()->SetRangeUser(0.0, 0.2);
//new TCanvas();
c5->cd(2);
TH2F *posXZr3 = new TH2F("posXZr3", "The distribution of neutrons escaped from the target geometry (XZ), reasonID==3", 401, -200.5, 200.5, 401, -800.5, -400.5);
if(posXZr3){posXZr3->GetXaxis()->SetTitle("X (mm)"); posXZr3->GetYaxis()->SetTitle("Z (mm)"); posXZr3->GetZaxis()->SetLabelSize(0.025);}
EverythingN->Draw("NPosZ:NPosX>>posXZr3", "reasonID==3","COLZ1");
posXZr3->Scale(1/NumberP);
//posXZr3->GetZaxis()->SetRangeUser(0.0, 0.2);
//new TCanvas();
c5->cd(3);
TH2F *posYZr3 = new TH2F("posYZr3", "The distribution of neutrons escaped from the target geometry (YZ), reasonID==3", 401, -200.5, 200.5, 401, -800.5, -400.5);
if(posYZr3){posYZr3->GetXaxis()->SetTitle("Y (mm)"); posYZr3->GetYaxis()->SetTitle("Z (mm)");posYZr3->GetZaxis()->SetLabelSize(0.025);}
EverythingN->Draw("NPosZ:NPosY>>posYZr3", "reasonID==3","COLZ1");
posYZr3->Scale(1/NumberP);
//posYZr3->GetZaxis()->SetRangeUser(0.0, 0.2);
c5->cd(4);
//new TCanvas();
TH2F *MonXYr3 = new TH2F("MomentumXYr3", "The momentum distribution of neutrons escaped from the target geometry (XY), reasonID==3", 201, -100.5, 100.5, 201, -100.5, 100.5);
if(MonXYr3){MonXYr3->GetXaxis()->SetTitle("P_{x}"); MonXYr3->GetYaxis()->SetTitle("P_{y}");}
EverythingN->Draw("mY:mX>>MomentumXYr3", "reasonID==3","");
}



TCanvas* c6 = new TCanvas("c6", "", 50, 50, 1000, 1000);
c6->Divide(2,2);
if(c6){

c6->cd(1);
TH2F *posXYdiff12 = new TH2F("posXYdiff12", "The distribution of neutrons escaped from the target geometry (XY), r1+2", 401, -200.5, 200.5, 401, -200.5, 200.5);
EverythingN->Draw("NPosY:NPosX>>posXYdiff12", "reasonID==1","COLZ1");
posXYdiff12->Add(posXYr2,1);
posXYdiff12->Draw("COLZ1");
TH2F *posXYdiff23 = new TH2F("posXYdiff23", "The distribution of neutrons escaped from the target geometry (XY), r12+3", 401, -200.5, 200.5, 401, -200.5, 200.5);
EverythingN->Draw("NPosY:NPosX>>posXYdiff23", "reasonID==3","COLZ1");
posXYdiff23->Add(posXYdiff12,1);
posXYdiff23->Draw("COLZ1");

TH2F *posXYdiffr0 = new TH2F("posXYdiffr0", "The distribution of neutrons r-1 - (r1+r2+r3) (XY)", 401, -200.5, 200.5, 401, -200.5, 200.5);
EverythingN->Draw("NPosY:NPosX>>posXYdiffr0", "reasonID==-1","COLZ1");
Double_t entry1 = posXYdiff12->GetEffectiveEntries();
Double_t entry2 = posXYdiff23->GetEffectiveEntries();
Double_t entry3 = posXYdiffr0->GetEffectiveEntries();
Double_t entry4 = entry3 - entry2;
//posXYdiffr0->SetEntries(entry2);
posXYdiffr0->Add(posXYdiff23,-1.0);
posXYdiffr0->Draw("COLZ1");

/*TH2F *h_out = (TH2F*)posXYdiffr0->Clone("");
int nx = h_out->GetXaxis()->GetNbins();
int ny = h_out->GetYaxis()->GetNbins();
for (int i=1;i<=nx;i++) {
for (int j=1;j<=ny;j++) {
double c1 =posXYdiffr0->GetBinContent(i,j);
double c2 =posXYdiff23->GetBinContent(i,j);
h_out->SetBinContent(i,j,c1-c2);
h_out->Draw("COLZ");
}
}*/
c6->cd(2);
TH2F *posXZdiff12 = new TH2F("posXZdiff12", "The distribution of neutrons escaped from the target geometry (XZ), r1+2",401, -200.5, 200.5, 401, -800.5, -400.5);
EverythingN->Draw("NPosZ:NPosX>>posXZdiff12", "reasonID==1","COLZ1");
posXZdiff12->Add(posXZr2,1.);
posXZdiff12->Draw("COLZ1");

TH2F *posXZdiff23 = new TH2F("posXZdiff23", "The distribution of neutrons escaped from the target geometry (XZ), r12+3", 401, -200.5, 200.5, 401, -800.5, -400.5);
EverythingN->Draw("NPosZ:NPosX>>posXZdiff23", "reasonID==3","COLZ1");
posXZdiff23->Add(posXZdiff12,1.);

TH2F *posXZdiffr0 = new TH2F("posXZdiffr0", "The distribution of neutrons r-1 - (r1+r2+r3) (XZ)", 401, -200.5, 200.5, 401, -800.5, -400.5);
EverythingN->Draw("NPosZ:NPosX>>posXZdiffr0", "reasonID==-1","COLZ1");

posXZdiffr0->Add(posXZdiff23,-1.0);
posXZdiffr0->Draw("COLZ1");


c6->cd(3);
TH2F *posYZdiff12 = new TH2F("posYZdiff12", "The distribution of neutrons escaped from the target geometry (YZ), r1+2",401, -200.5, 200.5, 401, -800.5, -400.5);
EverythingN->Draw("NPosZ:NPosY>>posYZdiff12", "reasonID==1","COLZ1");
posXZdiff12->Add(posYZr2,1.);
posXZdiff12->Draw("COLZ1");

TH2F *posYZdiff23 = new TH2F("posYZdiff23", "The distribution of neutrons escaped from the target geometry (YZ), r12+3", 401, -200.5, 200.5, 401, -800.5, -400.5);
EverythingN->Draw("NPosZ:NPosY>>posYZdiff23", "reasonID==3","COLZ1");
posYZdiff23->Add(posYZdiff12,1.);

TH2F *posYZdiffr0 = new TH2F("posYZdiffr0", "The distribution of neutrons r-1 - (r1+r2+r3) (YZ)", 401, -200.5, 200.5, 401, -800.5, -400.5);
EverythingN->Draw("NPosZ:NPosY>>posYZdiffr0", "reasonID==-1","COLZ1");

posYZdiffr0->Add(posYZdiff23,-1.0);
posYZdiffr0->Draw("COLZ1");



}




}