Add basic optical surface physics

[sethrj]

Next steps:

• [x] Boundary interaction result for classical surface physics
• [x] Apply interaction result to optical track
• [x] Extend the boundary action so we can have it in the stepping loop twice
• [x] Geometry surface specification: #1351 @sethrj
• [x] Multiple surface physics, surface<->physics mappings, with future extension to user physics
• [x] Importing surface properties from Geant4: @stognini

Surface physics types:

• reflection/refraction
• diffuse reflective/lambertian (LZ)
• polished crystal (calvision)
• user customizable
• (GLISUR, unified models?)

Boundary physics class:

• Map surface ID -> {physics ID, local physics surface index}
• Determines whether multiple bounces are necessary etc.
• Manages temporary data (sampled local normal, ...)

End of step if encountered boundary:

• Save previous volume
• Cross boundary into next volume
• Look up physics on the boundary (https://github.com/celeritas-project/celeritas/issues/1351 )
• Select boundary physics action based on current volume and (default: just crossing; )
• Sample physics (or whatever, may scatter or absorb), return interaction result
• Apply result using geometry normal: change geometry physical state + particle state (polarization), logical state if reflection
• Cross boundary again for particles reflecting (reuse boundary crossing kernel)

Boundary interaction result:

• No change
• Internal reflection
• Refraction with some polar angle change
• Absorption
• Absorption + detection (hit scoring: call back to geant4/host)?
• Custom?
• Isotropic?

Future optimizations:

• Peek at next volume rather than crossing in and out for reflection
• Definitely profile internal reflection in fiber optics: average reflections per volume before leaving

Future work:

• Special volume physics for total internal reflection

[hhollenb]

After reviewing some of the content of G4OpBoundaryProcess.cc, there seems to be a rough pattern between most of the models:

while dot(newMomentum, globalNormal) <= 0:
    calculate reflectivity and transmissivity
    sample reflection / absorption / transmission
    if reflection:
        choose reflection type
        calculate reflection
    else if absorption:
        do absorption based on material efficiency
    else if transmission:
        sample facet normal based on momentum and surface type
        sample new momentum direction from distribution
        determine new polarization from momentum proposal
        update old momentum and polarization with new quantities for next round of sampling

The main differences between the surface models enters for:

• How they calculate probabilities for transmission, absorption, reflection
• Possible reflection models (spike, lobe, back scattering, or Lambertian)
• Proposed momentum and polarization for transmission

The major complication is the while dot(newMomentum, globalNormal) <= 0 loop, which means we will have to make calls to model specific calculations perhaps multiple times. Here, globalNormal refers to the normal of the surface of the new volume, so the precondition for the boundary interaction would be dot(oldMomentum, globalNormal) < 0.

The question becomes how do we want to do actions with perhaps lots of branching for different models. Some of the model specific calculations do involve rejection sampling while others don't. I'm not sure what the tradeoffs look like yet nor how often we'll have photons undergoing boundary interactions each step, nor how many of them will be using the same models. I'll also do a cursory investigation to see if the loop is indeed necessary or if there's a way to ensure we don't sample invalid directions.

[hhollenb]

Going through the code a bit more, I think the looping is meant to simulate the chance of multiple reflections / refractions at a single surface crossing. If there's a rough surface, then reflecting off of a jagged edge (1st reflection) could mean that the photon's direction is still towards volume 2. Geant then keeps simulating these possible interactions until it has a momentum moving away from the surface.

I think with the way we do boundary crossings in Celeritas, we can just sample the interaction and update it's momentum and polarization at the boundary, updating the surface sense only if it has properly moved across or not.

[sethrj]

Ah excellent, yes you are quite right, we can support reentry like that: we should be able to handle a lot of direction changes "on" the surface without losing tracks or misbehaving.

[hhollenb]

Reference for the UNIFIED surface model in Geant4:

A. Levin and C. Moisan, "A more physical approach to model the surface treatment of scintillation counters and its implementation into DETECT," 1996 IEEE Nuclear Science Symposium. Conference Record, Anaheim, CA, USA, 1996, pp. 702-706 vol.2, doi: 10.1109/NSSMIC.1996.591410.

[hhollenb]

Here's my current notes on how the optical surface physics can be separated out. There's quite a few branching points for model choices, so I'm not sure if it'd be better to have them separated out into subactions or have just a few monolithic actions.

Determine local facet normal:

• Default: use global surface normal
• UNIFIED: sample facet normal from UNIFIED distribution based on sigma_alpha
• GliSurface: if polished, then sample based on a smear factor

Calculate reflectivity:

• Default: from reflectivity material property vector. If available should always be sampled first.
• Analytic: calculate from Fresnel equations based on index of refraction. Possible to use if reflectivity MPV is missing. For dielectric-metal interfaces, this should be called after the 1st iteration instead of the MPV calculation, and it will also determine polarization fractions for TE and TM modes.
• Dichroic: Lookup transmittance in dichroic material property vector. Only reflection and transmission, so reflectivity is one minus transmission.

Select interaction type:

Note that in Geant, reflectivity doesn't actually refer to the material's reflectivity, but rather one minus the absorption coefficient. It encompasses both reflection and refraction.

• Absorption: photon is absorbed on the surface. It's recorded in the sensitive detector if it passes an efficiency sampling.
• Transmission: photon is transmitted through the surface without any change. Transmittance is calculated only from the material transmittance MPV and doesn't correspond with the transmittance from Fresnel equations.
• Reflection: photon undergoes the surface model reflection / refraction

Snell's Laws:

If the surface allows transmission, use Snell's laws and Fresnel equations to determine if the photon is refracted or reflected.

• Refraction: Apply model specific refraction interactor
• Reflection: Select reflection mode

Select reflection mode:

• Specular spike: geometric reflection off of global normal with probability C_ss
• Specular lobe: model specific reflection off of facet normal with probability C_sl
• Diffuse / Lambertian: diffuse reflection following Lambert's law with probability C_dl
• Back scattering: directly back scattered with probability C_bs

Finishes and models

UNIFIED model finishes:

• Dielectric-Dielectric
  • Polished: supports refraction; always do lobe reflection with global surface normal
  • Ground: supports refraction; lobe reflection is spike off of local facet
  • (Polished/Ground) Front Painted: represents a wrapped surface. No refraction. Reflection is specular spike / Lambertian.
  • (Polished/Ground) Back Painted: implicitly cross into an air-gap first, then treat the air-gap as the corresponding front-painted.
• Dielectric-Metal
  • Polished / Ground: same as corresponding dielectric-dielectric front painted

Look Up Table Models:

• LUT: lookup table for reflection, no refraction
• LUT/DAVIS: lookup table for reflection and refraction

[sethrj]

Next step:

• Hardcode global kREFLECTIVITY for all surfaces (to be moved into and expanded to SurfacePhysicsParams)
• Implement surface interaction with reflection + transmission

Later work:

• Import G4 logical surfaces
• Import surface material tables
• SurfacePhysicsParams (includes reflectivity, finishes, models, etc)
• During importing require optical surfaces don't have finish (GetFinish) to avoid LUT stuff for now

Question:

• Does geometry handle mapping of (entering/exiting pv/lv) to surface?

[sethrj]

converted to mermaid:

flowchart TD
    A{{Determine Local Normal}} --> B1[Default: Global Normal]
    A --> B2[UNIFIED sigma_alpha]
    A --> B3[GliSurface]
    B1 --> C{{Calculate Reflectivity}}
    B2 --> C
    B3 --> C
    C --> D1[Default: From Reflectivity material property vector]
    C --> D2[Analytically from Fresnel equations]
    D1 --> E{{Select interaction type}}
    D2 --> E
    E --> F1[Absorption]
    E --> F2[Transmission: no change]
    E --> F3{{Snell’s law to determine reflection / refraction}}
    F3 --> G1[Model specific refraction]
    F3 --> H{{Select reflection mode}}
    H --> I1[Back scattering]
    H --> I2[Specular spike reflection]
    H --> I3[Diffuse Lambertian reflection]
    H --> I4[Model specific lobe reflection]

[hhollenb]

General surface properties to import:

• model
• reflectivity MPV
• transmittance MPV
• optical material 1 ID
• optical material 2 ID

GLISUR model:

• polish
• finish

UNIFIED model:

• sigma_alpha
• Spike reflection probability
• Lobe reflection probability
• Backscattering reflection probability
• finish

[sethrj]

From the discussion with @stognini and @hhollenb

Geometry surface mapping

@sethrj

Done via GeoParams:

• {pv, pv} -> border (interface) -> SurfaceId
• {lv} -> skin (boundary) -> SurfaceId

Optical surface physics

• SurfaceId -> {optical surf model, within-model index}

Import data

@stognini

Analogies

Instead of having an intermediate {pv,pv} -> surface -> surface properties (or optical surface properties), just map directly:

Volume: Volume ID -> Material ID == G4Material Surface: {skin or border via pv/lv IDs} -> Surface ID == G4SurfaceProperty/G4OpticalSurface

[sethrj]

From Friday morning's discussion:

• Support instantiating boundary crossing action twice
• Add outward surface normal
• Add a "paused" (nomenclature TBD) action ID for being on a surface
• Regular boundary crossing for optical photons (and maybe main loop) become pre-post
• Add post-post action for reflecting post-surface physics

[sethrj]

Notes from our surface discussion today (@hhollenb , @Rashika-Gupta , @amandalund , @stognini , @whokion ):

Independent of model:

• Reflectivity (depends on neighbor materials, user provided properties)
• Applying reflection/refraction/transmission

Surface physics model dependent:

• Surface normal calculation
• Surface interaction type calculation (absorption/refl/refra)
• Refraction model kernel (some overlap)
• Reflection model kernel (unified is a superset)

LUT: (energy, angle to surface) -> (exit angle PDF)

---

Surface physics state data:

• Pre/post material
• Surface ID from pre/post volume
• Surface physics model
• Sampled surface normal
• Calculated reflectivity (not random)
• Sampled interaction type

Surface physics for first step across a surface:

• Save pre-crossing logical/physical volume temporarily (and material ID)
• Cross boundary and relocate
• Use post-crossing logical/physical volume to find surface ID and material ID
• Use surface ID to select surface physics model
• Call EACH surface model to sample surface normal (masking based on particle's current surface model)
• Calculate reflectivity at surface (based on photon's direction and sampled surface normal); this is a single kernel, not per model: only depends on user-input grid or fresnel equation model (In the future, we could make this model-dependent to support dichroic surfaces?)
• Call EACH surface model to select and do the interaction type
  • Return absorption
  • Select, sample, apply reflection or refraction -> update dir, pol

Surface physics data:

• Reflectivity

Surface model functions:

• Calculate surface normal
• Select interaction type
• Sample reflection/refraction (uses helper classes to reduce overlap across models)

IGNORE TRANSMISSION FOR NOW: seems like no one uses it; lhcb has 1.0 for everything.

[sethrj]

Next steps:

1. Add surface/volume properties to core and optical params (@amandalund )
2. Extend "init boundary" to look up surface ID and kill (no need for SD check) if not found (@hhollenb )
3. Add skeleton physics params that includes surface normal, surface ID
4. Hardcoded smooth surface, scalar reflection probability, action for refract/reflect/transmit
5. Grid reflection (requires moving data from input to device, loading surface data)
6. Surface roughness: require careful "on the surface multiple reflection" logic
7. Multiple roughness models (polished + uniform smearing)
8. Multiple reflection models (?)
9. Interaction models
10. Layers? @hhollenb @amandalund @rahmans1 @amandalund @stognini @Rashika-Gupta @whokion

[hhollenb]

@stognini Here's a mapping of Geant4's models and finishes onto our layered surfaces.

Celeritas options:

Facet normal calculators:

• Polished: facet normal is same as surface normal (no sampling)
• Glisur: uniform smearing
• Unified: gaussian (sigma alpha)

Reflectivity calculators:

• Grid: grid-based lookup (kREFLECTIVITY)
• Analytic: from Fresnel equations based on index of refraction of incident and target volumes
  • Analytic (surface): uses surface material property table's index of refraction as the target (kREALRINDEX, kIMAGINARYRINDEX)

Interaction models:

• Only reflection
• Dielectric-Dielectric
• Dielectric-Metal

Reflection modes (4 real types):

• Specular spike: reflection wrt surface normal
• Specular lobe: reflection wrt facet normal
• Back-scattering: opposite momentum, polarization
• Lambertian: diffuse wrt Lambert's law

Geant4 input:

Model	Surface Type	Finish	Layer	Roughness	Reflectivity	Interaction	Reflection Mode
glisur	dielectric-dielectric	Polished	1	polished	grid / analytic	dielectric-dielectric	spike
		Ground	1	smear	grid / analytic	dielectric-dielectric	lobe
	dielectric-metal	Polished	1	polished	grid / analytic	dielectric-metal	spike
		Ground	1	smear	grid / analytic	dielectric-metal	lobe
unified	dielectric-dielectric	Polished	1	polished	grid / analytic	dielectric-dielectric	spike
		Ground	1	gaussian	grid / analytic	dielectric-dielectric	any
		PolishedFrontPainted	1	polished	grid / analytic (surface)	only reflect	spike
		GroundFrontPainted	1	gaussian	grid / analytic (surface)	only reflect	Lambertian
		PolishedBackPainted	1	gaussian	analytic (surface)	dielectric-dielectric	any
			2	polished	grid	only reflect	spike
		GroundBackPainted	1	gaussian	analytic (surface)	dielectric-dielectric	any
			2	polished	grid	only reflect	Lambertian
	dielectric-metal	Polished	1	polished	grid / analytic	dielectric-metal	spike
		Ground	1	ground	grid / analytic	dielectric-metal	any

Default to choosing grid for reflectivity if available (effectively a user override of the analytic calculations). We probably will append the surface r-index grids to the end of optical material grids, so all the analytic calculations just take two optical material IDs.

[sethrj]

That's amazing @hhollenb ! Thanks for that.

[sethrj]

1. Volume surface selector: map before/after geo state to surface ID (@hhollenb )
2. Surface model class (without currently inheriting from a step action!) that returns vector of surface layers (currently type aliased to SurfaceId) and surface physics params class that maps surface+layer -> model for model in {rough, refl, interact} (@sethrj ) [just CPU only mapping, no device data]
3. Input definitions: unfolding the "bundles" of Geant4 surface physics into the individual distribution types (don't worry yet about user-defined stuff)
4. Primary generator for testing (@amandalund)
5. Add GPU data to surface physics: surface normal, surface ID in state
6. Actually implement a surface model class: hardcoded smooth surface (extend SurfaceModel to inherit from step action)
7. (a, b, c) Add scalar reflection probability, action for refract/reflect/transmit
8. More physics (see previous list); note that "calculators" can be implemented asynchronously with previous steps for the most part
9. Add layers

[sethrj]

NAMING agreed upon by @sethrj @stognini @amandalund @Rashika-Gupta

Roughness types: polished, smear [used by GLISUR], gaussian [used by UNIFIED] Reflection types: grid [user defined], analytic [Fresnel equations] Interaction types: DielectricDielectric, DielectricMetal [both can use the same "interaction form" data struct]

[hhollenb]

Follow-up from our discussion today on possible interaction results:

InitBoundaryAction set as post-step action when the track is geometry step limited and has been already moved to the boundary. Geometry state says it is in the "pre"-volume. The executor is crosses the boundary and is responsible for initializing relevant boundary crossing state. Variable marked immutable will not be changed during a surface crossing since we'll need to know which volume would be re-entrant for efficient geometry traversal.

• Pre-volume instance ID [immutable]
• Post-volume instance ID [immutable]
• Surface ID [immutable]
• Global surface normal [immutable]: this is the geometric normal which local facet normals will be sampled from
• Physical sense: which volume the photon is physically in. This doesn't match the geometry sense to avoid re-traversing the geometry (When layers are supported later, this can be the layer index)

The layer surface normal will depend on the photon direction. For local calculations we use the convention that the facet normal $\hat{n}$ and the track momentum $\vec{p}$ satisfy

$$ \hat{n} \cdot \vec{p} < 0 $$

The global surface normal $\hat{n}_g$ is required to know which direction the layer was crossed (towards the pre-volume or towards the post-volume), so the layer normal $\hat{n}_l$ is just the global normal following the above convention:

$$ \hat{n}_l = \begin{cases} \hat{n}_g & \hat{n}_g \cdot \vec{p} < 0 \\ -\hat{n}_g & \text{else} \end{cases} $$

The local normal is then sampled as usual wrt the layer surface normal:

$$ \hat{n} = \text{roughness}(\hat{n}_l ; \sigma \dots) $$

The track will then run through its appropriate roughness, reflectivity, and interaction actions and get an interaction result.

• Absorbed: track is killed at next pre-step action (may want to invoke SD call in post post-step action?)
• Reflected: Physical sense is not changed (remains in the same layer)
• Refracted: Physical sense is changed (moves to different layer)
  • The relative motion between layers is determined by the pre-interaction momentum and the global normal:

$$ \text{layer} += \text{sgn}(-\hat{n}_g \cdot \vec{p}) $$

The new momentum is $\vec{p'}$. If it's in a pre-volume or post-volume and moving away from the surface, then the track is considered to be done with surface physics (cases 1 and 2) and it may resume normal tracking.

1. If in the pre-volume layer and $\hat{n}_g \cdot \vec{p'} > 0$ then reentrant into pre-volume
2. If in the post-volume layer and $\hat{n}_g \cdot \vec{p'} < 0$ then entrant into post-volume
3. Otherwise repeat boundary action

This decision can be done in the InteractionApplier or in a separate PostBoundaryAction.

Using traversal direction

To simplify some of the logic a bit further, introduce the traversal direction $D$:

$$ D = -\text{sgn}(\hat{n}_g \cdot \vec{p}) $$

An enum can then be used to describe the traversal direction of the track

• FORWARD: If $D = +1$ then the track is headed from the pre-volume towards the post-volume
• REVERSE: If $D = -1$ then the track is headed from the post-volume towards the pre-volume

The layer surface normal is then just $\hat{n}_l = D \hat{n}_g$, so the local normal is

$$ \hat{n} = \text{roughness}(D \hat{n}_g; \sigma \dots) $$

Similarly defining the new direction as $D_{new} = -\text{sgn}(\hat{n}_g \cdot \vec{p'})$ the post-step logic becomes

// traversal direction pre-interaction
int D = track.surface_physics().traversal_direction();
auto& layer = track.surface_physics().layer();
if (result.action == Absorbed)
{
    track.sim().status(TrackStatus::killed);
    return;
}
else if (result.action == Refracted)
{
    layer += D;
}
else
{
    // no layer change for reflected
}

track.geometry().set_dir(result.direction);
track.particle().polarization(result.polarization);

// traversal direction post-interaction
int D_new = track.surface_physics().traversal_direction();

if (layer == 0 && D_new == REVERSE)
{
    // re-entrant in pre-volume: switch to normal tracking
}
else if (layer == 1 && D_new == FORWARD)
{
    // entrant in post-volume: switch to normal tracking
    // can replace with (layer == num_layers) later when implemented
}
else
{
    // still performing boundary crossing, do nothing
}

[sethrj]

From today's planning: the surface physics actions will proceed through multiple stages.

1. Volume finding is part of the along-step (linear) kernel; it establishes the distance-to-boundary and moves the particle to the surface.
2. Surface crossing temporarily[^1] stores the pre-surface volume, crosses the surface (i.e. relocates into a new volume), and uses the pre- and post- volumes to determine the surface ID, which is then stored on the state. Even though the geometry says the particle is in the new volume, we initialize a supplemental "sense" (which I think will later be extended to be the layer index) that tells us we haven't yet done the surface crossing physics: it starts in the "pre-crossing" state.
3. The roughness/reflection/interaction kernels from all the models are applied in the correct order to get an updated sense and direction. Reflection preserves the sense (we haven't yet left the old material), while refraction changes it (we have entered the new material).
4. If our sense matches the current direction (i.e., if the old volume is left, on the left side of the surface, heading leftward) then we "commit" the surface change. If our sense is "pre-crossing" (left) then we do another relocation to update the volume to the one from step 1. Committing clears the surface ID (but it doesn't change any ORANGE surface information).

If we didn't complete the surface physics crossing in step 4, then the surface flag stays active and disables the rest of the stepping loop until the next time we try surface physics, effectively repeating step 3. [^2]

Stage	Action	Geometry volume	Interacting surface layer	Sense	Direction	Action
1	Find volume	Pre	—	—	Right	Cross and find
2	Cross and find surface	Post	0	Pre	Right	Post-step A
3a	Apply smooth reflection	Post	0	Pre	Left	Post-step B
4	"Commit" surface change	Pre	—	—	Left	Continue
3b	Apply rough reflection	Post	0	Pre	Right	Post-step A
3c	Refract	Post	0	Post	Right	Post-step B
4	"Commit" refraction	Post	—	—	Right	Continue

Here is the optical stepping loop as we've sketched it:

The next step is for @hhollenb to add a surface physics state that has the surface ID and sense (maybe we should go ahead and just call it layer ID??), add a kernel to "commit" (which for now just clears the on-surface ID) and mask the rest of the loop so that tracks on the surface are ignored.

[^1]: Once we have touchable-to-index save/restore, we could store the pre-surface crossing to reduce the cost of relocation when re-entering an old volume. [^2]: We might want to configure an option to repeat step 3 a couple times per physics step, since it'll be relatively cheap kernels compared to thrust algorithms and other actions that require synchronization.

[sethrj]

Lol jinx @hhollenb

[hhollenb]

I think it's best if we just use layer ID for the sense and have

• pre-volume: layer = 0
• post-volume: layer = 1

and then it's basically just the special case for 1 surface / 2 layers. I think this is basically what Orange does for surface sense anyways?

[sethrj]

Yes agreed, layer ID instead of trying to define a new sense. Forward/reverse for +1/-1 index. ORANGE has an inside/outside sense but it's universally defined by the quadric surface's coefficient, rather than before/after. And in fact we may need to think hard about how the traversal direction/layer index differ based on whether we're going from A->B versus B->A. (Does Geant4 behave rationally for skin surfaces with painting? Do users only define back/front painted with border surfaces?)

[hhollenb]

I believe yes, Geant4 does behave rationally for skin surfaces with painting. In the G4OpBoundaryProcess.cc most of that weird if statement chaining with polished/groundbackpainted is to handle whether the photon is entering on the gap side or on the wrapping side.

For ordering, let me hash out some more details. I'm thinking a simple way to approach it would be to have a per-track scratch array of layer IDs that's filled during the InitBoundaryAction s.t. layer_map[0] = pre_volume and layer_map[n] = post_volume.

Alternatively something like a LayerTraversalView that encapsulates all of the traversal logic for a surface and has some minimal +1/-1 directions for each trick might be more memory efficient.

[sethrj]

Alternatively something like a LayerTraversalView that encapsulates all of the traversal logic for a surface and has some minimal +1/-1 directions for each trick might be more memory efficient.

Ooh I like that idea a lot. The surface state will then have:

StateItems<SurfaceId> surface;
StateItems<LayerId> layer;
StateItems<Direction> direction; // Forward for pre-> post, reverse post->pre
StateItems<Direction> entry; // Forward if a one-directional boundary surface, or if exiting a volume with a skin; reverse if entering a volume with a skin

Thoughts:

• Maybe instead of putting the surface ID in a surface physics state, this should actually be a surface traversal state, and then the surface physics (which would store the geometry normal etc) would get the SurfaceLayer from it
• Maybe direction should be entering/exiting? Or have another enum for entry?
• We definitely want to have some view built on top of the traversal state, because we could potentially compact some of these variables to take up less memory (direction is each one bit of data but takes a byte by default; layer we can limit to a few bits since no one will have more than a couple of layers)

[hhollenb]

In classic physicist fashion, I think by describing them as coordinate frames we just need a single direction flag.

Surface Record Frame

The surface record frame is how the layers are directly stored in memory. For an interface surface from A to B they satisfy

• An array of N layer records (contains the specific physics of each layer)
• An array of N-1 interstitial optical materials (we'll implement this later)
• LayerId 0 always corresponds to volume A
• LayerId N always corresponds to volume B
• The layer physics record at LayerId x separates the layers x and x+1.

With the convention that forward = +1 and reverse = -1, the layer physics (in the surface record frame) is just given by

LayerId{x + static_cast<int>(d)};

For boundary surfaces, LayerId 0 is the logical volume and LayerId N is the arbitrary exterior volume.

Track Local Frame

Since tracks might be incident on different sides of a layer, we'll store their layer in a track local frame. For a given pre- and post-volume the track local data consists of:

• SurfaceId the track is currently on
• LayerId of in the track local frame
• Direction of the surface orientation (forward or reverse)
• LayerId 0 always corresponds to the pre-volume
• LayerId N always corresponds to the post-volume

We can then transform between these coordinates by:

LayerId surface_record_layer = this->layer();
if (surface_orientation == Direction::reverse)
{
    surface_record_layer = this->num_layers() - surface_record_layer;
}
return surface_record_layer;

Therefore we can just initialize the track local layer to 0, but we will also need to supply the surface orientation. This works for boundary surfaces by:

• boundary is pre-volume: surface orientation forward (inside out)
• boundary is post-volume: surface orientation reverse (outside in)

Tracking can then be done only in the track-local frame so our InitBoundaryAction and PostBoundaryAction can work directly with pre- and post-volumes. A track is leaving the surface if its layer record is out of bounds:

LayerId{this->layer() + static_cast<int>(direction)} >= this->num_layers();

Since the global normal $\hat{n}_g$ is in track local coordinates, it always points from post-volume into pre-volume. The direction given by the track's momentum is

$$ D' = -\text{sgn}(\hat{n}_g, \vec{p}) $$

which is in track local coordinates as well, and thus can be used to immediately test if moving off the surface or converted into the corresponding surface record layer.

[sethrj]

Sounds good to me! One point though is I think LayerId should be only the "global" reference frame corresponding to the data layout (i.e., a SurfacePhysicsRecord would have an ItemMap<LayerId, LayerRecordId>). The indexing from the local reference frame should be either a separate type (and your transformation would convert between them) or simply a size_type for "number of surfaces crossed from pre- to post".

Also we need to be careful about indexing and maybe naming "layer" something different? So if we have an interface surface from a volume with material A-> to material B defined with optional coating material C and wrapping W, would a photon going from left to right see the following "layers" :

Layer	Regular	Coating	Wrapping	Coated then wrapped
0	A	A	A	A
1	B	C	air?	C
2	—	B	W	air?
3	—	—	B	W
4	—	—	—	B

and "layer record 0" would have properties about the interface between 0 and 1? I'm worried that "layer" might be confusing here because in the simple case there really aren't layers (coatings/wrappings), just an interface. Maybe surface material or surface zone or something?

[sethrj]

Discussed yesterday:

• SubsurfaceMaterial : locally indexed 0 -> N while crossing a surface, includes volumetric materials at entry and exit (integers)
• SubsurfaceInterface: roughness and other parameters (half integers) defined by user
• User input: models/data for N - 1 interfaces

[sethrj]

Surface physics models will look like:

namespace optical
{
class SurfacePhysicsModel : virtual public SurfaceModel,
                            virtual public OpticalStepActionInterface
{
    using OpticalStepActionInterface::order;
    using OpticalStepActionInterface::step;
    using SurfaceModel::get_surfaces;
};

class SmearRoughnessModel final : public SurfacePhysicsModel,
                                  public StaticConcreteAction,
                                  public ParamsDataInterface<SmearParamsData>
{
  public:
    // Construct with action ID, name, and inp::SmearRoughnessModel

  private:
    std::vector<SurfaceId> surfaces_;
};
}  // namespace optical

[sethrj]

Addendum to last post:

• Surface physics models won't be concrete actions necessarily; they'll just have the same step interface and we'll embed them inside a "surface physics" action.
• Default surface physics will be determined by an empty SurfaceId{} in the input models, and changing that to an "extra" surface ID internally to the surface physics (see https://github.com/celeritas-project/celeritas/pull/1876/files/246d56b78e5ba582354737b2b73af8e57d6e8a51#r2254360609)

[sethrj]

@hhollenb Would you mark this issue as "complete" with your little fixes from validation? 😄