Problem 2113

Summary: Incorrect G4OpBoundaryProcess operation when G4ParallelWorldProcess enabled
Product: Geant4 Reporter: Wouter Deconinck <wdconinc>
Component: processes/opticalAssignee: asai
Status: CLOSED FIXED    
Severity: normal CC: asai, daren.sawkey, wdconinc
Priority: P4    
Version: 10.5   
Hardware: All   
OS: All   
See Also: https://bugzilla-geant4.kek.jp/show_bug.cgi?id=1953
https://github.com/JeffersonLab/remoll/issues/158
Attachments: patch to create ParallelWorld.hh and modify OpNovice2.cc
diff of parallel disabled and enabled with /tracking/verbose 3 (first 1k lines only)

Description Wouter Deconinck 2019-01-05 22:23:36 CET 
Created attachment 530 [details]
patch to create ParallelWorld.hh and modify OpNovice2.cc

We are experiencing incorrect OpBoundary processing for optical photons when parallel (read-out) worlds are enabled (without layered mass geometry). In particular, optical photons just seem to pass straight through optical boundaries when parallel worlds are enabled even if they reflect/refract correctly when the parallel world is disabled.

We have this problem with geant4.10.04.p02 and geant4.10.05 (but it likely occurs with earlier versions too; I just haven't gone back and tested explicitly).

Downstream report here in application code: https://github.com/JeffersonLab/remoll/issues/158 (though I don't expect you to debug there, of course, just to indicate this appears 'in the wild')

Possibly related: https://bugzilla-geant4.kek.jp/show_bug.cgi?id=1953 (and duplicate therein) which was marked as INVALID due to inability to reproduce. The referred to lines 179-181 of G4OpBoundaryProcess.cc are where the behavior diverges for this bug as well. In the code below both the parallel world geometry and the parallel world physics processes are loaded, and with the same parallel world name.

Demonstration on geant4.10.05 OpNovice2 example:

cp -r geant4.10.05/examples/extended/optical/OpNovice2 /tmp
cp ParallelWorld.patch /tmp/OpNovice2
cd /tmp/OpNovice2
patch -p0 < ParallelWorld.patch
mkdir build && cd build
cmake .. && make
./OpNovic2 surface.mac

1. Add include/ParallelWorld.hh which just sets up a box volume that is equal in size to the full material world (same happens with just an empty parallel world). Content is essentially (excluding #include lines, patch attached):

class ParallelWorld : public G4VUserParallelWorld
{
  public:
    ParallelWorld(G4String parallelWorldName)
    : G4VUserParallelWorld(parallelWorldName) { };
    virtual ~ParallelWorld() { };

    virtual void Construct() {
       G4VPhysicalVolume* ParallelWorldPhys = GetWorld();
       G4LogicalVolume* ParallelWorldLog = ParallelWorldPhys->GetLogicalVolume();
       G4VSolid* box = new G4Box("boxsolid",10.0*m,10.0*m,10.0*m);
       G4LogicalVolume* boxlog
         = new G4LogicalVolume(box,0,"boxlogical",0,0,0);
       //G4VPhysicalVolume* boxphys =
       new G4PVPlacement(0,G4ThreeVector(),boxlog,"boxphysical",
                        ParallelWorldLog,false,0);
    };
    virtual void ConstructSD() { };
};

2. Add ParallelWorld and G4ParallelWorldPhysics in OpNovice2.cc with lines like (patch attached):

G4String parallelWorldName = "ParallelWorld";
detector->RegisterParallelWorld(new ParallelWorld(parallelWorldName));
physicsList->RegisterPhysics(new G4ParallelWorldPhysics(parallelWorldName));

3. Run build/OpNovice2 surface.mac and observe differences.
Before adding parallel world (or when just commenting out the relevant RegisterPhysics() line):

    Run Summary
---------------------------------
Primary particle was: opticalphoton with energy 3 eV .
Material of world: G4_AIR
Material of tank:  G4_WATER

Average number of OpRayleigh per event:   0
Average number of OpAbsorption per event: 0

Surface events (on +X surface, maximum one per photon) this run:
# of primary particles:        100000
OpAbsorption before surface:       74
Total # of surface events:      99926
Unaccounted for:                    0

Surface events by process:
  Fresnel refraction:           45943
  Lambertian reflection:        49369
  Lobe reflection:                241
  Spike reflection:               504
  Backscattering:                2881
  Absorption:                     988
 Sum:                           99926
 Unaccounted for:                   0
---------------------------------


After adding parallel world:

    Run Summary
---------------------------------
Primary particle was: opticalphoton with energy 3 eV .
Material of world: G4_AIR
Material of tank:  G4_WATER

Average number of OpRayleigh per event:   0
Average number of OpAbsorption per event: 0

Surface events (on +X surface, maximum one per photon) this run:
# of primary particles:        100000
OpAbsorption before surface:       74
Total # of surface events:      99926
Unaccounted for:                    0

Surface events by process:
  Not at boundary:              99926
 Sum:                           99926
 Unaccounted for:                   0
---------------------------------

Expected behavior: would have assumed that with this way of adding parallel worlds (which is based on examples RE01, RE04, etc) there should be no difference in how optical photons are treated at optical boundaries.
Comment 1 Wouter Deconinck 2019-01-05 22:26:25 CET 
Created attachment 531 [details]
diff of parallel disabled and enabled with /tracking/verbose 3 (first 1k lines only)

Detailed diff output with /tracking/verbose 3 in our application code (I suspect similar in OpNovice2 code) is attached. The difference start with:

G4WT0 >     ++List of secondaries generated (x,y,z,kE,t,PID):   G4WT0 >     ++List of secondaries generated (x,y,z,kE,t,PID):
G4WT0 >  Photon at Boundary!                                  | G4WT0 >  *** NotAtBoundary *** 
G4WT0 >  thePrePV:  logicMotherVol_5open_PV                   <
G4WT0 >  thePostPV: quartzRec_50001                           <
G4WT0 >  Old Momentum Direction: (-0.01808152687973002,-0.020 <
G4WT0 >  Old Polarization:       (-0.6678976456559221,-0.7437 <
G4WT0 >  New Momentum Direction: (-0.01121623129885978,-0.021 <
G4WT0 >  New Polarization:       (0.258481474387094,-0.965845 <
G4WT0 >  *** FresnelRefraction ***                            <
G4WT0 >                                                         G4WT0 > 
G4WT0 >  **PostStepDoIt (after all invocations):                G4WT0 >  **PostStepDoIt (after all invocations):
G4WT0 >     ++List of invoked processes                         G4WT0 >     ++List of invoked processes 
G4WT0 >       1) CoupledTransportation                          G4WT0 >       1) CoupledTransportation
G4WT0 >       2) OpBoundary                                     G4WT0 >       2) OpBoundary
G4WT0 >                                                         G4WT0 > 

This indeed indicates divergence due to G4OpBoundaryProcess.cc near lines 179-181.
Comment 2 asai 2019-01-07 16:30:57 CET 
Hi Wouter,

Regardless of tracking world or parallel world, there must not be any daughter volume that is touching to the boundary of the world volume. This is a special restriction that is applied only to the world volume. Could you please make your box smaller than the size of the world volume and try to see the results?

Thanks,
Makoto
Comment 3 Wouter Deconinck 2019-01-07 16:34:19 CET 
The same problem occurs when the parallel world is empty (no placement at all). I will check with a smaller G4Box in the next two hours and report back.
Comment 4 Wouter Deconinck 2019-01-07 16:58:39 CET 
I can confirm that the same issue occurs when choosing a smaller G4Box size, but there are differences depending on the size. This may help in debugging the problem:

- G4Box of 2.0*m cubed (larger than the 1.0*m cubed tank in OpNovice2 material geometry):

    Run Summary
---------------------------------
Primary particle was: opticalphoton with energy 3 eV .
Material of world: G4_AIR
Material of tank:  G4_WATER

Average number of OpRayleigh per event:   0
Average number of OpAbsorption per event: 0

Surface events (on +X surface, maximum one per photon) this run:
# of primary particles:        100000
OpAbsorption before surface:       74
Total # of surface events:      99926
Unaccounted for:                    0

Surface events by process:
  Not at boundary:              99926
 Sum:                           99926
 Unaccounted for:                   0
---------------------------------


- G4Box of 0.5*m cubed (smaller than the 1.0*m cubed tank in OpNovice2 material geometry):

    Run Summary
---------------------------------
Primary particle was: opticalphoton with energy 3 eV .
Material of world: G4_AIR
Material of tank:  G4_WATER

Average number of OpRayleigh per event:   0
Average number of OpAbsorption per event: 0

Surface events (on +X surface, maximum one per photon) this run:
# of primary particles:        100000
OpAbsorption before surface:       74
Total # of surface events:      99926
Unaccounted for:                    0

Surface events by process:
  Same material:                99926
 Sum:                           99926
 Unaccounted for:                   0
---------------------------------


Note the difference between surface events by process, which went from "Not at boundary" to "Same material".
Comment 5 asai 2019-01-07 20:20:01 CET 
Thanks. Your results are very informative. I will look into it. 
Makoto
Comment 6 Daren Sawkey 2019-09-10 22:45:31 CEST 
Hi Wouter,

Could you please try to modify 
source/physics_lists/constructors/electromagnetic/src/G4OpticalPhysics.cc
as below.

It *seems* to work for OpNovice2. I'd be grateful for any feedback in applying this to your application code.

Thanks  --  Daren.


248   if( fScintillationProcess->IsApplicable(*particle) &&
249      fProcessUse[kScintillation] ) {
250        pManager->AddProcess(fScintillationProcess);
251        pManager->SetProcessOrderingToLast(fScintillationProcess,idxAtRest);
252        pManager->SetProcessOrderingToLast(fScintillationProcess,idxPostStep);
253      }
///////////////  add these lines  ////////////////////////
254    if (fBoundaryProcess->IsApplicable(*particle) && fProcessUse[kBoundary]) {
255       pManager->SetProcessOrderingToLast(fBoundaryProcess,idxPostStep);
256     }
257   }
Comment 7 Wouter Deconinck 2019-09-10 22:50:21 CEST 
Against which version should I apply that? 4.10.04.p02, 4.10.05 or 4.10.05.p01? Maybe it doesn't matter...
Comment 8 Daren Sawkey 2019-09-10 22:59:44 CEST 
Maybe it doesn't matter. I was testing against the development version. Whatever you try, please let me know the result.
Comment 9 Wouter Deconinck 2019-09-11 02:10:49 CEST 
Your suggested change is confirmed working in my application code with 4.10.05.p01.
Comment 10 Daren Sawkey 2019-09-11 18:17:01 CEST 
Great. I will update the code. Thanks for your patience.
Comment 11 asai 2019-12-10 17:14:59 CET 
The patch is released along version 10.6.
Thank you once again for your contribution.
Makoto