Support for materials in AdePT

[amadio]

Before we can have complex geometry and physics models added to AdePT, we need to add their dependencies. Material information is required both to describe detector geometry and to fetch and/or calculate cross sections for various processes. Therefore, we should look at what features are required to support at a minimum creating material data in memory with a simple API or when geometry is loaded from GDML, which would require at least support for describing composition based on elements and isotopes, and having density as auxiliary property. We should also identify what other properties are required for our initial prototype.

Additionally, we need to identify which features will be needed on the GPU, like what data needs to be moved over, what queries physics models make based on materials, etc.

It would be nice if we could have isotopes and elements built in using data from AME2016, IUPAC, NIST (linked below), rather than require these things to be defined in the input file (although customization should be possible).

In this issue, we should discuss:

• The in-memory data layout for materials
• The API we shall use to create/modify/access materials
• The connection (probably via a single numerical id) of materials and other components like geometry and physics models
• Built in support for elements and isotopes (maybe also reference materials) from various data sources
• Potential file formats for material data in AdePT

Design Goals:

• Creating/modifying/accessing material properties should be done without exposing the in-memory data layouts of materials and their properties to other parts of the code, like geometry and physics. This is important to make it possible to have different data layouts on CPU/GPU if necessary, and to let us move only the parts of the data for material properties that are required on the GPU. This also allows the data layout to be changed if needed without affecting other parts of the code.
• It should be possible to easily extend the system to handle new material properties, as new properties may be required by new physics models.
• The API should be reentrant (at least the part for accessing material properties, but ideally the whole of it).
• The API should not depend on publicly exposed variables or singletons.
• The code should be independently testable.

References

• GDML Manual
• Geant4 Materials
• Standard Atomic Weights
• The 2016 Atomic Mass Evaluation
• NIST Physical Reference Data
• G4HepEm (already has required support for materials on CPU/GPU for EM models)

[hahnjo]

My understanding was that material properties are precisely one of the things that we do not need during simulation: G4HepEm will calculate the cross sections during initialization (with data queried from Geant4) and that's it. Unless I'm missing something new in here, this was discussed several times (including yesterday).

[amadio]

Not all models calculate cross sections ahead of time and use just that. Some models have per element cross sections and the material composition and density information is still needed to compute the full cross section during the event loop. Some materials have, say, optical properties, that may be needed for models that handle surface interactions of optical photons, etc. I don't think that sending only cross sections and couple id to the device will be enough for everything, so it's worth having a discussion here about which models need what and use that as guidance to develop the required code in AdePT to handle materials. VecGeom doesn't handle material information in GDML files as far as I know, so we will need a solution that doesn't involve reading the input multiple times. We will need an in-memory geometry description that contains the material information (even if just a material id per physical volume), so the strucutre VecGeom currently provides is not sufficient.

[agheata]

@amadio in principle you are right, it is possible to do integration of differential cross sections on the fly, which can be rather expensive, and require material information. However, as @hahnjo mentioned, I also understood that this is handled in a different way in G4HepEm. It is not excluded that we might still need some material info, but we have to understand exactly in which context before integrating them. I think we should involve Mihaly (@mnovak42) in this discussion as he surely can bring extra clarification.

[WitekPokorski]

We should certainly not take into account the optical photons in our prototype. We don't aim to develop a general simulation toolkit, but a specialised EM calorimeter simulation code.

[amadio]

In the end I believe we do want to have a full simulation toolkit, but I'm not saying we need to add support for everything now. We just need to think about what we will need later on and not make it hard to support it by keeping things like extensibility in mind as we move along with the implementation of the prototype. For LHCb, for example, optical photons may be important to have working on GPU.

[WitekPokorski]

Geant4 already uses the NIST database. We should also not worry about "It should be possible to easily extend the system to handle new material properties, as new properties may be required by new physics models." for the time being. We want to demonstrate that EM shower can be simulated on GPUs with a set of hard-code physics processes. Once we demonstrated that, we will worry about adding all the flexibility.

