Add optical photon simulation

[whokion]

Requirements discussion

This is an initial draft for supporting optical photon simulation and a placeholder to collect and identify sub-tasks for implementing related codes and testing. Your additional input, refinement and contributions are welcome.

Major categories and sub-tasks:

Optical materials and properties

• [X] Determine common optical properties
• [ ] ~~Automatically populate common optical materials such as Water, Liquid Xeon and Argon~~
• [x] Add OpticalPropertyParams (constant properties, vectors of properties such as momentum dependent optical properties)
• [x] Construct OpticalPropertyView from material ID

Optical photon physics

• Add processes for generating optical photons in optical materials

• [x] Design consideration for interfaces to other physics actions (energy loss or deposition) and Photino (an object returned from processes of optical photon generation, which will be passed to the loop of optical photon propagation) — OpticalPhotonDistribution and related

• [x] Scintillation

• [x] Cerenkov

• [ ] ~~Transition Radiation~~

• Add processes for propagating optical photons in fiducial volumes/materials

• [x] Design consideration for step limiter (mfp) and interactor (final state sampling) and executing them on HPC systems

• [x] Absorption

• [x] Boundary and surface processes and associated models (GLISUR, UNIFIED)

• [x] Scattering (Rayleigh and Mie)

• [x] Wavelength shift or re-emission

I/O

• [x] Add input data structure, OpticalTrackView or OpticalPhotonInput for an optical action
• [ ] ~~Add default sensitive hits (detector Id, hit{position, energy, time}) or interfaces for user defined sensitive detectors~~

Testing and physics validation/verification

• [ ] ~~Verify physics results in the process level with respect to Geant4 - use the exiting physics validation suite~~
• [ ] Verify the hit position, energy, timing in sensitive detectors (PMT or Calorimeters) with respect to Geant4

Performance evaluation

• [ ] Evaluate performance (event throughput) with respect to the Geant4 optical photon simulation
• [ ] ~~Evaluate performance (event throughput) with respect to the Geant4+Opticks+OptiX with a test example~~
• [ ] Study performance and identify areas for improvement or optimization
• [ ] ~~Design consideration for efficient ray tracing on accelerators~~

Offloading interfaces and integration on experiment frameworks

• [ ] Design consideration for offloading interfaces (extension to the existing core and accel libraries or a separate library)
• [ ] Demonstration of integration (LZ, CalVision, LArG4)

Future work

• [ ] ~~GPU-optimized "particle guns" to determine number of internal reflections, etc.~~

Optical properties

• Partition Geant4 optical properties list into core material properties (e.g. diffraction index) and properties that are process-specific

Process	StepLimiter	Interactor
G4Scintillation		RESOLUTIONSCALE
		SCINTILLATIONCOMPONENT[1,2,3]
		SCINTILLATIONYIELD[1,2,3]
		SCINTILLATIONTIMECONSTANT[1,2,3]
		SCINTILLATIONRISETIME[1,2,3]
		[PTYPE]SCINTILLATIONYIELD[1,2,3]
G4Cerenkov		RINDEX
G4OpAbsorption	ABSLENGTH
G4OpBoundaryProcess		GROUPVEL
		RINDEX
		REALRINDEX
		IMAGINARYRINDEX
		REFLECTIVITY
		EFFICIENCY
		TRANSMITTANCE
		SURFACEROUGHNESS
(UNIFIED)		SPECULARLOBECONSTANT
		SPECULARSPIKECONSTANT
		BACKSCATTERCONSTANT
(Coated Dielectric)		COATEDRINDEX
		COATEDTHICKNESS
		COATEDFRUSTRATEDTRANSMISSION
G4OpMieHG	MIEHG	MIEHG_FORWARD_RATIO
		MIEHG_FORWARD
		MIEHG_BACKWARD
G4OpRayleigh	RAYLEIGH
	ISOTHERMAL_COMPRESSIBILITY*
	RINDEX*
	RS_SCALE_FACTOR*
G4OpWLS[2]	WLSABSLENGTH[2]	WLSCOMPONENT[2]
		WLSMEANNUMBERPHOTONS[2]
		WLSTIMECONSTANT[2]

