Using v4.10.3, in the function getMu (line 715 in G4ParticleHPThermalScattering.cc) the scattering angle of events is determined by randomly generating the cosine of the scattering angle, mu, within equal probability bins. The N average mu values given by the data library are used to define the limits of N-1 bins. The overall distribution at the edges of the mu range (within (-1,mu1) at the bottom, and (muN,1) at the top) must have the same cumulative probability as the other bins, and the division of events between the two edges is determined by a variable x. The variable x belongs to (0,1). A ratio is defined corresponding to the width of the low edge bin divided by the width of both edge bins; events are placed in the low edge bin if x <= ratio, and in the top edge bin if x>=ratio. This distributes events proportionately to the width of each edge bin. However, within the edge bins the mu value is generated using mu = ( mu_h - mu_l ) * x + mu_l where mu_h,l are the upper and lower bounds of the respective bin. The use of x to determine the placement of the event within the bin is problematic, because within each edge bin x is no longer distributed within (0,1); it belongs to (0,ratio) in the lower bin, and (ratio, 1) in the upper bin. Because of this, portions of each of these bins remain unsampled. This can be fixed either by replacing x with a new uniform random number, or by making the replacement [x -> x/ratio] for the bottom edge and [x -> (x-ratio)/(1-ratio)] for the top edge. Also, for consistency the ratio should be set equal to 0.5, rather than weighted by the bin widths. By using the average mu values as bin limits, the algorithm outlined above assumes that the average mu value within each bin is equivalent to the median. In this case the probability of an event being placed in the lower edge bin should be the same as that for the higher edge bin.
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.
