Get volume positioning

Hi,
I’m deciding at the PreUserTrackingAction, if a particle should be tracked or not by figuring out if it will hit a detector, given random initial position and direction, and loosely accounting for multiple scattering effects.
However, to do this I need the position of the detector in the world reference frame.
My understanding is that using the touchable history is not the solution because the particle’s tracking has not started yet… am I wrong ? How can I do it ?

Thanks,
Bernardo Tomé

You can do it by instantiating your own G4Navigator (maybe have a data member of your TrackingAction) and using that to get the global position of your detector volume.

That way, you’re not messing with (or messing up!) the real tracking navigator.

Another option (involving redesigning your simulation :frowning: ) might be to consider the “adjunct” or “reverse” MC method. If you’re only dealing with electromagnetic interactions (no high-energy nuclear or particle physics), Geant4 has a framework where you can set up your geometry, and launch “reverse particles” from your detector out to the world, and then it turns around and tracks just those particles to get the interactions and energy deposits.

https://geant4-userdoc.web.cern.ch/Doxygen/examples_doc/html/ExampleReverseMC01.html

Thanks for the prompt reply !
I’m not sure I understood the usage of another G4Navigator… wouldn’t that imply in any case tracking the particle through the geometry ?

Regarding using the reverse MC it’s a very nice feature I was not aware of. The catch here is that I’m simulating atmospheric muons, accounting for the energy vs direction correlations, so I might need to properly weight the events…

If just getting the absolute position of the detector in the world frame is not too complicated I think I would go to that solution.

Thanks again !

No, but I’m realizing that it wouldn’t do what you originally asked. G4Navigator takes a world-coordinate as input and returns a touchable for the daughter-most PV at that coordinate. You want the opposite: given a PV (or the unique name/copy-number of the PV?), get the world coordinate of its center (or equivalently, get the touchable and transform local (0,0,0) to world).

In our simulation framework, we borrowed code from the GAMOS experiment (@arce) that allows us to both interrogate the geometry stores by name, and also to build a touchable “from scratch”. This is non-trivial in general: with a highly nested geometry, you may have the same daughter PV and copy number within an LV that itself gets placed multiple times. You need to know the whole chain of PVs (that’s what the touchable encodes).

As I noted, in general it’s quite complex. Depending on your specific geometry, it may be simple enough to do by hand.

Right ! I’ll try to survive by keeping track of the chain of PVs by hand. It’s not that nested the geometry.
Thanks for the inputs.

Mike’s right. In general it’s complicated. Even if you know the name of the volume or have its pointer, its location depends on its mother and grandmother, etc. It may turn up in several different locations.

I’m guessing you have quite a simple setup in which there is no such ambiguity. Still, there’s no simple way to get its position. But you know where it is. So you could record it in your detector construction and pick it up from there.

In the vis system, we do need to know a volume’s position, and we have code that can be used. We have repurposed it to provide some handy commands. But these simply print out the properties - position, rotation, material,… If you want to have this information in your own stepping action or whatever you would have to look at how we do it.

To specify a volume you need to know its path - world, great-grandmother,…,mother, volume. Use

/vis/touchable/findPath <physical-volume-name>

This prints a list of paths (possible more than one) containing this volume. You can use the printout to supply arguments to

/vis/set/ouchable <path>

(see command guidance - use help). Then

/vis/touchable/dump

For example (example B1):

  /vis/touchable/findPath Shape2
 World 0 Envelope 0 Shape2 0 (mother logical volume: Envelope)
Use this to set a particular touchable with "/vis/set/touchable <path>"
or to see overlaps: "/vis/drawLogicalVolume <mother-logical-volume-name>"
  /vis/set/touchable World 0 Envelope 0 Shape2 0
  /vis/touchable/dump

Good luck
John

Forgot to say…a “touchable” is the name we give to a particular instance of a physical volume, depending its “path”. As I said, a physical volume may - is allowed to - turn up in multiple locations. Naturally, the vis system deals in touchables.

Many thanks for the guidance ! In fact I am setting the geometry through a gdml file, so don’t have explicitly the volume positioning in the code, but I guess I can fetch that somehow. And, as you anticipated, the setup is not very involved. In any case I will look at the code behind the commands you mentioned, which seem very useful even for future needs !

Thanks.

Yes, I was aware of it. Thanks.

Hi,
I managed to adapt the code in G4VisCommandsTouchable.cc to fit my needs.
It is working nicely. I can now fetch the position(s) of a given physical volume in the world ref frame.

Thanks again !
Bernardo

Well done!!!

Is there anything that might be useful to others? Could we, perhaps, provide a utility?

John

Indeed I was wondering if there shouldn’t be some utility in this context… The code is basically a cut&paste from the command class I mentioned. I hope I didn’t mess anything, but it is working for my needs. Here it is :

// The coode below was adapted from the code in G4VisCommandsTouchable.cc
// It is used to find a given physical volume in the world ref. frame.

std::vector
DetectorConstruction::GetPositionInWorld(const G4String& physName)
const
{
G4TransportationManager* transportationManager =
G4TransportationManager::GetTransportationManager ();

size_t nWorlds = transportationManager->GetNoWorlds();

std::vectorG4PhysicalVolumesSearchScene::Findings findingsVector;
std::vector<G4VPhysicalVolume*>::iterator iterWorld =
transportationManager->GetWorldsIterator();
for (size_t i = 0; i < nWorlds; ++i, ++iterWorld) {
G4PhysicalVolumeModel searchModel (*iterWorld); // Unlimited depth.
G4ModelingParameters mp; // Default - no culling.
searchModel.SetModelingParameters (&mp);
G4PhysicalVolumesSearchScene searchScene (&searchModel, physName);
searchModel.DescribeYourselfTo (searchScene); // Initiate search.

for (const auto& findings: searchScene.GetFindings()) {
  findingsVector.push_back(findings);
}

}

for (const auto& findings: findingsVector) {
G4cout
<< findings.fFoundBasePVPath
<< ’ ’ << findings.fpFoundPV->GetName()
<< ’ ’ << findings.fFoundPVCopyNo
<< " (mother logical volume: "
<< findings.fpFoundPV->GetMotherLogical()->GetName()
<< ‘)’
<< G4endl;

}

std::vector vPos;

if (!findingsVector.size()) {
G4cout << physName << " not found" << G4endl;
} else {
for (const auto& findings: findingsVector) {
vPos.push_back(findings.fFoundObjectTransformation.getTranslation());
}

for (const auto& pos: vPos) {
  G4cout << pos/m << G4endl;
}

}

return vPos;
}

Thanks, Bernado. I think that’s worth packaging up into a utility - I’ll do it for the next release.

Best wishes
John

Perfect !
All the best,
Bernardo

Hi, I have come across a similar issue where I am looking for a convenient way to find the global coordinates of a lot of highly nested volumes, and was glad to find that someone figured it out. I tried using your code, but am confused about how exactly to do so.
Where exactly would this snippet go? I tried putting it in the detector construction class of my project, but it keeps throwing this error on execution -

-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : modeling0012
      issued by : G4PhysicalVolumeModel::DescribeYourselfTo
No model.
*** Fatal Exception *** core dump ***
 **** Track information is not available at this moment
 **** Step information is not available at this moment

-------- EEEE -------- G4Exception-END --------- EEEE -------

It would be helpful if you could provide a complete example of how you utilized this function. I called it from inside a function that I am using to place all the volumes, after the volumes are placed. I think this is probably the mistake. Any help would be appreciated.

Thanks,
Advait