How to: Geant4-DNA chemistry
============================

This chapter explains how to run a Geant4-DNA chemistry simulation in GATE 10,
how to configure chemistry lists and track-structure EM physics, and how to
use chemistry actors.

.. warning::

   The chemistry functionality is currently in alpha stage. Interfaces,
   behavior, and outputs may still change, so results should be validated
   carefully before production use.

The current chemistry support is built around the Geant4-DNA workflow:

- activate track-structure EM in the region where DNA physics should run;
- select a Geant4 chemistry list;
- optionally customize that chemistry list;
- optionally define a chemistry world and scavenger reactions;
- add one or more chemistry-aware actors;
- read summary and counter results from actor outputs.


Minimal setup
-------------

A minimal chemistry simulation needs three ingredients:

1. a standard Geant4 physics list;
2. a Geant4-DNA chemistry list;
3. track-structure EM activated in the region of interest.

The smallest practical setup looks like this:

.. code-block:: python

    import opengate as gate

    sim = gate.Simulation()

    sim.physics_manager.physics_list_name = "G4EmStandardPhysics"
    sim.chemistry_manager.chemistry_list_name = "G4EmDNAChemistry_option3"
    sim.chemistry_manager.time_step_model = "IRT"

    water_box = sim.add_volume("Box", "water_box")
    water_box.material = "G4_WATER"
    water_box.size = [10 * gate.g4_units.um] * 3

    water_box.set_track_structure_em_physics("G4EmDNAPhysics_option2")

This does not yet add chemistry scoring, but it prepares the simulation so
that Geant4-DNA chemistry can run.


Track-structure EM physics
--------------------------

Chemistry requires track-structure EM physics in the relevant region. In GATE
10 this is configured with ``track_structure_em_physics``.

The recommended full Geant4 names are for example:

- ``G4EmDNAPhysics``
- ``G4EmDNAPhysics_option2``
- ``G4EmDNAPhysics_option4``
- ``G4EmDNAPhysics_option6``

The simplest way is to configure it directly on a volume:

.. code-block:: python

    target = sim.add_volume("Box", "target")
    target.material = "G4_WATER"
    target.size = [10 * gate.g4_units.um] * 3
    target.set_track_structure_em_physics("G4EmDNAPhysics_option2")

This creates or reuses the appropriate GATE region and activates Geant4-DNA
track-structure EM there.

You can also configure it explicitly through a region:

.. code-block:: python

    region = sim.physics_manager.add_region("dna_region")
    region.associate_volume(target)
    sim.physics_manager.set_track_structure_em_physics_in_region(
        "dna_region",
        "G4EmDNAPhysics_option4",
    )

If you use a chemistry actor, the actor may also request
``track_structure_em_physics`` for its attached volume. This is a convenience
feature, but configuring the region or volume explicitly is often clearer.


Choosing a chemistry list
-------------------------

The chemistry list is selected through ``ChemistryManager``:

.. code-block:: python

    sim.chemistry_manager.chemistry_list_name = "G4EmDNAChemistry_option3"

At the moment, one simulation must resolve to one coherent chemistry list.
This means:

- the manager can request a chemistry list;
- a chemistry actor can also request a chemistry list;
- all such requests must agree.

If multiple incompatible chemistry-list names are requested, GATE stops with a
fatal error.

Current built-in chemistry list names include:

- ``G4EmDNAChemistry``
- ``G4EmDNAChemistry_option1``
- ``G4EmDNAChemistry_option2``
- ``G4EmDNAChemistry_option3``


Chemistry time-step model
-------------------------

The chemistry time-step model is configured through:

.. code-block:: python

    sim.chemistry_manager.time_step_model = "IRT"

Current allowed values are:

- ``SBS``
- ``IRT``
- ``IRT_syn``

For many aqueous-radiolysis examples, ``IRT`` is a reasonable starting point.


Customizing the chemistry list
------------------------------

The active chemistry list is available as:

.. code-block:: python

    chem_list = sim.chemistry_manager.chemistry_list

You can extend the selected built-in Geant4 chemistry list with:

- additional chemical species;
- additional bimolecular reactions;
- additional dissociation channels.

