#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <sstream>
#include <cstdlib>
#include <math.h>

#include <vector>

#include <TApplication.h>
#include <TAxis.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <TH1F.h>
#include <TLatex.h>
#include <TLegend.h>
#include <TLine.h>
#include <TRandom3.h>
#include <TROOT.h>
#include <TString.h>
#include <TSystem.h>
#include <TStyle.h>
#include <TSpectrum.h>

using namespace std;

void simulation_analysis() {
  
  //depth fix - water box centred at -226 cm, half length = 10cm. 
  Double_t depthFix = 2360;
  

  //change this line if you modify the length and width of the box
  Float_t lengthBox = 200, widthBox = 200;
  
  // Create histograms (stopping distance and energy deposited)
  
  TH1D* stoppingDistance = new TH1D("stoppingDistance","62.5 MeV Proton Stopping Distance in Water", 100, 0., 65.);
  TH1D* waterEnergy = new TH1D ("waterDep", "Energy deposition in water box", 800, 0., 65.);
  TH1D* productionNeutrons = new TH1D ("neutronEnergy", "Kinetic energy of production neutrons", 800, 0., 65.);
  TH1D* SDEnergy = new TH1D ("SDDep", "Energy deposition of protons in user-defined sensitive detector", 800, 0., 65.);
  TH1D* kineticEnergySpectrum = new TH1D ("kinEn", "Kinetic energy of proton beam at specified z position along beamline", 800, 0., 65.);

  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  string data_file_name = "/unix/pbt/rstephens/ProtonPBFolder/ProtonPB_build/hits.out";

  ifstream resultsFile;
  resultsFile.open(data_file_name.c_str());
  
  if(resultsFile.is_open()){
    
    std::cout << "File selected to open: " << data_file_name << std::endl;
    
    string data;
    
    while (getline(resultsFile,data)) {

      stringstream dataReadin;
      dataReadin << data;
      
      Double_t distance = 0;

      dataReadin >> distance;
  
      stoppingDistance->Fill(distance+depthFix);
    }
    
    resultsFile.close();
  }
   
  //read in sensitive detector information - water box and user-defined sensitive detector

   TString stopperFile = "/unix/pbt/rstephens/ProtonPBFolder/ProtonPB_build/sensitiveDetectors.out";
   
  ifstream stopperRead;
  stopperRead.open(stopperFile);
  int counter = 0;
  if(stopperRead.is_open()){
    string dataStopper;
    std::cout << "File selected to open: " << stopperFile << std::endl;

    while (getline(stopperRead,dataStopper)) {

      stringstream dataStopperReadin;
      dataStopperReadin << dataStopper;
      // each sensitive detector has a unique ID
      int SD_ID;
      Double_t energyDep = 0;

      dataStopperReadin >> SD_ID >> energyDep;
      if(energyDep==0){
	counter++;
      }
      if(energyDep!=0 && SD_ID==0){
	waterEnergy->Fill(energyDep);
      }
      if(energyDep!=0 && SD_ID==1){
	SDEnergy->Fill(energyDep);
      }
    }
    stopperRead.close();
    std::cout<<"Number of zeroes"<< counter << std::endl;
  }

  TString neutronFile = "/unix/pbt/rstephens/ProtonPBFolder/ProtonPB_build/secondariesKE.out";
  
  ifstream neutronRead;
  neutronRead.open(neutronFile);
  if(neutronRead.is_open()){
    string neutronData;
    while(getline(neutronRead, neutronData)){
	stringstream neutronReadIn;
	neutronReadIn << neutronData;
	string particleName;
	Double_t secondaryEnergy = 0;
	neutronReadIn >> secondaryEnergy;
	productionNeutrons->Fill(secondaryEnergy);

    }
    neutronRead.close();
  }

  
  //read in kinetic energy info for specified z position along beamline
   TString energyFile = "/unix/pbt/rstephens/ProtonPBFolder/ProtonPB_build/protonKE.out";
   TNtuple *TNtuplePosition = new TNtuple("ntuplePosition","Position from simulation file", "xPos:yPos:zPos:energy");

  ifstream energyRead;
  energyRead.open(energyFile);
  int counter = 0;
  if(energyRead.is_open()){
    string dataEnergy;
    std::cout << "File selected to open: " << energyFile << std::endl;
    while (getline(energyRead,dataEnergy)) {

      stringstream dataEnergyReadin;
      dataEnergyReadin << dataEnergy;
      
      Double_t positionX;
      Double_t positionY;
      Double_t positionZ;
      Double_t KE = 0;

      dataEnergyReadin >> positionX >> positionY >> positionZ >> KE;
      kineticEnergySpectrum->Fill(KE);
      TNtuplePosition->Fill(positionX, positionY, positionZ, KE);

      
    }
    energyRead.close();
  }

  ///////////////////////////////////////////////////////////////////////////////////////////////////
 

     // Analysis of proton flux inside water box

    TString data_file_name4 = "/unix/pbt/rstephens/ProtonPBFolder/ProtonPB_build/ProtonFluxWaterLongitudinal.txt";

    ifstream resultsFile4;
    resultsFile4.open(data_file_name4);
    
    ofstream outfileFlux;
    TString fluxFileSim = "ProtonFluxWaterLongitudinal_Mod.txt";
    outfileFlux.open(fluxFileSim);
    
    string line;
    getline(resultsFile4,line);
    getline(resultsFile4,line);
    getline(resultsFile4,line);

    while( !resultsFile4.eof()){
	getline(resultsFile4,line, ',');
	outfileFlux << line << " ";
    }
    resultsFile4.close();
    outfileFlux.close();


    TNtuple *TNtupleSimFlux = new TNtuple("ntupleSimulationFlux","Flux from simulation file", "iX:jY:kZ:fluxSim");

    cout << "Reading file \" " << fluxFileSim << "\" ... ";
    Long64_t nlines = TNtupleSimFlux -> ReadFile(fluxFileSim, "iX:jY:kZ:fluxSim"); 
    if (nlines <=0){cout << "Error: Check file \"" << fluxFileSim << "\"\n"; return;}
    printf("%d Experimental points found\n", nlines);
        
    Float_t kZ, fluxSim;
    TNtupleSimFlux -> SetBranchAddress("kZ", &kZ);
    TNtupleSimFlux -> SetBranchAddress("fluxSim", &fluxSim);

     // Normalise and fill TNtupleSimFlux
    Int_t nentries = (Int_t)TNtupleSimFlux -> GetEntries();   
    TNtupleSimFlux -> GetEntry(0);
    Float_t maxDose = 0, depthSim = 0, slength = lengthBox/nentries;
    vector <Float_t> vec_dose, vec_kZ;

    for (Int_t l = 0; l<nentries; l++)
    {
      TNtupleSimFlux -> GetEntry(l);
      if (maxDose < fluxSim) { maxDose = fluxSim;
                               vec_dose.push_back(fluxSim);
			       vec_kZ.push_back(depthSim);                             
                              }
      else { vec_dose.push_back(fluxSim);
	     vec_kZ.push_back(depthSim);

            }
      depthSim += slength;
    }
    
    TNtupleSimFlux -> Reset(); 
    
    // normalisation to the maximum dose
    for (Int_t l = 0; l<vec_dose.size();l++)
    {
        kZ = vec_kZ[l];
	fluxSim = 100*(vec_dose[l]/maxDose);
        TNtupleSimFlux -> Fill(0,0,kZ,fluxSim);
    }
    
    TNtupleSimFlux -> Fill(0,0,kZ,fluxSim);

    //////////////////////////////////////////////////////////////////////////////////////////////////
    
  // Analysis of longitudinal proton energy deposition in water box


    TString data_file_name5 = "/unix/pbt/rstephens/ProtonPBFolder/ProtonPB_build/ProtonEnergyWaterLongitudinal.txt";

    ifstream resultsFile5;
    resultsFile5.open(data_file_name5);
    
    ofstream outfileEnergy;
    TString depFileSim = "ProtonEnergyWaterLongitudinal_Mod.txt";
    outfileEnergy.open(depFileSim);
    
    string line;
    getline(resultsFile5,line);
    getline(resultsFile5,line);
    getline(resultsFile5,line);

    while( !resultsFile5.eof()){
	getline(resultsFile5,line, ',');
	outfileEnergy << line << " ";
    }
    resultsFile5.close();
    outfileEnergy.close();


    TNtuple *TNtupleDep = new TNtuple("ntupleSimulationFlux","Flux from simulation file", "iX:jY:kZ:fluxSim");

    cout << "Reading file \" " << depFileSim << "\" ... ";
    Long64_t nlines = TNtupleDep -> ReadFile(depFileSim, "iX:jY:kZ:fluxSim"); 
    if (nlines <=0){cout << "Error: Check file \"" << depFileSim << "\"\n"; return;}
    printf("%d Experimental points found\n", nlines);
        
    Float_t kZ, fluxSim;
    TNtupleDep -> SetBranchAddress("kZ", &kZ);
    TNtupleDep -> SetBranchAddress("fluxSim", &fluxSim);
    
     // Normalise and fill TNtupleDep
    Int_t nentries = (Int_t)TNtupleDep -> GetEntries();   
    TNtupleDep -> GetEntry(0);
    Float_t maxDose = 0, depthSim = 0, slength = lengthBox/nentries;
    vector <Float_t> vec_dose, vec_kZ;

    for (Int_t l = 0; l<nentries; l++)
    {
      TNtupleDep -> GetEntry(l);
      if (maxDose < fluxSim) { maxDose = fluxSim;
                               vec_dose.push_back(fluxSim);
			       vec_kZ.push_back(depthSim);                             
                              }
      else { vec_dose.push_back(fluxSim);
	     vec_kZ.push_back(depthSim);

            }
      depthSim += slength;
    }
    
    TNtupleDep -> Reset(); 
     
    // normalisation to the maximum dose
    for (Int_t l = 0; l<vec_dose.size();l++)
    {
        kZ = vec_kZ[l];
	//	fluxSim = 100*(vec_dose[l]/maxDose);
	fluxSim = vec_dose[l];
        TNtupleDep -> Fill(0,0,kZ,fluxSim);
    }
   
    TNtupleDep -> Fill(0,0,kZ,fluxSim);
    
    //////////////////////////////////////////////////////////////////////////////////////////////////
    
  // Analysis of lateral energy deposition of protons


    TString data_file_name10 = "/unix/pbt/rstephens/ProtonPBFolder/ProtonPB_build/EnergyLateralMesh.txt";

    ifstream resultsFile10;
    resultsFile10.open(data_file_name10);
    
    ofstream outfileLateralEnergy;
    TString latEnFile = "EnergyLateralMesh_Mod.txt";
    outfileLateralEnergy.open(latEnFile);
    
    string line;
    getline(resultsFile10,line);
    getline(resultsFile10,line);
    getline(resultsFile10,line);

    while( !resultsFile10.eof()){
	getline(resultsFile10,line, ',');
	outfileLateralEnergy << line << " ";
    }
    resultsFile10.close();
    outfileLateralEnergy.close();


    TNtuple *TNtupleProtonLateral = new TNtuple("ntupleNeutronSimulationFlux","Neutron flux from simulation file", "iX:jY:kZ:fluxSim");

    cout << "Reading file \" " << latEnFile << "\" ... ";
    Long64_t nlines = TNtupleProtonLateral -> ReadFile(latEnFile, "iX:jY:kZ:fluxSim"); 
    if (nlines <=0){cout << "Error: Check file \"" << latEnFile << "\"\n"; return;}
    printf("%d Experimental points found\n", nlines);
        
    Float_t kZ, fluxSim;
    TNtupleProtonLateral -> SetBranchAddress("kZ", &kZ);
    TNtupleProtonLateral -> SetBranchAddress("fluxSim", &fluxSim);
    
     // Normalise and fill TNtupleProtonLateral
    Int_t nentries = (Int_t)TNtupleProtonLateral -> GetEntries();   
    TNtupleProtonLateral -> GetEntry(0);
    Float_t maxDose = 0, depthSim = 0, slength = lengthBox/nentries;
    vector <Float_t> vec_dose, vec_kZ;

    for (Int_t l = 0; l<nentries; l++)
    {
      TNtupleProtonLateral -> GetEntry(l);
      if (maxDose < fluxSim) { maxDose = fluxSim;
                               vec_dose.push_back(fluxSim);
			       vec_kZ.push_back(depthSim);                             
                              }
      else { vec_dose.push_back(fluxSim);
	     vec_kZ.push_back(depthSim);

            }
      depthSim += slength;
    }
    
    TNtupleProtonLateral -> Reset(); 
     
    // normalisation to the maximum dose
    for (Int_t l = 0; l<vec_dose.size();l++)
    {
        kZ = vec_kZ[l];
	//	fluxSim = 100*(vec_dose[l]/maxDose);
	fluxSim = vec_dose[l];
        TNtupleProtonLateral -> Fill(0,0,kZ,fluxSim);
    }
   
    TNtupleProtonLateral -> Fill(0,0,kZ,fluxSim);

    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


 // Analysis of neutron flux  along beamline


    TString data_file_name7 = "/unix/pbt/rstephens/ProtonPBFolder/ProtonPB_build/NeutronFluxBeamline.txt";

    ifstream resultsFile7;
    resultsFile7.open(data_file_name7);
    
    ofstream outfileNeutronFlux;
    TString neutronFluxFile = "NeutronFluxBeamline_Mod.txt";
    outfileNeutronFlux.open(neutronFluxFile);
    
    string line;
    getline(resultsFile7,line);
    getline(resultsFile7,line);
    getline(resultsFile7,line);

    while( !resultsFile7.eof()){
	getline(resultsFile7,line, ',');
	outfileNeutronFlux << line << " ";
    }
    resultsFile7.close();
    outfileNeutronFlux.close();


    TNtuple *TNtupleNeutronFlux = new TNtuple("ntupleNeutronSimulationFlux","Neutron flux from simulation file", "iX:jY:kZ:fluxSim");

    cout << "Reading file \" " << neutronFluxFile << "\" ... ";
    Long64_t nlines = TNtupleNeutronFlux -> ReadFile(neutronFluxFile, "iX:jY:kZ:fluxSim"); 
    if (nlines <=0){cout << "Error: Check file \"" << neutronFluxFile << "\"\n"; return;}
    printf("%d Experimental points found\n", nlines);
        
    Float_t kZ, fluxSim;
    TNtupleNeutronFlux -> SetBranchAddress("kZ", &kZ);
    TNtupleNeutronFlux -> SetBranchAddress("fluxSim", &fluxSim);
    
     // Normalise and fill TNtupleNeutronFlux
    Int_t nentries = (Int_t)TNtupleNeutronFlux -> GetEntries();   
    TNtupleNeutronFlux -> GetEntry(0);
    Float_t maxDose = 0, depthSim = 0, slength = lengthBox/nentries;
    vector <Float_t> vec_dose, vec_kZ;

    for (Int_t l = 0; l<nentries; l++)
    {
      TNtupleNeutronFlux -> GetEntry(l);
      if (maxDose < fluxSim) { maxDose = fluxSim;
                               vec_dose.push_back(fluxSim);
			       vec_kZ.push_back(depthSim);                             
                              }
      else { vec_dose.push_back(fluxSim);
	     vec_kZ.push_back(depthSim);

            }
      depthSim += slength;
    }
    
    TNtupleNeutronFlux -> Reset(); 
     
    // normalisation to the maximum dose
    for (Int_t l = 0; l<vec_dose.size();l++)
    {
        kZ = vec_kZ[l];
	//	fluxSim = 100*(vec_dose[l]/maxDose);
	fluxSim = vec_dose[l];
        TNtupleNeutronFlux -> Fill(0,0,kZ,fluxSim);
    }
   
    TNtupleNeutronFlux -> Fill(0,0,kZ,fluxSim);
   
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

 // Analysis of proton flux  along beamline


    TString data_file_name7 = "/unix/pbt/rstephens/ProtonPBFolder/ProtonPB_build/ProtonFluxBeamline.txt";

    ifstream resultsFile7;
    resultsFile7.open(data_file_name7);
    
    ofstream outfileProtonFlux;
    TString protonFluxFile = "ProtonFluxBeamline_Mod.txt";
    outfileProtonFlux.open(protonFluxFile);
    
    string line;
    getline(resultsFile7,line);
    getline(resultsFile7,line);
    getline(resultsFile7,line);

    while( !resultsFile7.eof()){
	getline(resultsFile7,line, ',');
	outfileProtonFlux << line << " ";
    }
    resultsFile7.close();
    outfileProtonFlux.close();


    TNtuple *TNtupleProtonFlux = new TNtuple("ntupleProtonSimulationFlux","Proton flux from simulation file", "iX:jY:kZ:fluxSim");

    cout << "Reading file \" " << protonFluxFile << "\" ... ";
    Long64_t nlines = TNtupleProtonFlux -> ReadFile(protonFluxFile, "iX:jY:kZ:fluxSim"); 
    if (nlines <=0){cout << "Error: Check file \"" << protonFluxFile << "\"\n"; return;}
    printf("%d Experimental points found\n", nlines);
        
    Float_t kZ, fluxSim;
    TNtupleProtonFlux -> SetBranchAddress("kZ", &kZ);
    TNtupleProtonFlux -> SetBranchAddress("fluxSim", &fluxSim);
    
     // Normalise and fill TNtupleProtonFlux
    Int_t nentries = (Int_t)TNtupleProtonFlux -> GetEntries();   
    TNtupleProtonFlux -> GetEntry(0);
    Float_t maxDose = 0, depthSim = 0, slength = lengthBox/nentries;
    vector <Float_t> vec_dose, vec_kZ;

    for (Int_t l = 0; l<nentries; l++)
    {
      TNtupleProtonFlux -> GetEntry(l);
      if (maxDose < fluxSim) { maxDose = fluxSim;
                               vec_dose.push_back(fluxSim);
			       vec_kZ.push_back(depthSim);                             
                              }
      else { vec_dose.push_back(fluxSim);
	     vec_kZ.push_back(depthSim);

            }
      depthSim += slength;
    }
    
    TNtupleProtonFlux -> Reset(); 
     
    // normalisation to the maximum dose
    for (Int_t l = 0; l<vec_dose.size();l++)
    {
        kZ = vec_kZ[l];
	//	fluxSim = 100*(vec_dose[l]/maxDose);
	fluxSim = vec_dose[l];
        TNtupleProtonFlux -> Fill(0,0,kZ,fluxSim);
    }
   
    TNtupleProtonFlux -> Fill(0,0,kZ,fluxSim);


 
    // oo000ooo000ooo000ooo000ooo00  PLOTTING HISTOGRAMS  00ooo000ooo000ooo000ooo000ooo000oo

   
  TCanvas* canvas = new TCanvas("Proton Simulation","Proton Simulation",0,0,1200,600);
  canvas->Divide(2,2); //split the canvas into 4 sections
  
  
  canvas->cd(1);
    waterEnergy->SetDirectory(0);
    waterEnergy->GetXaxis()->SetTitle("Energy, MeV");
    waterEnergy->GetYaxis()->SetTitle("Number of events");
    waterEnergy->GetXaxis()->SetNdivisions(506);
    waterEnergy->GetYaxis()->SetTitleOffset(2);
    waterEnergy->GetYaxis()->CenterTitle();
    waterEnergy->GetXaxis()->SetRangeUser(0, 64);
    TF1* gaussian = new TF1("gaussian","gaus");
    gaussian->SetLineColor(2);
    gaussian->SetLineWidth(2);
    waterEnergy->Fit("gaussian", "", "", 55, 65);
    waterEnergy->Draw();

  canvas->cd(2);
  
 
  stoppingDistance->SetDirectory(0);
  stoppingDistance->GetXaxis()->SetTitle("Depth, mm");
  stoppingDistance->GetYaxis()->SetTitle("Number of events");
  stoppingDistance->GetYaxis()->SetTitleOffset(2);
  stoppingDistance->GetYaxis()->CenterTitle();
  stoppingDistance->GetXaxis()->SetNdivisions(506);
  stoppingDistance->GetXaxis()->SetRangeUser(0, 300);
  stoppingDistance->Fit("gaussian","","",5, 34);
  stoppingDistance->Draw();

  
  canvas->cd(3);
     
    TNtupleDep-> SetMarkerStyle(20);
    TNtupleDep -> SetMarkerColor(4);
    TNtupleDep -> SetMarkerSize(0.5);
       
    TNtupleDep ->  Draw("fluxSim:kZ");
    TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
    TAxis *xaxis = htemp->GetXaxis();
    xaxis->SetRangeUser(0, 40);

    TAxis *yaxis = htemp->GetYaxis();
    xaxis->SetTitle("Depth [mm]");
    yaxis->SetTitle("Energy deposited in water box [MeV]");
    yaxis->SetTitleOffset(1.5);

    htemp -> SetTitle("Longitudinal energy deposition in water box");

  
  canvas->cd(4);

    TNtupleSimFlux-> SetMarkerStyle(20);
    TNtupleSimFlux -> SetMarkerColor(4);
    TNtupleSimFlux -> SetMarkerSize(0.5);
       
    TNtupleSimFlux ->  Draw("fluxSim:kZ");
    TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
    TAxis *xaxis = htemp->GetXaxis();
    xaxis->SetRangeUser(0, 40);

    TAxis *yaxis = htemp->GetYaxis();
    xaxis->SetTitle("Depth [mm]");
    yaxis->SetTitle("Percentage of initial proton flux, %");
    htemp -> SetTitle("Proton flux through water box");




    TCanvas* protonlat = new TCanvas("Lateral proton energy", "Lateral energy deposition of protons", 0, 0, 800, 400);

    TNtupleProtonLateral-> SetMarkerStyle(20);
    TNtupleProtonLateral -> SetMarkerColor(4);
    TNtupleProtonLateral -> SetMarkerSize(0.5);
       
    TNtupleProtonLateral ->  Draw("fluxSim:kZ");
    TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
    TAxis *xaxis = htemp->GetXaxis();
    xaxis->SetRangeUser(0, 200);

    TAxis *yaxis = htemp->GetYaxis();
    xaxis->SetTitle("Distance [mm]");
    yaxis->SetTitle("Energy deposition by protons [MeV]");
    htemp -> SetTitle("Lateral energy deposition of protons directly after nozzle");

    ///////////////////////////////////////////////////////

    TCanvas* nc = new TCanvas("Neutron Flux", "Neutron Flux along Beamline", 0, 0, 1800, 600);
    nc->Divide(2,1);

    nc->cd(1);
    TNtupleNeutronFlux-> SetMarkerStyle(20);
    TNtupleNeutronFlux -> SetMarkerColor(4);
    TNtupleNeutronFlux -> SetMarkerSize(0.5);
       
    TNtupleNeutronFlux ->  Draw("fluxSim:kZ");
    TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
    TAxis *xaxis = htemp->GetXaxis();
    xaxis->SetRangeUser(0, 200);

    TAxis *yaxis = htemp->GetYaxis();
    xaxis->SetTitle("Distance from proton source [cm]");
    yaxis->SetTitle("Neutron flux along beamline per cm^{2}");
    htemp -> SetTitle("Neutron flux along beamline");

    nc->cd(2);

    productionNeutrons->SetDirectory(0);
    productionNeutrons->GetXaxis()->SetTitle("Energy, MeV");
    productionNeutrons->GetYaxis()->SetTitle("Number of events");
    productionNeutrons->GetXaxis()->SetNdivisions(506);
    productionNeutrons->GetYaxis()->SetTitleOffset(2);
    productionNeutrons->GetYaxis()->CenterTitle();
    productionNeutrons->GetXaxis()->SetRangeUser(0, 64);
    productionNeutrons->Draw();

    ///////////////////////////////////////////////////////////////

 TCanvas* protonBeamlineFlux = new TCanvas("Proton Flux", "Proton Flux along Beamline", 0, 0, 1800, 600);

    TNtupleProtonFlux-> SetMarkerStyle(20);
    TNtupleProtonFlux -> SetMarkerColor(4);
    TNtupleProtonFlux -> SetMarkerSize(0.5);
       
    TNtupleProtonFlux ->  Draw("fluxSim:kZ");
    TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
    TAxis *xaxis = htemp->GetXaxis();
    xaxis->SetRangeUser(0, 200);

    TAxis *yaxis = htemp->GetYaxis();
    xaxis->SetTitle("Distance from proton source [cm]");
    yaxis->SetTitle("Proton flux along beamline per cm^{2}");
    htemp -> SetTitle("Proton flux along beamline");
    
    /////////////////////////////////////////////////////////////////////////////////////////////////////
  
    TCanvas* c1 = new TCanvas("SensitiveDetector","Energy deposition within user-specified beamline component", 0, 0, 600, 300);
    SDEnergy->SetDirectory(0);
    SDEnergy->GetXaxis()->SetTitle("Energy, MeV");
    SDEnergy->GetYaxis()->SetTitle("Number of events");
    SDEnergy->GetXaxis()->SetNdivisions(506);
    SDEnergy->GetYaxis()->SetTitleOffset(2);
    SDEnergy->GetYaxis()->CenterTitle();
    SDEnergy->GetXaxis()->SetRangeUser(0, 64);
    TF1* gaussian = new TF1("gaussian","gaus");
    gaussian->SetLineColor(2);
    gaussian->SetLineWidth(2);
    SDEnergy->Fit("gaussian", "", "", 60, 65);
    SDEnergy->Draw();

    ////////////////////////////////////////////////////////////////////////////////////////////////////

  TCanvas* canvas2 = new TCanvas("ProtonSimulation2","Distribution of kinetic energy and position of protons in beam",0,0,1200,600);
  canvas2->Divide(2,1);
  
  
  canvas2->cd(1);
    kineticEnergySpectrum->SetDirectory(0);
    kineticEnergySpectrum->GetXaxis()->SetTitle("Energy, MeV");
    kineticEnergySpectrum->GetYaxis()->SetTitle("Number of events");
    kineticEnergySpectrum->GetXaxis()->SetNdivisions(506);
    kineticEnergySpectrum->GetYaxis()->SetTitleOffset(2);
    kineticEnergySpectrum->GetYaxis()->CenterTitle();
    kineticEnergySpectrum->GetXaxis()->SetRangeUser(55, 65);
    TF1* gaussian = new TF1("gaussian","gaus");
    gaussian->SetLineColor(2);
    gaussian->SetLineWidth(2);
    kineticEnergySpectrum->Fit("gaussian", "", "", 55, 65);
    kineticEnergySpectrum->Draw();

  canvas2->cd(2);
     
    TNtuplePosition-> SetMarkerStyle(20);
    TNtuplePosition -> SetMarkerColor(4);
    TNtuplePosition -> SetMarkerSize(0.2);
       
    TNtuplePosition ->  Draw("energy:yPos:xPos");
    TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
    TAxis *xaxis = htemp->GetXaxis();
    TAxis *yaxis = htemp->GetYaxis();
    TAxis *zaxis = htemp->GetZaxis();
    xaxis->SetTitle("x [mm]");
    xaxis->SetTitleOffset(1.5);
    yaxis->SetTitle("y [mm]");
    yaxis->SetTitleOffset(1.5);
    zaxis->SetTitle("Kinetic energy [MeV]");
    zaxis->SetTitleOffset(1.5);
    xaxis->SetRangeUser(-25,25);
    yaxis->SetRangeUser(-25,25);
    htemp -> SetTitle("Spatial distribution of proton beam at specified z position along beamline");

    
 
    
   ///////////////////////////////////////////////////////////////////////////////////////////////////
 

}