Normal surface access

[froz1233]

How can I find the surface normal to a shape in my geometry in the World system of coordinates?
I am willing to get the cosinus between my particle entering or exiting a volume.
For this I have to compute the dot product between my particle direction ( in the world system of coordinates) with my shape normal at the crossing point (in the world system of coordinates).

Until now, I am trying to do this way:

G4StepPoint* aStepPoint = 0;
// OUT
if(aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary && aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "myvolume")
        aStepPoint = aStep->GetPostStepPoint();
//IN
else if(aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary && aStep->GetPostStepPoint()->GetPhysicalVolume()->GetName() == "myvolume")
        aStepPoint = aStep->GetPreStepPoint();
//NOT CROSSING VOLUME
else return;
G4ThreeVector pointInSurf = aStepPoint->GetPosition();
G4ThreeVector localPointInSurf = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(pointInSurf);
G4ThreeVector localNormal = aStep->GetPreStepPoint()->GetTouchableHandle()->GetSolid()->SurfaceNormal(localPointInSurf);
G4ThreeVector globalNormal = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().Inverse().TransformPoint(localNormal).unit();
G4double dot = globalNormal.dot(aStepPoint->GetMomentumDirection().unit());

I am not sure what I am doing wrong but it s quite easy to be wrong to decide if you need to use the post or pre step point, to change the coordinate system, to normalize the vector… I am a bit confuse with the local/global coordinate system and any help would be appreciated

[orru-nicolas]

I am also interested by this topic since I need the surface normal for particles crossing shape boundaries to obtain the crossing angle

[froz1233]

Dear Geant4 community,
Do you have a solution to my request or at least could you please review my code and let me know where I am mistaken?
Thanks in advance for your appreciated help.

[froz1233]

Dear @mkelsey ,
I have seen your answer in the thread Surface Normal in Complex Geometry - #2 by mkelsey
Would you have any input for my question?
Best regards

[mkelsey]

I find the global vs. local coordinate transforms confusing myself, and of course there are multiple conventions – is it “rotating” the vector, or is it rotating the coordinate system? Is the direct transform “from local to global” or is it “from global to local”?

My best recommendation (which is what I end up doing when I have to) is to add debugging output and look at what it’s doing. Print your position vector. Then print the computed “local position” vector. Does the latter correspond to what you think it should be, based on how you built your geometry? Print out the normal that was returned by the solid, and print out the global vector after the transform. Do they look right to you?

If you build G4 and your application with debugging enabled, you could run in GDB and look at the quantities without having to modify the source code.

One thing I can see that is wrong is that you are calling TransformPoint() but passing in what you intend to be a direction vector (the surface normal). You should be calling TransformAxis() for the surface normal.

[froz1233]

Thank you mkelsey, I will try your suggestion regarding TransformPoint() and TransformAxis().

I am already using gdb to test the quantities of interest but was just wondering how my code could be wrong. I performed some comparisons with other codes to assess surface fluence, having discrepancies while I get the correct volume fluence. Hence, it probably means there is a mistake in my way of assessing the surface normale, which is the tricky part.

Thank you.

[froz1233]

So, I have tried with a sphere and everything seems to be working well.
With a cylinder, there are some cases that do not work well, for instance, the following:
The cylinder is translated in X direction from the mother volume

We can see the track escaping the cylinder from one side but the normal vector found at the escape location is (0,0,1) as if it ran from the top of the cylinder. I do not understand well. Shouldn’t my code find the normal vector that should be (xn, yn, 0) ?
Only few tracks give such result with the cylinder, most of the tracks provide the correct normal vector.
Thanks for your help

[froz1233]

I am starting to think that GEANT4 will consider randomly the surface when the steppoint is close to multiple surfaces which is strange because I am using the steppoint that is on the surface I am interested in. In this case above, the steppoint is located on the circular shaped surface but the normal vector GEANT4 sends me back is linked with the plane surface on top of the cylinder.
Any idea ?

[froz1233]

Is there somebody able to respond?

[psaxioti]

I would suggest something like the following (based on your code from the original post).

G4StepPoint* aStepPoint = 0;
// OUT
if(aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary && aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "myvolume")
        aStepPoint = aStep->GetPreStepPoint();
//IN
else if(aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary && aStep->GetPostStepPoint()->GetPhysicalVolume()->GetName() == "myvolume")
        aStepPoint = aStep->GetPostStepPoint();
//NOT CROSSING VOLUME
else return;
G4ThreeVector pointInSurf = aStepPoint->GetPosition();
const G4VTouchable *aStepTouchable = aStepPoint->GetTouchable();
const G4AffineTransform *trans = aStepTouchable->GetHistory()->GetPtrTopTransform();
G4ThreeVector localPointInSurf = trans->TransformPoint(pointInSurf);
G4ThreeVector localNormal = aStepTouchable->GetVolume()->GetLogicalVolume()->GetSolid()->SurfaceNormal(localPointInSurf);
G4ThreeVector localDirection = trans->TransformAxis(aStepPoint->GetMomentumDirection())

G4double AngleToNormal = localDirection.angle(localNormal)

The main difference is the use of the aStepPoint for all the calculations (you were using PreStepPoint for localPointInSurf, localNormal and globalNormal).
Additionally, I exchanged the aStepPoint between IN and OUT in order to be the step point in your volume.

Furthermore, depending on what interests you I would say that in the case of OUT the localNormal should have a negative sign.

Hopefully, I am not thinking “bad” and the above works… Have a try!

[froz1233]

Thank you psaxioti, I just see your post now…
I did not have yet the time to test it but I plan to do a test during the Christmas holidays and I will submit the results in this thread.
Thanks again

[froz1233]

Thnak you @psaxioti. Using your corrections solved my issues.