Error occurs when an optical photon hits the edge of a cubic scintillator

Dear all,
I find a problem when I run the Geant4 example LXe in folder examples/extended/optical/LXe. I generate an optical photon from zero point with momentum direction (1,1,0), and it is supposed to hit the edge of the cubic scintillator. As a result, an error occurs and program terminates.
Here is my PrimaryGeneratorAction file, you can replace the original LXePrimaryGeneratorAction.cc with it.
LXePrimaryGeneratorAction.cc (3.8 KB)
And the error message is as follow:
Screenshot from 2022-09-14 12-49-15

For the last step, the optical boundary status becomes NoRindex and the track status is fStopAndKill, and I wonder why program is still terminated even though this “wrong” track had been killed already?
The probability is very small for this case, but it may happen when I run thousands of events, how can I solve/avoid this problem, could you give some advice?
Thank you so much!

1 Like

Hi,

There is an issue with boundary scattering, when the opticalphoton is reflected back into the volume it came from. When this happens, there are two G4Steps involved. One of them has zero distance and the boundary status is StepTooSmall. This occurs to set the pre and post step volumes correctly.

To do this, G4OpBoundaryProcess checks to see if the step is zero. But, sometimes the step is small but non-zero, I think because of geometrical tolerances. So G4OpBoundaryProcess checks to see if the step size is less than the geometrical tolerance, here:
https://geant4.kek.jp/lxr/source/processes/optical/src/G4OpBoundaryProcess.cc#L186

The LXe example has code to double-check this:
https://geant4.kek.jp/lxr/source/examples/extended/optical/LXe/src/LXeSteppingAction.cc#L168

However, it happens that sometimes this step has a step length of more than the geometrical tolerance (the value fCarTolerance in G4OpBoundaryProcess, which is 1e-9 mm). When this happens LXe gives an error.

There are several possible ways around this.

  1. It is a rare occurrence, so you could simply delete the FatalException in LXeSteppingAction (or change it to JustWarning). Still, it is wrong for that event.

  2. In G4OpBoundaryProcess, line 186 (link above), change the limit to a multiple of fCarTolerance. 10x seems to work. However, if there are steps that are actually this small length, G4OpBoundaryProcess will treat them incorrectly.

  3. Ideally we’d add a flag, rather than checking for small step lengths. I haven’t got it to work yet.

Any comments appreciated. The bug report is here:

https://bugzilla-geant4.kek.jp/show_bug.cgi?id=2510

Thanks for your detailed explanation! In the 2nd solution, you mentioned that the limit of fCarTolerance can be increased by 10 times. I changed it in source code, but how can this change be applied to program. Should I build and install Geant4 again?

Yes, it should just be a matter of running make and make install.

To be clear, don’t change fCarTolerance. Instead, multiply it by a factor in line 186:

  if(aTrack.GetStepLength() <= 10. * fCarTolerance)

This is because fCarTolerance is used elsewhere.

Also, as noted, I haven’t fully verified this solution. It seems to work with limited testing. So you should test your application as needed.

Hi, thanks for reply. I multiply the fCarTolerance by 10 in line 186, and also run make and make install, but it seems that the error still exists.

I am not sure if I make a mistake, but I think when a photon hits the edge, it is like a singularity (though impossible in real world), how does it reflect?

By the way, I simplify the geometry to only two parts (a cubic scintillator and a layer of reflector), and I find that when I change the optical surface from G4LogicalSkinSurface to G4LogicalBorderSurface (in LXeMainVolume.cc), the program runs without error.
Screenshot from 2022-09-19 10-55-04

Maybe I mis-diagnosed the problem. I see in your first post the step length is 898 um. Could you please run in debug mode and post the backtrace. Also, note you don’t have to modify PrimaryGeneratorAction–instead, you can modify the primary in the macro file.

The step length is 898 um, because the thickness of reflector is 500 um, the pre-step point is in the inner layer and the post-step point is in the outside.

I run example LXe by gdb ./LXe , is this debug mode? And the message is almost the same:

Adding task 0 to task-group...
G4WT7 > ### Run 0 starts on worker thread 7.
G4WT7 > 
-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : LXeExpl01
      issued by : LXeSteppingAction::UserSteppingAction()
LXeSteppingAction::UserSteppingAction(): No reallocation step after reflection!
Something is wrong with the surface normal or geometry