Add a reaction
~~~~~~~~~~~~~~

This is the most common customization pattern:

.. code-block:: python

    dm3_per_mole_s = (
        1e-3 * gate.g4_units.m3 / (gate.g4_units.mole * gate.g4_units.s)
    )

    chem_list.add_reaction(
        reactant_a="H",
        reactant_b="H",
        rate_constant=0.503e10 * dm3_per_mole_s,
        products=["H2"],
        reaction_type=0,
    )

Reactions are appended to the Geant4-DNA reaction table after the selected
built-in chemistry list has created its default content.

Add a species
~~~~~~~~~~~~~

You can add or modify a species by name:

.. code-block:: python

    chem_list.add_chemical_species(
        name="MySpecies",
        charge=0,
        diffusion_coefficient=1.0e-9 * gate.g4_units.m2 / gate.g4_units.s,
    )

If the molecular configuration already exists, GATE updates its configurable
properties. Otherwise it creates the necessary Geant4 molecule definition and
configuration.

Add a dissociation channel
~~~~~~~~~~~~~~~~~~~~~~~~~~

You can also add a dissociation channel:

.. code-block:: python

    chem_list.add_chemical_dissociation(
        parent="H2O2",
        products=["°OH", "°OH"],
        probability=1.0,
    )

Current model
~~~~~~~~~~~~~

The current model is additive:

- start from one built-in Geant4 chemistry list;
- then append species, reactions, and dissociations in GATE.

If you customize species, reactions, or dissociations, make sure a valid base
chemistry list is still selected.


Chemistry worlds and scavengers
-------------------------------

The optional ``ChemistryWorld`` object defines the chemistry-space box used by
Geant4 chemistry and can also carry background components and scavenger
reactions.

You usually create it from a volume:

.. code-block:: python

    chemistry_world = sim.chemistry_manager.create_chemistry_world(volume=target)
    chemistry_world.pH = 7

You can also create it explicitly from a box:

.. code-block:: python

    chemistry_world = sim.chemistry_manager.create_chemistry_world(
        translation=[0, 0, 0],
        half_size=[10 * gate.g4_units.um] * 3,
    )

Once the chemistry world exists, you can define background components:

.. code-block:: python

    chemistry_world.add_component("H2O", 55.5)
    chemistry_world.add_component("O2", 2.5e-2 * gate.g4_units.mole / gate.g4_units.liter)

and then add scavenger reactions:

.. code-block:: python

    dm3_per_mole_s = (
        1e-3 * gate.g4_units.m3 / (gate.g4_units.mole * gate.g4_units.s)
    )

    chemistry_world.add_scavenger_reaction(
        tracked_molecule="e_aq",
        scavenger="O2",
        products=["O2m"],
        rate_constant=1.74e10 * dm3_per_mole_s,
    )

The ``scavenger`` species must already be present as a chemistry-world
component. Scavenger products and tracked molecules must resolve to known
chemistry species in the active Geant4-DNA chemistry setup.

``pH`` is optional. At the moment, Geant4's stock scavenger material expects an
integer pH value.

This chemistry-world configuration is independent from
``confine_chemistry_to_volume``:

- ``ChemistryWorld`` defines the chemistry-space box and optional scavenger
  reservoir;
- ``confine_chemistry_to_volume`` defines a GATE-side policy that kills
  chemistry tracks starting outside a chosen volume subtree.


Adding a chemistry actor
------------------------

The current reference chemistry actor is ``ChemicalCountingActor``.
Further chemistry actors can be developed and made available in the future.

It is intended as a passive chemistry-scoring actor. Simulation-wide chemistry
control stays on ``ChemistryManager``.

It provides:

- chem6-like energy-loss and LET bookkeeping;
- chemistry-time species sampling;
- reaction counting;
- built-in and configured chemistry counter outputs.

A typical setup looks like this:

.. code-block:: python

    chem_actor = sim.add_actor("ChemicalCountingActor", "chem_actor")
    chem_actor.attached_to = target
    chem_actor.number_of_time_bins = 50

The actor must be attached to a volume. This volume defines the local region
for the actor logic.

