Custom Muon Angular Distribution

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


Hi all. I have hopefully an easy question to resolve. I have a custom anisotropic underground muon distribution that I am trying to use as an injection distribution to simulate the muon-induced secondaries underground. I am sampling the zenith and azimuth angles (as a pair) from my CDF. Currently, I generate a random spatial position on a hemisphere (just inside of the World Volume), and then I use my sampled angles to generate the momentum direction of the muon. I would greatly appreciate if someone could help me understand if this sampling method would properly reproduce the custom muon distribution for the secondary production.

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) {
    std::pair<double, double> sampled_angles = SampleAnglesFromCDF();
    double sampled_zenith_deg = sampled_angles.first;
    double sampled_azimuth_deg = sampled_angles.second;
    
    double particle_energy = FindClosestEnergy(sampled_zenith_deg, sampled_azimuth_deg);
    fParticleGun->SetParticleEnergy(particle_energy);
    
    double sampled_zenith = sampled_zenith_deg * CLHEP::pi / 180.0;
    double sampled_azimuth = sampled_azimuth_deg * CLHEP::pi / 180.0;

    const G4double sphereRadius = 12.851 * m;

    G4double theta = std::acos(G4Random::getTheEngine()->flat()); // Angle from Z-axis
    G4double phi = 2.0 * CLHEP::pi * G4Random::getTheEngine()->flat();

    // Calculate position
    G4double x = sphereRadius * std::sin(theta) * std::cos(phi);
    G4double y = sphereRadius * std::sin(theta) * std::sin(phi);
    G4double z = sphereRadius * std::cos(theta);

    G4ThreeVector position(x, y, z);

    G4cout << "Generated position magnitude: " << position.mag() / m << G4endl;
    //G4cout << "Generated position (x, y, z): " << position.x() / m << " " << position.y() / m << " " << position.z() / m << G4endl;

    fParticleGun->SetParticlePosition(position);

    // Set momentum direction
    double cosTheta2 = std::cos(sampled_zenith);
    double sinTheta2 = std::sin(sampled_zenith);
    double cosPhi2 = std::cos(sampled_azimuth);
    double sinPhi2 = std::sin(sampled_azimuth);

    G4ThreeVector direction(sinTheta2 * cosPhi2, sinTheta2 * sinPhi2, -cosTheta2);
    fParticleGun->SetParticleMomentumDirection(direction);

    fParticleGun->GeneratePrimaryVertex(anEvent);
}

I have a suggestion: you can keep one additional scorer (same as the world’s material) to get the distribution.