*** Fatal Exception *** core dump ***
G4WT7 > G4Track (0x7fffb0833cd0) - track ID = 1, parent ID = 0
G4WT7 >  Particle type : opticalphoton - creator process : not available
G4WT7 >  Kinetic energy : 7.07 eV  - Momentum direction : (-0.707107,0.707107,0)
G4WT7 >  Step length : 898.026 um  - total energy deposit : 7.07 eV 
G4WT7 >  Pre-step point : (89,89,0) - Physical volume : housing (Al)
G4WT7 >  - defined by : Transportation - step status : 1
G4WT7 >  Post-step point : (88.365,89.635,0) - Physical volume : expHall (Vacuum)
G4WT7 >  - defined by : Transportation - step status : 1
G4WT7 >  *** Note: Step information might not be properly updated.
G4WT7 > 
-------- EEEE -------- G4Exception-END --------- EEEE -------

G4WT7 > 
QObject::setParent: Cannot set parent, new parent is in a different thread

Now I guess the problem is due to G4OpticalSurface. When the reflector is replaced by air along edges, photons will reflect correctly, and there are no errors.
Screenshot from 2022-09-20 21-41-07

I am sorry I mentioned the usage of G4LogicalBorderSurface in my last post, I checked it and found that I reversed the order of physical volumes, so photons are totally absorbed.

The message “No reallocation step…”, I believe, only happens with the error I described earlier. Try a few of things.

Set /track/verbose 1, send 1 opticalphoton to hit the surface not at the edge, and see 2 steps at the surface. One will have 0 step length. Then let the optical photon hit the edge and compare the printouts.

Print out the boundary status each step (it’s the variable ‘boundaryStatus’ a few lines above).

Print out what UserSteppingAction thinks the step length is (right before, or in, the FatalException). Use theTrack->GetStepLength().

I set /tracking/verbose 1 and /process/optical/boundary/verbose 2.

For the photon hits the surface, the boundary status is always “×××Reflection” followed by “StepTooSmall”, and the “StepTooSmall” step has 0 step length as you said. I list the last 3 steps (it reflects many times and is absorbed by scintillator in the last step):

G4WT2 > Step#    X(mm)    Y(mm)    Z(mm) KinE(MeV)  dE(MeV) StepLeng TrackLeng  NextVolume ProcName
G4WT2 >  *** SpikeReflection ***
G4WT2 >    13       50    -42.5        0   2.8e-06        0     22.7       477     housing Transportation
G4WT2 > G4VDiscreteProcess::PostStepGetPhysicalInteractionLength() - [ OpBoundary]
G4WT2 >  Particle type - opticalphoton
G4WT2 >    mass:        0[GeV]
G4WT2 >    charge:      0[e]
G4WT2 >    Direction x: -0.9438583563660174, y: 0.3303504247281061, z: 0
G4WT2 >    Total Momentum = 2.8e-09[GeV]
G4WT2 >    Momentum: -2.642803397824849e-09[GeV], y: 9.249811892386971e-10[GeV], z: 0[GeV]
G4WT2 >    Total Energy   = 2.8e-09[GeV]
G4WT2 >    Kinetic Energy = 2.8e-09[GeV]
G4WT2 >  MagneticMoment  [MeV/T]: 0
G4WT2 >    ProperTime     = 0[ns]
G4WT2 >  in Material  Al
G4WT2 > InteractionLength= 1.797693134862316e+307[cm] 
G4WT2 >  Photon at Boundary! 
G4WT2 >  thePrePV:  housing
G4WT2 >  thePostPV: scintillator
G4WT2 >  *** StepTooSmall ***
G4WT2 >    14       50    -42.5        0   2.8e-06        0        0       477 scintillator Transportation
G4WT2 > G4VDiscreteProcess::PostStepGetPhysicalInteractionLength() - [ OpBoundary]
G4WT2 >  Particle type - opticalphoton
G4WT2 >    mass:        0[GeV]
G4WT2 >    charge:      0[e]
G4WT2 >    Direction x: -0.9438583563660174, y: 0.3303504247281061, z: 0
G4WT2 >    Total Momentum = 2.8e-09[GeV]
G4WT2 >    Momentum: -2.642803397824849e-09[GeV], y: 9.249811892386971e-10[GeV], z: 0[GeV]
G4WT2 >    Total Energy   = 2.8e-09[GeV]
G4WT2 >    Kinetic Energy = 2.8e-09[GeV]
G4WT2 >  MagneticMoment  [MeV/T]: 0
G4WT2 >    ProperTime     = 0[ns]
G4WT2 >  in Material  LYSO
G4WT2 > InteractionLength= 1.797693134862316e+307[cm] 
G4WT2 >    15     20.7    -32.2        0   2.8e-06  2.8e-06       31       508 scintillator OpAbsorption

