LocateGlobalPointAndSetup in volume with G4PVReplica

Geometry, Fields and Transportation
1

Dear all,

I would like to manually step through my geometry to get a hand on the material properties for studies and tracking etc.

My idea is to use the G4Navigator with the LocateGlobalPointAndSetup and ComputeStep method to step though the single materials. Starting from a given point the ComputeStep method is called to obtain the stepLength towards the next boundary in the respective direction, then with a tiny push one can cross the boundary and is in the next volume,

This method works nicely with all my geometry (CAD - tessellated, reflections, etc.) but unfortunately not with a volume containing a G4PVReplica → segmentation violation in LocateGlobalPointAndSetup.

Also in works with no problems using the normal G4 simulation. Therefore I assume that the geometry is fine and I am just missing something in my LocateGlobalPointAndSetup method

Example for Geometry creation with G4PVReplica:

// 1. Mother volume - sensor 
G4Box* sensor = new G4Box( "Sensor", (pixelSizeX * numberPixelX) / 2.0, (pixelSizeY * numberPixelY) / 2.0,pixelTickness / 2.0 );
G4LogicalVolume* sensorLV = new G4LogicalVolume(sensor, helium, sensor->GetName());

// 2. Replica volume - row
G4Box* row = new G4Box( "row", (pixelSizeX * numberPixelX) / 2.0, pixelSizeY / 2.0, pixelThickness / 2.0);
G4LogicalVolume* rowLV = new G4LogicalVolume(rowSolid,helium,rowSolid->GetName());

// 3. Create replica row volume in mother sensor
new G4PVReplica(rowLV->GetName(), rowLV, sensorLV, kYAxis, numberPixelY, pixelSizeY);

Example for manual stepping using LocateGlobalPointAndSetup and ComputeStep:

G4Navigator* g4Navigator = new G4Navigator();
g4Navigator->SetWorldVolume( worldConstruction->getWorldVolume() );

const G4ThreeVector direction(0,0,1); // Simply go along z-axis
const double proposedPhysicsStep = 10 * CLHEP::m; 
while( currentZ < endZ ) {
    G4ThreeVector point( 0 * CLHEP::mm, 0 * CLHEP::mm, currentZ);
    G4VPhysicalVolume *volume = g4Navigator->LocateGlobalPointAndSetup( point, nullptr, true, true );

    const double radLen = volume->GetLogicalVolume()->GetMaterial()->GetRadlen();
    const double density = volume->GetLogicalVolume()->GetMaterial()->GetDensity();

    double newSafety;
    double stepLength = g4Navigator->ComputeStep( point, direction.unit(), proposedPhysicsStep, newSafety );
    if( stepLength > proposedPhysicsStep  ) {
        stepLength = proposedPhysicsStep;
    }
    stepLength += 0.1 * CLHEP::um; // Cross boundary - flexible
    currentZ += stepLength; // Proceed along z
}

Maybe there is also a smarter way to do this - but fromt the principle point of view it works - as long as one does not use G4PVReplica volumes.

G4 Version: 10.5.0

Thank you very much for any feedback on this.

2

Update:

After some investigation I discovered that especially in the G4PVReplica the voxel headers where missing - those are created during the optimization of the geometry. This does not happen if one uses the code snippets mentioned above. Normally it is done via the runManager->BeamOn() function.

Nevertheless, if one wants to only use the geometry itself my solution in this case is to manually close the geometry before studying it - and opening it at the end to have it properly deleted.

Basically one has to wrap the above code in:

// First create your geometry
// Then close it
G4GeometryManager* geomManager = G4GeometryManager::GetInstance();
if (!geomManager->IsGeometryClosed()) {
    geomManager->CloseGeometry(true);
}

// Do fancy stuff here
// ...........

// At the very end: open geometry again for cleanup
geomManager->OpenGeometry();

I am not 100% sure whether this is the way it is suppose to be - but it seems to work. I hope I do no miss anything, but please let me know if I did!

Cheers!

3

Yes, the optimisation structure gets generated at the time the geometry is closed and is necessary to perform efficient navigation. The optimisation structure is not applied to replica volumes, as not needed.