Add distance to nearest surface field to GermaniumOutputScheme

[EricMEsch]

Fixes #154

For some graphical illustration i made a Google Presentation: https://docs.google.com/presentation/d/1gId7s7pBYLD5PBw2EEabf4iIHPTsUEJl24yW60s9TQU/edit#slide=id.p1 (Sorry, in PowerPoint this would look way better)

TLDR:

• The DistanceToOut() and DistanceToIn() functions are used by G4Navigator itself. This means they should be very reliable (If we start to question the navigation, we start to question the whole simulation). This also means they are computed for every single step, meaning some additional calls in Sensitive Detectors should not slow down the simulation. (I thought about just using the safety variable of the G4Navigator, but it also does additional stuff that we dont want, like ignoring recently exited volumes)
• DistanceToOut() and DistanceToIn() functions are purely virtual in G4VSolid and are implemented differently for each solid. This means we have to check the implementation of each solid we want to use. I checked the most common:
  • For G4Box it is accurate (I checked the source code and some simple test simulations). Should be accurate for all cartesian solids.
  • For G4GenericPolycone (which are all detectors i think?) it should compute the exact geometrical distance (check the presentation for how) according to the source code. The only exception is if the calculated distance is < 0.5 * kCarTolerance (which is the safety tolerance of Geant4 to account for floating-point inaccuracies) it returns 0.

Although i did not expect the performance to be severely affected, i still did some tests using the L200.gdml created with the pygeom-l200 package, registered every single detector and simulated runs of 1e6 10 MeV gammas isotropically in the center of the Strings. This means not always a detector was hit, but in ~40% of the cases at least one was hit.

• ReMaGe before any changes processed: 1084 events/s (averaged over 3 runs, taken at different times) (On my local pc)
• ReMaGe after the changes processed: 1076 events/s (averaged over 3 runs, taken at different times) (On my local pc)

Open Questions:

• I checked the output for G4GenericPolycones by eye and they looked reasonable. But i did not do any quantitative cross-checks.
• The only way i can think of for doing some checks is to calculate the distance myself using Python (basically implementing the exact code Geant4 already has, that i am pretty certain is correct) and compare with the Geant4 output. Which i am not sure is necessary?
• Do we want to have this as an option hidden behind a macro, or should this be permanent part of the OutputScheme (as it is currently)?

[EricMEsch]

Ah, right, the tests would also have to be adjusted to the ntuple size. After we have decided if this is supposed to be a permanent part of the tuple or only added via macro

[ManuelHu]

The only way i can think of for doing some checks is to calculate the distance myself using Python (basically implementing the exact code Geant4 already has, that i am pretty certain is correct) and compare with the Geant4 output. Which i am not sure is necessary?

I have another idea for a more definitive check:

1. Create a simple geometry with one Polycone detector at the origin (i.e. wih some python code and legend-pygeom-hpges). You might be interested in some WIP docs on this
2. Confine the gps inside the detector and simulate some 1 MeV (or maybe even less) gammas; record the output.
3. using python, read the output and convert x,y,z to r,z.
4. plot profiles of the hits, cutting on the SurfaceDist_in_m in the output, and see if you have visual outliers.

edit: I mean some plots like this in the WIP reboost docs Of course without the difference in surface types

[EricMEsch]

Ah, that looks like a very good idea! I will have a look at the docs next week, thanks for the references!

[ManuelHu]

Apart from those comments, thank you very much for the test. We already did assume this, but had not looked deeper in it.

The overall performance looks also quite promising. It might be also worth comparing the performance again with the geometry with only a single detector. (this "microscopic" check would be - sort of - an upper bound on simulation slowdon)

Do we want to have this as an option hidden behind a macro, or should this be permanent part of the OutputScheme (as it is currently)?

what do you think, @tdixon97?

For G4GenericPolycone (which are all detectors i think?) [...]

yes and no. The detectors are polycones, but might have subtractions of other volumes (i.e. boxes). But like you, I would also assume that Geant4 is pretty well validated in the geometry part.