For the photon hits the edge, there are only 2 steps (“SpikeReflection” and “NoRindex”):

G4WT5 > Step#    X(mm)    Y(mm)    Z(mm) KinE(MeV)  dE(MeV) StepLeng TrackLeng  NextVolume ProcName
G4WT5 >  *** SpikeReflection ***
G4WT5 >     1       50       50        0   2.8e-06        0     70.7      70.7     housing Transportation
G4WT5 > G4VDiscreteProcess::PostStepGetPhysicalInteractionLength() - [ OpBoundary]
G4WT5 >  Particle type - opticalphoton
G4WT5 >    mass:        0[GeV]
G4WT5 >    charge:      0[e]
G4WT5 >    Direction x: -0.7071067811865476, y: 0.7071067811865476, z: 0
G4WT5 >    Total Momentum = 2.8e-09[GeV]
G4WT5 >    Momentum: -1.979898987322333e-09[GeV], y: 1.979898987322333e-09[GeV], z: 0[GeV]
G4WT5 >    Total Energy   = 2.8e-09[GeV]
G4WT5 >    Kinetic Energy = 2.8e-09[GeV]
G4WT5 >  MagneticMoment  [MeV/T]: 0
G4WT5 >    ProperTime     = 0[ns]
G4WT5 >  in Material  Al
G4WT5 > InteractionLength= 1.797693134862316e+307[cm] 
G4WT5 >  Photon at Boundary! 
G4WT5 >  thePrePV:  housing
G4WT5 >  thePostPV: expHall
G4WT5 >  Old Momentum Direction: (-0.7071067811865476,0.7071067811865476,0)
G4WT5 >  Old Polarization:       (0.4496521907555314,0.4496521907555314,0.7717679798342909)
G4WT5 >  *** NoRINDEX ***
G4WT5 >     2       49       51        0   2.8e-06  2.8e-06     1.41      72.1     expHall Transportation

As for getting step length in UserSteppingAction, the result is the same as /tracking/verbose 1 prints out, and the only difference is numerical precision. For example when photon hits edge, the two steps’ length is:

G4WT5 > step number: 1 step length: 70.71067811865476
G4WT5 > step number: 2 step length: 1.414213562373095

Here is backtrace:

