//
// ********************************************************************
// * 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 "G4Scintillation.hh"
#include "G4Step.hh"
#include "G4Event.hh"
#include "G4RunManager.hh"
//#include "g4analysis.hh"
#include "G4AnalysisManager.hh"
#include "G4LogicalVolume.hh"
#include "G4Track.hh"
#include "G4ParticleMomentum.hh"
#include "G4SteppingManager.hh"
#include "G4ios.hh"
#include "G4SteppingManager.hh"
#include "G4Step.hh"
#include "G4StepPoint.hh"
#include "G4VPhysicalVolume.hh"
#include "G4ios.hh"

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

B1SteppingAction::B1SteppingAction(B1EventAction* eventAction, G4bool PrimaryPositronElectronExitFlag, G4bool SecondaryPhotonsExitFlag, G4bool SecondaryMuonsExitFlag, G4bool SecondaryPosEleExitFlag, G4bool SecondaryNeutrinosExitFlag, G4bool PrimaryPositronElectronEnteringFlag, G4bool RegeneratedPositronFlag, G4bool KillNoGammaFlag, G4bool ParticlePlastScintFlag, G4int TargExt, G4bool SecondaryPionExitFlag, G4bool SecondaryMuonsInDummyFlag, G4int DummyIn, G4bool ScintillatorShield)
: G4UserSteppingAction(),
  fEventAction(eventAction),
  fPrimaryPositronElectronExitFlag(PrimaryPositronElectronExitFlag),
  fSecondaryPhotonsExitFlag(SecondaryPhotonsExitFlag),
  fSecondaryMuonsExitFlag(SecondaryMuonsExitFlag),
  fSecondaryPosEleExitFlag(SecondaryPosEleExitFlag),
  fSecondaryNeutrinosExitFlag(SecondaryNeutrinosExitFlag),
  fPrimaryPositronElectronEnteringFlag(PrimaryPositronElectronEnteringFlag),
   fRegeneratedPositronFlag(RegeneratedPositronFlag),
   fKillNoGammaFlag(KillNoGammaFlag),
   fParticlePlastScintFlag(ParticlePlastScintFlag),
  fTargExt(TargExt),  
   fSecondaryPionExitFlag(SecondaryPionExitFlag), 
    fSecondaryMuonsInDummyFlag(SecondaryMuonsInDummyFlag),  
	fDummyIn(DummyIn),
  fScoringVolume(0),  
  fScoringVolume2(0),
   fScoringVolume3(0),
    fScoringVolume4(0),
     fScoringVolume5(0),
   fScoringVolumeReg(0),
   fScoringVolumeScint(0), 
   fScintillatorShield(ScintillatorShield)
{}

//....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();
    fScoringVolume3 = detectorConstruction->GetScoringVolume3();
    fScoringVolume4 = detectorConstruction->GetScoringVolume4();
    fScoringVolume5 = detectorConstruction->GetScoringVolume5();
    fScoringVolumeReg = detectorConstruction->GetScoringVolumeReg();
    fScoringVolumeScint = detectorConstruction->GetScoringVolumeScint();
}

	//Event
	const G4Event* evt = G4RunManager::GetRunManager()->GetCurrentEvent();

  // 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);
}

	if (volume == fScoringVolume3) {
  G4double edepStep3 = step->GetTotalEnergyDeposit();
  fEventAction->AddEdep3(edepStep3);
}

	if (volume == fScoringVolume4) {
  G4double edepStep4 = step->GetTotalEnergyDeposit();
  fEventAction->AddEdep4(edepStep4);
}

	if (volume == fScoringVolume5) {
  G4double edepStep5 = step->GetTotalEnergyDeposit();
  fEventAction->AddEdep5(edepStep5);
}

	if (volume == fScoringVolumeReg) {
  G4double edepStepReg = step->GetTotalEnergyDeposit();
  fEventAction->AddEdepReg(edepStepReg);
}
	if (volume == fScoringVolumeScint) {
  G4double edepStepScint = step->GetTotalEnergyDeposit();
  fEventAction->AddEdepScint(edepStepScint);
}

 //
 G4String TargetExt;
 if (fTargExt==1) TargetExt="Target1";
 else if (fTargExt==2) TargetExt="Target2";
 else if (fTargExt==3) TargetExt="Target3";
 else if (fTargExt==4) TargetExt="Target4";
 else if (fTargExt==5) TargetExt="Target5";
  G4String DummyInString;
 if (fDummyIn==1) DummyInString="Dummy1";
 else if (fDummyIn==2) DummyInString="Dummy2";
 
 // Retrieve the secondary particles
  	G4SteppingManager*  steppingManager = fpSteppingManager;
    fSecondary = steppingManager -> GetfSecondary();
    G4String secondaryParticleName;
    G4String process;
    G4int fSectrackID;
    G4int fSsecparID;
    G4String volumeName;
    
    //Definitions
    G4Track* track = step->GetTrack();
    
    	
    