I would be slightly more worried about performance with the subtractions: It might also be interesting to repeat the "microscopic" performance measurement just with P00664b. (The visual check of distance does not work with this detector, as the radial symmetry is broken.)

[gipert]

Thanks a lot @EricMEsch! This looks very promising.

For G4GenericPolycone (which are all detectors i think?) it should compute the exact geometrical distance (check the presentation for how) according to the source code.

If you are right we don't even need to calculate distances in post-processing.

Although i did not expect the performance to be severely affected, i still did some tests using the L200.gdml created with the pygeom-l200 package, registered every single detector and simulated runs of 1e6 10 MeV gammas isotropically in the center of the Strings. This means not always a detector was hit, but in ~40% of the cases at least one was hit.

I guess that a better figure for the performance penalty would be shooting electrons uniformly in germanium.

As mentioned by @ManuelHu, it would be good to check the results with the few detectors made with subtractions.

[EricMEsch]

I will do the "microscopic" part mentioned by @ManuelHu and can also change gammas to electrons and shoot them at the detector as you mentioned @gipert. I can also make additional tests for the subtraction solids.

I would not expect there to be any issues with subtractions, as the DistanceToOut() function of G4SubtractionSolid is just a few lines of code calling to the original volumes

G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const 
{
  G4double dist=0.0;

if( Inside(p) == kOutside )
  { 
   ... // (there is a warning print here)
  }
  else
  {
     dist= std::min(fPtrSolidA->DistanceToOut(p),
                      fPtrSolidB->DistanceToIn(p) ) ; 
  }
  return dist; 
}

[tdixon97]

This is great, I agree the test with electrons is more meaningfull. Unfortunately, this cannot replace totally computation of distances in post processing since you need distances to only some surfaces. But probbaly, having this precomputed still enables a dramatic speed up.

…