Backtrace:
[PID=11032, TID=6][ 0/21]> /lib/x86_64-linux-gnu/libc.so.6(gsignal+0xcb) [0x7fd346cc800b]
[PID=11032, TID=6][ 1/21]> /lib/x86_64-linux-gnu/libc.so.6(abort+0x12b) [0x7fd346ca7859]
[PID=11032, TID=6][ 2/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4global.so(_Z11G4ExceptionPKcS0_19G4ExceptionSeverityS0_+0x1117) [0x7fd3472f3e87]
[PID=11032, TID=6][ 3/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4global.so(_Z11G4ExceptionPKcS0_19G4ExceptionSeverityRNSt7__cxx1119basic_ostringstreamIcSt11char_traitsIcESaIcEEE+0xa5) [0x7fd3472f4025]
[PID=11032, TID=6][ 4/21]> ./MLC(+0x208b9) [0x5624938b68b9]
[PID=11032, TID=6][ 5/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4tracking.so(_ZN17G4SteppingManager8SteppingEv+0x586) [0x7fd349ab6936]
[PID=11032, TID=6][ 6/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4tracking.so(_ZN17G4TrackingManager15ProcessOneTrackEP7G4Track+0x118) [0x7fd349acc488]
[PID=11032, TID=6][ 7/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4event.so(_ZN14G4EventManager12DoProcessingEP7G4Event+0x893) [0x7fd349b0d303]
[PID=11032, TID=6][ 8/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4tasking.so(_ZN22G4WorkerTaskRunManager15ProcessOneEventEi+0x3c) [0x7fd34a3cbc1c]
[PID=11032, TID=6][ 9/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4tasking.so(_ZN22G4WorkerTaskRunManager11DoEventLoopEiPKci+0x17e) [0x7fd34a3cbb5e]
[PID=11032, TID=6][10/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4tasking.so(_ZN22G4WorkerTaskRunManager6DoWorkEv+0x15d) [0x7fd34a3cbead]
[PID=11032, TID=6][11/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4tasking.so(_ZN22G4TaskRunManagerKernel17ExecuteWorkerTaskEv+0x6a) [0x7fd34a3c8d3a]
[PID=11032, TID=6][12/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4tasking.so(+0x3a2d6) [0x7fd34a3b82d6]
[PID=11032, TID=6][13/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4tasking.so(_ZNSt13__future_base13_State_baseV29_M_do_setEPSt8functionIFSt10unique_ptrINS_12_Result_baseENS3_8_DeleterEEvEEPb+0x2d) [0x7fd34a3c038d]
[PID=11032, TID=6][14/21]> /lib/x86_64-linux-gnu/libpthread.so.0(+0x114df) [0x7fd34598e4df]
[PID=11032, TID=6][15/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4tasking.so(_ZN3PTL4TaskIvJEEclEv+0x12c) [0x7fd34a3c050c]
[PID=11032, TID=6][16/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4ptl.so.0(_ZN3PTL10ThreadPool14execute_threadEPNS_14VUserTaskQueueE+0x351) [0x7fd3471ed6d1]
[PID=11032, TID=6][17/21]> /home/wangjunhao/geant4/geant4-v11.0.2-install/lib/libG4ptl.so.0(_ZN3PTL10ThreadPool12start_threadEPS0_PSt6vectorISt10shared_ptrINS_10ThreadDataEESaIS5_EEl+0x1a3) [0x7fd3471ede33]
[PID=11032, TID=6][18/21]> /lib/x86_64-linux-gnu/libstdc++.so.6(+0xd6de4) [0x7fd3470b7de4]
[PID=11032, TID=6][19/21]> /lib/x86_64-linux-gnu/libpthread.so.0(+0x8609) [0x7fd345985609]
[PID=11032, TID=6][20/21]> /lib/x86_64-linux-gnu/libc.so.6(clone+0x43) [0x7fd346da4133]

: Aborted (Signal sent by tkill() 11032 1000)
Aborted (core dumped)

Thanks for doing these tests. I see what is happening. Usually, when an opticalphoton in volume A is incident on the surface with volume B, if it is reflected, it is reflected back into volume A. However, in this case, the opticalphoton hits the edge of the volume, and is reflected into volume C. You can see this be comparing thePostPV in the output above.

Volume C is expHall, which has no RINDEX defined. The boundary status should be both StepTooSmall and NoRINDEX, but it can only be one. It gets set to NoRINDEX, and the check in LXESteppingAction fails.

This is a literal edge case!

What to do?

  • you could verify the above by defining RINDEX for the expHall. Maybe also give expHall a short ABSORPTIONLENGTH so that optical photons are not propagated in it.
  • without an RINDEX, the correct behaviour is for the opticalphoton to be killed. I believe this would happen correctly, except the SteppingAction check is failing first. You could
    • modify the SteppingAction FatalException to allow NoRINDEX (but I don’t know if there are other cases that should correctly fail)
    • change the FatalException to JustWarning (we’ve made this change for future Geant4 versions), and make sure they are few, and/or understood
  • the probability of this happening must be pretty low (maybe two floats have to be equal). Don’t shoot photons at the edge :slight_smile:

I’d consider this a bug, so you can file a bug report if you like.

Thank you for your suggestions, I tried all of them.

First, the expHall has a defined RINDEX property in LXeDetectorConstruction.cc:

  G4MaterialPropertiesTable* vacuum_mt = new G4MaterialPropertiesTable();
  vacuum_mt->AddProperty("RINDEX", "Air");
  fVacuum->SetMaterialPropertiesTable(vacuum_mt);

  fExperimentalHall_log =
    new G4LogicalVolume(fExperimentalHall_box, fVacuum, "expHall_log", 0, 0, 0);

And I give expHall 1 um ABSORPTIONLENGTH, the result does not changed. However, it does work when change the FatalException to JustWarning!

In fact, I found this problem when I run my own simulation rewritten from LXe which requires a lot of events, so the “hit edge” event is unavoidable even though the probability is low.

And thank you again for your help, I think my problem is solved!