Dear all, in addition to what was reported in https://bugzilla-geant4.kek.jp/show_bug.cgi?id=2352 there seem to be more (most?/all?) photons generated from G4NDL/Capture/FSMF6 data with way too high energies. This is much more prominent in 10.6 and higher compared to 10.1 since the G4NDL 4.6 has simply a lot more isotopes with energy-angle correlation data in G4NDL/Capture/FSMF6 compared to G4NDL 4.5. One particular obvious case is thermal neutrons on Fe56 in Hadr03 with any _HP physics list. This gives in each event N mono-energetic photons all having as energy the Q-value of the reaction (7.6457 MeV) with the probability of N falling exponentially: 48% for N=1, 25.4% for N=2, 16.6% for N=3, 6.6% for N=4, 2.6% for N=5, 0.7% for N=6, 0.1% for N=7. The reported Q-value is then 14.75 MeV in 10.6 - almost twice the correct number. I have not tested all FSMF6 isotopes but Cu63 is another one with dramatic energy increase from photons from neutron capture. Since both copper and iron are common shielding materials this impacts radiation background estimates for detectors in high-energy physics experiments significantly. At the same time for these the detailed angular distributions might not be as important as the energy ... So I am wondering if a) the FSMF6 data (or its treatment) could be fixed in 10.7 and 11 and/or b) if the FS capture data could be made available in addition (and with a switch to use as environment variable) in G4NDL 4.6 to at least get the energies correct without having to completely bypass the HP treatment by setting the environment variable G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION ... Thanks & Cheers, --Sven
Dear all, I believe I found one problem contributing to the wrong photon energies. It is this line: https://geant4.kek.jp/lxr/source/processes/hadronic/models/particle_hp/src/G4ParticleHPContAngularPar.cc?v=10.6.p3#L196 the delta is calculated using theAngular[i] but the loop is actually incrementing j not i! The index i is initialized to 0 and hence all discrete lines receive the probability of the first and most energetic line. While this is clearly a bug it does not alone explain why all photon energies are the same value for thermal neutron capture on Fe56 in 10.6 - it would result just in equal probability for all of the discrete lines ... This bug was also present in 10.1 and likely explains some discrepancies observed in that version. Cheers, --Sven
Dear all, digging a bit deeper now it becomes clear that the interpolation approach in G4ParticleHPContAngularPar is broken. The discrete energies are ignored and the partial code that is available is broken (even if it would run). Picking again version 10.6 (but all available versions are equally broken) the bad code is here: https://geant4.kek.jp/lxr/source/processes/hadronic/models/particle_hp/src/G4ParticleHPContAngularPar.cc?v=10.6.p3#L736 The nDiscreteEnergies from here: https://geant4.kek.jp/lxr/source/processes/hadronic/models/particle_hp/src/G4ParticleHPContAngularPar.cc?v=10.6.p3#L745 is 0 since the set to store the interpolation in is not yet filled. Even if it would be available the following lines make no sense. The indices ie1 and ie2 are set to -1 if the discrete energy to interpolate the angular data for is found and to 0 if it is *not* found in the sets to average. They should be -1 if not found and set to the actual respective index of the discrete energy if found ... Even if this would be correct - directly after this part whatever would have been filled in theAngular is deleted again here: https://geant4.kek.jp/lxr/source/processes/hadronic/models/particle_hp/src/G4ParticleHPContAngularPar.cc?v=10.6.p3#L790 So any FSMF6 dataset with discrete energies is currently ignored. Since Fe56 has discrete energies only for thermal neutron energies - this ends up with just the max energy for the produced gammas. In 10.1 the NeutronHP was still in use instead of ParticleHP. The interpolation as described above seems only present in ParticleHP and hence in 10.1 the discrete energies are not ignored in 10.1 ... but NeutronHP also suffered from the [i] instead of [j] bug mentioned in the previous comment. Maybe it is possible to revert to the NeutronHP interpolation scheme? Thanks & Cheers, --Sven
Created attachment 749 [details] patch for particle_hp to fix the FSMF6 treatment Dear all, please find attached a proposed patch for the described problem. There were several issues: 1) the index issue (´i´ instead of ´j´) in building the cumulative probabilities 2) only for a new interaction the interpolation should be calculated from scratch - not for every photon produced in the same interaction 3) building the interpolation table was completely broken for discrete energies 4) ClearHistory() should not be called after the Sample() call because Sample() is invoked for each photon and the remaining energy needs to be known 5) the ClearHistory() needs to be called on the cached angular parameters when the parent class' ClearHistory() is invoked This seems to work in 10.7.3 for thermal neutrons on Fe56 with G4NDL4.6. I have not tried other targets. Cheers, --Sven
Thank you for submitting this fix. It has been implemented and, after testing, will appear in the next Geant4 patch.
The code failed in run testing with a heap corruption error. This is apparently due to code in G4ParticleHPContAngularPar::BuildByInterpolation where the interpolation table is re-built. I copied your segments of code directly and they compile and run correctly in some cases but not others. I have not been able to pinpoint the error. One possible problem is that you are working in Geant4 10.6, while I am working in Geant4 11.0 patch 01. Please try updating your Geant4 version, re-building and re-running your code. If you will send your versions of the entire files G4ParticleHPContAngularPar.hh and G4ParticleHPContAngularPar.cc, I will try to spot the problem there.
Hi Denis, the patch (as the name of the patch suggests) is for 10.7.3. I have not checked the code vs. 10.6 or 11.X. If you can confirm that there is a problem with the patch in 10.7.3 I can try to find out more (or send you the full patched file from 10.7.3 should that still be needed ...). Cheers, --Sven
Geant4 does not test against old releases, so even if I did build and test 10.7 with the fix, it wouldn't be accepted for the current release without building and testing in 11.0.
(In reply to dennis.herbert.wright from comment #8) > Geant4 does not test against old releases, so even if I did build and test > 10.7 with the fix, it wouldn't be accepted for the current release without > building and testing in 11.0. o.k. but if you run the test that complained in 11.0 in 10.7.3 what do you get? I'd need a bit more info on what failed in what condition ... Thanks, --Sven
I don't have 10.7.3 set up. I would have to download and rebuild it, then rebuild and rerun the tests. I can do this but it will take me a while and I will not have much time for this in the near future. In the meantime, I'm asking for the fix to be removed from Geant4 testing so it doesn't hold things up there.
(In reply to dennis.herbert.wright from comment #10) > I don't have 10.7.3 set up. I would have to download and rebuild it, then > rebuild and rerun the tests. I can do this but it will take me a while and > I will not have much time for this in the near future. In the meantime, I'm > asking for the fix to be removed from Geant4 testing so it doesn't hold > things up there. Hi Dennis, then at least - what test failed for you in 11.0? Could you send me a log? ATLAS does not use Geant4 11.0 - so I'd have to download and build it as well ... but even then I'd need to know what test failed ... Thanks, --Sven
Created attachment 755 [details] patch for particle_hp in Geant4.11.0.1 to fix the FSMF6 treatment Hi Dennis, all, I downloaded 11.0.1 and found one difference w.r.t. 10.7.3 that impacts the previous patch: The remaining energy is updated after each sampled photon only if adjustResult is set to true. This in turn is set by the DoNotAdjustFinalState flag and that used to be false in my 10.7.3 tests but in 11.0.1 is true by default. The bug is that this is coupled to the adjustResult flag at all. Regardless of that switch photons should be generated only if there is enough remaining energy to produce them. So the new patch supersedes the previous one and removes that dependency. The Q-value looks better with DoNotAdjustFinalState set to false - but not totally crazy for true (with the new patch). I don't know if this fixes the test-jobs in 11.0 but at least for Fe56 neutron capture it works in 11.0.1 ... Cheers, --Sven
I'll test it with test04.
I tried the new code with Geant4 test04. It still fails with the following: test04(7449,0x70000530c000) malloc: Heap corruption detected, free list is damaged at 0x600021a49b50 *** Incorrect guard value: 4567829261830659808 test04(7449,0x70000530c000) malloc: *** set a breakpoint in malloc_error_break to debug ### CAUGHT SIGNAL: 6 ### address: 0x7ff816b3d00e, signal = SIGABRT, value = 6, description = abort program (formerly SIGIOT). test04(7449,0x700005289000) malloc: Heap corruption detected, free list is damaged at 0x600021a4cd30 *** Incorrect guard value: 4567474613092648320 test04(7449,0x700005289000) malloc: *** set a breakpoint in malloc_error_break to debug At this point, I think it's best if you send me both of the complete .cc files. That way we could at least rule out the possibility of me making a mistake inputting the code to my version of Geant4.
Created attachment 757 [details] patched G4ParticleHPContAngularPar.cc for 11.0.1
Created attachment 758 [details] patched G4ParticleHPContEnergyAngular.cc for 11.0.1
Hi Dennis, I've just attached both patched .cc files suitable for 11.0.1 Cheers, --Sven
Created attachment 760 [details] test04
Created attachment 761 [details] test04
These files failed as well. Test04 runs successfully for a while before failing and if I make small, trivial code changes, it fails in different places. I've attached a tar file of test04 so you can build and test it. I'll keep looking, but it seems like a nasty problem. By the way, your G4ParticleHPContEnergyAngular.cc by itself does not cause any problem. G4ParticleHPContAngularPar.cc is also OK if I only include your changes up to and including the swap of indices from i to j in theAngular. Beyond that, it fails.
Created attachment 763 [details] patch for already patched G4ParticleHPContAngularPar.cc removing extra max energy point Hi Dennis, thanks for the test cases! I finally managed to get a segv too. And I think I found one problem. I was wondering before, why there is an extra energy point added to the list of energies ... Removing that artificially added endpoint runs now for 5e6 neutron captures (it crashed before after ~1e5). Probably the random number normalization is to blame ... Could you please try to patch the already patched G4ParticleHPContAngularPar.cc with the latest patch (or simply comment out line 751) and test again? Thanks & Cheers, --Sven
It got much further this time, but failed while doing the Pb tests. G4WT1 > /mydet/SelectMaterial G4_Pb test04(52759,0x700002d2e000) malloc: Heap corruption detected, free list is damaged at 0x6000219a7200 *** Incorrect guard value: 4567829261830659808 test04(52759,0x700002d2e000) malloc: *** set a breakpoint in malloc_error_break to debug ### CAUGHT SIGNAL: 6 ### address: 0x7ff816b3d00e, signal = SIGABRT, value = 6, description = abort program (formerly SIGIOT). test04(52759,0x700002c28000) malloc: Heap corruption detected, free list is damaged at 0x6000219a05c0 *** Incorrect guard value: 4567474613092648320 test04(52759,0x700002c28000) malloc: *** set a breakpoint in malloc_error_break to debug
Created attachment 764 [details] patch for particle_hp to fix the FSMF6 treatment for G4.11.0.1 Hi Dennis, thanks for the next round of testing! Indeed, I can confirm to see the same on Pb207 with neutrons at 1.1 keV. The ultimate reason I believe was the vector of G4ParticleHPList objects vAngular. Making that a vector of pointers to such objects is much more robust. I also added some other measures to increase robustness - a copy constructor of the angular parameters since in the interpolation potentially one of the maps being interpolated is also the target. I've attached a new patch for all files needed relative to the unpatched G4.11.0.1 sources. So this one obsoletes all previous patches/codes. Could you please test once more with this? Thanks & Cheers, --Sven
Would you send the two .cc files, please? That would leave much less chance of my making a mistake in doing the patch.
Created attachment 767 [details] gziped tar archive with the 3 modified files (one .hh, two .cc) for G4.11.0.1 Hi Dennis, I've attached a .tgz file containing all 3 files that changed: include/G4ParticleHPContAngularPar.hh src/G4ParticleHPContAngularPar.cc src/G4ParticleHPContEnergyAngular.cc Please note that it includes a header file as well - so there will be O(20) files re-compiled in particle_hp. All three are from particle_hp ... Cheers, --Sven
I downloaded your zip file, but for some reason it was highly garbled and I could not open it. Please try sending it to my cern e-mail account, dennis.herbert.wright@cern.ch. I hope that works.
(In reply to dennis.herbert.wright from comment #26) > I downloaded your zip file, but for some reason it was highly garbled and I > could not open it. Please try sending it to my cern e-mail account, > dennis.herbert.wright@cern.ch. I hope that works. Hi Dennis, it's not a zip file ... could you please try tar ztvf G4ParticleHPCont-files.tgz on the downloaded file? That should list the content. To extract the content in the local directory do: tar zxvf G4ParticleHPCont-files.tgz If that does not work I'll send you the files directly per e-mail ... Cheers, --Sven
It passed test4 for me and in Geant4 testing, but now it's failing test11 in Geant4 testing. I'll build test11 and see what's going on.
Created attachment 771 [details] tar file of test11 code
The code passed test04 but failed test11 in Geant4 testing. I have attached a tar file of test11 so you can build it. I did some checking and it looks like the trouble happens in G4ParticleHPContAngularPar::BuildByInterpolation at the loop beginning after the comment // add the ones in set2 but not in set1 The code crashes in the loop, but does not seem to have run beyond the number elements in the vector. Partial stack trace: ### CAUGHT SIGNAL: 11 ### address: 0x0, signal = SIGSEGV, value = 11, description = segmentation violation. Address not mapped to object. Backtrace: [PID=40352, TID=0][ 0/22]> 0 ??? 0x0000000000000003 0x0 + 3 [PID=40352, TID=0][ 1/22]> 1 libG4processes.dylib 0x0000000107119761 _ZN29G4ParticleHPContEnergyAngular6SampleEddd + 3409 [PID=40352, TID=0][ 2/22]> 2 libG4processes.dylib 0x00000001070f974d _ZN19G4ParticleHPProduct6SampleEdi + 669 [PID=40352, TID=0][ 3/22]> 3 libG4processes.dylib 0x00000001071266f7 _ZN28G4ParticleHPEnAngCorrelation6SampleEd + 2327 [PID=40352, TID=0][ 4/22]> 4 libG4processes.dylib 0x000000010713fabf _ZN27G4ParticleHPInelasticCompFS14CompositeApplyERK15G4HadProjectileP20G4ParticleDefinition + 1839 [PID=40352, TID=0][ 5/22]> 5 libG4processes.dylib 0x0000000107159c74 _ZN24G4ParticleHPNInelasticFS13ApplyYourselfERK15G4HadProjectile + 36 [PID=40352, TID=0][ 6/22]> 6 libG4processes.dylib 0x000000010710cee6 _ZN23G4ParticleHPChannelList13ApplyYourselfEPK9G4ElementRK15G4HadProjectile + 1574 [PID=40352, TID=0][ 7/22]> 7 libG4processes.dylib 0x00000001071381a8 _ZN21G4ParticleHPInelastic13ApplyYourselfERK15G4HadProjectileR9G4Nucleus + 7
Hi Dennis, thanks for the test11 code. It appears that the crash is for 10 MeV neutrons on Pm148m1 which uses Inelastic/F04/61_148m1_Promethium.z The problematic code is now the continuous spectrum (a part I did not touch at all ...). And here the problem is that PrepareTableInterpolation() and BuildByInterpolation() assume different things about the number of energies to calculate. The former basically adds all energies from 2 sets (minus first and last) while the latter assumes the 2nd set defines the total number of energies to populate. This causes a segv when the 2nd set has less energies than the combined set in the prepare step. There seem to be cases where the prepare step is not run - so one needs to check for that. I will try a few options and report back here - but clearly this entire code is extremely fragile ... Cheers, --Sven
Hi Dennis, I found quite a number of issues with the continuous code but even with that fixed the test11 kept crashing in Pm148m1 but at a different file (G4NDL4.6/Inelastic/F01/61_148m1_Promethium.z) Ultimately the reason was that discrete lines at 1.481000e+05 eV appear twice in that file. The code assumes that discrete lines are unique for the interpolation - leaving a hole with a NULL pointer. I don't know if the duplicate energy is simply a bug in F01/61_148m1_Promethium.z that can be fixed or if the code should be capable of dealing with identical discrete lines. And for the latter case the question is how to interpolate then? What if one set has 2 and the other 1 or 3 lines with same energy? Judging from the sequence of discrete energies in that region steps are accurate to 0.1 keV and range between 3 and 5.5 keV but the gap to the next lower energy for this duplicate line is 6.8 keV - so the 2nd one might in fact have been closer to 145 keV ... How to proceed? Thanks, --Sven
Created attachment 772 [details] patch for particle_hp in Geant4.11.0.1 to fix the FS treatment Hi Dennis, here's a new patch (will attach gziped tar archive with the 4 files next). This touches now 4 files: include/G4ParticleHPContAngularPar.hh include/G4ParticleHPContEnergyAngular.hh src/G4ParticleHPContAngularPar.cc src/G4ParticleHPContEnergyAngular.cc in particle_hp and is relative to the unchanged 11.0.1 source. New are the following changes: 1) PrepareTableInterpolation() runs without argument since the continuous spectra calculation in this step was wrong and has to be re-done in BuildByInterpolation anyway 2) Because some data files (by mistake? by rounding?) contain discrete energy lines with identical energies the insertion into the unique map of energies for lookup needs to tweak the lookup table in PrepareTableInterpolation(). This does not re-assign energies - the tweaked entries are just used for lookup. So the data files still should be fixed ... 3) BuildByInterpolation(): the calculation of the continuous spectrum has been fixed. This also required some changes on the prior changed discrete line code. For example the min and max range of the interpolated set and also the total number of energies in the interpolated set are known at later points only. This runs in test11 for me. Could you please try too? Thanks & Cheers, --Sven
Created attachment 773 [details] gziped tar archive with the 4 modified files (two .hh, two .cc) for G4.11.0.1 gziped tar archive corresponding to the 4 files as in the patch (id 772)
Created attachment 774 [details] patch for particle_hp in Geant4.11.0.1 to fix the FS treatment small improvement for the duplicate lines - using the tweaked energy at all places for lookup to build interpolated table.
Created attachment 775 [details] gziped tar archive with the 4 modified files (two .hh, two .cc) for G4.11.0.1 gziped tar archive corresponding to id 774
The fix finally passed testing and has been merged with the Geant4 code. There is more cleanup required, but I will take care of that. I have asked our G4NDL expert about the duplicate energy lines and am waiting to hear from him. Thanks for your work on solving this.
