Hi all.
I have a problem to deal with G4Step::Settotalenergydeposit() function. I attach the code below…
I assume this problem comes from the “const”… but I do not sure it is because of it or not.
Also, I wonder is that correct that change the total energy deposit by set total energy deposit function if I do want to change the value that EDEPMAP takes. In my program, it uses EDEPMAP for saving total energy deposit. but I want to give weighting factor depending on the particle type like Deposit energy * weighting factor(e.g. 20). To do this, I though that I have to change the total energy deposit in advance using Settotalenergydeposit function!
Thank you for in advance!
below is G4steppingaction.hh
#ifndef TETSteppingAction_h
#define TETSteppingAction_h 1
#include “G4UserSteppingAction.hh”
#include “globals.hh”
#include “G4ThreeVector.hh”
#include “G4Step.hh”
#include “G4RunManager.hh”
#include “G4UnitsTable.hh”
class G4LogicalVolume;
// *********************************************************************
// With very low probability, because of the internal bug of G4TET, the
// particles can be stuck in the vertices of tetrahedrons. This
// UserSteppingAction class was written to slightly move these stuck
// particles.
// – UserSteppingAction: Slightly move the stuck particles.
// *********************************************************************
class TETSteppingAction : public G4UserSteppingAction
{
public:
TETSteppingAction();
virtual ~TETSteppingAction();
virtual void UserSteppingAction(const G4Step*);
private:
G4double kCarTolerance;
G4int stepCounter;
G4bool checkFlag;
// G4double radiationfactor;
// G4double energys;
// G4double equivalent_dose;
G4Step* G4_not_const_step;
G4double count;
G4RunManager* stRun;
};
#endif
Below is G4Steppingaction.cc
#include “TETSteppingAction.hh”
#include “G4Event.hh”
#include “G4MTRunManager.hh”
#include “G4SystemOfUnits.hh”
TETSteppingAction::TETSteppingAction()
: G4UserSteppingAction(), kCarTolerance(1.0000000000000002e-07), stepCounter(0),
checkFlag(0)
{}
TETSteppingAction::~TETSteppingAction()
{}
void TETSteppingAction::UserSteppingAction(const G4Step* step)
{
// G4cout<<“stepping action cc”<<G4endl;
// G4double radiationfactor;
// G4double energys;
// G4double equivalent_dose;
// Slightly move the particle when the step length of five continuous steps is
// shorter than the tolerance (0.1 nm)
//
G4Track* theTrack = step->GetTrack();
G4bool CheckingLength = (step->GetStepLength() < kCarTolerance);
G4double radiationfactor;
G4double equivalent_dose;
G4double energys;
// G4double cn = step->GetPreStepPoint()->GetPhysicalVolume()->GetCopyNo();
// G4String namee = step->GetPreStepPoint()->GetPhysicalVolume()->GetName();
G4String name = step->GetTrack()->GetParticleDefinition()->GetParticleName();
G4double delta_ke = step->GetDeltaEnergy();
energys = step->GetTotalEnergyDeposit();
// G4String post_pv_name = step->GetPostStepPoint()->GetPhysicalVolume()->GetName();
G4String post_pv_name = step->GetPreStepPoint()->GetPhysicalVolume()->GetName();
// if (post_pv_name == "wholePhantom"){
//// G4cout<<name<<G4endl;
// if (name == "neutron"){
// G4cout<<step->GetTotalEnergyDeposit()<<G4endl;
// }
// }
// if (post_pv_name =="wholePhantom" && delta_ke>0.0*MeV){
// if (name == "neutron"){
// G4cout<<name<<G4endl;}}
// G4cout<<energys<<G4endl;
if (post_pv_name =="wholePhantom" && energys>0.0*MeV){
// if (name == "neutron"){
// G4cout<<name<<G4endl;}
if (name=="gamma"){
radiationfactor = 1.0;
// G4cout<<"gamma parts"<<G4endl;
}
else if (name=="e-"){
radiationfactor = 1.0;
// G4cout<<"e- parts"<<G4endl;
}
//
// // else if (name=="alpha"){
// // radiationfactor = 20.0;
// // // G4cout<<"e- parts"<<G4endl;
// // }
else if (name=="neutron"){
count +=1;
if (energys< 1*MeV){
radiationfactor = 2.5+18.2*exp(-1*(pow(log(energys),2))/6);
// G4cout<<"neutron 1"<<G4endl;
}
else if (energys>= 1*MeV && energys <=50*MeV){
radiationfactor = 5.0+17.0*exp(-1*(pow(log(2*energys),2))/6);
// G4cout<<"neutron 2"<<G4endl;
}
// G4cout<<radiationfactor<<G4endl;
}
else {
theTrack->SetTrackStatus(fStopAndKill);
radiationfactor=0;
}
if (count!=0){
// G4cout<<"neutron detected"<<G4endl;
}
// G4cout<<radiationfactor<<G4endl;
// G4cout<<radiationfactor<<std::setw(5)<<energys<<std::setw(5)<<radiationfactor<<G4endl;
equivalent_dose = energys*radiationfactor;
// G4cout<<equivalent_dose<<G4endl;
// G4cout<<name<<std::setw(10)<<radiationfactor<<std::setw(5)<<energys<<std::setw(5)<<equivalent_dose<<G4endl;
G4cout << "before"<<step->GetTotalEnergyDeposit()/MeV<<G4endl;
// step->ResetTotalEnergyDeposit();
step->SetTotalEnergyDeposit(equivalent_dose);
// step->SetTotalEnergyDeposit(equivalent_dose);
G4cout << “after”<GetTotalEnergyDeposit()/MeV<<G4endl;
}
// G4double equivalent_dose = step->SetTotalEnergyDeposit(energys*radiationfactor);
// // G4double energy_ke = step->GetPreStepPoint()->GetTotalEnergy();
// //
// G4int ptlID = step->GetTrack()->GetParticleDefinition()->GetParticleDefinitionID();
// G4cout<<name<<G4endl;
// if (name != "gamma" || name != "neutron"){
// theTrack->SetTrackStatus(fStopAndKill);
//// G4cout << name <<G4endl;
// }
// G4cout << name <<G4endl;
// else{};
// if (namee =="wholePhantom"){ //
// if (energys >0.0){
// G4cout<< cn <<G4endl;
// // G4cout << name<<" " <<energys<<G4endl;
// }}
if(CheckingLength)
{
++stepCounter;
if( checkFlag && stepCounter>=5 )
{
// kill the track if the particle is stuck even after the slight move
// (this hardly occurs)
theTrack->SetTrackStatus(fStopAndKill);
stepCounter=0;
checkFlag=0;
}
else if( stepCounter>=5 )
{
// if a particle is at the same position (step length < 0.1 nm) for five consecutive steps,
// slightly move (0.1 nm) the stuck particle in the direction of momentum
theTrack->SetPosition(theTrack->GetPosition() + theTrack->GetMomentumDirection()*kCarTolerance);
checkFlag=1;
}
}
else stepCounter=0;
}