We found that the behavior of Geant4 changes depending on whether the USolids/VecGeom option is enabled or not. In particular, generating surface points via calling GetPointOnSurfaces() for G4Tubs solids is affected by this. The following two issues have been seen:
VecGeom provides its own implementation of SamplePointOnSurface(), which Geant4 calls into from GetPointOnSurface() when USolids is enabled. It turns out that the VecGeom version of G4Tubs chooses points by using r = random(rmin, rmax). This produces a flat distribution in r and not in the area which goes as r². As a consequence the density of random surface points is larger for small radii which causes a bias towards the center of a cylinder’s face. In Geant4’s definition of G4Tubs this issue has been resolved long ago (version G4-9.6). We confirmed the difference for both G4-10.06.p03 and G4-10.07.p03 with and without USolids. The plots below illustrate the difference for G4-10.07.p03.
VecGeom does not use the CLHEP random engine, but instead calls the standard C++ rand() function. This is done without doing any seed initialization, so every job will get exactly the same sequence of random points from VecGeom calls. This also means that the Geant4 interface for initializing the engine or getting the engine state does not affect VecGeom at all. There are alternative random engines in VecGeom (e.g. by making use of CUDA), but the issue of overriding Geant4’s default random engine remains the same.
Attached is a simple example demonstrating the call of GetPointOnSurface() to create N random points on the surface of a G4Tubs solid. The example was used to visualize the plots of 10k surface points displayed above. When executed multiple times with Geant4 using USolids, the same sequence of random points is created ignoring the initialization of the Geant4 random engine. surfacePoints.cc (1.6 KB)
We are concerned that these issues, combined with the fact that the main LCG installation of Geant4 make use of the “experimental” USolids feature, could be affecting a broader community of Geant4 users. We have also submitted a bug report to the VecGeom developers here.
Regarding having a separate generator for this kind of computation, in general, it does not represent a problem, as it does NOT affect tracking in any way and therefore event reproducibility… you certainly have a peculiar use-case here in having the selection of the starting positions depending on the random positions on the surface…
It is true that the separate VecGeom random generator does not affect event reproducibility, but it does affect event uniqueness. If we run multiple jobs to get high statistics, then the sequence of VecGeom randoms will be identical across all of those jobs, even if we use use the Geant4 “setSeeds” or equivalent methodology to initialize the Geant4 random engine.
As to the unusual use case. We (and I imagine other low-background experiments) use GetPointOnSurface() to generate events which model surface contaminants (such as radon adsorption or dust) on our apparatus components.
The function accesses and modifies internal state objects, which may cause data races with concurrent calls to rand or srand.
Some libraries provide an alternative function that explicitly avoids this kind of data race: rand_r (non-portable).
C++ library implementations are allowed to guarantee no data races for calling this function.
I guess that we would have to investigate the different compilers (GCC, LLVM, Intel, etc.) individually to determine whether rand is thread safe or not.
Use of GetPointOnSurface() is only foreseen for static generation of points for use in detecting overlaps or computing static quantities. This is the case also in vanilla Geant4 solids, as the generation of points is NOT happening through the static random generator, but using a dedicated fast algorithm.
For your specific use-case, you need to implement your own sampling of points on the surface for setting the initial position.
For all solids it is also possible to generate pseudo-random points lying on their surfaces, by invoking the method
G4ThreeVector GetPointOnSurface() const
which returns the generated point in local coordinates relative to the solid. To be noted that this function is not meant to provide a uniform distribution of points on the surfaces of the solids.
If GetPointOnSurface() is not intended to be used by user code such as source generation, it would be good to say that explicitly. As it is, it appears to be the obvious method for simulating surface radiocontamination.
Yes, right. And this was made on purpose to make more efficient such computation for detecting overlaps in complex setups. The generation of points on surface was never meant for use associated with tracking.
The documentation will be updated to explicitly mention this.
Sigh. So Geant4 no longer provides a mechanism for generating sources to model surface contamination. I suppose that we could copy over all of the past implementation of GetPointOnSurface() (which includes optimizations for many of the basic volume types) in order to recover the performance of G4 10.6 moving forward.
inline G4double G4QuickRand()
{
static const G4double f = 1. / 4294967296.; // 2^-32
// Algorithm "xor" from p.4 of G.Marsaglia, "Xorshift RNGs"
static G4ThreadLocal uint32_t y = 2463534242;
uint32_t x = y;
x ^= x << 13;
x ^= x >> 17;
x ^= x << 5;
y = x;
return x * f;
}
As a temporary easy solution I would suggest to reimplement G4QuickRand() in your local version of Geant4:
Thanks, Gabriele. As we are a distributed experiment (CDMS), modifying “our local version of Geant4” would involve multiple sites, and would not deal with the problem of centralized builds such as the LCG distribution which is used (via CVMFS) by many of our groups.
I see that in 10.07, G4QuickRand() is inlined, which means we also can’t override it in our own library.
It is not at all obvious how we can solve this critical blockage to our simulations. We have a very complex apparatus geometry (not as complex as ATLAS, of course), using many different kinds of solids, with many instances of complex booleans. It is not feasible, nor do I think it is reasonable to expect, that experiments should be expected to redesign and reimplement the complex code already provided in all of the classes’ GetPointOnSurface() implementations.
We also cannot simply move to using USolids and modify our unique seeding with a call to srand(). We have already found three different bugs (or badly designed features) in USolids. How many more will there be?
Geant4 has been reliable and consistent for us for a decade. It is not obivious yet if there is a way for us to do the simulations we need to do to model backgrounds and determine our detector sensitivity.
In case where there is no method to set a seed, it’s always possible to emulate seeding, for example:
// Set seed for default random number engine
...
// Make seeding for G4QuickRand()
const uint32_t nloop = G4UniformRand()*4294967295; // nloop < 2^32
G4double rnd = 0.;
for (uint32_t i = 0; i < nloop; ++i) { rnd = G4QuickRand(); }
if (rnd > 0.5) rnd = 1. - rnd;
I hope that a not very long period of G4QuickRand() (2^32-1) is not a problem.
We have avoided these problems, both G4QuickRand and USolids, by implementing the general (shape-independent) algorithm described by Detweiler et al. (IEEE Trans.Nucl.Sci.55:2329-2333,2008, [0802.2960] A Generic Surface Sampler for Monte Carlo Simulations). This was very straightfoward, using strictly member functions of G4VSolid and vector arithmetic. It does eat multiple random values per selected point, but that’s true of the existing implementations as well.