where PTYPE = [PROTON|DEUTERON|TRITON|ALPHA|IONS|ELECTRON] if GetScintillationYieldByParticleType is used, and * in G4OpRayleigh are needed for the on-the-fly m.f.p calculation if the property RAYLEIGH is not provided.

[sethrj]

From today's discussion:

We'll almost certainly want to have a separate photon event loop that operates either in parallel or serially with the main tracking loop, immediately after along-step.

graph TB
  subgraph generate
    scintillation
    cerenkov
    initializers["buffered secondaries"]
  end
  subgraph physics
    pre-step --> propagate
    propagate --> boundary
    boundary --> reflection["reflection*"]
    boundary --> cross["cross boundary"]
    cross --> refraction["refraction*"]
    cross --> absorb
    propagate --> volume-interaction
    volume-interaction --> absorption
    volume-interaction --> scatter["scatter*"]
    volume-interaction --> reemission["reemission*"]
  end
  steps --> generate
  generate --> physics
  absorb --> hits
  absorption --> kill
  hits --> kill

Asterisks denote multiple models per action

Much of the loop should look almost identical to the celeritas one.

Note that like Geant4, generating optical photons violates conservation of energy. If we want detailed energy balance, we'd have to subtract the generated photon energy from the local dE/dx, and tally photon absorption inside SD regions to create hits. This detailed balance is not needed for known applications.

Step properties

Two cases:

• Geant4 (hadronic or EM) to Celeritas optical
• Celeritas EM to Celeritas optical

Input to photon event loop can be translated and copied to GPU (like acceleritas offloading), or transferred natively on device memory, and needs the following properties:

• Start/End position
• Start direction
• Step length
• Volume ID
• Mean (global) time
• Particle ID
• Start and end value of beta (speed/c)
• Energy loss over the step

Possibly for GPU we can pass in the geometry and a track ID so that we can copy the initial state.

The number of photons generated per step (the yield) depends on dE/dx, the particle type, and the "yield factor".

Volume Interaction

• absorption/scattering state change
• direction
• polarization
• updated energy
• occasional secondaries
• time delay

Surface crossing

Surface processes (e.g. diffraction) require materials on both sides of a surface, so we need to be able to query the new material and safely perform actions on the boundary, and potentially ignore the boundary crossing in favor of reflection, or at least be able to restore the state to the other surface side.

[sethrj]

First step:

• Generate optical photons
• Raytrace loop until we hit a PMT
• Send hits back to main loop

Optical track data:

• Geometry state
• Particle state (energy or wavelength)
• Time

Comparison problem:

• Turn off all optical surface/volume physics in geant4

graph TB
    gun["gun or external"]
  subgraph main-celeritas-loop
    scintillation
    cerenkov
  end
  subgraph photon-input-data
    scintillation-dist
    cerenkov-dist
  end

    scintillation-dist
    cerenkov-dist

    scintillation --> scintillation-dist
    cerenkov --> cerenkov-dist
    gun --> photon-primary

  subgraph photon-generator
    scintillation-dist --> scintillation-generator
    cerenkov-dist  --> cerenkov-generator
 end

  scintillation-generator --> photon-primary
  cerenkov-generator --> photon-primary

Initial tasks

• PhotonPrimary: raw single photon data as if from a gun or sampled from a generator, used to initialize a photon track state
  • Position
  • Direction
  • Energy
  • Time
  • Volume ID
  • (parent track ID? event ID?)
• Distributions are input data for sampling many photons in the "generator" (number of photons, start/stop position, step length, ...)
  • ScintillationDist @whokion
  • CerenkovDist @amandalund
• Pre-generators gather pre/post step info as part of main event loop to create "distributions" to send to GPU
  • ScintillationPreGenerator @stognini
  • CerenkovPreGenerator @stognini
• Generators populate the photon primary buffer on gpu (or as part of special loop) with individual photons sampled via models from the distributions
  • ScintillationGenerator @whokion
  • CerenkovGenerator @amandalund

Later:

• Patch generators into "models" (or some other kind of action) but these will just be action kernels responsible for offloading photons to special loop
• Buffering and flushing of distributions
• Buffering and flushing photon primaries
• Optical photon event loop
• Hit callback (should be able to reuse DetectorStepOutput etc)

Much later:

• Surface processes
• Volume processes

[sethrj]

OpticalInitStateData will be shared between:

