#include "SteppingAction.hh"
#include "EventAction.hh"
#include "DetectorConstruction.hh"
#include "PrimaryGeneratorAction.hh"
#include "RunAction.hh"
#include "HistoManager.hh"
#include "G4GeneralParticleSource.hh"
#include "G4VPhysicalVolume.hh"
#include "G4ParticleTypes.hh"
#include "G4Step.hh"
#include "G4Event.hh" 
#include "G4RunManager.hh"
#include "G4LogicalVolume.hh"
#include "g4root.hh"
#include "G4SteppingManager.hh"
#include "G4StepStatus.hh"

SteppingAction::SteppingAction(RunAction* runaction, EventAction* eventaction, PrimaryGeneratorAction* primary, HistoManager* histo, DetectorConstruction* det): G4UserSteppingAction(), fRunAction(runaction), fEventAction(eventaction), fPrimary(primary), fHistoManager(histo),fDetector(det),
energy_mu(0.), dirx_mu(0.), diry_mu(0.), dirz_mu(0.), pX_mu(0.), pY_mu(0.), pZ_mu(0.), p_mu(0.), 
posX_mu(0.), posY_mu(0.), posZ_mu(0.), Eta_mu(0.), Phi_mu(0.), Theta_mu(0.), PDG_mu(0.), 
Charge_mu(0.), Spin_mu(0.), Mass_mu(0.), 
energy_n(0.), dirx_n(0.), diry_n(0.), dirz_n(0.), pX_n(0.), pY_n(0.), pZ_n(0.), p_n(0.), 
posX_n(0.), posY_n(0.), posZ_n(0.), Eta_n(0.), Phi_n(0.), Theta_n(0.), PDG_n(0.), 
Charge_n(0.), Spin_n(0.), Mass_n(0.)
{ }

SteppingAction::~SteppingAction() { }

void SteppingAction::BeginOfSteppingAction(const G4Step*)
{
  energy_mu = 0;
  dirx_mu = 0;
  diry_mu = 0;
  dirz_mu = 0;
  pX_mu = 0;
  pY_mu = 0;
  pZ_mu = 0;
  p_mu = 0; 
  posX_mu = 0;
  posY_mu = 0;
  posZ_mu = 0;
  Eta_mu = 0;
  Phi_mu = 0;
  Theta_mu = 0;
  PDG_mu = 0; 
  Charge_mu = 0;
  Spin_mu = 0;
  Mass_mu = 0;

  energy_n = 0;
  dirx_n = 0;
  diry_n = 0;
  dirz_n = 0;
  pX_n = 0;
  pY_n = 0;
  pZ_n = 0;
  p_n = 0; 
  posX_n = 0;
  posY_n = 0;
  posZ_n = 0;
  Eta_n = 0;
  Phi_n = 0;
  Theta_n = 0;
  PDG_n = 0; 
  Charge_n = 0;
  Spin_n = 0;
  Mass_n = 0;

}