//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 positrons-electrons
if (fSecondaryPosEleExitFlag){
	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 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);  	
		}
}

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

//get regenerated positrons
if (fRegeneratedPositronFlag){ 

if (fKillNoGammaFlag){ 
//Kill no photons entering in the regenerator
	if(NextVol && ThisVol->GetName()=="World" && NextVol->GetName()=="Regenerator" && step->GetTrack()->GetDynamicParticle()->GetPDGcode()!=22 ) { //if you want kill the no-gamma that could enter in the regenerator even-if you setted a magnetic field
	step->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
	//		G4cout<<"Particle Killed"<<G4endl;
}
	else if(NextVol && ThisVol->GetName()=="World" && NextVol->GetName()=="Regenerator" && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==22 ) { //get photon entering in the regenerator
			 		    int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(4,0, pdg);
			 			double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(4,1, kinEnergy);
			 			double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(4,2, MomDirx);
			 			double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(4,3, MomDiry);
			 			double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(4,4, MomDirz);
			 			double Vertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(4,5, Vertx);
						double Verty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(4,6, Verty);
						double Vertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(4,7, Vertz);
						double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(4,8, Mom);
						double Dirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(4,9, Dirxsp); 
			 			double Dirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(4,10, Dirysp); 
			 			double Dirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(4,11, Dirzsp); 
			 			double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(4,12, Momx);
						double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(4,13, Momy);
						double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(4,14, Momz);
						double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(4,15, xprime);
						double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(4,16, yprime);
						double ParID=step->GetTrack()->GetParentID();
						analysisManager->FillNtupleDColumn(4,17, ParID);
			   			analysisManager->AddNtupleRow(4);  	
			}
		}
	
	
 	else if (!fKillNoGammaFlag){ 
		if(NextVol && ThisVol->GetName()=="World" && NextVol->GetName()=="Regenerator"  ) { //get all particles entering in the regenerator
			 		    int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(3,0, pdg);
			 			double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(3,1, kinEnergy);
			 			double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(3,2, MomDirx);
			 			double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(3,3, MomDiry);
			 			double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(3,4, MomDirz);
			 			double Vertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(3,5, Vertx);
						double Verty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(3,6, Verty);
						double Vertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(3,7, Vertz);
						double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(3,8, Mom);
						double Dirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(3,9, Dirxsp); 
			 			double Dirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(3,10, Dirysp); 
			 			double Dirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(3,11, Dirzsp); 
			 			double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(3,12, Momx);
						double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(3,13, Momy);
						double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(3,14, Momz);
						double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(3,15, xprime);
						double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(3,16, yprime);
						double ParID=step->GetTrack()->GetParentID();
						analysisManager->FillNtupleDColumn(3,17, ParID);
			   			analysisManager->AddNtupleRow(3);  	
		}
		
		if(NextVol && ThisVol->GetName()=="World" && NextVol->GetName()=="Regenerator" && step->GetTrack()->GetParentID()!=00 && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==22 ) { //get only photons entering in the regenerator
			 		    int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(4,0, pdg);
			 			double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(4,1, kinEnergy);
			 			double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(4,2, MomDirx);
			 			double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(4,3, MomDiry);
			 			double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(4,4, MomDirz);
			 			double Vertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(4,5, Vertx);
						double Verty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(4,6, Verty);
						double Vertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(4,7, Vertz);
						double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(4,8, Mom);
						double Dirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(4,9, Dirxsp); 
			 			double Dirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(4,10, Dirysp); 
			 			double Dirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(4,11, Dirzsp); 
			 			double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(4,12, Momx);
						double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(4,13, Momy);
						double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(4,14, Momz);
						double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(4,15, xprime);
						double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(4,16, yprime);
						double ParID=step->GetTrack()->GetParentID();
						analysisManager->FillNtupleDColumn(4,17, ParID);
			   			analysisManager->AddNtupleRow(4);  	
		}
	
	}
		
		 for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++){
     secondaryParticleName =  (*fSecondary)[lp1]->GetDefinition() -> GetParticleName();  // name of the secondary
     process = (*fSecondary)[lp1]-> GetCreatorProcess()-> GetProcessName();   // process creating it
     fSectrackID = (*fSecondary)[lp1]->GetTrackID(); //track id
     fSsecparID = (*fSecondary)[lp1]->GetParentID();   //parent id
     volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName(); // volume where secondary was generated 
    // if (process=="conv" && secondaryParticleName =="e+" && volumeName == "Regenerator" && ThisVol->GetName()=="Regenerator"){

    //G4cerr << G4endl;
     //G4cerr << "Event n°: " << lp1 << G4endl;
     //G4cerr << "secondaryParticleName: " << (*fSecondary)[lp1]->GetDefinition() -> GetParticleName() << G4endl;
     //G4cerr << "process: " << (*fSecondary)[lp1]-> GetCreatorProcess()-> GetProcessName() << G4endl;
     //G4cerr << "fSectrackID: " << (*fSecondary)[lp1]->GetTrackID() << G4endl;
     //G4cerr << "fSsecparID: " << (*fSecondary)[lp1]->GetParentID() << G4endl;
     //G4cerr << "volumeName: " << (*fSecondary)[lp1] -> GetVolume() -> GetName() << G4endl;
 //    G4cerr << G4endl;
 	}
	 	
	if(NextVol && ThisVol->GetName()=="Regenerator" && NextVol->GetName()=="World" && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-11 && volumeName == "Regenerator" ) { //get positrons outgoing the regenerato
					//    G4cout << "Process " << process << G4endl; 
						int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(5,0, pdg);
			 			double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(5,1, kinEnergy);
			 			double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(5,2, MomDirx);
			 			double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(5,3, MomDiry);
			 			double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(5,4, MomDirz);
			 			double Vertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(5,5, Vertx);
						double Verty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(5,6, Verty);
						double Vertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(5,7, Vertz);
						double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(5,8, Mom);
						double Dirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(5,9, Dirxsp); 
			 			double Dirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(5,10, Dirysp); 
			 			double Dirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(5,11, Dirzsp); 
			 			double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(5,12, Momx);
						double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(5,13, Momy);
						double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(5,14, Momz);
						double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(5,15, xprime);
						double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(5,16, yprime);
						double ParID=step->GetTrack()->GetParentID();
						analysisManager->FillNtupleDColumn(5,17, ParID);
			   			analysisManager->AddNtupleRow(5);  	
			}
		}
	
	//get particles hitting the scintillator