________________________________ From: Eric Esch ***@***.***> Sent: Saturday, December 7, 2024 5:09 PM To: legend-exp/remage ***@***.***> Cc: Dixon, Toby ***@***.***>; Mention ***@***.***> Subject: Re: [legend-exp/remage] Add distance to nearest surface field to GermaniumOutputScheme (PR #182) ⚠ Caution: External sender I will do the "microscopic" part mentioned by @ManuelHu<https://github.com/ManuelHu> and can also change gammas to electrons and shoot them at the detector as you mentioned @gipert<https://github.com/gipert>. I can also make additional tests for the subtraction solids. I would not expect there to be any issues with subtractions, as the DistanceToOut() function of G4SubtractionSolid is just a few lines of code calling to the original volumes G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const { G4double dist=0.0; if( Inside(p) == kOutside ) { ... // (there is a warning print here) } else { dist= std::min(fPtrSolidA->DistanceToOut(p), fPtrSolidB->DistanceToIn(p) ) ; } return dist; } — Reply to this email directly, view it on GitHub<#182 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ANSET46CWGQTWMPBBFHTNM32EMMTDAVCNFSM6AAAAABTGHTTESVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKMRVGIZDQMRWGE>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[EricMEsch]

I have done the check @ManuelHu proposed, for Detector V02160A (The only subtraction solid i think?) (left) and detector B00035B (right). Looks promising:

I am not 100% sure about the subtraction solid, as detector one would expect a small dent at the bottom of one side, which is not shown here. ( And also not represented in the scatter cause i did just randomly assign a +- like in reboost)

I will probably benchmark tomorrow, as for that i want to restart my pc to have a clean environment.

Question: Where should the python code for generating this plot go? As i guess it is pretty important to have the opportunity to validate this at any time, but also this can not be part of a Pipeline, as it is a visual confirmation. Also the code depends on legendhpges , pyg4ometry and legendmeta, which is not supposed to be part of remage

[gipert]

P00664B and V02160A are the two only cylindrical asymmetric detectors. If you made the plots with legend-pygeom-hpges the asymmetric features will not be shown... But if you are plotting simulated steps then I would expect to see some faint deficit at one corner (the "dent" you mentioned).

Question: Where should the python code for generating this plot go? As i guess it is pretty important to have the opportunity to validate this at any time, but also this can not be part of a Pipeline, as it is a visual confirmation. Also the code depends on legendhpges , pyg4ometry and legendmeta, which is not supposed to be part of remage

This is a good question, the thing is essentially tracked in #179. I will try to find some time soon to figure out how to integrate some Python tools.

[tdixon97]

@EricMEsch are you sure the option was not set to treat the subtraction solid as a polycone ?
Otherwise, this all looks very good to me, another nice check might be comparing the distances from Geant4 and legend-hpges (this can be done with the distance_to_surface method of the HPGe class. This would also be a nice check on our python implementation.

[EricMEsch]

are you sure the option was not set to treat the subtraction solid as a polycone ?

Pretty sure, Pyg4ometry shows the "dent", the .gdml lists a subtraction solid and Geant4 also shows the subtraction.

I used the distance_to_surface method to compare just like you proposed @tdixon97, altough this seems to not be implemented for subtractions yet, so i can only compare results for B00035B. I basicly did

print(np.sum(np.abs(Dist_G4 - Dist_HPGe) > 1e-6)
> 0
print(np.sum(np.abs(Dist_G4 - Dist_HPGe) > 1e-7)
> 65744112

So the difference is smaller than 1 nm (1e-6 mm). (It might be relevant to note that the original distance in Geant4 was saved in m)

Regarding the performance: The ratio of Events_per_sec_post/Events_per_sec_pre (with post meaning after applying the changes in this PR) is given as: 91% for the subtraction solid (V02160A) and 93% for B00035B

Some more information on that:

• V02160A:
  • For 1 MeV gammas sampled inside the detector: 17241 events/s post vs 18882 events/s pre
  • For 200 keV electrons sampled inside the detector: 42792 events/s post vs 47169 events/s pre
• B00035B:
  • For 1 MeV gammas sampled inside the detector: 25746 events/s post vs 27266 events/s pre
  • For 200 keV electrons sampled inside the detector: 47117 events/s post vs 50933 events/s pre

[gipert]

I guess the effect of the "dent" could be hard to see in that plot, mixed with the rest of the polar angle. Could you try to only select steps close to that region?

[gipert]

A request from the point of view of coding conventions: remage tries to adopt snake casing (snake_case) as much as possible (despite Geant4 predominantly using camel case (snakeCase)). Could you update this PR?

[EricMEsch]

I am not sure ReMaGe is really consistent considering camel case or snake case (e.q. fHitsCollection), and i guess it is also pretty confusing when all of the Geant4 stuff is in camel case 😅

[gipert]

I am not sure ReMaGe is really consistent considering camel case or snake case (e.q. fHitsCollection), and i guess it is also pretty confusing when all of the Geant4 stuff is in camel case 😅

Yes you are right, it's not consistent but there is some logic... for example, class members still follow the fMember convention like in G4, otherwise it would be very wierd...

[ManuelHu]

so, to summarize, I think there are still two open issued with this:

• where is the dent in the geant4 derived distance distributon? -> probably need to look at slices at a specific angle in r,z,phi coordinates?
• should we have a flag to disable the calculation and output field?

or is there still anything else?

[EricMEsch]

Yes, regarding the dent i have made a new plot, where i did a cut on (np.abs(xPos_det0) > 34) & (np.abs(yPos_det0) < 10) (to only get locations close to the dent) and displayed the side based on the +- sign of the xPos, rather than randomly assigning it. I have also changed the colormap, to allow for better differentiation between "0 distance entries" and "no entries" and reduced the markersize to be more accurate.

this leaves the remaining question

should we have a flag to disable the calculation and output field?

[gipert]

should we have a flag to disable the calculation and output field?

Maybe not really needed, since we get basically no performance penalty.

I think this is good to merge. OK?