How to get individual segment length crossing voxels

Hello everybody,

I am working on a particle simulation program using ROOT (geometry package, histograms, etc) and I need to know what is the segment length inside some voxel of a straight line during the transportation of a particle. Is there any easy/fast way to do it using ROOT?

Thanks, Luís

Assuming
Double_t x,y,z,xp,yp,zp;
the current point and the direction cosinus in a given volume,
simply compute the distance to the next boundary with
Double_t dist = 9999999;
gGeoManager->FindNextBoundary(dist);

see Users Guide for more details.

Rene

Hello Rene,

Some time ago I’ve started to use the way suggested by you, but I got some troubles those I can’t remenber now. Anyway, I’ll try it again and if I get troubles I’ll bring it here.

Thanks, Luís

Hello Rene,

I have changed my code in order to use the TGeoManager navigation methods and it is working good without voxelization. My geometry is pretty simple: a cylinder phantom inside a world (wich is also a cylinder). When I divide my target in slices I get some strange results from the queries of TGeoManager. For example, when a primary particle is entering in one of cylinder’s face, the gGeoManager->FindNextBoundary (99999.) gives me these results in each step:

sliceZ_1
sliceZ_50
sliceZ_50
sliceZ_50
sliceZ_50
sliceZ_50
phantom_phys_1
phantom_phys_1
phantom_phys_1
phantom_phys_1
phantom_phys_1
phantom_phys_1
phantom_phys_1

The node slicez_1 is the first voxel and the sliceZ_50 the last one. The strange point is that all my voxels are inside the phantom_phys_1 :open_mouth: . The voxelization is done like this:

	sliceZ = fTargetVolume->Divide ("sliceZ", 3,50, 0, 0);
	fTargetVolume->Voxelize("");

How to unmix these things?

Another problem is that part of my code needs to know the SafeDistance and the Step for the phantom ignoring the voxelization, and the other part needs to know the Step inside the voxelized geometry. Both at the same time! Do you have any sugestion?

Thanks, Luís

Hi Luis,

You mix voxelisation with dividing a volume. They are both used for improving navigation performance, but in slightly different ways and cannot be mixed. When dividing a volume you create a set of daughters of the volume by enforcing a pattern. No extra voxelisation of daughters is needed. TGeoVolume::Voxelize is called AUTOMATICALLY when closing the geometry. Your call to Voxelize() after Divide() has no effect - you should not call Voxelize() at all…

The sequence of boundaries that you get when tracking your primary highly depends on the way you cross boundaries. For instance in many cases it is not enough to propagate from the last track position EXACTLY with the distance given by FindNextBoundary. There is another method that deals with that:
TGeoNode *TGeoManager::FindNextBoundaryAndStep(Double_t stepmax, Bool_t comp_safe)
The difference is that this method makes the navigator to effectively cross the next found boundary. If you need the safety to be computed, just activate “comp_safe”. You can retreive the step and safety after calling either FindNextBoundary or FindNextBoundaryAndStep:

Double_t step = gGeoManager->GetStep(); //step to next boundary
Double_t safety = gGeoManager->GetSafeDistance(); // sfety to closest boundary

If you still have problems, attach a small macro and I will try to see what you are doing wrong.

Regards,

Hi Andrei,

Thank you for your help, I have implemented your sugestion about TGeoManager::Voxelize() and TGeoManager::FindNextBoundaryAndStep().

I have done a simple macro where I try to understand some of features of tracking from TGeoManager. I have tried 3 ways to get the Step/SafeDistance using a voxelized phantom and a non-voxelized phantom.

1st case - gGeoManager->FindNextBoundary (99999.):
With a non-voxelized phantom this method gives me the distance to the next boundary like I need. In a voxelized phantom I get as the first step 4.8 (entrance face of the last voxel) instead of 5.0.

2nd case - gGeoManager->FindNextBoundaryAndStep(.5, kTRUE) :
In non-voxelized phantom I get as proposed step 0.5, but actualy the step taken is 1.0 :exclamation: What is wrong? In a voxelized phantom I get the same results.

3rd case - gGeoManager->FindNextDaughterBoundary (currpos, currdir, nodeid); Well, the step chosen depends on the last test! This is my fault, but I don’t know where is the bug.

I appreciate any comments and suggestions.