template<Ownership W, MemSpace M>
struct OpticalInitStateData
{
    // Buffer to be filled by main loop/step collectors
    Items<OpticalDistributionData> distributions;
    // Buffer of actual optical primaries
    Items<OpticalPrimary> primaries;
// Stack allocator for initializers?
};

For each pre-generator: Input: StepState<MemSpace> step Output: push onto initializers

TODO: do we want StepStateTrackView??

class CerenkovPregenerator
{
  constructor:
 - step state
 - track ID (ASSUME we have step data from this track ID)
 - general optical material data
 - cerenkov-specific params??
 - particle track view
 - stack allocator<distributions>

size_type operator()(RNG& )
{
// Stack allocator to push onto distributions
// return number of distributions created, either zero or one?
}
};

[sethrj]

Plan for this week's optical physics hackathon: basic optical tracking loop

• collection of optical distributions
• optical core data
• initialize optical primaries/tracks from distributions
• no secondaries for this week but we'll eventually need them
• boundary crossing
• maybe absorption?

[sethrj]

The full list

Everything we need for a nicely integrated complete optical physics package

Main "core" transport loop

• [ ] define import processes for optical physics (@stognini )
• [ ] test for optical data in "import data" and create optical data structures
• [x] step collectors for capturing step data (@amandalund ); these will be added to action registry if optical physics is enabled (#1184, #1173 )
• [ ] associating optical data with the core params/state (these should be shared with the step collectors)
• [x] "pre-generators" for converting step data into optical distributions (@stognini) completed with #1153; these will be called by step collectors. Effectively the same as the start of G4Cerenkov and G4Scintillation
• [ ] buffer checking (sum of photons in queued optical distributions) to determine whether to execute the optical stepping loop

Geant4 integration

Possible integration from least to most difficult

• [ ] "offload" logic for optical photons: no optical distributions, no celeritas cernekov, but celeritas optical processes (same approach as CATS apparently)
• [ ] optical photons only from celeritas EM tracks: distributions and such from celeritas tracks, import optical physics from geant4, but geant4 does e.g. optical photons from muons
• [ ] replace Geant4 cerenkov+scintillation physics after importing data from it, so we can directly create "distributions" from the cerenkov/scintillation: no actual optical tracks ever created in geant4 (@whokion )

Note: Importing material properties (from the Geant4 material property table) and building optical data by celeritas optical processes or generators are independent from Geant4 optical processes. As we builds optical property data or params (such as angle_integral, Rayleigh_spectrum data, etc. need for optical generators or interactors) from the imported material properties directly, no extra functionality or control is needed as the material property table is built by the detector construction and independently available once the geometry is initialized - i.e., Celeritas does not import optical physics data from Geant4, but build its own optical data or params with imported material properties which are independent from Geant4 physics processes. For a benchmark application, we just need an option to toggle between Geant4 and Celeritas optical physics processes inside a user physics list.

• [ ] if using Celeritas for EM, combination of the last two

Optical stepping loop

• [x] Optical "core" state and params data (mostly, missing some bits)
• [ ] Host-only params/state wrappers (analogous to CoreParams, CoreState<M>), loading/creating some data as needed (e.g. volume -> optical material ID mapping)
• [x] Use action interface to manage kernel launching for optical stepping loop
• [ ] Initializing tracks from distributions
• [ ] Initializing tracks from primaries (for testing or for one-by-one optical track integration with Geant4)
• [ ] Initializing tracks from photon secondaries (needed only after WLS optical process is integrated)
• [ ] Optical physics process management and creating cross section data tables and optical models aka explicit actions (@hhollenb )
• [ ] Step length calculation (distance to interaction) via cross section calculation for the sum of the optical physics processes (@hhollenb )
• [ ] Propagation (easy, linear) and things like time update
• [ ] Surface physics (reflection, refraction) (@whokion nominates @amandalund )
• [ ] Discrete interactions
• [ ] Geant4 hit reconstruction (?) and/or other diagnostics

Minimal working example

load up geant4 physics, just EM tracking with celer-sim on EM particles, create but don't use cerenkov/scintillation; CPU only to make it even easier

• [x] Import cerenkov/scintillation (C/S) data
• [x] Add C/S step collector to create distributions
• [ ] Immediately launch optical tracking loop
• [ ] Initialization of tracks from distributions
• [ ] No optical physics at all: no cross sections, no distance to interaction, and no negativity of any kind... and also no consistency
• [ ] Tracking across boundaries
• [ ] Updating time
• [ ] If we hit a sensitive region, instantantly "absorb" and send hit back to geant4

[sethrj]

Tagging #1162, #1163, #1173, #1182

[sethrj]

Separate tasks from today for most minimal case:

• [x] Add optical data construction in celer-sim/Runner from imported and core_params_ @stognini
• [x] OpticalCollector will manage state/params data, constructing these from CoreParams plus optical data @sethrj
• [x] New action to perform optical stepping loop (would own optical action sequence, see #1160) (#1386 #1417 #1424)
• [ ] Optical stepping actions:
  • [x] Initialize primaries (same purpose as "track initializers" in main stepping loop) from distributions @amandalund (#1379 #1399 #1414)
  • [x] Fill empty track slots from primaries and initialize the tracks @amandalund (#1438)
  • [ ] Propagate to geometry boundary @sethrj
  • [x] Clear killed tracks, find vacancies, count remaining photons @amandalund (#1441)

Later:

• @stognini is working on optical process importing, @hhollenb has a draft ImportOpticalProcess as well
• @hhollenb to split up #1162, first goal being to add absorption

[sethrj]

TODO:

• Unify mfp/absorption length "import" accessors for different optical processes

How to treat missing data?

• Differently on a per-model basis?
• Warn, error at startup, error at runtime, opaque, transparent

Model/physics construction ordering:

• User imports optical data (from geant importer)
• Core params constructed (?)
• Create empty optical action registry
• Create optical properties/material params
• Create optical models from imported data and add to registry sequentially
• Create optical physics params from vector of optical models and other optical data
  • Ensure that optical models are consecutive and the most recently added (to the action registry)
  • Create virtual actions after the last optical model ID
  • Query optical models for absorption lengths to build cross section tables
• Create optical loop management
  • Inputs: core params, optical action registry, optical properties
  • Builds pre-generators
  • Builds generators
  • Integrates into core action loop

Physics implicit actions:

• Don't need integral rejection since we're not slowing down
• Don't include failure action since (for now) we can just require a secondary stack twice as big as the number of tracks
• Obvs. don't need msc or eloss

[sethrj]

Pseudo-code from hacking with @hhollenb :

struct OpticalModelXs
{
    ItemMap<OpticalModelId, GridId> per_model_mfp;
};

struct OpticalPhysicsData
{
    // ItemMap<OpticalMaterialId, OpticalModelXs> per_material_modelxs;

    Collection<OpticalModelXs, W, M, OpticalMaterialId> model_xs;
    Collection<GenericGridData>;
    // ....
};


void calc_step_limit()
{
    auto models = params[opt_mat];

    real_type total_xs{0};
    for (ModelId mod_id : models)
    {
        auto xs = 1/phys.calc_mfp(opt_mat, mod_id);
        pstep.per_model_xs(mod_id) = xs;
        total_xs += xs;
    }

    // Normalize per-process xs

    // Calculate step limit by multiplying into sampled MFP
    limit.step = physics.interaction_mfp() / total_macro_xs;
}

OpticalPhys()
{
    HostVal data;
    OpticalPhysicsBuilder builder{&data};
    for (mat : materials)
    {
        for (mod : models)
        {
            OpticalXsBuilder local_builder(&builder ,mat, mod_id);
            mod.build(local_builder);
            // Assert that the correct data was inserted?
        }
    }

}

class OpticalPhysicsBuilder
{
  public:
    OpticalProperties const& props();

    OpticalModelId action_to_model(ActionId) const;

    GridId add_xs(...);

  private:
    DedupeCollectionBuilder<real_type> reals_;
};

class OpticalXsBuilder
{
  public:
    OpticalXsBuilder(OpticalPhysicsBuilder* phys,
                     OpticalMaterialId,
                     OpticalModelId);

    void add_xs(...)
    {
        return phys->insert_grid(...);
    }

  private:
};


// Call once *per model*
void Absorption::build(OpticalXsBuilder& builder)
{
    CELER_EXPECT(builder.opt_mat_id() < import_.size());
    builder.add_xs(import_[opt_mat_id.get()].absorption_length);
}

[sethrj]

See also #1302 and #1314