Please fill out the following information to help in answering your question, and also see tips for posting code snippets. If you don’t provide this information it will take more time to help with your problem!
_Geant4 Version: geant4v11.2.1
_Operating System: ubuntu22.04
_CMake Version: 3.22.1
Hello,everyone!
I have a problem from rdecay02 example.The SteppingAction.cc script provides information such as energy deposition,trackID,weight and so on.I simulate the radioactive decay from Ac-227.When I print the information ,the energy deposition looks well.However,when I fill the results in histogram.Some information is lost.
#
# Macro file for "rdecay02.cc"
#
/control/verbose 1
###/control/cout/ignoreThreadsExcept 0
###/run/numberOfThreads 2
/run/verbose 1
#
/process/em/verbose 0
/process/had/verbose 0
#
/run/initialize
#
# Set a very high time threshold to allow all decays to happen
/process/had/rdm/thresholdForVeryLongDecayTime 1.0e+60 year
#
/process/had/rdm/selectVolume Target
/gun/particle ion
/gun/ion 89 227
#/gun/ion 57 138
#
/analysis/setFileName alpha
/analysis/h1/set 0 1000 0. 10. MeV
/analysis/h1/set 1 1000 0. 10. MeV
/analysis/h1/set 2 1000 0. 10. MeV
/analysis/h1/set 3 1000 0. 10. MeV
/analysis/h1/set 4 1000 0. 10. MeV
/analysis/h1/set 5 1000 0. 10. MeV
/analysis/h1/set 6 1000 0. 10. MeV
/analysis/h1/set 7 1000 0. 10. MeV
/analysis/h1/set 8 1000 0. 3. MeV
/analysis/h1/set 9 1000 0. 3. MeV
#
/run/printProgress 100
/run/beamOn 1000
Below is the code snippet and result of my output information.
G4int iVol = 0;
if (lVolume == fDetector->GetLogicTarget()) iVol = 1;
if (lVolume == fDetector->GetLogicDetector()) iVol = 2;
// count processes
//
const G4StepPoint* endPoint = aStep->GetPostStepPoint();
const G4VProcess* process = endPoint->GetProcessDefinedStep();
run->CountProcesses(process, iVol);
G4String processName = aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
//if(iVol == 1) G4cout << " step in laBr : "<< processName <<G4endl;
// else if(iVol == 2) G4cout << " step in CsI : "<< processName <<G4endl;
// energy deposit
//
G4double edepStep = aStep->GetTotalEnergyDeposit(); // Get Energy
if (edepStep <= 0.) return;
G4double time = aStep->GetPreStepPoint()->GetGlobalTime();
G4double weight = aStep->GetPreStepPoint()->GetWeight();
G4Track* trackStep = aStep->GetTrack();
G4int trackNum = trackStep->GetTrackID();
G4String particleID = aStep->GetTrack()->GetParticleDefinition()->GetParticleName();
if(iVol == 1) {
G4cout << " particle ID : " << particleID << "\t"
<< " trackID : " << trackNum << "\t"
<< " energy Deposition : " << edepStep <<"\t"
<< " weight : "<<weight <<"\t"
<< G4endl;
}
fEventAction->AddEdep(iVol, edepStep, time, weight);
G4WT11 > particle ID : alpha trackID : 74 energy Deposition : 6.81905 weight : 1
G4WT11 > particle ID : Po215 trackID : 73 energy Deposition : 0.127057 weight : 1
G4WT11 > particle ID : alpha trackID : 76 energy Deposition : 7.386 weight : 1
G4WT11 > particle ID : Pb211 trackID : 75 energy Deposition : 0.140247 weight : 1
G4WT11 > particle ID : e- trackID : 79 energy Deposition : 0.0695136 weight : 1
G4WT11 > particle ID : Bi211[831.960] trackID : 77 energy Deposition : 6.7096e-07 weight : 1
G4WT11 > particle ID : e- trackID : 83 energy Deposition : 0.0123528 weight : 1
G4WT11 > particle ID : e- trackID : 82 energy Deposition : 0.0362911 weight : 1
G4WT11 > particle ID : Bi211[404.866] trackID : 80 energy Deposition : 4.64148e-07 weight : 1
G4WT11 > particle ID : gamma trackID : 88 energy Deposition : 0.001556 weight : 1
G4WT11 > particle ID : e- trackID : 90 energy Deposition : 0.218489 weight : 1
G4WT11 > particle ID : e- trackID : 89 energy Deposition : 0.184821 weight : 1
G4WT11 > particle ID : Bi211 trackID : 87 energy Deposition : 4.17116e-07 weight : 1
G4WT11 > particle ID : alpha trackID : 92 energy Deposition : 6.62216 weight : 1
Below is the code snippet from EventAction.cc.
void EventAction::AddEdep(G4int iVol, G4double edep,
G4double time, G4double weight)
{
// initialize t0
if (fTime0 < 0.) fTime0 = time;
// out of time window ?
const G4double TimeWindow (100*nanosecond);
if (std::fabs(time - fTime0) > TimeWindow) return;
if (iVol == 1) { fEdep1 += edep; fWeight1 += edep*weight;}
if (iVol == 2) { fEdep2 += edep; fWeight2 += edep*weight;}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventAction::EndOfEventAction(const G4Event* event)
{
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
G4double Etot = fEdep1 + fEdep2;
G4double Wtot = (fWeight1 + fWeight2)/Etot;
//G4cout << "fEdep1 : " <<fEdep1 << "\t" << " fWeight1 : " <<fWeight1 <<G4endl;
// pulse height in target
//
if (fEdep1 > 0.) {
fWeight1 /= fEdep1;
analysisManager->FillH1(0, fEdep1,fWeight1);
analysisManager->FillNtupleDColumn(1,fEdep1);
}
// pulse height in detector
//
if (fEdep2 > 0.) {
fWeight2 /= fEdep2;
analysisManager->FillH1(2, fEdep2, fWeight2);
analysisManager->FillNtupleDColumn(2,fEdep2);
}
// total
//
analysisManager->FillH1(2, Etot, Wtot);
analysisManager->FillNtupleDColumn(3,Etot);
analysisManager->AddNtupleRow();// must have
// threshold in target and detector
const G4double Threshold1(10*keV), Threshold2(10*keV);
//coincidence, anti-coincidences
//
G4bool coincidence = ((fEdep1 >= Threshold1) && (fEdep2 >= Threshold2));
G4bool anti_coincidence1 = ((fEdep1 >= Threshold1) && (fEdep2 < Threshold2));
G4bool anti_coincidence2 = ((fEdep1 < Threshold1) && (fEdep2 >= Threshold2));
if (coincidence) analysisManager->FillH1(3, fEdep2, fWeight2);
if (anti_coincidence1) analysisManager->FillH1(4, fEdep1, fWeight1);
if (anti_coincidence2) analysisManager->FillH1(5, fEdep2, fWeight2);
// pass energies to Run
//
Run* run = static_cast<Run*>(
G4RunManager::GetRunManager()->GetNonConstCurrentRun());
run->AddEdep (fEdep1, fEdep2);
}
}
