// Run with program with 
//
//		root DrawHistograms.cpp



#include "TFile.h"
#include "TH1F.h"
#include "TVector3.h"
#include "TTreeReader.h"
#include "TTreeReaderValue.h"

#include "TCanvas.h"
#include "TGraph.h"
#include "TAxis.h"
#include "TPad.h"
#include "TSpectrum.h"
#include "TRandom.h"

#include <TString.h>
#include "TSystem.h"

#include <fstream>
#include <iostream>
using namespace std;

void DrawIntegral() {


// PRESETTINGS

	// Set path to root file
	std::string myPath = "../Integral/";

	// Set name of root file
	std::string myFileName = "Out_run_0.root";

	// Open the file containing the tree
	TFile *myFile = TFile::Open(std::string(myPath+myFileName).c_str());

	// Create a TTreeReader for the tree by passing the TTree's name and the TFile it is in
	TTreeReader myReader("Secondaries", myFile);

	// Set variables to access the data in the tree
	
	//The branch "ParticleID" contains a vector<int>; access them as vID
   TTreeReaderValue<vector<int>> vPDGcode (myReader, "ParticleID"); //New ADD
   
   // A histogram for Particle ID ranging from 1 to 100  with 200 bins
	TH1F* hPDGcode = new TH1F("hPDGcode", ";Number;#entries", 2500, -30, 30); //New ADD
	

	// The branch "ParticleEnergy" contains a vector<double>; access them as vEnergies
	TTreeReaderValue<vector<double>> vEnergies(myReader, "ParticleEnergy");
	
	// Create histograms for the values we read

	// A histogram for the energies ranging from 0 MeV to 20 MeV with 200 bins
	TH1F* hEnergyGamma = new TH1F("hEnergyGamma", ";E in MeV;#entries", 600, 0,10);
	
	TH1F* hEnergyAll = new TH1F("hEnergyAll", ";E in MeV;#entries", 600, 0, 10);


// READ THE TREE
  
	// Loop over all entries of the TTree.
    
    while (myReader.Next()) {
    int energySize = vEnergies->size();
		
		for(int i = 0; i<energySize; ++i){
			
			double myEnergy= vEnergies->at(i);
			int myPDG= vPDGcode->at(i);
		
			hEnergyAll->Fill(myEnergy);
			if (myPDG==22){
				hEnergyGamma->Fill(myEnergy);
				}
				
				
			}
    } 

     TAxis  *axis = hEnergyGamma->GetXaxis();
     double binx1 = axis->FindBin(4.3);
     double binx2=  axis->FindBin(4.6);

     double integral = hEnergyGamma->Integral(binx1, binx2);
         
     
     //std::cout << "#####$$$$$  Area under the Peak(4.4) = " << integral << std::endl;
     

/*// DRAW THE HISTOGRAMS


	// Energy histogram

	// Create a canvas to draw into
	TCanvas* cEnergy = new TCanvas();
	// Change into this canvas (now it is the active one and everything drawn will end up in this canvas)
	cEnergy->cd();
	// Draw the histogram
	hEnergyGamma->Draw();
	// Force an update of the graphical view
	cEnergy->Update();
	// Save the image as pdf
	cEnergy->SaveAs("../images/EnergySpectrum.pdf");*/

	
}