#include "mgSteppingAction.hh"
#include "mgEventAction.hh"
#include "mgDetectorConstruction.hh"
#include "mgPrimaryGeneratorAction.hh"

#include "G4Step.hh"
#include "G4Event.hh"
#include "G4RunManager.hh"
#include "G4LogicalVolume.hh"
#include "G4UnitsTable.hh"
#include "G4SystemOfUnits.hh"

mgSteppingAction::mgSteppingAction(mgEventAction* eventAction) 
: G4UserSteppingAction(),
  fEventAction(eventAction),
  fShell1Volume(0), fShell2Volume(0), fShell3Volume(0)
{}

mgSteppingAction::~mgSteppingAction() 
{}

void mgSteppingAction::UserSteppingAction(const G4Step* step)
{

  //  G4cout << "Step Made: part 1" << G4endl;


  if (!fShell1Volume) {
    const mgDetectorConstruction* detectorConstruction
      = static_cast<const mgDetectorConstruction*>
      (G4RunManager::GetRunManager()->GetUserDetectorConstruction());
    fShell1Volume = detectorConstruction->GetShell1Volume();
    fShell2Volume = detectorConstruction->GetShell2Volume();
    fShell3Volume = detectorConstruction->GetShell3Volume();
    //    fShell0Volume = detectorConstruction->GetShell0Volume();

  }
  //  G4cout << "Step Made: part 2" << G4endl;

  G4String name1 = fShell1Volume->GetName();
  G4String name2 = fShell2Volume->GetName();
  G4String name3 = fShell3Volume->GetName();

  // G4cout << "Shell name: " << name << G4endl;

  G4LogicalVolume* volume = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetLogicalVolume();
  
  G4String testname = volume->GetName();
  //  G4cout << name1 << "  " << name2 << "  " << name3 << "  "<< testname <<G4endl;

  //  G4cout << "Step Made: volume defined" << G4endl;
  G4String particlename = step->GetTrack()->GetParticleDefinition()->GetParticleName();
  G4double px=0;
  G4double py=0;
  G4double pz=0;
  G4double xx=0;
  G4double xy=0;
  G4double xz=0;
  G4double kinEnergy=0;
  G4ThreeVector momentum = G4ThreeVector(0,0,0);
  G4ThreeVector position = G4ThreeVector(0,0,0);
  G4double radius = 0;

  if (particlename == "proton") {
    step->GetTrack()->SetTrackStatus(fStopAndKill);
    return;
  }

  if (particlename == "electron") {
    step->GetTrack()->SetTrackStatus(fStopAndKill);
    return;
  }

  if (particlename == "gamma") {
    step->GetTrack()->SetTrackStatus(fStopAndKill);
    return;
  }

  //check if neutron is entering the detector shell and print results is it is. also kill the neutron.
  if (testname == name1 && particlename=="neutron") {
    kinEnergy = step->GetTrack()->GetKineticEnergy();
    if (kinEnergy < 2*MeV) {
      step->GetTrack()->SetTrackStatus(fStopAndKill);
      return;
    }   
    momentum = step->GetTrack()->GetMomentumDirection();
    px = momentum.x();
    py = momentum.y();
    pz = momentum.z();
    position = step->GetTrack()->GetPosition();
    xx = position.x();
    xy = position.y();
    xz = position.z();
    radius = pow((pow(xx,2)+pow(xy,2)+pow(xz,2)),(0.5));
    
    G4cout << "Sd" << "  " << particlename << "  " << kinEnergy/MeV << "  " << px << "  " << py << "  " << pz << "  " << xx/cm << "  " << xy/cm << "  " << xz/cm <<"  "<< radius/cm<< G4endl;
    step->GetTrack()->SetTrackStatus(fStopAndKill);
    return;
  }

  //check if particle is a neutron in the second (buffer-water) shell and output a line if it is
  if ((testname==name2 && particlename=="neutron")) {
    kinEnergy = step->GetTrack()->GetKineticEnergy();
    if (kinEnergy < 2*MeV) {
      step->GetTrack()->SetTrackStatus(fStopAndKill);
      return;
    }
    momentum = step->GetTrack()->GetMomentumDirection();
    px = momentum.x();
    py = momentum.y();
    pz = momentum.z();
    position = step->GetTrack()->GetPosition();
    xx = position.x();
    xy = position.y();
    xz = position.z();
    
    G4cout << "Sb" << "  " << particlename << "  " << kinEnergy/MeV << "  " << px << "  " << py << "  " << pz << "  " << xx/cm << "  " << xy/cm << "  " << xz/cm << G4endl;
    return;
  }

  //check if particle is a neutron in the third (outer) shell and print results if it is.
  if ((testname==name3 && particlename=="neutron")) {
    kinEnergy = step->GetTrack()->GetKineticEnergy();
    if (kinEnergy < 2*MeV) {
      step->GetTrack()->SetTrackStatus(fStopAndKill);
      return;
    }    
    momentum = step->GetTrack()->GetMomentumDirection();
    px = momentum.x();
    py = momentum.y();
    pz = momentum.z();
    position = step->GetTrack()->GetPosition();
    xx = position.x();
    xy = position.y();
    xz = position.z();
    
    G4cout << "S3" << "  " << particlename << "  " << kinEnergy/MeV << "  " << px << "  " << py << "  " << pz << "  " << xx/cm << "  " << xy/cm << "  " << xz/cm << G4endl;
    step->GetTrack()->SetTrackStatus(fStopAndKill);
    return;
  }
}