Created attachment 656 [details] Histogram containing the exit angle of neutrons in a water phantom for Geant4 10.6 (and later) yields wrong neutron fluences, as compared to Geant4 10.5 and earlier. The reason is a regression: the momentum change should be randomly rotated around the Z-axis, but this was removed. The issue had been reported here: https://geant4-forum.web.cern.ch/t/thermal-neutron-physics-drastic-change-from-10-5-to-10-6/4011 To further test the issue, I created a histogram with the output angle of neutrons with respect to the Z-axis in elastic hadronic processes (divided by sin(theta)). In 10.5 and the patched 10.6 runs, the curves are almost flat (except for the expected bias towards theta=0 since the neutrons are shot in the +Z direction). In 10.6, almost all neutrons exit in the plane perpendicular to the Z-direction. This figure is attached. The system is a water phantom, the source consists of 10 keV neutrons in the +Z direction. Please ask if further details are needed.
Created attachment 657 [details] Patch adding the removed randomization around the Z-axis The patch consists of selected lines from Geant4 10.5.1.
This has already been fixed according to bug report 2290. It will be included in the first patch to Geant4 10.7
Good to know. I will add a comment in the discussion.
I have updated my geant4 to geant4.10.7.2. But I haven't seen the bug fixed! Is there something wrong? Here is the part of the code in greant4.10.7.2: // Check the result for catastrophic energy non-conservation // cannot be applied because is not guranteed that recoil // nucleus is created // result = CheckResult(theProj, targNucleus, result); // directions G4ThreeVector indir = track.GetMomentumDirection(); G4ThreeVector outdir = result->GetMomentumChange(); if(verboseLevel>1) { G4cout << "Efin= " << result->GetEnergyChange() << " de= " << result->GetLocalEnergyDeposit() << " nsec= " << result->GetNumberOfSecondaries() << " dir= " << outdir << G4endl; } // energies G4double edep = std::max(result->GetLocalEnergyDeposit(), 0.0); G4double efinal = std::max(result->GetEnergyChange(), 0.0); // primary change theTotalResult->ProposeEnergy(efinal); if(efinal > 0.0) { outdir.rotateUz(indir); theTotalResult->ProposeMomentumDirection(outdir); } else { if(part->GetProcessManager()->GetAtRestProcessVector()->size() > 0) { status = fStopButAlive; } else { status = fStopAndKill; } theTotalResult->ProposeTrackStatus(status); }
