// ExR04Histograms.cc
// 29-Sep-2000 Bill Seligman
// This is an example class for booking histograms of hit information
// in GEANT4.

#include "ExR04Histograms.hh"

// Include files for ROOT.
#include "Rtypes.h"
#include "TROOT.h"
#include "TFile.h"
#include "TH1.h"

// Needed to convert numbers to strings
#include <strstream>
#include <string>

// Include files for the G4 classes
#include "globals.hh"
#include "G4Run.hh"
#include "G4Event.hh"

// Classes needed for access to hit information.
#include "G4SDManager.hh"
#include "G4HCofThisEvent.hh"
#include "ExN02TrackerHit.hh"

ExR04Histograms::ExR04Histograms() {

  // Since we're not creating the TROOT object in the same scope as we're
  // creating our histogram objects, we must declare it to be static.
  // Otherwise, the TROOT object will go away when the constructor is over.

  // Note: an alternative would be to create rootBase in the main program
  // (exampleR04.cc).  But if we did that, the main program would "know" that
  // ROOT is being used to make histograms (as opposed to HBOOK, LHC++, etc.).
  // I'd rather keep that information internal to this class.
  static TROOT rootBase("simple","Test of histogramming and I/O");

  // Open a new ROOT file.
  m_histFile = new TFile("exampleR04.root","RECREATE","Demo ROOT file with histograms"); 
}

ExR04Histograms::~ExR04Histograms() {

  // Close the file. Note that this is automatically done when you leave
  // the application.
  m_histFile->Close();

  // We created the file with "new", so we have to destroy it with "delete":
  delete m_histFile;
}


void ExR04Histograms::RecordBeginOfRun(const G4Run* a_run)
{
  // Create new histograms for this run.

  // Note the use of ROOT data types as arguments to the ROOT method.
  // Need "char *" because ROOT doesn't know about STL strings, only const char *.
  char * histTitle;
  Int_t numBins;
  Axis_t zLow, zHigh;

  // Create a histogram name that changes with the run number.
  string histName;
  ostrstream convert2strA;
  // "ends" is the ostrsteam equivalent of "endl"; that is, it ends the string.
  convert2strA << "Hits" << a_run->GetRunID() << std::ends;
  // The str() method of the ostrstream class converts a (ostrstream) into a (string).
  histName = convert2strA.str();
  // The following obscure statement "thaws" the ostrstream so that it is deleted
  // when this method finishes.
  convert2strA.rdbuf()->freeze(0);

  // The c_str() method of the string class converts a (string) into a (const char *).
  m_eHitHist = 
    new TH1D(histName.c_str(),
	     histTitle = "z-distribution of energy deposited by hits",
	     numBins = 240,
	     zLow  = 0.,
	     zHigh = 600.);
  m_eHitHist->SetXTitle("z-position [cm]");
  m_eHitHist->SetYTitle("Energy [GeV]");

  ostrstream convert2strB;
  convert2strB << "Track" << a_run->GetRunID() << std::ends;
  histName = convert2strB.str();
  convert2strB.rdbuf()->freeze(0);
  m_eTrackHist = 
    new TH1D(histName.c_str(),
	     histTitle = "z-distribution of energy deposited in Aluminum",
	     numBins = 240,
	     zLow  = 0.,
	     zHigh = 600.);
  m_eTrackHist->SetXTitle("z-position [cm]");
  m_eTrackHist->SetYTitle("Energy [GeV]");

  ostrstream convert2strC;
  convert2strC << "Tot" << a_run->GetRunID() << std::ends;
  histName = convert2strC.str();
  convert2strC.rdbuf()->freeze(0);
  m_eTotHist = 
    new TH1D(histName.c_str(),
	     histTitle = "z-distribution of total energy deposited in the detector",
	     numBins = 240,
	     zLow  = 0.,
	     zHigh = 600.);
  m_eTotHist->SetXTitle("z-position [cm]");
  m_eTotHist->SetYTitle("Energy [GeV]");

  // Side note: If you use the ExN02DetectorConstruction class from
  // the original exampleN02, then all three of the above histograms
  // will be the same.  This is because the detector consists of
  // aluminum plates separated by air, and the "sensitive detector"
  // used to record hits are those same plates. Therefore, the only
  // energy deposited in the detector is in the plates in all three
  // histograms.
}


void ExR04Histograms::RecordEndOfRun(const G4Run* a_run)
{
  // Save the histogram, then delete it.
  m_histFile->Write();

  delete m_eHitHist;
  delete m_eTrackHist;
  delete m_eTotHist;
}


void ExR04Histograms::RecordEndOfEvent(const G4Event* a_event)
{  
  // Access the Hit Collection of this event
  // (from section 3.12.4.2 of the Geant4 Users' Guide).

  // First, get the pointer to the singleton Sensitive Detector Manager.
  G4SDManager* fSDM = G4SDManager::GetSDMpointer();

  // Get the collection ID number for the particular type of hit
  // we're interested in.  "trackerCollection" was defined to be the
  // name of a hit collection in ExN02TrackerSD.cc.  If there were
  // many types of hits (e.g., calorimeter hits, chamber hits, etc.)
  // then the following code could be repeated as needed for each
  // type.
  G4int collectionID = fSDM->GetCollectionID("trackerCollection");

  // Get the pointer to the hit collections of this event.
  G4HCofThisEvent* HCofEvent = a_event->GetHCofThisEvent();

  // Get the hit collection that corresponds to the particular collection ID number.
  // Note the type cast that converts the "G4VHitsCollection*" returned by GetHC
  // into a pointer to our custom type of hit collection.
  ExN02TrackerHitsCollection* trackerHC = 
    (ExN02TrackerHitsCollection*) (HCofEvent->GetHC(collectionID));

  // For this example, we want to loop through all the tracker hits.
  // Get the number of hits in this hit collection.
  G4int numberHits = trackerHC->entries();

  for (G4int i = 0; i < numberHits; i++)
    {
      // Get the "i-th" hit in the collection.
      ExN02TrackerHit* aHit = (*trackerHC)[i];

      // Again, note the use of ROOT data types as arguments to the ROOT method.
      // Also note the unit conversion before we fill the histogram.
      Stat_t edep = aHit->GetEdep() / GeV;

      // GetPos returns a G4ThreeVector; the "z()" method returns the z-component.
      Axis_t zpos = ( aHit->GetPos() ).z() / cm;

      // Add the hit to the histogram. 
      m_eHitHist->Fill(zpos,edep);
    }
}



void ExR04Histograms::RecordStep(const G4Step* a_step)
{
  // Energy deposited during this step.  Again (see RecordEndOfEvent
  // above) note the use of ROOT data types, and the unit conversion.
  Stat_t edep = a_step->GetTotalEnergyDeposit() / GeV;

  // z-position of this G4Step.
  Axis_t zpos = ( a_step->GetPreStepPoint()->GetPosition() ).z() / cm;

  // Add this step to the total energy histogram.
  m_eTotHist->Fill(zpos,edep);

  // Now histogram only the energy that's deposited in the tracker.

  string material = a_step->GetPreStepPoint()->GetMaterial()->GetName();
  // The tracker is made of aluminum
  if (material == "Al") 
    {
      m_eTrackHist->Fill(zpos,edep);
    }
}