Getting informations from a GDML file using a root script

Hi, I have a .gdml file and I want to be able to read informations of every volume in the geometry using a root script, without parsing the gdml code.

The informations I would like to obtain are the dimensions and position of each volume, and ideally it should not depend on the knowledge of the content of the gdml.

An example of what i would like to do is:

TGeoManager *geom = TGeoManager::Import("geometry.gdml");  //correct
int n_volumes = geom->GetListOfVolumes()->GetEntries();            //correct

// the following is invented
for(int i=0; i<n_volumes; i++){
    // variables declared before
    x = geom->GetVolume(i)->GetPositionX();
    y = geom->GetVolume(i)->GetPositionY();
    name = geom->GetVolume(i)->GetVolumeName();
    shape = geom->GetVolume(i)->GetShape();
    // and maybe save all these in a TTree or cout them or whatever

I tried to look for examples or documentation but i could not find much.

Thank you for the help,


This is too specific to be written in the docs or examples. Your snippet is fine, you can do the loop as:

TGeoVolume *vol;
TIter next(geom->GetListOfVolumes());
while ((vol = (TGeoVolume*)next())) {
  TString name = vol->GetName();
  auto shape = vol->GetShape();
  // There is no such thing as position of a volume, they are all logical entities positioned relative to their mother - see the geometry documentation.
 // The TGeoShape interface is also described in the documentation

Hi, thank you for your answer.
After a bit of more research and trial and error I managed to obtain something similar to what you have suggested but using nodes instead of volumes

TGeoManager *geom = TGeoManager::Import(filename);
  TGeoVolume *top = gGeoManager->GetTopVolume();

  TGeoIterator next(top); 
  TGeoNode *node;
  while ((node=next())){

Passing via the volumes as you said is of course more straightforward.
I understand that in GDML like in Geant4 positions are relative to the Mother volume, but also that information can be really useful as I should be able to reconstruct where the element is in the Top volume.

Also, I tried to trigger the TGeoShape::Contains() but I failed. Indeed I could simply go element by element using the while loop showed above and check if a given particle is inside or outside that volume.

I tried:

  double point[3];
  point[0]=0.05;   //x  in meters wrt center of the volume
  point[1]=0.05;   //y  in meters wrt center of the volume
  point[2]=0.00;   //z  in meters wrt center of the volume
  while ((vol = (TGeoVolume*)next())){

It fails, but I suspect it can be a beginner error… The Contains() function could be the most promising for what I have to do actually, as it would allow me to simply check whether a particle is inside or outside a given volume, regardless of its position in the main volume.

Hi, it doesn’t work like this, if you need to enumerate all physical instances of volumes in world coordinates, you should use a TGeoIterator which is iterating the physical tree from top to leaves, giving you the global transformation matrix. A physical node (real object in the geometry) is represented by a path made of logical nodes: /TOP_1/A_1/B_2 where the global transformation is the multiplication of the matrices for TOP_1, A_1 and B_2. You can do:

TGeoIterator next(geom->GetTopVolume());
TGeoNode*current; TGeoVolume *vol;
TString path;
while ((current=next())) {
  vol = current->GetVolume();
  TGeoMatrix *global = next.GetCurrentMatrix();
  // if you want to see where is the center (origin) of this object in TOP coordinates:
  const double *local_bbox_orig = ((TGeoBBox*)vol->GetShape())->GetOrigin(); 
  double master_orig[3];
  global->LocalToMaster(local_bbox_orig, master_orig);
  std::cout << path << ": volume=" << vol->GetName() << "  position: (" << master_orig[0] <<   ", " << master_orig[1] << ", " << master_orig[2] << ")\n";

Sorry for the typos if any, I haven’t tested this snippet. Note that if you have a million physical volumes, you should buy a big screen before running this :slight_smile:

To check where a point is in your setup, rather do:

double global[3] = {...}; // these should be coordinates in the world reference.
TGeoNode * node = gGeoManager->FindNode(global);
if (!node) std::cout << "point is outside the setup\n";
  std::cout << "point is inside " << node->GetName() << "  in path: " << gGeoManager->GetPath() << "\n";

Thank you very much, your last two answers explain exactly what I was looking for!

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.