If you want the actor to request DNA EM on its attached volume itself, you can
do:

.. code-block:: python

    chem_actor.track_structure_em_physics = "G4EmDNAPhysics_option2"

If you have already configured track-structure EM explicitly on the volume or
region, you can leave this as ``None``.

The actor also supports:

- ``track_only_primary``
- ``primary_pdg_code``
- ``energy_loss_min``
- ``energy_loss_max``
- ``min_kinetic_energy``
- ``let_cutoff``
- ``times_to_record``
- ``number_of_time_bins``

Example:

.. code-block:: python

    chem_actor.track_only_primary = True
    chem_actor.primary_pdg_code = 11
    chem_actor.number_of_time_bins = 50


Chemistry confinement
---------------------

Chemistry confinement is configured globally through ``ChemistryManager``, not
per actor.

If you want chemistry tracks starting outside a given volume subtree to be
killed before chemistry processing continues, configure:

.. code-block:: python

    sim.chemistry_manager.confine_chemistry_to_volume = target

where ``target`` can be either a volume object or a volume name.

This is a simulation-wide chemistry-control policy. It is intentionally kept
outside chemistry actors so that chemistry actors can remain focused on
probing and scoring.


Chemistry actor outputs
-----------------------

``ChemicalCountingActor`` currently exposes five outputs:

- ``results``
- ``molecule_counter``
- ``reaction_counter``
- ``configured_reaction_counter``
- ``configured_species_counter``

Summary output
~~~~~~~~~~~~~~

The ``results`` output contains actor-level summary quantities such as:

- ``recorded_events``
- ``chemistry_starts``
- ``chemistry_stages``
- ``pre_time_step_calls``
- ``post_time_step_calls``
- ``reaction_count``
- ``killed_particles``
- ``aborted_events``
- ``accumulated_primary_energy_loss``
- ``total_energy_deposit``
- ``mean_restricted_let``
- ``std_restricted_let``
- ``species``
- ``times_to_record``

You can read it like this:

.. code-block:: python

    sim.run()
    results = chem_actor.results.get_data()
    print(results.reaction_count)
    print(results.mean_restricted_let)

Counter outputs
~~~~~~~~~~~~~~~

The detailed chemistry histories live on the dedicated counter outputs:

- ``chem_actor.molecule_counter``
- ``chem_actor.reaction_counter``
- ``chem_actor.configured_reaction_counter``
- ``chem_actor.configured_species_counter``

All of these detailed counter outputs follow the same in-memory structure:

- a dictionary;
- keys are molecule or reaction labels;
- each value is a structured NumPy array with fields ``time`` and ``count``.

Example:

.. code-block:: python

    molecule_data = chem_actor.molecule_counter.get_data()
    reaction_data = chem_actor.reaction_counter.get_data()

    for molecule_name, series in molecule_data.items():
        print(molecule_name)
        print(series["time"])
        print(series["count"])

The counts are cumulative.

Configured counters
~~~~~~~~~~~~~~~~~~~

The configured counters let you track selected observables without relying only
on the stock Geant4 built-in counter classes.

``configured_reaction_counter`` tracks cumulative counts of selected chemistry
reaction signatures versus chemistry time.

.. code-block:: python

    chem_actor.counters.configured_reaction_counter.tracked_reactions = [
        {
            "reactant_a": "H",
            "reactant_b": "H",
            "products": ["H2"],
        }
    ]
    chem_actor.counters.configured_reaction_counter.active = True

``configured_species_counter`` tracks cumulative appearances of selected
species versus chemistry time.

.. code-block:: python

    chem_actor.counters.configured_species_counter.tracked_species = [
        "e_aq",
        "O2m",
    ]
    chem_actor.counters.configured_species_counter.active = True

This can be especially useful in scavenger workflows, where counting the
appearance of signature product species is sometimes more robust than trying to
count a direct reaction callback.

Naming conventions
~~~~~~~~~~~~~~~~~~

Chemistry names appear at several layers, and the exact Geant4-facing labels
should be treated as examples rather than a stable long-term naming contract:

