Problem 2442

Summary:	Wrong velocity after an optical transition with non-polished surface
Product:	Geant4	Reporter:	Pierre <pgorel>
Component:	processes/optical	Assignee:	Daren Sawkey <daren.sawkey>
Status:	ASSIGNED ---
Severity:	minor
Priority:	P4
Version:	10.7
Hardware:	All
OS:	All
Attachments:	Mac file to reproduce the bug Log file with tracking verbose 3

Description Pierre 2021-11-02 15:04:20 CET

Good day,
after applying the propose patch for problem #2438 (works great), I found another issue with the velocity during transition between two dielectrics.

I have two volumes (scintillator and WLS) sharing a non-polished surface. Optical photons are emitted by the WLS, and reach the surface. The index of refraction of the scintillator is lower than the index of refraction of the WLS.  G4OpBoundaryProcess::DielectricDielectric() is eventually called. Because the surface is not polished, it is possible that we go through the "do/while (!done)" loop more than once. 

When we enter DielectricDielectric, Material1=WLS and Material2=Scint (with respective Rindex)

If we have have the following sequence:
- FresnelRefraction but Done = false (again, it is a non-polished material so using a FacetNormal) then we go through the loop again.
- "Through" is now true so we swap Materials, Rindex and invert theGlobalNormal.
- FresnelReflection and then exit the loop.

Problem: when we get back into PostStepDoIt, the materials have been swapped, so
Material1=Scint and Material2=WLS. Status at this time is "FresnelReflection" since it is the last process that has happened.

We have call proposeVelocity at the beginning of PostStepDoIt, but at that time (coming from aTrack) it was the velocity of the WLS. No new velocity is proposed so we are going to keep it. It is the wrong velocity as we are now in Scint.

If we had exited after the FresnelRefraction (without swapping), the velocity of Material2 (Scint) would have been proposed, correctly.

The way I see it, it enough to propose the velocity of the Material1 if theStatus is neither FresnelRefraction nor Transmission (and remove the proposeVelocity at the beginning).


I propose that at the end of PostStepDoIt, after the "if(theStatus == FresnelRefraction || theStatus == Transmission)" we add a else statement proposing the velocity of Material1.

Comment 1 Daren Sawkey 2022-04-28 23:29:26 CEST

I'm not able to reproduce this (yet). Using OpNovice2, with the SteppingVerbose updated to print out the velocity each step. In G4OpBoundaryProcess I added printouts when the particle reflects or refracts, and when velocity is changed. In the situation described, where the G4OpBoundaryProcess first refracts then reflects, the velocity does change.

Refractive index in tank is 1.4. World is 1.1.

Can you please give a set of parameters, and ideally /tracking/verbose 3 output, that shows this bug. (verbose 3 is very verbose. Let me know if you'd like the updated OpNovice2 SteppingVerbose.)

********************************************************************************************************* 
* G4Track Information:   Particle = opticalphoton,   Track ID = 1,   Parent ID = 0 
********************************************************************************************************* 
 
Step#      X         Y         Z        KineE     Velocity   dEStep  StepLeng  TrakLeng       Volume     Process 
    0     90 cm      0 fm      0 fm      5 eV   21.2 cm/ns     0 eV      0 fm      0 fm         Tank    initStep 
Initial velocity: 212.288 
PartChange::ProposeVelocity: 212.288 
ProposeVelocity1: 212.288 
refract 
reflect 
    1      1 m    1.73 cm   20.8 cm      5 eV   21.2 cm/ns     0 eV   23.1 cm   23.1 cm         Tank  Transportation 
Initial velocity: 296.824 
PartChange::ProposeVelocity: 296.824 
ProposeVelocity1: 296.824 
    2   2.87 m    69.5 cm   59.2 cm      5 eV   21.2 cm/ns     0 eV   2.02 m    2.26 m         World  OpRayleigh 
Initial velocity: 296.824 
PartChange::ProposeVelocity: 296.824 
ProposeVelocity1: 296.824 
    3   3.72 m    1.05 m    2.15 m       5 eV   29.7 cm/ns     0 eV   1.81 m    4.07 m         World  OpRayleigh 
Initial velocity: 296.824 
PartChange::ProposeVelocity: 296.824 
ProposeVelocity1: 296.824


OpNovice2 macro:
/control/verbose 2
/tracking/verbose 2
/run/verbose 1
/process/optical/verbose 1
/run/numberOfThreads 1
#/control/cout/ignoreThreadsExcept 0

/opnovice2/boxProperty RINDEX    0.000002 1.4     0.000008 1.41
/opnovice2/boxProperty ABSLENGTH 0.000002 100 0.000005 200 0.000008 300
/opnovice2/boxProperty RAYLEIGH  0.000002 1000   0.000008 1000

/opnovice2/worldProperty RINDEX    0.000002 1.01    0.000008 1.01
/opnovice2/worldProperty ABSLENGTH 0.000002 1000000 0.000005 2000000 0.000008 3000000
/opnovice2/worldProperty RAYLEIGH  0.000002 1000   0.000008 1000

/opnovice2/surfaceModel unified
/opnovice2/surfaceType dielectric_dielectric

/opnovice2/surfaceFinish ground
/opnovice2/surfaceSigmaAlpha 0.5

#/opnovice2/surfaceProperty REFLECTIVITY  0.000002 1.0 0.000008 1.0

/run/initialize
#
/gun/particle opticalphoton
/gun/energy 5 eV
/gun/position 90. 0 0 cm
/gun/direction 1 0 0
/opnovice2/gun/optPhotonPolar 0. deg
/opnovice2/gun/randomDirection true
#
/run/printProgress 100000
/run/beamOn 10000

Comment 2 Pierre 2022-11-15 00:34:41 CET

Created attachment 790 [details]
Mac file to reproduce the bug

Comment 3 Pierre 2022-11-15 03:26:12 CET

Created attachment 791 [details]
Log file with tracking verbose 3

I modified G4StepVerbose to display the velocity (as calculated by the steplength/steptime).
Also added a cout in G4ParticleChange::ProposeVelocity
and few cout in G4OpBoundaryProcess::DielectricDielectric to track the changes for "theStatus", "Through" and "Swap"

Comment 4 Pierre 2022-11-15 03:31:19 CET

Good day,
I managed to reproduce the bug with OpNovice2. Please find attached the mac file and the log. Look for "FresnelReflection" in the log, you will find that the next step will be "World" with the velocity of "Tank" (line 38672 and line 152779).
Sorry it took so long to get back to you.

Regards

Comment 5 Daren Sawkey 2023-10-04 17:35:00 CEST

Thanks for reporting this. I am able to reproduce it.

I believe the problem is that optical photons are erroneously reflecting into the next volume. According to (my understanding of) the Unified paper, if reflection from a facet leaves the particle momentum directed into the next volume, another facet is to be chosen. However, there are cases where a photon in volume A is reflected into volume B.

The code expects that reflected particles remain in the volume they started in, so the velocity is not set by the boundary process when reflection occurs.

Note that the velocity is set for every step. The incorrect velocity is for one step.

Also note that your original observation used the Glisur model, and see the effect with FresnelReflection. I don't know what the expectation for reflection from a facet is in the Glisur model. However, the Unified paper is clear, and it is the same code doing the reflection for both models. Using unified, the apparently incorrect reflection is LambertianReflection.

It is possible that the inverse problem occurs, namely a photon refracting back into the original volume, but I haven't yet seen it.

I believe the fix is choose another facet and do another interaction, until the photon either refracts into the next volume or reflects back into the first volume. However the code logic is pretty opaque so this is not an easy task.

Please let me know if you disagree or have comments.