Isotropic lambertian sampling from planar surface

_Geant4 Version: 4-11.2
_Operating System: Windows 11
_Compiler/Version: Visual Studio 17.9
_CMake Version: 3.29.2

Hello. I am hoping that someone can assist with a problem I am having with sampling an isotropic distribution. I am attempting to propagate geantinos toward a coincidence detector setup, by randomly sampling positions at the top planar surface of my world volume, 12.2m x 12.2m area (the area is MUCH larger than my detector area), and sample the momentum direction isotropically towards the detector. I have put my current code below. I defined the theta sampling based in this paper: []

I have ran my simulation with the world volume 1m above my top detector and 1cm above it, and I would have expected the outputs to be extremely similar if the sampling was indeed isotropic. However, they differ by around 18%, which is a lot. I am essentially trying to determine to effective area of the detector, which can only be done is the flux is isotropic. Any assistance would be greatly appreciated.

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) {
    G4double x0 = (G4UniformRand()-0.5)*12.2*m;
    G4double y0 = (G4UniformRand()-0.5)*12.2*m;
    G4double z0 = 0.555 * m;

    fParticleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0));

    G4double phi = 2 * CLHEP::pi * G4UniformRand();
    G4double sinTheta = sqrt(G4UniformRand());
    G4double cosTheta = sqrt(1 - sinTheta * sinTheta);

    G4double px = sinTheta*std::cos(phi);
    G4double py = sinTheta*std::sin(phi);
    G4double pz = -cosTheta;

    fParticleGun->SetParticleMomentumDirection(G4ThreeVector(px, py, pz));

Try this:

G4double cosTheta = G4UniformRand();
G4double sinTheta = sqrt(1 - cosTheta * cosTheta);

BTW, to specify a uniformly distributed direction within a hemisphere you can use G4RandomDirection(cosTheta) defined in G4RandomDirection.hh