Problem 1971

Indeed, G4ParticleHPThermalScattering::getMu is not sampling properly the angular distributions. I just histogrammed the sampled values of mu from a single distribution (using the same binning), and noticed that the lower and upper bins are about half filled than the rest. This is consistent with the Katy's description of the problem.

The implementation of this funciton is, however, unnecessarily complicated. I propose instead the following one, which is, in my opinion, more simple, clear and computationally efficient:

G4double G4ParticleHPThermalScattering::getMu(const E_isoAng *anEPM) const
	{
	const G4int nEdges = anEPM->n; // The stored value of n is actually the number of edges, not the number of bins, i.e., the size of vector isoAngle.

	const G4double x = G4UniformRand()*(nEdges+1);
	if ( x == nEdges+1 ) return 1.; // Just in case, the output of G4UniformRand() happens to be 1.

	const G4int i = G4int(x);

	const G4double MuL = ( i == 0 ) ? -1. : anEPM->isoAngle[i-1];
	const G4double MuH = ( i < nEdges ) ? anEPM->isoAngle[i] : 1.;

	return MuL+(x-i)*(MuH-MuL);
	}

The same test with the new implementation showed that the sampled values are uniformly distribuited over the entrire [-1,1] range.

The algorithm proposed by Aczel above assumes that the mu values provided in the library are bin edges. However, this is not the case for libraries processed with NJOY (such as G4NDL). In these libraries the mu values are bin averages, and as a result this algorithm will treat a library with N bins as though there were N+1 bins.

Hi Katy and Aczel,

Thank you for reproting this problem.
I've just commited the fix and it will be included future release of Geant4.

Tatsumi

Hi Tatsumi,

Katy and I have been taking a look at this. After a deep review of the issue I agree totally with her. The algorithm I proposed would be suitable if the values in the library were bin edges. This is what the Physics Reference Manual says and what G4ParticleHPThermalScattering assumes when reading the library. However, the current implementation of G4ParticleHPThermalScattering::getMu is consistent with the interpretation of the values as bin averages, which is what NJOY calculates as Katy said. The function getMu has the flaws that Katy pointed out. I corrected it and got a slight better result when compared the thermal flux with equivalent simulations with MCNPX. Here is my correction:

G4double G4ParticleHPThermalScattering::getMu(const E_isoAng *anEPM) const
	{
	G4double Mu = 0.;

	const G4int n = anEPM->n;

	const G4double x1 = G4UniformRand()*n;
	const G4int i = G4int(x1);

	// Sampling within the n-1 intervals between contiguous average mu values uniformly, 7/8 of the total probability.
	if ( i > 0 )
		{
		const G4double MuL = anEPM->isoAngle[i-1];
		const G4double MuH = anEPM->isoAngle[i];

		Mu = MuL+(x1-i)*(MuH-MuL);
		}
	// Sampling within the interval [-1; mu[0]] half the 1/8 remaining fraction of the total probability, and the other half within [mu[n-1];1.].
	else
		{
		const G4double x2 = G4UniformRand(); // Sampling a new random number within [0.;1] to uncorrelate the selection of the edge interval from the sampling of mu within the selected interval.

		if ( x1 <= 0.5 )
			{
			const G4double MuL = -1.;
			const G4double MuH = anEPM->isoAngle[0];

			Mu = MuL+x2*(MuH-MuL);
			}
		else
			{
			const G4double MuL = anEPM->isoAngle[n-1];
			const G4double MuH = 1.;

			Mu = MuL+x2*(MuH-MuL);
			}
		}

	return Mu;
	}

Either the manual is changed and both the library and the reading in G4ParticleHPThermalScattering are fixed to make them consistent with bin averages, or the library is generated with bin edges, which can be calculated from the averages. Katy and I agree it should be the first, as it is the direct output of NJOY.