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

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

void SameHisto()
{
   // PRESETTINGS
   // Set path to root file
   std::string myPath = "../analysis/";

   // Set name of root file
   std::vector<std::string> myFileName;
   myFileName.push_back("Out_run_0.root");
   myFileName.push_back("Out_run_1.root");

   // Create histograms for the values we read
   // A histogram for the energies ranging from 0 MeV to 20 MeV with 200 bins
   std::vector<TH1F *> EnergyAll;
   std::vector<TH1F *> EnergyGamma;
   EnergyAll.push_back(new TH1F("EnergyAll[1]", ";E in MeV;#entries", 20, 0, 20));
   EnergyAll.push_back(new TH1F("EnergyAll[2]", ";E in MeV;#entries", 20, 0, 20));
   EnergyGamma.push_back(new TH1F("EnergyGamma[1]", ";E in MeV;#entries", 20, 0, 20));
   EnergyGamma.push_back(new TH1F("EnergyGamma[2]", ";E in MeV;#entries", 20, 0, 20));

   for (int j=0; j<myFileName.size(); ++j) {
      // Open the file containing the tree
      TFile *myFile = TFile::Open(std::string(myPath+myFileName[j]).c_str());
      // Create a TTreeReader for the tree by passing the TTree's name and the TFile it is in
      TTreeReader myReader("Secondaries", myFile);
      // The branch "ParticleEnergy" contains a vector<double>; access them as vEnergies
      TTreeReaderValue<vector<double>> vEnergies(myReader, "ParticleEnergy");
      //The branch "ParticleID" contains a vector<int>; access them as vID
      TTreeReaderValue<vector<int>> vPDGcode (myReader, "ParticleID");

      // 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);
            EnergyAll[j]->Fill(myEnergy);
            if (myPDG==22) {
               EnergyGamma[j]->Fill(myEnergy);
            }
         }
      }
      myFile->Close();
      delete myFile;
   }

   // DRAW THE HISTOGRAM
   // Energy Histogram
   TCanvas *cEnergy = new TCanvas();
   cEnergy->cd();
   THStack *hs = new THStack("hs", ";E in MeV;#entries");

   EnergyGamma[0]->SetFillColor(kBlue);
   EnergyGamma[0]->SetFillStyle(3004);

   EnergyGamma[1]->SetFillColor(kRed);
   EnergyGamma[1]->SetFillStyle(3005);

   hs->Add(EnergyGamma[0]);
   hs->Add(EnergyGamma[1]);
   hs->Draw("nostack");
   //Adding Legend
   TLegend *legend = new TLegend(0.7, 0.7, 0.95, 0.95);
   legend->SetHeader("Energy Spectrum");
   legend->AddEntry(EnergyGamma[0], "PMMAPhantom EnergySpectrum", "l");
   legend->SetTextColor(kBlue);
   legend->AddEntry(EnergyGamma[1]s, "CTDICOM EnergySpectrum", "l");
   legend->SetTextColor(kRed);
   legend->Draw();
   cEnergy->Update();
   cEnergy->SaveAs("../images/EnergySpectrum.pdf");
}