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.
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.
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
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.
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".
Thanks. Your results are very informative. I will look into it. Makoto
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 }
Against which version should I apply that? 4.10.04.p02, 4.10.05 or 4.10.05.p01? Maybe it doesn't matter...
Maybe it doesn't matter. I was testing against the development version. Whatever you try, please let me know the result.
Your suggested change is confirmed working in my application code with 4.10.05.p01.
Great. I will update the code. Thanks for your patience.
The patch is released along version 10.6. Thank you once again for your contribution. Makoto