void SteppingAction::UserSteppingAction(const G4Step* step)
{
  G4int trackID = step->GetTrack()->GetTrackID();
  G4double Edep = step->GetTotalEnergyDeposit(); 
  if(Edep > 0.) fEventAction->SumEnergyDeposited(trackID, Edep);

  if (trackID > 1) return;

  //count processes
  
  const G4StepPoint* prePoint = step->GetPreStepPoint();
  const G4StepPoint* endPoint = step->GetPostStepPoint();
  const G4VProcess* process = endPoint->GetProcessDefinedStep();
  G4String procName = process->GetProcessName();
  G4int subtype = process->GetProcessSubType();
  G4int nbsec = step->GetNumberOfSecondariesInCurrentStep();
  if ((subtype==2)&&(nbsec==0)) procName = "Edep alone";
  fRunAction->CountProcesses(procName);
  
  //step size and track length

  G4double stepSize = step->GetStepLength();
  fRunAction->TrackLength(stepSize);
  fHistoManager->FillHisto(1,stepSize);

  if (nbsec == 0) return;   //no secondary particles

  //energy transferred to secondary particles

  const std::vector<const G4Track*>* secondaries = step->GetSecondaryInCurrentStep();

  G4double Etransfer = 0.;
  for (G4int i=0; i<nbsec; i++)
    {
      const G4Track* track =(*secondaries)[i];
      const G4ParticleDefinition* particle = track->GetParticleDefinition();
      G4String name = particle->GetParticleName();
      G4double energy = track->GetKineticEnergy();
      fRunAction->EnergySpectrumOfSecondaries(name, energy);
      G4int ih = 0;
      if (particle == G4Gamma::Gamma()) ih = 11;
      else if (particle == G4Electron::Electron()) ih = 12;
      //else if (particle == G4Positron::Positron()) ih = 13;
      if (ih > 0) fHistoManager->FillHisto(ih, energy);
      if(subtype == 4) energy = track->GetTotalEnergy();
      Etransfer += energy;
    }

  fEventAction->SumEnergyTransfered(process, Etransfer);

  G4VPhysicalVolume* volume = prePoint->GetTouchableHandle()->GetVolume();  

  //G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
  G4ParticleDefinition* prim_particle = fPrimary->GetParticleGun()->GetParticleDefinition();
  G4String partName = prim_particle->GetParticleName();
  //G4StepPoint* point1 = step->GetPreStepPoint();
  //G4StepPoint* boundarypoint = step->GetPostStepPoint();
  //G4TouchableHandle touch1 = point1->GetTouchableHandle();
  //G4VPhysicalVolume* volume = touch1->GetVolume();
  //G4TouchableHandle touch = boundarypoint->GetTouchableHandle();
  //G4VPhysicalVolume* volume = touch->GetVolume();
  //G4String name = volume->GetName();
  //G4cout << "Volume name : "<< name << G4endl;
  //G4bool condition = step->IsLastStepInVolume();
  //G4cout << "condition is " << condition << G4endl;

  //if ((prim_particle == G4MuonMinus::MuonMinus()) /*&& (boundarypoint->GetStepStatus() == fGeomBoundary) && condition == true step->IsLastStepInVolume()*/){

  if (endPoint->GetStepStatus() == fGeomBoundary) {
  //if (step->IsLastStepInVolume()) {
        G4double length = fDetector->GetSize();
	G4cout << "Size of geometry is " << length << G4endl;
	G4cout << "Step is the last step in the volume " << volume->GetName()<< "\n" <<  G4endl;
	G4Material* nextMaterial = endPoint->GetMaterial();
	G4cout << "Next Material is " << nextMaterial << G4endl;
	energy_mu =  endPoint->GetKineticEnergy();
	G4cout << "Muon energy on exit " << energy_mu << G4endl; 
	G4ThreeVector direction = endPoint->GetMomentumDirection();
	dirx_mu = direction.x();
	diry_mu = direction.y();
	dirz_mu = direction.z();

	//G4double momX_mu = step->GetDeltaMomentum().x();
	//G4double momY_mu = step->GetDeltaMomentum().y();
	//G4double momZ_mu = step->GetDeltaMomentum().z();

	pX_mu = endPoint->GetMomentum().x();
	pY_mu = endPoint->GetMomentum().y();
	pZ_mu = endPoint->GetMomentum().z();
        p_mu = sqrt(pX_mu*pX_mu + pY_mu*pY_mu + pZ_mu*pZ_mu);

	posX_mu = endPoint->GetPosition().x();
	posY_mu = endPoint->GetPosition().y();
	posZ_mu = endPoint->GetPosition().z();
	Eta_mu = endPoint->GetPosition().eta();
        Phi_mu = endPoint->GetPosition().phi();
	Theta_mu = endPoint->GetPosition().theta();
	PDG_mu = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
	Charge_mu = step->GetTrack()->GetParticleDefinition()->GetPDGCharge();
	Spin_mu = step->GetTrack()->GetParticleDefinition()->GetPDGSpin();
	Mass_mu = step->GetTrack()->GetParticleDefinition()->GetPDGMass();
	//G4cout << "energy is "<<  energy_mu << G4endl;

	//G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
	fHistoManager->FillNtuple_step_muon(energy_mu, dirx_mu, diry_mu, dirz_mu, 
		pX_mu, pY_mu, pZ_mu, p_mu, posX_mu, posY_mu, posZ_mu, Eta_mu, Phi_mu, 
		Theta_mu, PDG_mu, Charge_mu, Spin_mu, Mass_mu);
        }

	for (G4int i=0; i<nbsec; i++)
        {
	      const G4Track* track =(*secondaries)[i];
	      const G4ParticleDefinition* particle = track->GetParticleDefinition();
	      G4String particle_name = particle->GetParticleName();
	  if (particle == G4Neutron::Neutron() /* && step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary*/) 			{

		//G4cout << "fGeomBoundary in volume " << volume->GetName() << "\n " << G4endl;
		energy_n = endPoint->GetKineticEnergy(); 
		G4ThreeVector direction1 = endPoint->GetMomentumDirection();
		dirx_n = direction1.x();
		diry_n = direction1.y();
		dirz_n = direction1.z();

		//G4double momX_n = step->GetDeltaMomentum().x();
		//G4double momY_n = step->GetDeltaMomentum().y();
		//G4double momZ_n = step->GetDeltaMomentum().z();

		pX_n = endPoint->GetMomentum().x();
		pY_n = endPoint->GetMomentum().y();
		pZ_n = endPoint->GetMomentum().z();
                p_n = sqrt(pX_n*pX_n+pY_n*pY_n+pZ_n*pZ_n);
		posX_n = endPoint->GetPosition().x();
		posY_n = endPoint->GetPosition().y();
		posZ_n = endPoint->GetPosition().z();
		Eta_n = endPoint->GetPosition().eta();
		Phi_n = endPoint->GetPosition().phi();
		Theta_n = endPoint->GetPosition().theta();
		PDG_n = track->GetParticleDefinition()->GetPDGEncoding();
		Charge_n = track->GetParticleDefinition()->GetPDGCharge();
		Spin_n = track->GetParticleDefinition()->GetPDGSpin();
		Mass_n = track->GetParticleDefinition()->GetPDGMass();

	    	//G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
		fHistoManager->FillNtuple_step_neutron(energy_n, dirx_n, diry_n, dirz_n, 
		pX_n, pY_n, pZ_n, p_n, posX_n, posY_n, posZ_n, Eta_n, Phi_n, 
		Theta_n, PDG_n, Charge_n, Spin_n, Mass_n);
	}

}

}