For the life of me I cannot find the cause of a series of strange spikes in
the histogram of the energy deposited as a function of depth in water for
proton irradiation.
I have included the stepping action code below. If anyone has any suggestions
please post.
What I am doing here is simply a "virtual" voxelization of a water phantom.
By pulling the particle position at each step, using the preStep information,
I bin the deposited energy. Everything works fine and I have checked this
against using a full 3d voxelization and then integrating in the y-z plane
at fixed depth (x-axis). However, the problem arises with respect to the
production cut. If I make the electron and proton production cut greater
than 0.1mm I see a series of spikes in the binned distribution of deposited
energy with a spacing and intensity between spikes that decreases with
depth in water. I have checked the histogram method by writing the output
at each step and performing the histogram with a commercial analysis
package. The spikes are still there. At this point I need to see if
there is someone willing to take a look to rule out a misuse of the
steppingAction methodology.
Thanks...
void protonTherapySteppingAction::UserSteppingAction(const G4Step* aStep)
{
G4int ix, iy, iz;
G4double x_ave, y_ave, z_ave;
G4double hphantsize = phtparameters[1]/2.;
FILE *fptr;
G4double EdepStep = aStep->GetTotalEnergyDeposit();
G4StepPoint* get_prepoint = aStep->GetPreStepPoint();
G4StepPoint* get_postpoint = aStep->GetPostStepPoint();
G4VPhysicalVolume* phys_vol = get_prepoint->GetPhysicalVolume();
G4String phys_vol_name = phys_vol->GetName();
const G4Material* currentMat = get_prepoint->GetMaterial();
const G4Track* currentTrack = aStep->GetTrack();
if(phys_vol_name == "water_phantom")
{
const G4ParticleDefinition* currentParticle = currentTrack->GetDefinition();
G4String particle_name = currentParticle->GetParticleName();
G4double Ekin_PreStep = get_prepoint->GetKineticEnergy();
G4double Ekin_PostStep = get_postpoint->GetKineticEnergy();
// G4double Table_Eloss = G4EnergyLossTables::GetDEDX(currentParticle, Ekin_PreStep, currentMat);
const G4ThreeVector &position_pre = get_prepoint->GetPosition();
const G4ThreeVector &position_post = get_postpoint->GetPosition();
const G4VProcess* process = get_postpoint->GetProcessDefinedStep();
G4String process_name = process->GetProcessName();
//Now get pre- and post- step positions
G4double x_pre = position_pre.x();
G4double y_pre = position_pre.y();
G4double z_pre = position_pre.z();
// G4double x_pre=position_pre.getX();
// G4double y_pre=position_pre.getY();
// G4double z_pre=position_pre.getZ();
G4double x_post = position_post.x();
G4double y_post = position_post.y();
G4double z_post = position_post.z();
if((Ekin_PreStep/keV) > 10.)
{
int ipos = (int)((x_pre/mm)/0.5 + 0.5);
int ieng = (int)((Ekin_PreStep/MeV)/0.001 + 0.5);
if(process_name == "hLowEIoni") {nioniz[ipos] += 1;}
else if(process_name == "inelastic") {ninelastic[ipos] += 1;}
else if(process_name == "LElastic") {nLElastic[ipos] += 1;}
else if(process_name == "LowEnergyCompton") {nLECompton[ipos] += 1;}
else if(process_name == "LowEnBrem") {nLowEnBrem[ipos] += 1;}
else if(process_name == "LowEnPhotoElec") {nLowEnPhotoElec[ipos] += 1;}
else if(process_name == "LowEnergyIoni") {nLowEnergyIoni[ipos] += 1;}
else if(process_name == "msc") {nscat[ipos] += 1;}
//else if(process_name == "eIoni") {neioniz[ipos] += 1;}
//Standard hadronic ionization is called "hIoni".
double Energyloss_per_Step = (EdepStep/MeV);
pEngDepo[ipos] += Energyloss_per_Step;
fprintf(file_ptr, "%f %f %f\n", x_pre/mm, x_post/mm, Energyloss_per_Step);
if(particle_name == "proton")
{
//Output range here, based upon an energy threshold of 100keV
if(((Ekin_PostStep/keV) <= 100.) && (ptr_EvtAction->flag != false))
{
int xint = (int)(x_pre/0.5 + 0.5);
int yint = (int)(y_pre/0.5 + 0.5);
int zint = (int)(z_pre/0.5 + 0.5);
prangex[xint] += 1;
prangey[yint] += 1;
prangez[zint] += 1;
ptr_EvtAction->flag = false;
}
// Look to see if the proton is a secondary or primary.
// Primary has a track ID, fTrackID = 1.
G4int fTrackID = currentTrack->GetTrackID();
if(fTrackID != 1) // secondary protons
{
nspro[ipos] += 1;
psprotDose[ipos] += Energyloss_per_Step;
}
else // primary protons
{
pprotDose[ipos] += Energyloss_per_Step;
preStep_Eng[ipos] += Ekin_PreStep/MeV;
postStep_Eng[ipos] += Ekin_PostStep/MeV;
counter[ipos] += 1.;
}
}
else if(particle_name == "e-")
{
nelec[ipos] += 1;
peDose[ipos] += Energyloss_per_Step;
}
else if(particle_name == "e+")
{
npos[ipos] += 1;
pposDose[ipos] += Energyloss_per_Step;
}
else if(particle_name == "gamma")
{
ngam[ipos] += 1;
pgDose[ipos] += Energyloss_per_Step;
}
else if(particle_name == "neutron")
{
nneu[ipos] += 1;
pnDose[ipos] += Energyloss_per_Step;
}
else if(particle_name == "alpha")
{
nalpha[ipos] += 1;
paDose[ipos] += Energyloss_per_Step;
}
else if(particle_name == "deuteron")
{
ndeut[ipos] += 1;
pdDose[ipos] += Energyloss_per_Step;
}
else if(particle_name == "triton")
{
ntriton[ipos] += 1;
ptDose[ipos] += Energyloss_per_Step;
}
x_ave = x_pre/mm;
y_ave = y_pre/mm;
z_ave = z_pre/mm;
y_ave += hphantsize;
z_ave += hphantsize;
// G4cout << "Stepping Action!" << G4endl;
// GetSocialStatus(&x_ave, &y_ave, &z_ave, &ix, &iy, &iz, &phtsize, &num_cell);
// pindex = (int)(1 + ix + (iy - 1)*400 + (iz - 1)*160000);
// eng_Density[ix][iy][iz] += Energyloss_per_Step;
}
}
}
|
UserSteppingAction histogram method 
Feedback
]
icon
and there is no indentation. 
