_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);
}