Problem 1314 - G4VEmProcess::PostStepGetPhysicalInteractionLength bias manager check on current step number
Summary: G4VEmProcess::PostStepGetPhysicalInteractionLength bias manager check on curr...
Status: RESOLVED FIXED
Alias: None
Product: Geant4
Classification: Unclassified
Component: processes/electromagnetic (show other problems)
Version: 9.5
Hardware: All All
: P5 normal
Assignee: Vladimir.Ivantchenko
URL:
Depends on:
Blocks:
 
Reported: 2012-05-18 13:13 CEST by Iwan Cornelius
Modified: 2012-08-01 16:02 CEST (History)
1 user (show)

See Also:


Attachments
screenshot of test case (319.38 KB, image/png)
2012-05-18 13:48 CEST, Iwan Cornelius
Details
patch for G4VEmProcess.cc (1.66 KB, patch)
2012-06-08 07:57 CEST, Iwan Cornelius
Details | Diff
Patch for G4EmBiasingManager.cc (832 bytes, text/plain)
2012-06-08 08:04 CEST, Iwan Cornelius
Details
Geant4 vs EGSnrc dose kernel E=1.25MeV - repositioning algorithm (73.55 KB, image/png)
2012-06-08 08:21 CEST, Iwan Cornelius
Details
Geant4 vs EGSnrc dose kernel E=1.25MeV - forced interactions - PRE patch (66.50 KB, image/png)
2012-06-08 08:23 CEST, Iwan Cornelius
Details
Geant4 vs EGSnrc dose kernel E=1.25MeV - forced interactions - POST patch (66.24 KB, image/png)
2012-06-08 08:23 CEST, Iwan Cornelius
Details

Note You need to log in before you can comment on or make changes to this problem.
Description Iwan Cornelius 2012-05-18 13:13:14 CEST
I have activated forced interactions for gammas with G4ComptonScattering in my Physics List, like so: 
 
G4ComptonScattering* cs   = new G4ComptonScattering;
            cs->ActivateForcedInteraction(1*nanometer,"forced", true);
            cs->SetVerboseLevel(2);

Nb. my primary photons are being generated from within the region "forced". 

In G4VEmProcess::PostStepGetPhysicalInteractionLength a check is made to determine if the biasManager pointer is defined, if the parentID is 0, and whether the current step number is zero. 

Unfortunately, the track number is incremented prior to the G4VEmProcess::PostStepGetPhysicalInteractionLength being called. As such, the current step number is never zero, and the biasFlag is never set (as below).  

if(0 == track.GetCurrentStepNumber()) {
        biasFlag = true; 
	biasManager->ResetForcedInteraction(); 
      }

changing 0 to 1 rectifies this problem. 

I have attached an image of the code snippet from the method in question, as well as the console output  (nb. the comment ***DEBUG increment step number ***  is printed in the relevant method of G4Track). 

I hope this helps. If you have any questions, please don't hesitate to contact me. 

Best Regards, 
Iwan Cornelius
Comment 1 Iwan Cornelius 2012-05-18 13:48:53 CEST
Created attachment 163 [details]
screenshot of test case
Comment 2 Iwan Cornelius 2012-06-08 07:57:30 CEST
Created attachment 166 [details]
patch for G4VEmProcess.cc

Attached is a patch for G4VEmProcess

A change is required in the method:
G4double G4VEmProcess::PostStepGetPhysicalInteractionLength

in the revised approach, the interaction path lengths are derived as per usual

if the step is a primary particle, and it is taking its first step in the forced region, the steplimit enforced by the biasmanager is used to scale the PIL.

Doing this ensures that the relative probability of each interaction is maintained, whilst ensuring that an interaction occurs within the steplimit. 

*See associated patch for G4EMBiasManager
Comment 3 Iwan Cornelius 2012-06-08 08:04:47 CEST
Created attachment 167 [details]
Patch for G4EmBiasingManager.cc

removed line in G4double G4EmBiasingManager::GetStepLimit(G4int coupleIdx, 
                                          G4double previousStep)

whereby the steplimit is randomly sampled for each process. 

