Compute the global transform for all instances of a volume

I am trying to solve this problem: given a G4VPhysicalVolume*, find all instances in the global geometry and compute the net transform. The following code works if all volumes are singly placed (i.e., with G4PVPlacement):

std::vector<G4AffineTransform> LocateAllVolumes(const G4VPhysicalVolume* pv,
                                                G4NavigationHistory* geomtree)
{
  //this function is recursive. If geomtree is non-null, this is a recursive call
  // if it is null, we need to build it first
  bool cleanup = false;
  if(!geomtree){
    geomtree = new G4NavigationHistory;
    geomtree->SetFirstEntry(world);
    cleanup = true;
  }

  if(geomtree->GetTopVolume() == pv){
    return std::vector<G4AffineTransform>{geomtree->GetTopTransform().Inverse()};
  }

  //otherwise we need to loop through children
  std::vector<G4AffineTransform> result;

  //now loop through children and check if we belong
  G4LogicalVolume* currentLog = geomtree->GetTopVolume()->GetLogicalVolume();
  for(size_t i=0; i<currentLog->GetNoDaughters(); ++i){
    G4VPhysicalVolume* daughter = currentLog->GetDaughter(i);
    if(daughter == pv || daughter->GetLogicalVolume()->IsAncestor(pv)){
        geomtree->NewLevel(daughter);
        std::vector<G4AffineTransform> subs = LocateAllVolumes(pv, geomtree);
        result.insert(result.end(), subs.begin(), subs.end());
        geomtree->BackLevel();
    }
  }

  if(cleanup){
    delete geomtree;
  }

  return result;
}

My problem is, how to extend this or replace with the proper way to do things when there are replicas and or parameterized volumes? It seems to correctly construct a touchable for a replica, I have to set the transform of the PVReplica before calling G4NavigationHistory::NewLevel. And it looks like the correct code to calculate that exists only in G4ReplicaNavigation::ComputeTransform, i.e., nowhere I can call it from. But all the G4Navigator interactions are for computing transforms of points, and I don’t have any coordinates to work with, just a volume pointer.

What am I missing?

It is not an easy problem. I solved it long time ago and you can find the code in GAMOS: source/GamosCore/GamosGeometry/src/GmTouchable. You can download GAMOS and use it freely, even just copy the file (of course, as usual mentioning it would be nice). Look for GAMOS in the Geant4 Applications web page

If I’ve understood correctly, this is similar to a problem faced by the vis system. When a user types /vis/drawVolume Gap (I’m thinking of example B4a), we need to find all occurrences of Gap. You will find the relevant code under /vis/scene/add/volume. Here is a digest that would fit at the end of the detector construction:

#include "G4PhysicalVolumesSearchScene.hh"
...
  G4String name = "Gap";
  G4int copyNo = -1;  // All copy numbers
  G4PhysicalVolumeModel searchModel (worldPV);  // Unlimited depth.
  G4ModelingParameters mp;  // Default - no culling.
  searchModel.SetModelingParameters (&mp);
  G4PhysicalVolumesSearchScene searchScene (&searchModel, name, copyNo);
  searchModel.DescribeYourselfTo (searchScene);  // Initiate search.
  const auto& findings = searchScene.GetFindings();

and you can get at the information like this

  for (const auto& finding: findings) {
    G4Scale3D scale;
    G4Translate3D translation;
    G4Rotate3D rotation;
    finding.fFoundObjectTransformation.getDecomposition(scale,rotation,translation);
    G4cout
    << finding.fpFoundPV->GetName()
    << ':' << finding.fFoundPVCopyNo
    << " found at depth " << finding.fFoundDepth
    << " with translation ("
    << translation.dx() << ',' << translation.dy() << ','<< translation.dz() << ')'
    << G4endl;
  }

Hope. this helps.

The code which is used in ParameterisedNavigation to create a new level for a touchable is at line 647:

  // Obtain solid (as it can vary) and obtain its parameters
  //
  pSolid = IdentifyAndPlaceSolid( replicaNo, pPhysical, pParam ); 

  // Setup history
  //
  history.NewLevel(pPhysical, kParameterised, replicaNo);

You can copy or better use the 3 line method IdentifyAndPlaceSolid by creating your object of this class.

At a more general level do you want to create something that works in one application, or that can be used with a wide variety of applications ? In the latter case it could also need to work with Nested Parameterisation, for which the full path from the world to the current (parent) volume must be available in a touchable (passed to the parameterisation).

I don’t think this solves my immediate problem, because replicated volumes don’t have an associated paramerisation. However I do also want to solve the problem generally, including nested parameterisations as you mention.

The end goal is to be able to distribute primaries uniformly in all copies of a volume without knowing any details of the geometry

I suggest to copy an extracted part of

G4ReplicaNavigation::LevelLocate( G4NavigationHistory& history,
const G4VPhysicalVolume* blockedVol,
const G4int blockedNum,
const G4ThreeVector&, // globalPoint
const G4ThreeVector*, // globalDirection
const G4bool, // pLocatedOnEdge
G4ThreeVector& localPoint )

You may be able just to call this method. Else create an object and then use something like:

G4VPhysicalVolume *motherPhysical, *pPhysical;
G4LogicalVolume *motherLogical;
G4SmartVoxelHeader *motherVoxelHeader;
G4int nodeNo;

motherPhysical = history.GetTopVolume();
motherLogical = motherPhysical->GetLogicalVolume();
motherVoxelHeader = motherLogical->GetVoxelHeader();
pPhysical = motherLogical->GetDaughter(0);

nodeNo = VoxelLocate(motherVoxelHeader, localPoint);

replicaNav->ComputeTransformation(nodeNo, pPhysical, localPoint);
history.NewLevel(pPhysical, kReplica, nodeNo);
pPhysical->SetCopyNo(nodeNo);