//
// This is to allow profiling of the energy deposition 
// in the iron abosrbers. It obviously is not a detector.
// Infomation is written directly to ROOT histograms.
//
// For this to get called, the appropriate lines in
// CPPDetectorConstruction.cc must be uncommented.
//


#ifndef _CPPSensitiveDetectorIronAbsorber_
#define _CPPSensitiveDetectorIronAbsorber_

#include <map>
#include <mutex>
using namespace std;

#include <TFile.h>
#include <TH2.h>

#include "G4VSensitiveDetector.hh"
#include "G4SystemOfUnits.hh"
#include "G4VProcess.hh"


class CPPSensitiveDetectorIronAbsorber: public G4VSensitiveDetector
{
public:
	CPPSensitiveDetectorIronAbsorber(const G4String& name):G4VSensitiveDetector(name){

		lock_guard<mutex> lck(hmutex);	
		if(!rootfile           )           rootfile = new TFile("CPPsim_IronAbsorber.root", "RECREATE");
		if(!dE_vs_z            )            dE_vs_z = new TH1D("dE_vs_z" ,           "Energy deposited in Iron Abosrber", 9000, 900.0, 1200.0);
		if(!dE_vs_xy           )           dE_vs_xy = new TH2D("dE_vs_xy",           "Energy deposited in Iron Abosrber", 170, -85.0, 85.0, 170, -85.0, 85.0);
		if(!Ncharged_vs_z      )      Ncharged_vs_z = new TH1D("Ncharged_vs_z",      "Number of charged particles in Iron Abosrber", 3000, 900.0, 1200.0);
		if(!KE_leaving_vs_z    )   KE_leaving_vs_z  = new TH1D("KE_leaving_vs_z",    "Kinetic energy of particles leaving the Iron Abosrber", 3000, 900.0, 1200.0);
		if(!z_initial          )          z_initial = new TH1D("z_initial",          "First interaction point of pion in the Iron Abosrber", 3000, 900.0, 1200.0);;
		if(!Ncharged_vs_deltaz ) Ncharged_vs_deltaz = new TH1D("Ncharged_vs_deltaz", "Number of charged particles in Iron Abosrber with z relative to first interaction point", 3000, -100.0, 200.0);;
		if(!dE_vs_deltaz       )       dE_vs_deltaz = new TH1D("dE_vs_deltaz",       "Energy deposited in the Iron Abosrber with z relative to first interaction point", 3000, -100.0, 200.0);;
	}
	
	virtual ~CPPSensitiveDetectorIronAbsorber(){
		
		lock_guard<mutex> lck(hmutex);
		if(rootfile){
			rootfile->Write();
			rootfile->Close();
			delete rootfile;
			rootfile = NULL;
		}
	}
	
	
	//------------------------------
	// Initialize
	void Initialize(G4HCofThisEvent* hitCollection)
	{
		if(hitCollection == (void*)0x1) cout<< "Huh?"; // avoid compiler warning
		
		z_first_interaction = -1000.0;
	}
	
	//------------------------------
	// ProcessHits
	G4bool ProcessHits(G4Step* step, G4TouchableHistory* history)
	{
		if(history == (void*)0x1) cout<< "Huh?"; // avoid compiler warning
		
		G4double dEsum = step->GetTotalEnergyDeposit() / MeV;
		
//		if (dEsum==0.) return false;
		
		// Get pre-step and post-step points
		G4StepPoint *sp_pre = step->GetPreStepPoint();
		G4StepPoint *sp_post = step->GetPostStepPoint();
		const G4TouchableHandle &hand = sp_pre->GetTouchableHandle();

		bool is_first_interaction_point = false;
		if(z_first_interaction < 0.0){
			const G4VProcess *proc = sp_post->GetProcessDefinedStep();
			if(proc->GetProcessType() > 2) {
				z_first_interaction = sp_post->GetPosition().z()/cm;
				is_first_interaction_point = true;
			}
		}

		// Consider point where the particle is at as the mean between the
		// pre-step and post-step		
		G4ThreeVector pos = 0.5*(sp_pre->GetPosition() + sp_post->GetPosition());
		G4ThreeVector mom = 0.5*(sp_pre->GetMomentum() + sp_post->GetMomentum());
		//G4double t =  0.5*(sp_pre->GetGlobalTime() + sp_post->GetGlobalTime()) / ns;

		// To count charged particles as a function of z without double counting
		// look to see if step spansa 1mm boundary in z.
		double z_pre = sp_pre->GetPosition().z()/mm;
		double z_post = sp_post->GetPosition().z()/mm;
		bool crossed_boundary = floor(z_pre) != floor(z_post);
		double q = step->GetTrack()->GetDynamicParticle()->GetCharge();

		lock_guard<mutex> lck(hmutex);
		if(dEsum>0.0) dE_vs_z->Fill( pos.z()/cm , dEsum);
		if(dEsum>0.0) dE_vs_xy->Fill( pos.x()/cm, pos.y()/cm , dEsum);
		if(q!=0.0 && crossed_boundary) Ncharged_vs_z->Fill(pos.z()/cm);
		if (sp_post->GetStepStatus() == fGeomBoundary) KE_leaving_vs_z->Fill(sp_post->GetPosition().z()/cm, sp_post->GetKineticEnergy());
		if(is_first_interaction_point) z_initial->Fill(z_first_interaction);
		if(z_first_interaction>=0.0){
			if(q!=0.0 && crossed_boundary) Ncharged_vs_deltaz->Fill(pos.z()/cm - z_first_interaction);
			if(dEsum>0.0) dE_vs_deltaz->Fill( pos.z()/cm - z_first_interaction , dEsum );
		}
		return false;
	}
	
	//------------------------------
	// EndOfEvent
	void EndOfEvent(G4HCofThisEvent* hitCollection){}
	
private:

	double z_first_interaction;

	static mutex hmutex;
	static TFile *rootfile;
	static TH1D *dE_vs_z;
	static TH2D *dE_vs_xy;
	static TH1D *Ncharged_vs_z;
	static TH1D *KE_leaving_vs_z;
	static TH1D *z_initial;
	static TH1D *Ncharged_vs_deltaz;
	static TH1D *dE_vs_deltaz;

};

mutex CPPSensitiveDetectorIronAbsorber::hmutex;
TFile *CPPSensitiveDetectorIronAbsorber::rootfile  = NULL;
TH1D* CPPSensitiveDetectorIronAbsorber::dE_vs_z  = NULL;
TH2D* CPPSensitiveDetectorIronAbsorber::dE_vs_xy = NULL;
TH1D* CPPSensitiveDetectorIronAbsorber::Ncharged_vs_z  = NULL;
TH1D* CPPSensitiveDetectorIronAbsorber::KE_leaving_vs_z  = NULL;
TH1D* CPPSensitiveDetectorIronAbsorber::z_initial = NULL;
TH1D* CPPSensitiveDetectorIronAbsorber::Ncharged_vs_deltaz = NULL;
TH1D* CPPSensitiveDetectorIronAbsorber::dE_vs_deltaz = NULL;


#endif // _CPPSensitiveDetectorIronAbsorber_