Extract voxelized properties from geometry

Geometry, Fields and Transportation
1

Given an arbitrary geometry defined in Geant4, is it possible to obtain a voxelized map of some arbitrary property of the geometry, for example material density (at centre of voxel, or averaged over the voxel), for (regular) user-defined voxel sizes?

Alternatively, is it possible to query a geometry for some property at an arbitrary point within the geometry?

2

That used to be possible, and maybe still is. Once upon a time, Geant4 had a built-in tool where it would shoot geantinos and generate a plot for you of the interaction length or absorption length or something at points along the trajectories, basically giving you a “particle’s eye view” along the beam direction.

This is definitely possible, and fairly easy to do, though not clearly written down. The geometry has to be “closed” (i.e., after /run/initialize), so that the voxelization has been completed. Then you can instantiate your own G4Navigator instance (do not use the tracking navigator!), initialize it with the world volume, and query it with a position. You’ll get back a G4Touchable for the PV at that location, which gives you full access to everything you can know about the geometry.

3

BTW, here’s a simplified version of the code in our experiment’s framework that I use for this purpose:

// Create non-tracking Navigator for use with position finding below

G4Navigator* G4CMP::GetNavigator() {
  static G4ThreadLocal G4Navigator* theNavigator = 0;
  if (!theNavigator) theNavigator = new G4Navigator;

  // Make sure current world volume is the one in use  
  G4VPhysicalVolume* theWorld =
    G4TransportationManager::GetTransportationManager()->
      GetNavigatorForTracking()->GetWorldVolume();

  if (theNavigator->GetWorldVolume() != theWorld)
    theNavigator->SetWorldVolume(theWorld);

  return theNavigator;
}

G4VTouchable* G4CMP::CreateTouchableAtPoint(const G4ThreeVector& pos) {
  G4VTouchable* touchable = new G4TouchableHistory;
  GetNavigator()->LocateGlobalPointAndUpdateTouchable(pos, touchable, false);
  return touchable;
}
4

Hehe … I was wondering whether I might end up having to write something like this.

Great, navigators weren’t really on my radar, so that’s a massive clue.

Brilliant, that should save be a lot of time bumbling around. Many thanks!

5

Yeah. It’s something that’s typically treated as “internal” to G4 – it’s one of the main drivers for the whole tracking system, as you can imagine – so it really isn’t(*) well documented in the Application Developer’s Guide. But we use it extensively for exactly your purpose – mapping positions back to the geometry, and using information from that.

EDITED TO ADD: (*) “isn’t” the last time I was trying to learn about it a few years ago. That may have changed more recently, but since I wrote the code I haven’t had to look again.

6

BTW, you can find the relevant information here in the Application Developers Guide:
https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Detector/Geometry/geomNav.html

and specifically, about using the navigator (alternative or not) to locate points etc…:
https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Detector/Geometry/geomNav.html#using-the-navigator-to-locate-points

If something is not well documented, please just let us now!

1 Like
7

Thank you, it works.

However, I’m finding the queries surprisingly slow. For example, looking up just the density, I’m getting around 550 results per second. (I need to build up fine-grained 3D maps of the geometry, so this seems prohibitively expensive.)

Does this sort of rate sound correct, or am I doing something daft?

8

I suspect that you’re telling your G4Navigator to start over from scratch and find each point individually. The code is heavily optimized to be really fast for tracking, where the points come sequentially and are generally close to one another. So the G4Navigator will cache it’s previously found point, and can start the next search from there. Unless you tell it to start from the beginning each time.

9

Yes, I guess that’s what I’m doing.

What would starting over from scratch look like? Repeated calls to LocateGlobalPointAndUpdateTouchable()?

What would NOT restarting from scratch look like? Calling something like ComputeStep() repeatedly, after an initial LocateGlobalPointAndUpdateTouchable()?

10

Using LocateGlobalPointAndUpdateTouchable() in relative mode (last argument = true) gives me 4 orders of magnitude speedup … but the results stop making sense. Maybe I’m failing to get the first point correctly in absolute mode.

11

