#include "SteppingAction.hh"
#include "EventAction.hh"
#include "G4Step.hh"
#include "G4StepPoint.hh"
#include "DetectorConstruction.hh"
#include "RunAction.hh"
#include "Randomize.hh"
#include "G4RunManager.hh"
#include "G4SteppingManager.hh"
#include "G4EmSaturation.hh"
#include "G4LossTableManager.hh"
#include "g4root.hh"
#include "G4SystemOfUnits.hh"
#include <map>
#include <iostream>

using namespace std;

SteppingAction::SteppingAction(DetectorConstruction* det, EventAction* evt)
  :G4UserSteppingAction(),fDetector(det), eventaction(evt)
{ 
  filename = "protonKE.out"; //kinetic energy info output file
  filename2 = "secondariesKE.out"; //production neutron kinetic energy info output file
  
  componentMap.insert( std::pair<std::string,double>("Collimator1", 60.));
  componentMap.insert( std::pair<std::string,double>("ScatteringFoil1",80.025));
  componentMap.insert( std::pair<std::string,double>("BrassStopper",306.6));
  componentMap.insert( std::pair<std::string,double>("ScatteringFoil2",306.625));
  componentMap.insert( std::pair<std::string,double>("KaptonWindow",356.05));
  componentMap.insert( std::pair<std::string,double>("ScatterCollimator1",1135.));  
  componentMap.insert( std::pair<std::string,double>("NozzleBefore",1692.));
  componentMap.insert( std::pair<std::string,double>("NozzleAfter",1760.)); 
  componentMap.insert( std::pair<std::string,double>("WaterSurface",1840.));

}

SteppingAction::~SteppingAction()
{ }

void SteppingAction::UserSteppingAction(const G4Step* step)
{
  // G4AnalysisManager* man = G4AnalysisManager::Instance();
  
  G4VPhysicalVolume* volume = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
  G4String particleName = step->GetTrack()->GetDefinition()->GetParticleName();
  
  //filter for protons
  if(step->GetTrack()->GetDefinition()->GetParticleName()== "proton"){ 
  
    // Get z position of proton and add posShift to find distance from particle source (zpos)
    double zposRef = step->GetPreStepPoint()->GetPosition().z();
    double posShift = 4200.000; //distance from centre of room to source (mm)
    double zpos = zposRef + posShift; //obtain a value for z position with the proton source at 0mm

    //retrieve position value for a particular component from map
    double cPos = (double)componentMap.find("NozzleAfter")->second;
   
    
 //if step is at a certain z position write proton kinetic energy to file
     if(fabs(zpos-cPos) < 1.0e-8){
    // if(fabs(zpos-1770) < 1.0e-8){
         std::ofstream hitsfile(filename, std::ios::app);
      if(hitsfile.is_open()){
	G4double ke =  step->GetPreStepPoint()->GetKineticEnergy();
	G4double xPos = step->GetPreStepPoint()->GetPosition().x();
	G4double yPos = step->GetPreStepPoint()->GetPosition().y();

	hitsfile << xPos << " "
		 << yPos << " "
		 << zpos << " "
		 << ke <<G4endl;
	hitsfile.close();
      }
     }
  }
  // record kinetic energy of neutrons produced in a specific volume
  if(volume->GetName()=="WaterBox"){ // modify this line to change volume to observe neutron production
    if(particleName == "neutron" && step->GetTrack()->GetCurrentStepNumber()==1){
      double productionNeutronEnergy = step->GetTrack()->GetVertexKineticEnergy();
      std::ofstream hitsfile(filename2, std::ios::app);
      if(hitsfile.is_open()){
	hitsfile << productionNeutronEnergy <<G4endl;
	hitsfile.close();
      }
    }

  }
  
}