Thanks, Luís
geomtracktest.C (2.34 KB)

Hi Luis,

I attached your macro with some fixes/comments. The starting point is on the boundary between 2 division “voxels”, so the first step may depend on what the navigator decides it is the starting node. You should not use FindNextDaughterBoundary() since this is looking just inside the current volume (it is an utility that should become protected). When using FindNextBoundaryAndStep() (the method that you should use) you do not have to propagate the point using Step() since this is done automatically. The Step() method tries just blindly to cross the boundary by adding an extra step of 1E-6.

See the comments in the code.
The output of your macro looks now like below, which is correct.

Let me know if you have more problems.
Cheers,
[size=75]===========> Event: 0
Start: sliceZ_25
On boundary fStep = 1e-10 NodeID = sliceZ_26 Next point = (0, 0, 1e-10)
On boundary fStep = 0.2 NodeID = sliceZ_27 Next point = (0, 0, 0.2)
On boundary fStep = 0.2 NodeID = sliceZ_28 Next point = (0, 0, 0.4)
On boundary fStep = 0.2 NodeID = sliceZ_29 Next point = (0, 0, 0.6)
On boundary fStep = 0.2 NodeID = sliceZ_30 Next point = (0, 0, 0.8)
On boundary fStep = 0.2 NodeID = sliceZ_31 Next point = (0, 0, 1)
On boundary fStep = 0.2 NodeID = sliceZ_32 Next point = (0, 0, 1.2)
On boundary fStep = 0.2 NodeID = sliceZ_33 Next point = (0, 0, 1.4)
On boundary fStep = 0.2 NodeID = sliceZ_34 Next point = (0, 0, 1.6)
On boundary fStep = 0.2 NodeID = sliceZ_35 Next point = (0, 0, 1.8)
On boundary fStep = 0.2 NodeID = sliceZ_36 Next point = (0, 0, 2)
On boundary fStep = 0.2 NodeID = sliceZ_37 Next point = (0, 0, 2.2)
On boundary fStep = 0.2 NodeID = sliceZ_38 Next point = (0, 0, 2.4)
On boundary fStep = 0.2 NodeID = sliceZ_39 Next point = (0, 0, 2.6)
On boundary fStep = 0.2 NodeID = sliceZ_40 Next point = (0, 0, 2.8)
On boundary fStep = 0.2 NodeID = sliceZ_41 Next point = (0, 0, 3)
On boundary fStep = 0.2 NodeID = sliceZ_42 Next point = (0, 0, 3.2)
On boundary fStep = 0.2 NodeID = sliceZ_43 Next point = (0, 0, 3.4)
On boundary fStep = 0.2 NodeID = sliceZ_44 Next point = (0, 0, 3.6)
On boundary fStep = 0.2 NodeID = sliceZ_45 Next point = (0, 0, 3.8)
On boundary fStep = 0.2 NodeID = sliceZ_46 Next point = (0, 0, 4)
On boundary fStep = 0.2 NodeID = sliceZ_47 Next point = (0, 0, 4.2)
On boundary fStep = 0.2 NodeID = sliceZ_48 Next point = (0, 0, 4.4)
On boundary fStep = 0.2 NodeID = sliceZ_49 Next point = (0, 0, 4.6)
On boundary fStep = 0.2 NodeID = sliceZ_50 Next point = (0, 0, 4.8)
On boundary fStep = 0.2 NodeID = world_1 Next point = (0, 0, 5)
SafeDistance = 0 fStep = 0.5 NodeID = world_1 Next point = (0, 0, 5.5)
SafeDistance = 0.5 fStep = 0.5 NodeID = world_1 Next point = (0, 0, 6)
SafeDistance = 1 fStep = 0.5 NodeID = world_1 Next point = (0, 0, 6.5)
SafeDistance = 1.5 fStep = 0.5 NodeID = world_1 Next point = (0, 0, 7)
SafeDistance = 2 fStep = 0.5 NodeID = world_1 Next point = (0, 0, 7.5)
SafeDistance = 2.5 fStep = 0.5 NodeID = world_1 Next point = (0, 0, 8)
SafeDistance = 2 fStep = 0.5 NodeID = world_1 Next point = (0, 0, 8.5)
SafeDistance = 1.5 fStep = 0.5 NodeID = world_1 Next point = (0, 0, 9)
SafeDistance = 1 fStep = 0.5 NodeID = world_1 Next point = (0, 0, 9.5)
On boundary fStep = 0.5 Next point = (0, 0, 10)
[/size]
geomtracktest.C (3.3 KB)

