Generic surface sampler

[tdixon97]

This implements the generic surface sampling following the algorithm of @jasondet et. al's paper [link].
No geantino are used instead only the methods G4VSolid::DistanceToIn(G4ThreeVector &p,G4ThreeVector &v) and DistanceToOut are needed. These essentially track the distance to the next surface along a line and so can give us the number of intersections.
I started on testing.

• unit tests (of the C++ code) for the basic methods (this could be expanded to more types of volume),
• some visualisation (it was tricky to get this to work well), for different simple shapes,
• also some plots in python of the vertices and a check of selecting points on each surface and checking the fraction is equal to the ratio of surface areas,
  The last test is implemented but for some reason fails for the Union of 2 G4Box, in addition I find some generated points not on the surface, I will investigate.
  Update : This was related to the origin of the bounding sphere which wasn't accounted for - is now fixed.

Once we have this sorted we should think how to speed it up for difficult to sample solids, one idea is to check whether its impossible for a direction vector to intersect with the solid (eg. it points the wrong way) and then avoid tracking.

[gipert]

Should we start adding some docstrings to the important functions? Can you give it a try? https://www.doxygen.nl/manual/docblocks.html

Please also add a link to the paper.

[tdixon97]

Should we start adding some docstrings to the important functions? Can you give it a try? https://www.doxygen.nl/manual/docblocks.html

Please also add a link to the paper.

I agree but there seems to be many options, maybe lets pick one ?

[tdixon97]

The number of PDF produced by the tests is now large. I think we should consider how to organise better the various confinement tests but also perhaps a strategy to merge the figures into one document?

[jasondet]

Yes, point in the +Z direction, then rotate about y by theta then about Z by phi to get a vector in the direction theta phi.

…

------------------------------------- Jason Detwiler Department of Physics, Box 351560 University of Washington Seattle, WA 98195 (206) 543-4054

On Mon, Jan 13, 2025 at 2:13 AM Toby Dixon ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In src/RMGVertexConfinement.cc <https://urldefense.com/v3/__https://github.com/legend-exp/remage/pull/226*discussion_r1912952060__;Iw!!K-Hz7m0Vt54!g2nRwyep0alWPCTbEBvCf7Z6zpV6nuymOClJ84_Yp9zvDyxvbnEsJ3hjo_8rA-xl9ZAFDddM1ITDfR7hEDtkJxM$> : > + + // start on z-axis, pointing down. + pos = G4ThreeVector(0.0, 0.0, bounding_radius); + dir = G4ThreeVector(0.0, 0.0, -1.0); + + // push in rho direction by some impact parameter + G4double diskPhi = 2.0 * CLHEP::pi * G4UniformRand(); + G4double diskR = sqrt(G4UniformRand()) * bounding_radius; + + pos += G4ThreeVector(cos(diskPhi) * diskR, sin(diskPhi) * diskR, 0); + + // now rotate pos and dir by some random direction + G4double theta = acos(2.0 * G4UniformRand() - 1.0); + G4double phi = 2.0 * CLHEP::pi * G4UniformRand(); + + pos.rotateY(theta); I took this from Mage we should check if it works well and also test it — Reply to this email directly, view it on GitHub <https://urldefense.com/v3/__https://github.com/legend-exp/remage/pull/226*discussion_r1912952060__;Iw!!K-Hz7m0Vt54!g2nRwyep0alWPCTbEBvCf7Z6zpV6nuymOClJ84_Yp9zvDyxvbnEsJ3hjo_8rA-xl9ZAFDddM1ITDfR7hEDtkJxM$>, or unsubscribe <https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/AARCTXUS22UGIH3FLIXDVZL2KOGTPAVCNFSM6AAAAABVBREJ4OVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDKNBWGIYDMNBYGQ__;!!K-Hz7m0Vt54!g2nRwyep0alWPCTbEBvCf7Z6zpV6nuymOClJ84_Yp9zvDyxvbnEsJ3hjo_8rA-xl9ZAFDddM1ITDfR7h1NZaiw0$> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[gipert]

The number of PDF produced by the tests is now large. I think we should consider how to organise better the various confinement tests but also perhaps a strategy to merge the figures into one document?

Yes, please feel free to reorganize things. If you come up with a proposal on how to merge things, I can implement something quick

[gipert]

I agree but there seems to be many options, maybe lets pick one ?

I like the first one:

/**
 *
 */

[tdixon97]

Well I think it would also be good to have some headings. And maybe also some text / captions could be included (eg. for geant4 vis where I havent been able to add text directly)

[tdixon97]

@jasondet Your algorithm generates N intersections and then picks one (or zero).
I understand that formally this keeps the sampled points independent. But I wonder if its really necessary for our applications.
In particular, this is a bit challenging since the maximum number of intersections has to be known in advance.

[tdixon97]

My geant4 visualisations don't work properly in CI, fine locally...

[jasondet]

@tdixon97 it is indeed probably not necessary to pick one of the N crossings at random -- at high stats you still get distribution that is essentially uniform over any surface, the only price you pay is that some consecutive points have some hint of a correlation (they lie along the same ray). That kind of correlation is probably irrelevant for our application (rad sims).

However inputting N_max is not difficult. We would just set it to something that seemed reasonable (like 16 or 32 or something) and then bumped it up by a factor of 2 and re-ran if we ever saw the error message that it was too small. Then we didn't have to worry about arguing whether the faint correlation was okay or not. So we opted to do that.

[tdixon97]

@tdixon97 it is indeed probably not necessary to pick one of the N crossings at random -- at high stats you still get distribution that is essentially uniform over any surface, the only price you pay is that some consecutive points have some hint of a correlation (they lie along the same ray). That kind of correlation is probably irrelevant for our application (rad sims).

However inputting N_max is not difficult. We would just set it to something that seemed reasonable (like 16 or 32 or something) and then bumped it up by a factor of 2 and re-ran if we ever saw the error message that it was too small. Then we didn't have to worry about arguing whether the faint correlation was okay or not. So we opted to do that.

I wrote our implementation so only one G4VSolid is ever sampled at a time Nmax can probably be usually something around 6 in our case (so less bad). However, since we weight the intersections by N/Nmax this still means we though away a lot of points and so are also using a lot of calls to GetRandom(). I guess that the vertex confinement is usually so fast this doesn't matter. But there are cases where it might...
I agree the correlation is clearly there but probably minimal.

[tdixon97]

Should now be ready @gipert