Message: UserSteppingAction histogram method Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

None UserSteppingAction histogram method 

Forum: Event and Track Management
Date: 25 Apr, 2006
From: Dan Fry <Dan Fry>

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

Inline Depth:
 1 1
 All All
Outline Depth:
 1 1
 2 2
 All All
Add message: (add)

1 None: Re: UserSteppingAction histogram method   (michel maire - 26 Apr, 2006)
 Add Message Add Message
to: "UserSteppingAction histogram method"

 Subscribe Subscribe

This site runs SLAC HyperNews version 1.11-slac-98, derived from the original HyperNews


[ Geant 4 Home | Geant 4 HyperNews | Search | Request New Forum | Feedback ]