Hello Andrei,

Thank you very much for your helpfull hints! I think now I know how to implement my tranposrt code. If I have more problems I’ll bring it here.

Cheers, Luís

Hello Andrei,

I still have two points to mention:

  1. In my code I transport a particle through several voxels in just one step, but never crossing a boundary between diferent shapes with diferent materials. Each shape must be divided into voxels in order to do the energy deposition. So I need to know at least the SafeDistance to the next boundary of a voxelized shape (not taking account the voxels) AND the Step through the voxels for compute the energy deposition. Do you know how can I do it?

  2. It seams that the newest ROOT (5.11/02) have a bug because my results using your code are different! But with ROOT 5.08/00 gives me the results exatly the same as you.
    See the result of ROOT 5.11/02:

 ===========> Event: 0
Start: phantom_phys_1
SafeDistance = 2.5      fStep = 0.5     NodeID = phantom_phys_1 Next point = (0, 0, 0.5)
SafeDistance = 2.5      fStep = 0.5     NodeID = phantom_phys_1 Next point = (0, 0, 1)
SafeDistance = 2.5      fStep = 0.5     NodeID = phantom_phys_1 Next point = (0, 0, 1.5)
SafeDistance = 2.5      fStep = 0.5     NodeID = phantom_phys_1 Next point = (0, 0, 2)
SafeDistance = 2.5      fStep = 0.5     NodeID = phantom_phys_1 Next point = (0, 0, 2.5)
SafeDistance = 2.3      fStep = 0.5     NodeID = phantom_phys_1 Next point = (0, 0, 3)
SafeDistance = 1.8      fStep = 0.5     NodeID = phantom_phys_1 Next point = (0, 0, 3.5)
SafeDistance = 1.3      fStep = 0.5     NodeID = phantom_phys_1 Next point = (0, 0, 4)
SafeDistance = 0.8      fStep = 0.5     NodeID = phantom_phys_1 Next point = (0, 0, 4.5)
On boundary     fStep = 0.3     NodeID = sliceZ_50      Next point = (0, 0, 4.8)
On boundary     fStep = 0.2     NodeID = world_1        Next point = (0, 0, 5)
SafeDistance = 0        fStep = 0.5     NodeID = world_1        Next point = (0, 0, 5.5)
SafeDistance = 0.5      fStep = 0.5     NodeID = world_1        Next point = (0, 0, 6)
SafeDistance = 1        fStep = 0.5     NodeID = world_1        Next point = (0, 0, 6.5)
SafeDistance = 1.5      fStep = 0.5     NodeID = world_1        Next point = (0, 0, 7)
SafeDistance = 2        fStep = 0.5     NodeID = world_1        Next point = (0, 0, 7.5)
SafeDistance = 2.5      fStep = 0.5     NodeID = world_1        Next point = (0, 0, 8)
SafeDistance = 2        fStep = 0.5     NodeID = world_1        Next point = (0, 0, 8.5)
SafeDistance = 1.5      fStep = 0.5     NodeID = world_1        Next point = (0, 0, 9)
SafeDistance = 1        fStep = 0.5     NodeID = world_1        Next point = (0, 0, 9.5)
On boundary     fStep = 0.5     Next point = (0, 0, 10)

Cheers, Luís

Hi Luis, see below:

If you want to compute the safety ignoring the voxels, you have to:

  • check if the current node is a voxel (meaning the particle is also in the “phantom” volume)
  • get a pointer to the global transformation matrix of the phantom volume
  • convert the global point in the local coordinates of the phantom
  • call the Safety() method for the shape of the phantom volume providing the local coordinates.
    This is a bit more tricky but I have attached the modified macro that does this.
  1. I have the CVS version of /geom and I cannot reproduce your second output. It looks like your phantom_phys is not divided at all, which can be a result in your macro if you call geomBuild(v=kFALSE)

Cheers,
geomtracktest.C (4.18 KB)

Hi Andrei,

The ROOT version from CVS is working fine, that problem occured somewhre after 5.08/00 until 5.11/02.
Thanks for all your advices.

Cheers,