Dear all, I am having difficulties calculating the primary particles in my structure. My code represents a molecular cloud with 13 concentric layers, where I calculate the incident energy, the angle of incidence, and the energies deposited in each layer. However, I need help calculating the distance that the primary particles have traveled (in parsecs). When I try to perform the calculation in SteppingAction.cc, I encounter interference with the calculation of the deposited energy. Please take a look:
EventAction.cc
#include "EventAction.hh"
#include "Run.hh"
#include "HistoManager.hh"
#include "G4Event.hh"
#include "G4EventManager.hh"
#include "G4RunManager.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventAction::EventAction(DetectorConstruction* det)
:G4UserEventAction(), fDetector(det)
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventAction::~EventAction()
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventAction::BeginOfEventAction(const G4Event*)
{
//energy deposited per event
for (G4int k=0; k<kMaxAbsor; k++) { fEdepAbsor[k] = 0.0; }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
std::ofstream saida1("Events1.dat");
void EventAction::EndOfEventAction(const G4Event* anEvent)
{
//get Run
Run* run = static_cast<Run*>(
G4RunManager::GetRunManager()->GetNonConstCurrentRun());
//plot energy deposited per event
G4double TotalEdep(0.);
for (G4int k=1; k<=fDetector->GetNbOfAbsor(); k++) {
if (fEdepAbsor[k] > 0.) {
run->AddEdep(k,fEdepAbsor[k]);
G4AnalysisManager::Instance()->FillH1(k, fEdepAbsor[k]);
TotalEdep += fEdepAbsor[k];
}
}
G4double PrimaryEnergy = anEvent->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy();
G4ThreeVector direction;
direction = anEvent->GetPrimaryVertex()->GetPrimary()->GetMomentumDirection();
G4double Theta;
G4double cosTheta =direction.z();
Theta = std::acos(-cosTheta);
saida1 << PrimaryEnergy/CLHEP::eV << " " << Theta / CLHEP::deg << " " << fEdepAbsor[1]/CLHEP::eV << " " << fEdepAbsor[2]/CLHEP::eV << " " << fEdepAbsor[3]/CLHEP::eV << " " << fEdepAbsor[4]/CLHEP::eV << " " << fEdepAbsor[5]/CLHEP::eV << " " << fEdepAbsor[6]/CLHEP::eV << " " << fEdepAbsor[7]/CLHEP::eV << " " << fEdepAbsor[8]/CLHEP::eV << " " << fEdepAbsor[9]/CLHEP::eV << " " << fEdepAbsor[10]/CLHEP::eV << " " << fEdepAbsor[11]/CLHEP::eV << " " << fEdepAbsor[12]/CLHEP::eV << " " << fEdepAbsor[13]/CLHEP::eV << " " << TotalEdep/CLHEP::eV << G4endl;
if (TotalEdep > 0.) {
run->AddTotEdep(TotalEdep);
}
std::ofstream particleCountFile("ParticleCount.dat", std::ios::app);
run->PrintParticleData(particleCountFile);
particleCountFile.close();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
and SteppingAcion.cc
#include "SteppingAction.hh"
#include "EventAction.hh"
#include "G4Step.hh"
#include "G4Track.hh"
#include "G4Event.hh"
#include "G4SystemOfUnits.hh"
#include "G4ParticleDefinition.hh"
#include <fstream>
#include <string>
SteppingAction::SteppingAction(DetectorConstruction* detector, EventAction* eventAction)
: fDetector(detector), fEventAction(eventAction),
fSumOfStepLength(0.), fDepositedEnergy(0.), fNumberOfSteps(0), fPrimaryPDGCode(0)
{
fOutputFile.open("ParticulasRastreadas.dat", std::ios::out | std::ios::trunc);
}
SteppingAction::~SteppingAction()
{
if (fOutputFile.is_open()) {
fOutputFile.close();
}
}
void SteppingAction::UserSteppingAction(const G4Step* step)
{
const G4Track* track = step->GetTrack();
if (track->GetParentID() == 0) {
fPrimaryPDGCode = track->GetDefinition()->GetPDGEncoding();
if (track->GetTrackID() == 1) {
fSumOfStepLength += step->GetStepLength();
fDepositedEnergy += step->GetTotalEnergyDeposit();
fNumberOfSteps += 1;
fCurrentProcess = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
G4double primaryEnergy = track->GetKineticEnergy() / MeV;
if (fOutputFile.is_open()) {
fOutputFile << std::setw(10) << fPrimaryPDGCode
<< std::setw(15) << fSumOfStepLength / pc
<< std::setw(20) << primaryEnergy
<< std::setw(25) << fDepositedEnergy / MeV
<< std::setw(10) << fNumberOfSteps
<< std::setw(30) << fCurrentProcess << std::endl;
}
}
}
}
_Geant4 Version: 11.1.2
_Operating System: Ubuntu 20.04
_Compiler/Version g++ --version
g++ (Ubuntu 9.4.0-1ubuntu1~20.04.2) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
_CMake Version: 3.16.3