if (fParticlePlastScintFlag){
	if (!fScintillatorShield){
	if(NextVol && ThisVol->GetName()=="World" && NextVol->GetName()=="PlasticScintillator") { 
						int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(6,0, pdg);
			 			double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(6,1, kinEnergy);
			 			double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(6,2, MomDirx);
			 			double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(6,3, MomDiry);
			 			double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(6,4, MomDirz);
			 			double Vertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(6,5, Vertx);
						double Verty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(6,6, Verty);
						double Vertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(6,7, Vertz);
						double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(6,8, Mom);
						double Dirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(6,9, Dirxsp); 
			 			double Dirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(6,10, Dirysp); 
			 			double Dirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(6,11, Dirzsp); 
			 			double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(6,12, Momx);
						double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(6,13, Momy);
						double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(6,14, Momz);
						double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(6,15, xprime);
						double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(6,16, yprime);
			   			analysisManager->AddNtupleRow(6);  
					}
		}
			else if (fScintillatorShield){
	if(NextVol && ThisVol->GetName()=="VacuumScintShield" && NextVol->GetName()=="PlasticScintillator") { 
						 int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(6,0, pdg);
			 			double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(6,1, kinEnergy);
			 			double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(6,2, MomDirx);
			 			double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(6,3, MomDiry);
			 			double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(6,4, MomDirz);
			 			double Vertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(6,5, Vertx);
						double Verty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(6,6, Verty);
						double Vertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(6,7, Vertz);
						double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(6,8, Mom);
						double Dirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(6,9, Dirxsp); 
			 			double Dirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(6,10, Dirysp); 
			 			double Dirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(6,11, Dirzsp); 
			 			double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(6,12, Momx);
						double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(6,13, Momy);
						double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(6,14, Momz);
						double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(6,15, xprime);
						double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(6,16, yprime);
			   			analysisManager->AddNtupleRow(6);  
					}
		}
}

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


	
		/*auto proc_man = step->GetTrack()->GetParticleDefinition()->GetProcessManager();
						G4String proc_name = step->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();
					//	 G4cout << "Process " << proc_name <<  G4endl;
    					G4int n_scint = 0;
						 if(proc_name == "Scintillation"){
						 	G4Scintillation* scint = new G4Scintillation("Scintillation");
						 	n_scint    = scint->GetNumPhotons();
					       	analysisManager->FillNtupleDColumn(8,0, n_scint);
					       		analysisManager->AddNtupleRow(8);  	
				            }*/
				            
				            
				G4int photons = 0;
				G4SteppingManager* fpSteppingManager = 
  				G4EventManager::GetEventManager()->GetTrackingManager()->GetSteppingManager();
				G4StepStatus stepStatus = fpSteppingManager->GetfStepStatus();
				if (stepStatus != fAtRestDoItProc) {
  					G4ProcessVector* procPost = fpSteppingManager->GetfPostStepDoItVector(); 
  					size_t MAXofPostStepLoops = fpSteppingManager->GetMAXofPostStepLoops(); 
  					for (size_t i3 = 0; i3 < MAXofPostStepLoops; i3++) {
    					if ((*procPost)[i3]->GetProcessName() == "Scintillation") {
      						G4Scintillation* proc1 = (G4Scintillation*) (*procPost)[i3];
      						photons += proc1->GetNumPhotons();
    					}
  					}	
				}	
      			if(photons > 0){ 
			  		G4cout << "In this step, "<< photons << " scintillation photons were produced." << G4endl;
			  			analysisManager->FillNtupleDColumn(8,0, photons);
					       		analysisManager->AddNtupleRow(8);  
      			}
      			
      			
    if(NextVol && ThisVol->GetName()=="PlasticScintillator" && NextVol->GetName()=="PlasticScintillator") { 
		int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
		analysisManager->FillNtupleDColumn(9,0, pdg);
		int trackID = step->GetTrack()->GetTrackID();
		analysisManager->FillNtupleDColumn(9,1, trackID);
		int evtID = evt->GetEventID();
		analysisManager->FillNtupleDColumn(9,2, evtID);
		double partEdep = step->GetTotalEnergyDeposit();
		analysisManager->FillNtupleDColumn(9,3, partEdep);
		analysisManager->AddNtupleRow(9);  
	}
		

  if (volume != fScoringVolume || volume != fScoringVolume2 || volume != fScoringVolume3 || volume != fScoringVolume4 || volume != fScoringVolume5 || volume !=fScoringVolumeReg || volume !=fScoringVolumeScint) return;
  
}

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