if the step limit is randomly sampled for each process (eg. PE, Compt, Conv), then each process is equally likely, resulting in incorrect results. Instead, the step limit for this region is returned, and is used to scale the PIL in order to ensure the correct interaction occurs within the forced interaction length (see patch for G4VEmProcess.cc

if startTracking is false, a zero step limit is returned, as we don't wish to consider forced interactions if it is not the first step in the forced interction region.
Comment 4 Iwan Cornelius 2012-06-08 08:16:47 CEST
In addition to the problem described from 2012-05-18 13:13:14 CEST. I have uncovered another problem with the forced interaction functionality of Geant4.9.5. 

My applicatoin: I am using Geant4 to generate dose kernels for radiotherapy. I am considering a monoenergetic photon beam originating at the origin, and in the (0,0,1) direction. The world volume is 3m^3 and comprised of water. 

I wish to force the beam to interact at the origin, and tally the energy deposition events in the water volume. I am doing this using my own spherical coord system scorer which scores total energy deposition per unit mass per primary particle.

http://irs.inms.nrc.ca/software/egsnrc/documentation/pirs702/node94.html

In order to force the interactions I am using 2 approaches: 
1. using the forced interaction functionality of the G4Vemprocesses. 
process->ActivateForcedInteraction( 1*nanometer,"forced");
where "forced" is a 2x2x2 nanometer volume centred on the origin
2. by recording the position of the first interaction in my water volume, then for each non zero edep event, calculating the displacement vector relative to this, and using this to create my dose distribution. 

I am comparing results to those obtained using EGS4 user code EDKNRC. 
http://irs.inms.nrc.ca/software/egsnrc/documentation/pirs702/node97.html

I was getting good agreement for approach #2, but not with approach #1. I narrowed it down to the fact that the biasingmanager was returning a randomly sampled steplimit for each process. As such PE, COMPT, RAY, CONV were occurring with equal probability. The resulting dose kernel was almost isotropic (compared to the expected forward-peaked dose distribution). After implmenting the attached patches, I have achieved good results with approach #1 (see attached image). 


I hope this helps. Please don't hesitate to contact me should you require further info, including my application.

Best Regards, 
Iwan Cornelius
Comment 5 Iwan Cornelius 2012-06-08 08:21:22 CEST
Created attachment 168 [details]
Geant4 vs EGSnrc dose kernel E=1.25MeV - repositioning algorithm

Figure attached shows dose kernels, D(r,theta), as predicted by Geant4 (left) and EGS4 (right). This approach does not use the forced interaction functionality of Geant4.9.5
Comment 6 Iwan Cornelius 2012-06-08 08:23:12 CEST
Created attachment 169 [details]
Geant4 vs EGSnrc dose kernel E=1.25MeV - forced interactions - PRE patch

Shows the Geant4 (left) and EGS4 dose kernels using existing forced interactions functionality available in geant4.9.5 

geant4 results clearly demonstarted an almost isotropic dose distribution
Comment 7 Iwan Cornelius 2012-06-08 08:23:55 CEST
Created attachment 170 [details]
Geant4 vs EGSnrc dose kernel E=1.25MeV - forced interactions - POST patch

Shows the Geant4 (left) and EGS4 dose kernels using revised forced interactions functionality 

geant4 results show good agreement to EGS4 predictions
Comment 8 Daren Sawkey 2012-06-11 19:27:27 CEST
Hi Iwan,

Thank you for the report and patch. I've changed G4VEmProcess and G4VEnergyLossProcess to check for step number == 1. 

As for the second part, forcing the interactions to occur at at point, I think this would be a new feature. I agree it would be a useful feature, and I'll look at adding it.

Regards  --  Daren.
Comment 9 Vladimir.Ivantchenko 2012-08-01 11:11:46 CEST
Hello Iwan,

Thanks you very much for the detailed report on Geant4 problems. We have addressed reported issues in Geant4 public version 9.6beta. Would it be possible to test it? We also add some improvements for Geant4 biasing on top of 9.6beta, if you have an interest to try out the latest version we can privately provide one in mid-August.

Vladimir
Comment 10 Iwan Cornelius 2012-08-01 16:02:31 CEST
Hi Vladimir, 

You're very welcome. I will be more than happy to test out this application with version 9.6b in mid August. Please forward to iwan.cornelius@qut.edu.au when the time comes. 

Best Regards, 
Iwan