//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
/// \file B1SteppingAction.cc
/// \brief Implementation of the B1SteppingAction class

#include "B1SteppingAction.hh"
#include "B1EventAction.hh"
#include "B1DetectorConstruction.hh"
#include "G4DynamicParticle.hh"
#include "G4Step.hh"
#include "G4Event.hh"
#include "G4RunManager.hh"
#include "g4analysis.hh"
#include "G4LogicalVolume.hh"
#include "G4Track.hh"
#include "G4ParticleMomentum.hh"

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

B1SteppingAction::B1SteppingAction(B1EventAction* eventAction, G4bool PrimaryPositronElectronExitFlag, G4bool SecondaryPhotonsExitFlag, G4bool SecondaryMuonsExitFlag, G4int TargExt)
: G4UserSteppingAction(),
  fEventAction(eventAction),
  fPrimaryPositronElectronExitFlag(PrimaryPositronElectronExitFlag),
  fSecondaryPhotonsExitFlag(SecondaryPhotonsExitFlag),
  fSecondaryMuonsExitFlag(SecondaryMuonsExitFlag),
  fTargExt(TargExt),  
  fScoringVolume(0),  
  fScoringVolume2(0)
{}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

B1SteppingAction::~B1SteppingAction()
{}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void B1SteppingAction::UserSteppingAction(const G4Step* step)
{
	
   auto analysisManager = G4AnalysisManager::Instance();
  if (!fScoringVolume) { 
    const B1DetectorConstruction* detectorConstruction
      = static_cast<const B1DetectorConstruction*>
        (G4RunManager::GetRunManager()->GetUserDetectorConstruction());
    fScoringVolume = detectorConstruction->GetScoringVolume();   
    fScoringVolume2 = detectorConstruction->GetScoringVolume2();
}

  // get volume of the current step
  G4LogicalVolume* volume 
    = step->GetPreStepPoint()->GetTouchableHandle()
      ->GetVolume()->GetLogicalVolume();
      
      G4VPhysicalVolume* ThisVol = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
	G4VPhysicalVolume* NextVol = step->GetPostStepPoint()->GetTouchableHandle()->GetVolume();
  
	if (volume == fScoringVolume) {
  G4double edepStep = step->GetTotalEnergyDeposit();
  fEventAction->AddEdep(edepStep);
	}
	
	if (volume == fScoringVolume2) {
  G4double edepStep2 = step->GetTotalEnergyDeposit();
  fEventAction->AddEdep2(edepStep2);
}

 //
 G4String TargetExt;
 if (fTargExt==1) TargetExt="Target1";
 else if (fTargExt==2) TargetExt="Target2";
 
 
//get primary particles
if (fPrimaryPositronElectronExitFlag){
	if(NextVol && ThisVol->GetName()==TargetExt && NextVol->GetName()=="World" && step->GetTrack()->GetParentID()==00 && (step->GetTrack()->GetDynamicParticle()->GetPDGcode()==11 || step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-11)) { 
			 		    int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(1,0, pdg);
			 			double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(1,1, kinEnergy);
			 			double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(1,2, MomDirx);
			 			double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(1,3, MomDiry);
			 			double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(1,4, MomDirz);
			 			double Vertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(1,5, Vertx);
						double Verty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(1,6, Verty);
						double Vertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(1,7, Vertz);
						double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(1,8, Mom);
						double Dirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(1,9, Dirxsp); 
			 			double Dirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(1,10, Dirysp); 
			 			double Dirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(1,11, Dirzsp); 
			 			double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(1,12, Momx);
						double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(1,13, Momy);
						double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(1,14, Momz);
						double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(1,15, xprime);
						double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(1,16, yprime);
			   			analysisManager->AddNtupleRow(1);  	
		}
}

//get secondary muons
if (fSecondaryMuonsExitFlag){
	if(NextVol && ThisVol->GetName()==TargetExt && NextVol->GetName()=="World" && step->GetTrack()->GetParentID()!=00 && (step->GetTrack()->GetDynamicParticle()->GetPDGcode()==13 || step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-13)) { 
						 int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(1,0, pdg);
			 			double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(1,1, kinEnergy);
			 			double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(1,2, MomDirx);
			 			double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(1,3, MomDiry);
			 			double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(1,4, MomDirz);
			 			double Vertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(1,5, Vertx);
						double Verty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(1,6, Verty);
						double Vertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(1,7, Vertz);
						double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(1,8, Mom);
						double Dirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(1,9, Dirxsp); 
			 			double Dirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(1,10, Dirysp); 
			 			double Dirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(1,11, Dirzsp); 
			 			double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(1,12, Momx);
						double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(1,13, Momy);
						double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(1,14, Momz);
						double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(1,15, xprime);
						double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(1,16, yprime);
			   			analysisManager->AddNtupleRow(1);  	
		}
}

//get secondary photons
if (fSecondaryPhotonsExitFlag){
		if(!step->GetTrack()->GetNextVolume() && step->GetTrack()->GetParentID()!=00 && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==22 ) { 
			 		    int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(1,0, pdg);
			 			double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(1,1, kinEnergy);
			 			double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(1,2, MomDirx);
			 			double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(1,3, MomDiry);
			 			double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(1,4, MomDirz);
			 			double Vertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(1,5, Vertx);
						double Verty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(1,6, Verty);
						double Vertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(1,7, Vertz);
						double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(1,8, Mom);
						double Dirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(1,9, Dirxsp); 
			 			double Dirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(1,10, Dirysp); 
			 			double Dirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(1,11, Dirzsp); 
			 			double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(1,12, Momx);
						double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(1,13, Momy);
						double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(1,14, Momz);
						double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(1,15, xprime);
						double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(1,16, yprime);
			   			analysisManager->AddNtupleRow(1);  	
		}
}
		

  if (volume != fScoringVolume || volume != fScoringVolume2) return;
  
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......