Problem 1480

Summary:	rdecay02 produces negative weight in the root tree produced by the macro, u2380-b.mac
Product:	Examples/Extended	Reporter:	John Watts <john.w.watts>
Component:	radioactivedecay/exrdm	Assignee:	dennis.herbert.wright
Status:	RESOLVED FIXED
Severity:	normal	CC:	dennis.herbert.wright
Priority:	P5
Version:	9.6
Hardware:	All
OS:	Linux

Description John Watts 2013-06-14 21:24:37 CEST

When I run the example rdecay02 with the macro, u238-b.mac and look at the weights in the root file created some of them are negative. Below is a scan of the first 100 events in the tree, tT3. Trees, tT1 and tRDecayProducts, also have some negative weights. This was motivated by trying to look at the decay of Americium 241 in variance reduction mode where I also saw negative weights. Is this a bug?

root [0] TBrowser b; root [1] (class TFile*)0x2247fb0 File name :
u238-b.root

************************************************
*    Row   *    Energy *      Time *    Weight *
************************************************
*        0 * 2.7853682 * 10000.702 * 2.749e-26 *
*        1 * 2.5189876 * 10000.702 * 2.749e-26 *
*        2 * 2.9832148 * 10000.702 * 7.036e-32 *
*        3 * 1.5333689 * 10000.702 * 7.036e-32 *
*        4 * 2.8262794 * 10000.702 * 2.746e-23 *
*        5 * 2.4780764 * 10000.702 * 2.746e-23 *
*        6 * 2.7881465 * 10000.702 * 3.521e-28 *
*        7 * 1.7284371 * 10000.702 * 3.521e-28 *
*        8 * 3.1026215 * 10000.702 * 1.676e-24 *
*        9 * 1.4139622 * 10000.702 * 1.676e-24 *
*       10 * 2.6442267 * 10000.702 * 2.174e-29 *
*       11 * 2.0375437 * 10000.702 * 2.174e-29 *
*       12 * 2.7544729 * 10000.702 * 2.285e-29 *
*       13 * 1.9272977 * 10000.702 * 2.285e-29 *
*       14 * 2.9806013 * 10000.702 * 1.632e-25 *
*       15 * 1.6629138 * 10000.702 * 1.632e-25 *
*       16 * 3.1151557 * 10000.702 * 7.835e-34 *
*       17 * 0.6045311 * 10000.702 * 7.835e-34 *
*       18 * 2.9392805 * 10000.702 * 1.643e-34 *
*       19 * 0.7804062 * 10000.702 * 1.643e-34 *
*       20 * 2.5905015 * 10000.702 * 4.289e-30 *
*       21 * 3.0821621 * 10000.702 * 4.289e-30 *
*       22 * 1.2294195 * 10000.702 * 4.289e-30 *
*       23 * 2.4770472 * 10000.702 * 4.118e-23 *
*       24 * 2.8341496 * 10000.702 * 4.118e-23 *
*       26 * 2.5544788 * 10000.702 * 2.060e-19 *
*       27 * 2.8283758 * 10000.702 * 2.060e-19 *
*       28 * 2.3039419 * 10000.702 * 2.060e-19 *
*       29 * 2.6172049 * 10000.702 * 3.491e-29 *
*       30 * 2.9899339 * 10000.702 * 3.491e-29 *
*       31 * 0.9238945 * 10000.702 * 3.491e-29 *
*       32 * 2.5979037 * 10000.702 * 1.031e-24 *
*       33 * 2.2287290 * 10000.702 * 1.031e-24 *
*       34 * 2.7688245 * 10000.702 * 3.245e-19 *
*       35 * 2.0155024 * 10000.702 * 3.245e-19 *
*       36 * 2.5915708 * 10000.702 * 6.009e-24 *
*       37 * 1.7794442 * 10000.702 * 6.009e-24 *
*       38 * 2.7426013 * 10000.702 * 2.943e-19 *
*       39 * 2.0327587 * 10000.702 * 2.943e-19 *
*       40 * 2.6744175 * 10000.702 * 4.677e-17 *
*       41 * 1.5276826 * 10000.702 * 4.677e-17 *
*       42 *         0 *         0 *         0 *
*       43 * 2.8376867 * 1000.3080 * 2.749e-26 *
*       44 * 2.4666688 * 1000.3080 * 2.749e-26 *
*       45 * 2.7498765 * 1000.3080 * 7.036e-32 *
*       46 * 1.7667073 * 1000.3080 * 7.036e-32 *
*       47 * 2.9209964 * 1000.3080 * 2.746e-23 *
*       48 * 2.3833591 * 1000.3080 * 2.746e-23 *
*       49 * 2.8206787 * 1000.3080 * 2.886e-23 *
*       51 * 3.0944950 * 1000.3080 * 1.676e-24 *
*       52 * 1.4220887 * 1000.3080 * 1.676e-24 *
*       53 * 2.8102359 * 1000.3080 * 3.266e-32 *
*       54 * 1.8332790 * 1000.3080 * 3.266e-32 *
*       55 * 3.0982220 * 1000.3080 * 4.568e-33 *
*       56 * 1.5835487 * 1000.3080 * 4.568e-33 *
*       57 * 2.8220512 * 1000.3080 * 2.174e-29 *
*       58 * 1.8597195 * 1000.3080 * 2.174e-29 *
*       59 * 2.8450727 * 1000.3080 * 2.285e-29 *
*       60 * 1.8366979 * 1000.3080 * 2.285e-29 *
*       61 * 2.7993373 * 1000.3080 * 7.826e-31 *
*       62 * 0.9203493 * 1000.3080 * 7.826e-31 *
*       63 * 2.9851815 * 1000.3080 * 8.226e-31 *
*       64 * 0.7345051 * 1000.3080 * 8.226e-31 *
*       65 * 2.6072919 * 1000.3080 * 4.289e-30 *
*       66 * 2.7685637 * 1000.3080 * 4.289e-30 *
*       67 * 1.5262278 * 1000.3080 * 4.289e-30 *
*       68 * 2.5885696 * 1000.3080 * 2.471e-29 *
*       69 * 2.9763240 * 1000.3080 * 2.471e-29 *
*       70 * 1.0447734 * 1000.3080 * 2.471e-29 *
*       71 * 2.6180970 * 1000.3080 * 2.143e-23 *
*       72 * 2.8751776 * 1000.3080 * 2.143e-23 *
*       73 * 1.4088084 * 1000.3080 * 2.143e-23 *
*       74 * 2.7525107 * 1000.3080 * 3.491e-29 *
*       76 * 0.4739412 * 1000.3080 * 3.491e-29 *
*       77 * 2.6603517 * 1000.3080 * 4.944e-23 *
*       78 * 2.7694582 * 1000.3080 * 4.944e-23 *
*       79 * 1.2563291 * 1000.3080 * 4.944e-23 *
*       80 * 2.8865802 * 1000.3080 * 2.060e-19 *
*       81 * 2.6028671 * 1000.3080 * 2.060e-19 *
*       82 * 2.6710798 * 1000.3080 * 2.233e-23 *
*       83 * 1.6728264 * 1000.3080 * 2.233e-23 *
*       84 * 3.1670901 * 1000.3080 * 8.673e-27 *
*       85 * 0.6613824 * 1000.3080 * 8.673e-27 *
*       86 * 3.1573665 * 1000.3080 * 2.886e-26 *
*       87 * 0.9519953 * 1000.3080 * 2.886e-26 *
*       88 * 2.7728865 * 1000.3080 * 4.677e-17 *
*       89 * 1.4292136 * 1000.3080 * 4.677e-17 *
*       90 *         0 *         0 *         0 *
*       91 * 2.8279478 * 10000.651 * -5.92e-26 *
*       92 *  2.476408 * 10000.651 * -5.92e-26 *
*       93 * 2.7955148 * 10000.651 * -1.51e-31 *
*       94 * 1.7210689 * 10000.651 * -1.51e-31 *
*       95 * 2.8671462 * 10000.651 * -5.91e-23 *
*       96 * 2.4372093 * 10000.651 * -5.91e-23 *
*       97 * 2.7106945 * 10000.651 * -6.21e-23 *
*       98 * 2.5936610 * 10000.651 * -6.21e-23 *
*       99 * 2.9626374 * 10000.651 * -2.96e-19 *

Comment 1 dennis.herbert.wright 2017-12-12 03:05:10 CET

Still working on this.  Problem narrowed down to generalized application of Bateman equations.

Comment 2 dennis.herbert.wright 2018-04-06 18:10:37 CEST

This is now fixed.  The main problem was the accumulation of cancellation errors which are especially apparent in long decay chains.  In the weight calculation there are sums of products of factors like tau1/(tau1 - tau2) which often alternate sign.  Each term in these sums is then multiplied by factors of (1 - exp(...) ).  The final result, which should always be positive, is sometimes a small negative number due to loss of precision.

The main part of the solution is to use the C++11 function std::expm1 which is exp - 1.  This function avoids the large cancellation errors which occur when using 1 - std::exp for arguments near zero.  The result is that the magnitude of negative weights is reduced by around ten orders of magnitude, for example 1e-20 down to 1e-30.  This is an indication that the decay chain math is correct and is attempting to produce a positive number, but that some small cancellation errors remain.  In this case the code now sets negative values to zero since the error of doing so is quite small.

The fix has been added to the Geant4 repository and will be available in the Geant4 10.5 beta release in June.