[WitekPokorski]

For optical photons, there is the Opticks library, which is dedicated to simulated those on GPUs.

[amadio]

it is possible to do integration of differential cross sections on the fly, which can be rather expensive, and require material informatio

I'm not referring to differential cross sections, just cross sections that depend on Z, and for a given material, even the total cross section depends on the material information (i.e. for PbWO₄, photon absorption is different for each element in the material). Maybe EM models just pre-compute things for the material as whole, but some models might need split information to generate secondaries by choosing an element in the material for the interaction. For example, for ionization, the binding energy of the electrons for each element may be required.

[amadio]

For optical photons, there is the Opticks library, which is dedicated to simulated those on GPUs.

Do you mean that we will not consider supporting optical photons even later on? I do not mean we should support them now, but we should not forget about them if we intend to have support later on.

[mnovak42]

There are several complete material description libraries (the Geant4 one or the one I developed for GeantV). These are complete: includes isotope data bases, element material compositions including all data and functionalities required for material composition, description and compute all correction factors required by physics models. These complete set of information is needed only in the case of initialisations in HEP EM physics of e-/e+ and gamma while only a handful information, regarding the materials and their element compositions, is needed at run time. G4HepEm contains and provides material-cuts, material and element related data structures containing the bare minimum of data needed at run time. I would like to invite you to check the corresponding part of the documentation (especially G4HepEmMatCutData , G4HepEmMaterialData and G4HepEmElementData).

So the material-cuts couple index, material index are the only information that G4HepEm needs to obtain any information that is needed at run time form outside. It stores all other information regarding materials, cross sections, other integrated quantities and even much more.

[amadio]

We want to demonstrate that EM shower can be simulated on GPUs with a set of hard-code physics processes. Once we demonstrated that, we will worry about adding all the flexibility.

If we don't worry about keeping things flexible now, when we decide that we need it later on it may be too late. That is, we may write code now that works, but will be very hard to change later, like it's hard to move Geant4 to GPU now.

[WitekPokorski]

I think the way forward to succeed on GPUs (and to speed up things on CPUs) is to have specialised library, which deal only with relatively narrow domain. This is what Mihaly's is doing with G4HepEm and this is what I believe we should do in AdePT. We know that all the simulation applications that succeeded on GPUs were limited to one domain (optical photons, low-E neutrons, medical applications, etc). Optical photons, in particular, are quite a separate domain, so I doubt we will ever manage to integrate them. But again, we still need to prove that we can simulate at least something on GPUs, before discussing possible extensions. I agree we should keep in mind the potential functionality we will require, but it should be the minimum needed to have a self-contain module we can plug into the current full simulation (the way fast simulation models are plugged into Geant4 now). The main reason why we can't move Geant4 to GPUs now is because it is too general and not because it is not flexible enough.

[amadio]

These complete set of information is needed only in the case of initialisations in HEP EM physics of e-/e+ and gamma while only a handful information, regarding the materials and their element compositions, is needed at run time.

Thank you for the pointers to documentation, I will check it out. Is this "handful of information" needed at runtime moved to the GPU by G4HepEm? The idea with this issue is to identify what work is still necessary to have eveything we need on the GPU side. If I recall correctly from previous discussions, there were still a few parts missing on the GPU.

[mnovak42]

All those material related information are already available on device as well in G4HepEm. So you don't need to worry about the material-cuts, material, element and related informations that are needed at run-time: given the index of the material-cuts or material (i.e. a single integer) everything can be accessed directly even on the device (although these will be be accessed mainly through G4HepEm functionalities even on the device).

Everything is available on the GPU in G4HepEm (except one interval sampling table that I developed for one of the berm models but it will come). When the development keeps going, further components will be transferred as well (e.g. gamma related cross sections).

[amadio]