The expected protocol for calling these methods is:

  • For the first step call either LocateGlobalPointAndSetup() or LocateGlobalPointAndUpdateTouchable() which calls the first;
  • Call ComputeStep() for the same point
  • Either a) call the method that informs that you are taking the full step and then one of the above Locate methods if you ‘took’ the full computed step and thus arrived at a boundary; or b) call LocateGlobalPointWithinVolume if the step did not.
  • Call ComputeStep() again - and so on.
12

In the case when you are only locating points, you can repeatedly call one of the first two Locate methods, with a new point each time.

If you are getting only 550 results per second, it seems that either: your geometry is incredibly complicated, or, more likely, you did not close the geometry enabling the optimisation voxels to be built - which speed up calls to Locate and ComputeStep.

The documentation is clear that the geometry must be closed! Else each call to Locate must ask all daughter volumes at each level until it find one that contains the relevant point…

13

Here is the code I have been using, which led to my previous observation. It contains some ‘superfluous’ features in order to enable exploring alternatives, which I’ll explain below.

  run_manager -> Initialize();
  auto navigator = std::make_unique<G4Navigator>();
  auto world = G4TransportationManager::GetTransportationManager()
    -> GetNavigatorForTracking()
    -> GetWorldVolume();
  navigator -> SetWorldVolume(world);
  auto touchable = std::make_unique<G4TouchableHistory>();

  G4ThreeVector old_point{0,0,0}; \\ [1]
  auto relative_search = false;   \\ [2]
  for (auto iz=0; iz<nz; ++iz) {
    for (auto iy=0; iy<ny; ++iy) {
      for (auto ix=0; ix<nx; ++ix) {
        auto x = (dx-DX) / 2 + ix * dx;
        auto y = (dy-DY) / 2 + iy * dy;
        auto z = (dz-DZ) / 2 + iz * dz;
        G4ThreeVector new_point{x,y,z};
        auto old_to_new_vecj = new_point - old_point; \\ [3]
        old_point = new_point;
        navigator -> LocateGlobalPointAndUpdateTouchable(new_point, \\ [4]
                                                         touchable.get(),
                                                         relative_search=false); \\ [5]
        relative_search = true; // Only the first one needs to be absolute  [6]
        auto density = touchable
          -> GetVolume()
          -> GetLogicalVolume()
          -> GetMaterial()
          -> GetDensity() / (kg/m3);
        write << static_cast<float>(density);
      }
    }
  }

  • [5]

    • Keeping relative_search=false gives correct results at a rate of mere hundreds of results per second.

    • Switching to relative_search=true speeds up the search to a rate of over a million results per second, but gives incorrect results.

    • Switching to relative_search (set in [2] and [6]) also gives a much faster rate, but the results are still wrong.

      According to the documentation

      A ‘relative’ search may be used for the first call of the function which will result in the search defaulting to a search from the root node of the hierarchy.

      which I interpret to mean that this explicit treatment of the first search as a special case, is implicitly taken care of by the function in the relative_search=true version. So this last variation should be equivalent to the previous (and yet my results appear to be slightly different).

  • [1,3] The functions have LocateGlobalPoint in their names, so I would expect to have to pass an absolute position. Just in case they require a displacement relative to the previous point, when used in relative search mode, here is support for calculating such displacements. To use it, change [4] to old_to_new_vec. This also gives fast but incorrect results.

  1. The region of the geometry being inspected, is almost trivial. (Elsewhere there are very many small elements, but the hierarchy is never very deep.)
  2. If calling G4RunManager::Initialize() closes the geometry, then the geometry is closed.
  3. Switching from absolute to relative mode in LocateGlobalPointAndUpdateTouchable results in a speedup which gives over a million results per second … but the results cease to make sense.

In summary, at present I appear to have the choice between excruciatingly slow correct results, or blazingly fast generation of garbage. I’m guessing that I don’t understand the relative search mode of the LocateGlobalPoint* functions.

14

Hi Jacek,

Can you share a GDML file of your setup, or at least the volume in question. (Privately by email is fine.)

Discussing your results in a vacuum is challenging.

Cheers,
John

15

Hi John,

There was a bug in my geometry (which I still don’t fully understand). When I isolate the part of the geometry that needs to be scanned, thereby taking the bug out of the picture, everything works as expected (over 2 million correct results per second, both in relative and non-relative mode).

Thanks for your help.