- Python chemistry configuration usually uses names such as ``e_aq``, ``O2``,
  or ``O2m``;
- molecule-counter outputs often expose Geant4 runtime names such as
  ``e_aq^-1``, ``O_2^-1``, or ``HO_2°^0``.

Configured counters accept chemistry-species names that resolve through the
Geant4 molecule table, while returned output keys may follow the runtime naming
used by the underlying counter. As development evolves, these runtime labels may
change. In practice, inspect the keys returned by ``get_data()`` rather than
assuming one naming convention everywhere.

Activating and deactivating outputs
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

By default, the predefined counter outputs are active. You can deactivate them
through the normal actor-output interface:

.. code-block:: python

    chem_actor.molecule_counter.active = False
    chem_actor.reaction_counter.active = False
    chem_actor.configured_reaction_counter.active = False
    chem_actor.configured_species_counter.active = False

The same can also be controlled through the actor-owned counter objects:

.. code-block:: python

    chem_actor.counters.molecule_counter.active = False
    chem_actor.counters.reaction_counter.active = False
    chem_actor.counters.configured_reaction_counter.active = False
    chem_actor.counters.configured_species_counter.active = False

For most user scripts, using the actor output interface directly is the
clearest choice.


Complete example
----------------

This combines the main pieces in one small script:

.. code-block:: python

    import opengate as gate

    sim = gate.Simulation()
    um = gate.g4_units.um
    keV = gate.g4_units.keV

    sim.world.material = "G4_WATER"
    sim.world.size = [1 * gate.g4_units.km] * 3

    sim.physics_manager.physics_list_name = "G4EmStandardPhysics"
    sim.chemistry_manager.chemistry_list_name = "G4EmDNAChemistry_option3"
    sim.chemistry_manager.time_step_model = "IRT"

    target = sim.add_volume("Box", "chem_box")
    target.material = "G4_WATER"
    target.size = [10 * um, 10 * um, 10 * um]
    target.set_track_structure_em_physics("G4EmDNAPhysics_option2")
    sim.chemistry_manager.confine_chemistry_to_volume = target

    source = sim.add_source("GenericSource", "source")
    source.particle = "e-"
    source.energy.mono = 2 * keV
    source.position.type = "point"
    source.position.translation = [0, 0, 0]
    source.direction.type = "momentum"
    source.direction.momentum = [0, 0, 1]
    source.n = 1

    chem_actor = sim.add_actor("ChemicalCountingActor", "chem_actor")
    chem_actor.attached_to = target
    chem_actor.number_of_time_bins = 50

    sim.run()

    results = chem_actor.results.get_data()
    species = chem_actor.molecule_counter.get_data()
    reactions = chem_actor.reaction_counter.get_data()

    print(results.reaction_count)
    print(species.keys())
    print(reactions.keys())


Current limitations
-------------------

The chemistry support on this branch is already usable, but a few constraints
are worth keeping in mind:

- one simulation currently needs one coherent chemistry list;
- one simulation currently has at most one global GATE-side chemistry
  confinement policy;
- chemistry counter writing to disk is not implemented yet on the python side;
- chemistry counter merged data currently uses a simple successive-run merge;
- scavenger-related observables may require configured species counting rather
  than direct reaction counting, depending on which Geant4 callback path the
  process uses;
- some advanced reaction-counter scenarios may still be fragile;
- ``ChemicalCountingActor`` currently assumes at most one molecule counter for its
  built-in species-sampling path.

These are implementation limits, not usage mistakes.


Practical tips
--------------

- Start with ``G4EmStandardPhysics`` plus regional
  ``track_structure_em_physics`` rather than trying to replace the whole
  physics list with a DNA-only list.
- Configure DNA EM only where you need it. Track-structure physics is much more
  expensive than standard condensed-history EM.
- Keep the chemistry list request consistent across the simulation and all
  chemistry actors.
- Use ``sim.chemistry_manager.confine_chemistry_to_volume`` if you want to
  restrict chemistry spatially at the GATE level.
- If you only need summary chemistry information, ``chem_actor.results`` is the
  easiest entry point.
- If you need time-resolved species or reaction histories, use the dedicated
  counter outputs.