Ok, sounds good. Then we can probably forget about this particular issue for now and come back to it later or just use it to discuss the connection between geometry and material data, as VecGeom does not read that from the GDML file yet. and we will need a solution for that even for an EM-only prototype.

[amadio]

Just updating the issue for what has been discussed this morning about geometry <-> materials connection on the GPU. The options seem to be one of:

1. Read the geometry twice, once with Geant4 and once with VecGeom, then somehow map volumes from one to the other and copy over the material index used by Mihaly's code onto the VecGeom geometry before loading onto the GPU
2. Use the native Geant4 geometry in memory to create an equivalent VecGeom geometry containing all required information and upload that to the GPU
3. Extend VecGeom's GDML reader to also read in materials and update Mihaly's code to be able to work with that rather than the geometry from Geant4

It seems that for now option 1 is the preferred choice, but it's not yet clear where in the current VecGeom geometry this information should be stored (or if the map of volumes -> materials) will remain as a separate table, possibly like how navigation indices are handled now. As long as a material id can be quickly fetched given a physical and/or logical volume id, it should be ok.

[drbenmorgan]

Asking a question on Mattermost in relation to #87 turned into a discussion, so summarising here for documentation and next steps:

• In G4HepEmData materials are only indexed by int, with a map between Geant4 and G4HepEm indices for each material
• In GDML, materials are only indexed, or rather identified by a string name/ref (NB: this can be either the plain name, e.g. G4_Air, or qualified by pointer, e.g. G4_Air0x...

• In GDML, each logical volume has a materialref, which is the string name name above, e.g.

  <material name="G4_Air">  
  ...
  </material>
  
  <volume name="MyVol">
    <materialref ref="G4_Air"/>
    ...
  </volume>

In exporting/importing geometry/material data in example7 we lose the information needed to map a Geant4 material name to its index, and hence we cannot easily map a VecGeom logical volume index to its material index. To ensure this link can be made (which is essentially @amadio's third option above) we would need:

1. Add a new array (or similar structure) somewhere in G4HepEmData whose entries would be the Geant4 material names. Basically, it's there to provide the mapping G4materialindex <-> G4materialname that we need.
   • It should be pretty lightweight, but need to think about qualifying names with pointers like GDML!
2. Upgrade VecGeom's vgdml parser and maybe core to provide a map between logical volumes (names or indices) and the material name from the GDML.
   • vgdml does handle parsing of materials so the data is there, just needs exposing appropriately.

With those in place, we'd have all the info to create an AdePT internal map between G4HepEmData material indices and the VecGeom logical volume indices using them.

First things first though, we should discuss these new features in G4HepEM and VecGeom.

[drbenmorgan]

Also noted in the SFT meeting by @mnovak42 - what matters as far as physics goes is mapping between logical volume and the materials cut couple. If we want to explore the volume-MCC connection, example7 only has a single region, but the DataImportExport test in G4HepEm:

https://github.com/mnovak42/g4hepem/blob/master/testing/DataImportExport/src/SimpleFakeG4Setup.cc

does exercise the "same material in two regions with different MCC" use case if we need an example setup.

So probably the first thing to fully clarify is what interfaces of G4HepEm we're using (or want to use). @hahnjo are TestEm3's electrons.cu and gammas.cu our primary implementation at the moment?

I'll also start looking in detail at what's in the VecGeom GDML parser both for materials and to use/support extensions.

[mnovak42]

For using G4HepEm outside of its Geant4 environment (that was used for its initialisation), not the material index but the material-cut couple index is required. This is required in each simulation step. These material-cut couple indices (actually the objects) are always known in a G4Track during the tracking and depend not only the material but also the detector region.

Note, that while there is a single e.g. Silicon material object, there are as many Silicon material-cut couple objects as appearances of this Silicon material in different detector regions. This is because secondary production thresholds can be set per detector regions, leading to different restricted integrated (macroscopic cross section, stopping power, range, etc.) data for e-/e+ in the different detector regions.

G4HepEm identifies the physics data, need to be used at a given simulation step, based on this material-cut couple index and this is what needs to be known also in an AdePT at each simulation step.

[hahnjo]

@hahnjo are TestEm3's electrons.cu and gammas.cu our primary implementation at the moment?

With respect to physics, TestEm3 is the same as example9. What differs is how they load / interact with geometry: example9 loads a GDML and only uses a single material and material-cut couple. In TestEm3, I construct and place the shapes by hand and have three materials with one material-cut couple each. The code to map it is a very simple array between the ID of the physical volume (that I know, by construction, is only used once) and the MC index: https://github.com/apt-sim/AdePT/blob/95e6309088c49d159e95b123a8aa404d25c710af/examples/TestEm3/TestEm3.cpp#L195-L204

Then during simulation, I can get this back with a simple lookup: https://github.com/apt-sim/AdePT/blob/95e6309088c49d159e95b123a8aa404d25c710af/examples/TestEm3/electrons.cu#L43-L45 The code in gammas.cu works exactly the same. No guarantees that I'm not messing up between Geant4 and G4HepEm indices at some point!

[agheata]

@drbenmorgan there is a merge request in VecGeom dealing with this: https://gitlab.cern.ch/VecGeom/VecGeom/-/merge_requests/805

[drbenmorgan]

Thanks to @agheata's MR, we should be able to read GDML <auxiliary> and <userinfo> tags. Just to confirm the data we'll write/read:

1. On the Geant4 side: map of Geant4 logical volume name to cuts couple index, i.e.

   G4LogicalVolume* lvStore = G4LogicalVolumeStore::GetInstance();
   for (const G4LogicalVolume* lv : *lvStore) {
     std::string lvName = lv->GetName();
     int lvCutCoupleIndex = lv->GetMaterialCutsCouple()->GetIndex();
   }

   G4HepEm provides the mapping between the G4 MCC index and its own through G4HepEmData::fTheMatCutData which provides the G4HepEmMatCutData::fG4MCIndexToHepEmMCIndex array.

2. Write this to the GDML file in the <userinfo> block, e.g.

   <userinfo>
     <auxiliary auxtype="VolumeCutsCoupleMap" auxvalue="UNUSED">
       <auxiliary auxtype="ALVName" auxvalue="0"/>
       <auxiliary auxtype="BLVName" auxvalue="0"/>
       <auxiliary auxtype="CLVName" auxvalue="1"/>
       ...
     </auxiliary>
   </userinfo>

3. On the AdePT/VecGeom side, we read in that GDML and the G4HepEm data files. Andrei's MR will provide the interface read the <userinfo> and hence together with the VecGeom volume stores and read G4HepEmMatCutData object, we can connect up volumes->MCC like Jonas has done by construction.

[agheata]

👍

[agheata]

Now we are also able to map the material-cut couple indices build by G4HepEm based on transient Geant4 logical volumes, with the corresponding VecGeom logical volumes (keeping corresponding Geant4 and VecGeom geometry in parallel in CPU memory). These indices can be used during stepping to talk to G4HepEm physics for sampling the step and final state. The procedure is demonstrated starting with example12:

• reading G4 geometry from a GDML file, fill G4 material-cut couple table based on the defined materials and regions, and using them for G4HepEm physics initialization
• reading VecGeom from the same GDML file using vgdml parser
• matching every vecgeom::LogicalVolume with the corresponding G4LogicalVolume and creating a LUT LogicalVolume::id->MCC_index (MCC_index is the material-cut couple index of G4HepEm for the given volume)
• copying the MCC LUT to device and using it in the physics kernels to access the correct MCC index during stepping

This approach is generic and complementary to the one described above by Ben, having the Geant4+G4HepEm initialization step in the same application with the AdePT simulation, rather than decoupling them using a persistent layer. Probably the combination of these 2 approaches is enough for the material functionality needed by the first prototype. Closing for now.