Observable validator

The validator is designed to take in a userfile (such as those made in Getting started) and a which_piece or generic process (see section Available processes and which_piece) and perform various checks of the measurement function to see if the observables it implements are physical.

The measurement function acts as a perfect observer of the system. This can have some draw backs: for example two indistinguishable particles can be distinguished by the measurement function and can lead to results that are physically impossible to obtain in real life.

The design of the validator is such that installing McMule in full on the computer running the tests is not required. Instead, pymule provides the mocks of the functions made available by mcmule.mod (cf. Section libmcmule user guide). Should the user function require more McMule internals, a full installation is required and needs to be manually initialised.

The checks results are sorted into three different categories: passed, warnings, and failed. Failing any test means that the measurement function is most certainly unphysical and should be fixed (if you believe a test is failed incorrectly, please open an issue with the full logs). Warnings are merely recommendations on how the user function might be modified to make it more universally applicable.

Tests that may result in an error include:
  • Soft photon limit: taking a soft photon’s energy to zero should have the same result as the photon not existing

  • Swap symmetry: swapping two indistinguishable particles should not change the observables

  • Soft-limit swap-symmetry check: a combination of swapping indistinguishable photons, one of which is soft

  • Out-of-bounds check: all events fit inside the histogram bounds

Tests that may result in a warning include:
  • Rotation check: it should not matter if the system decays in the Z axis or any other direction

  • Lorentz boost check: ideally a frame should be chosen and boosted into inside the measurement function

  • Flip check: symmetrically flipping the system in any axis should have no impact on the observables

Further usage and interpretation of each test can be found below.

Using the validator

When using the validator, we can either compile the user file ourselves

$ gfortran -I/path/to/mcmule.mod --shared -o user.so -fPIC user.f95 -fdefault-real-8
$ pymule validator <arguments> user.so

or use pymule to compile

$ pymule validator <arguments> user.f95

A full list of arguments is given in below.

A measurement function is always tested against a which_piece or collection of pieces. To test eg. m2enng0 we type

$ pymule validator -l INFO -w m2enng0 user.so

where we specify the logging to be INFO so we see all the tests.

The output here can be split into 3 sections

  1. Loading the required libraries

At the logging level of INFO we see each of the shared libraries which must be loaded for McMule to work correctly, including the userfile.

INFO:pymule.validator.importer:Loading shared libraries
INFO:pymule.validator.importer:Loading library: libcollier.so
INFO:pymule.validator.importer:Loading library: libolcommon.so
INFO:pymule.validator.importer:Loading library: libtrred.so
INFO:pymule.validator.importer:Loading library: liboneloop.so
INFO:pymule.validator.importer:Loading library: libcuttools.so
INFO:pymule.validator.importer:Loading library: libopenloops.so
INFO:pymule.validator.importer:Loading library: libqedc19.so
INFO:pymule.validator.importer:Loading library: libmcmule.so
INFO:pymule.validator.importer:Loading library: ./user.so
  1. Loading the measurement function and variables

Next the functions and variables required to test the userfile are found and loaded, including running inituser() if it can be found (in this case simply printing This is an empty user file).

INFO:pymule.validator.importer:Function 'QUANT' found with name: __user_MOD_quant
INFO:pymule.validator.importer:Function 'inituser' found with name: __user_MOD_inituser
INFO:pymule.validator.importer:Variable 'filenamesuffixLen' found with name: __user_MOD_filenamesuffixlen
This is a empty user file
INFO:pymule.validator.importer:Variable 'nrq' found with name: __user_MOD_nq
INFO:pymule.validator.importer:Variable 'nrbins' found with name: __user_MOD_nbins
  1. Testing of the measurement function

Then the actually testing of the measurement function can begin. The validator will automatically find and run the appropriate tests for a each which_piece, in this case:

INFO:pymule.validator.validator:Testing whichpiece: m2enng0 with 6 tests
INFO:pymule.validator.validator:Test 1: SwapSymmetryValidator(2, 3)... passed
INFO:pymule.validator.validator:Test 2: RotationValidator(0.08644597231507764, 0, 0)... passed
INFO:pymule.validator.validator:Test 3: RotationValidator(5.104560253717271, 3.868168781137662, 5.213091066969582)... passed
INFO:pymule.validator.validator:Test 4: FlipValidator(x)... passed
INFO:pymule.validator.validator:Test 5: FlipValidator(y)... passed
INFO:pymule.validator.validator:Test 6: FlipValidator(z)... passed
INFO:pymule.validator.validator:Overall pass rate: 6 / 6, 100.0%

For this userfile all the appropriate tests ran correctly and passed. Different which_pieces may run different tests with different results.

If a Non-fatal test fails then, instead of reporting at the INFO level, it reports at WARNING level.

WARNING:pymule.validator.validator:Test 4: FlipValidator(x)... failed,
...
WARNING:pymule.validator.validator:Overall pass rate: 5 / 6, 83.33%

If a Fatal test fails it reports at the ERROR level

ERROR:pymule.validator.validator:Test 1: SwapSymmetryValidator(2, 3)... failed,
...
ERROR:pymule.validator.validator:Overall pass rate: 5 / 6 with 1 fatal failure! 83.33%

Similarly a generic piece can have all the which_pieces under it be checked in one go.

For example testing all orders of the generic piece ee2ee is as simple as typing

$ pymule validate -g ee2ee -l INFO user.so
...
INFO:pymule.validator.validator:Testing at order LO:
INFO:pymule.validator.validator:Testing whichpiece: ee2ee0 with 8 tests
...
INFO:pymule.validator.validator:Testing at order NLO:
INFO:pymule.validator.validator:Testing whichpiece: ee2eeF with 8 tests
INFO:pymule.validator.validator:Testing whichpiece: ee2eeA with 8 tests
...
INFO:pymule.validator.validator:Testing at order NNLO:
INFO:pymule.validator.validator:Testing whichpiece: ee2eeFF with 8 tests
...

If we only want to test a single order (such as NLO) we could write

$ pymule validate -g ee2ee -o NLO -l INFO user.so
...
INFO:pymule.validator.validator:Testing at order NLO:
INFO:pymule.validator.validator:Testing whichpiece: ee2eeF with 8 tests
INFO:pymule.validator.validator:Testing whichpiece: ee2eeA with 8 tests
...

Logging Levels

Table 2 Various logging levels for the validator

Integer

Name

Effect

-1

CRITICAL

No output of anykind from the validator, mainly used for internal code testing

0

ERROR

Only Fatal failures output, suppresses graphing

1

WARNING

Fatal and Non-fatal failures are output. Allows graphing

2

INFO

All tests being run are output, even if they pass. Allows graphing

3

DEBUG

All output possible, including forcing graphs to be generated if applicable

Measurement function tests in detail

Each test checks some aspect of the measurement function. This section will go through example cases of where each test may fail and how to prevent that failure.

Out-of-bounds validator

The out-of-bounds validator generates a fixed number of Monte Carlo events with uniform weight and calculates the binning of the events that pass the cuts. If any event is ends up outside the histogram bounds, this test fails. As this test is fairly slow, it needs to be activated by specifying the number of points to be used. After all pieces have been considered, the validator will print a list of recommended selection cuts, eg.

Observable    Selection       min_val    act. min     act. max    max_val
------------  -----------  ----------  ----------  -----------  ---------
p+             51: 80        0.45        0.450094     0.522873    1
p-             51: 80        0.45        0.450053     0.522873    1
thav            1:600        1           1.00002      2.14154     2.14159
dphi            1:600       -0           0            0.149923    0.15
xi              1:394       -0           0            0.164099    0.25
mxx           252:185      299.167     387.413      422.634     700
cth+          126:471       -0.646667   -0.605093     0.593247    0.64
cth-          156:446       -0.646667   -0.540288     0.540256    0.64
p+c           inf:-inf       0.448     inf         -inf           1
p-c           inf:-inf       0.448     inf         -inf           1

The columns of this table indicate

  • observable: the name of the observable

  • selection: the recommended selection that guarantees the histogram bounds are as tight as possible without decreasing the bin width

  • min_val and max_val: the user-provided min_val and max_val

  • actual minimum and maximum: the actual bounds of the histogram as found using Monte Carlo

As McMule requires a constant number of histogram bins, one might artificially broaden the min_valmax_val range. This makes the out-of-bounds validator less powerful. It is hence recommended to include a selection hint as part of the observable specification. For example

integer, parameter :: nrq = 10
integer, parameter :: nrbins = 600
real(kind=prec), parameter :: &
   min_val(nrq) = (/ 0.40, 0.40, 1.0, 0.0, 0.00, 200., -1.0, -1.0, 0.40, 0.40 /)
real(kind=prec), parameter :: &
   max_val(nrq) = (/ 1.00, 1.00, pi-1, 0.15, 0.25, 700., 1.0, 1.0, 1.00, 1.00 /)

integer :: selection_hint(2*nrq) = (/ &
  51, 51, 1, 1, 1, 120, 107, 107, 49, 49, &
  -1, -1,-1,-1,-1,  -1, 492, 492, -1, -1 /)

The selection_hint is an array with 2N entries listing first the indices of the first bins one wishes to use and then the indices of the last bins where -1 defaults to the last bin. In the above example, the first histogram is defined to for the range [0.40, 1.00]. However, the cuts will mean that the first 50 bins [0.40, 0.45] are empty. Providing the selection_hint ensures pymule knows about this.

Swap symmetry validator

The swap symmetry validator checks that two indistinguishable particles can be swapped without the observables changing. In this example we look at the ee2ee0 which_piece which has 2 distinguishable electrons in and 2 indistinguishable electrons out. Taking \(S\) to be our measurement function we want to check

\[S(e_1, e_2) = S(e_2, e_1),\]

i.e. that the measurement function is the same after swapping the indistinguishable electrons.

An example measurement function which fails this test is given below, along with the validator being called on the code,

Listing 11 An example measurement function which is not swap symmetry safe
58! This should fail swap symmetry !
59FUNCTION QUANT(q1,q2,q3,q4,q5,q6,q7)
60    ! q1 -> e_in_1, q2 -> e_in_2 | q3 -> e_out_1, q4 -> e_out_2
61    real (kind=prec), intent(in) :: q1(4),q2(4),q3(4),q4(4), q5(4),q6(4),q7(4)
62    real (kind=prec) :: q1_rf(4), q2_rf(4), q3_rf(4), q4_rf(4)
63    real (kind=prec) :: quant(nr_q)
64    !! ==== keep the line below in any case ==== !!
65    call fix_mu
66
67    q1_rf = BOOST_RF(q1, q1)
68    q2_rf = BOOST_rf(q1, q2)
69    q3_rf = BOOST_rf(q1, q3)
70    q4_rf = BOOST_rf(q1, q4)
71
72    names(1) = 'Ee1'
73    quant(1) = q3_rf(4)
74
75    names(2) = 'Ee2'
76    quant(2) = q4_rf(4)
77
78END FUNCTION QUANT
$ pymule validate -w ee2ee0 -l INFO FailsSwapSymmetry.so
...
INFO:pymule.validator.validator:Testing whichpiece: ee2ee0 with 8 tests
ERROR:pymule.validator.validator:Test 1: SwapSymmetryValidator(2, 3)... failed fatally,
INFO:pymule.validator.validator:Test 2: LorentzBoostValidator(0)... passed
INFO:pymule.validator.validator:Test 3: LorentzBoostValidator(1)... passed
INFO:pymule.validator.validator:Test 4: RotationValidator(4.224052347507555, 0, 0)... passed
INFO:pymule.validator.validator:Test 5: RotationValidator(1.7262535645602621, 3.5353046863936015, 5.066130454491858)... passed
INFO:pymule.validator.validator:Test 6: FlipValidator(x)... passed
INFO:pymule.validator.validator:Test 7: FlipValidator(y)... passed
INFO:pymule.validator.validator:Test 8: FlipValidator(z)... passed
ERROR:pymule.validator.validator:Overall pass rate: 7 / 8 with 1 fatal failure! 87.5%

The measurement function correctly chooses a frame (the restframe of q1) so passes the LorentzBoostValidators (Tests 2 and 3). However it directly looks at the energies of the two electrons and assigns them to electron 1 and electron 2. This is something which could not be done in the real world since they are indistinguishable. The proper way of doing this would be looking at the higher energy electron and lower energy electron:

Listing 12 An example measurement function which is swap symmetry safe
58! This should fail swap symmetry !
59FUNCTION QUANT(q1,q2,q3,q4,q5,q6,q7)
60    ! q1 -> e_in_1, q2 -> e_in_2 | q3 -> e_out_1, q4 -> e_out_2
61    real (kind=prec), intent(in) :: q1(4),q2(4),q3(4),q4(4), q5(4),q6(4),q7(4)
62    real (kind=prec) :: q1_rf(4), q2_rf(4), q3_rf(4), q4_rf(4), qHigh(4), qLow(4)
63    real (kind=prec) :: quant(nr_q)
64    !! ==== keep the line below in any case ==== !!
65    call fix_mu
66
67    q1_rf = BOOST_RF(q1, q1)
68    q2_rf = BOOST_RF(q1, q2)
69    q3_rf = BOOST_RF(q1, q3)
70    q4_rf = BOOST_RF(q1, q4)
71
72    if (q3_rf(4) > q4_rf(4)) then
73        qHigh = q3_rf; qLow = q4_rf
74    else
75        qHigh = q4_rf; qLow = q3_rf
76    endif
77
78    names(1) = 'High'
79    quant(1) = qHigh(4)
80
81    names(2) = 'Low'
82    quant(2) = qLow(4)
83
84END FUNCTION QUANT
$ pymule validate -w ee2ee0 -l INFO PassesSwapSymmetry.so
...
INFO:pymule.validator.validator:Testing whichpiece: ee2ee0 with 8 tests
INFO:pymule.validator.validator:Test 1: SwapSymmetryValidator(2, 3)... passed
INFO:pymule.validator.validator:Test 2: LorentzBoostValidator(0)... passed
INFO:pymule.validator.validator:Test 3: LorentzBoostValidator(1)... passed
INFO:pymule.validator.validator:Test 4: RotationValidator(2.052666305609015, 0, 0)... passed
INFO:pymule.validator.validator:Test 5: RotationValidator(1.1268264348854389, 5.94661892034897, 1.2312470761319962)... passed
INFO:pymule.validator.validator:Test 6: FlipValidator(x)... passed
INFO:pymule.validator.validator:Test 7: FlipValidator(y)... passed
INFO:pymule.validator.validator:Test 8: FlipValidator(z)... passed
INFO:pymule.validator.validator:Overall pass rate: 8 / 8, 100.0%

Looking at the higher and lower energy electrons (instead of electron 1 and electron 2) means that swapping them has no effect on the observable, so the swap symmetry validator passes.

Lorentz boost validator

In the previous section the measurement function very deliberately chose to boost into the rest frame of the first incoming electron. This means that no matter how the system was boosted before this point, a frame is chosen and boosted into. This overwrites any previous boosts and gives a unique frame.

Note

LorentzBoostValidator(0) refers to testing the system by boosting into the rest frame of the zeroth particle (i.e. q1 in the measurement function)

Listing 13 An example measurement function which is not lorentz boost safe
58! This should fail Lorentz Boost !
59FUNCTION QUANT(q1,q2,q3,q4,q5,q6,q7)
60    ! q1 -> e_in_1, q2 -> e_in_2 | q3 -> e_out_1, q4 -> e_out_2
61    real (kind=prec), intent(in) :: q1(4),q2(4),q3(4),q4(4), q5(4),q6(4),q7(4)
62    real (kind=prec) :: qHigh(4), qLow(4)
63    real (kind=prec) :: quant(nr_q)
64    !! ==== keep the line below in any case ==== !!
65    call fix_mu
66
67    if (q3(4) > q4(4)) then
68        qHigh = q3; qLow = q4
69    else
70        qHigh = q4; qLow = q3
71    endif 
72
73    names(1) = 'High'
74    quant(1) = qHigh(4)
75
76    names(2) = 'Low'
77    quant(2) = qLow(4)
78
79END FUNCTION QUANT
$ pymule validate -w ee2ee0 -l INFO FailsLorentzBoost.so
...
INFO:pymule.validator.validator:Testing whichpiece: ee2ee0 with 8 tests
INFO:pymule.validator.validator:Test 1: SwapSymmetryValidator(2, 3)... passed
WARNING:pymule.validator.validator:Test 2: LorentzBoostValidator(0)... failed, assuming the center of mass frame is allowed but not advisable
WARNING:pymule.validator.validator:Test 3: LorentzBoostValidator(1)... failed, assuming the center of mass frame is allowed but not advisable
INFO:pymule.validator.validator:Test 4: RotationValidator(4.354392616889299, 0, 0)... passed
INFO:pymule.validator.validator:Test 5: RotationValidator(4.167936819187568, 5.129621429511378, 1.185060876837443)... passed
INFO:pymule.validator.validator:Test 6: FlipValidator(x)... passed
INFO:pymule.validator.validator:Test 7: FlipValidator(y)... passed
INFO:pymule.validator.validator:Test 8: FlipValidator(z)... passed
WARNING:pymule.validator.validator:Overall pass rate: 6 / 8, 75.0%

Choosing a frame to boost into (like in the swap symmetry example) means that the test no longer fails.

Listing 14 An example measurement function which is Lorentz boost safe
58        ! This should fail swap symmetry !
59        FUNCTION QUANT(q1,q2,q3,q4,q5,q6,q7)
60            ! q1 -> e_in_1, q2 -> e_in_2 | q3 -> e_out_1, q4 -> e_out_2
61            real (kind=prec), intent(in) :: q1(4),q2(4),q3(4),q4(4), q5(4),q6(4),q7(4)
62            real (kind=prec) :: q1_rf(4), q2_rf(4), q3_rf(4), q4_rf(4), qHigh(4), qLow(4)
63            real (kind=prec) :: quant(nr_q)
64            !! ==== keep the line below in any case ==== !!
65            call fix_mu
66        
67            q1_rf = BOOST_RF(q1, q1)
68            q2_rf = BOOST_RF(q1, q2)
69            q3_rf = BOOST_RF(q1, q3)
70            q4_rf = BOOST_RF(q1, q4)
71
72            if (q3_rf(4) > q4_rf(4)) then
73                qHigh = q3_rf; qLow = q4_rf
74            else
75                qHigh = q4_rf; qLow = q3_rf
76            endif
77            
78            names(1) = 'High'
79            quant(1) = qHigh(4)
80        
81            names(2) = 'Low'
82            quant(2) = qLow(4)
83        
84        END FUNCTION QUANT
$ pymule validate -w ee2ee0 -l INFO FailsLorentzBoost.so
...
INFO:pymule.validator.validator:Testing whichpiece: ee2ee0 with 8 tests
INFO:pymule.validator.validator:Test 1: SwapSymmetryValidator(2, 3)... passed
INFO:pymule.validator.validator:Test 2: LorentzBoostValidator(0)... passed
INFO:pymule.validator.validator:Test 3: LorentzBoostValidator(1)... passed
INFO:pymule.validator.validator:Test 4: RotationValidator(0.20776147014583715, 0, 0)... passed
INFO:pymule.validator.validator:Test 5: RotationValidator(0.21958279258314367, 1.470413694055576, 1.8065232010826624)... passed
INFO:pymule.validator.validator:Test 6: FlipValidator(x)... passed
INFO:pymule.validator.validator:Test 7: FlipValidator(y)... passed
INFO:pymule.validator.validator:Test 8: FlipValidator(z)... passed
INFO:pymule.validator.validator:Overall pass rate: 8 / 8, 100.0%

Soft limit validator and soft-swap validator

A soft photon is one whose energy can be taken to zero. For example in the which_piece m2enngR (\(\mu^- \to e^- \bar{\nu_e} \nu_{\mu} \gamma \gamma\)) there is one soft photon (represented by the R) and one always real photon (represented by g).

Mathematically, in the soft limit validator we are testing,

\[\lim_{\gamma \to 0} S(q, \gamma) = S(q),\]

and in the case of the soft-swap validator we are testing,

\[\lim_{\gamma_1 \to 0} S(q, \gamma_1, \gamma_2) = \lim_{\gamma_1 \to 0} S(q, \gamma_2, \gamma_1) = S(q, \gamma_2)\]

where \(S\) is our measurement function.

An example of failing this test is looking directly at number of photons (or experimentally number of photons detected)

If we say a photon exists if it has any energy at all then we could write

Listing 15 An example measurement function which is not soft limit safe
58! This should fail soft limit !
59FUNCTION QUANT(q1,q2,q3,q4,q5,q6,q7)
60    ! m(q1) -> e(q2) n(q3) n(q4) g(q5) R/0(q6)
61    real (kind=prec), intent(in) :: q1(4),q2(4),q3(4),q4(4), q5(4),q6(4),q7(4)
62    real (kind=prec) :: quant(nr_q)
63    integer :: numphotons
64    !! ==== keep the line below in any case ==== !!
65    call fix_mu
66
67    numphotons = 0
68
69    if (q5(4) > 0) then
70        numphotons = numphotons + 1
71    endif
72    if (q6(4) > 0) then
73        numphotons = numphotons + 1
74    endif
75
76    names(1) = 'numpho'
77    quant(1) = numphotons
78
79    names(2) = 'const'
80    quant(2) = 5
81
82END FUNCTION QUANT
$ pymule validate -w m2enngR -l INFO FailsSoftLimit.so
...
INFO:pymule.validator.validator:Testing whichpiece: m2enngR with 9 tests
# Graph appears!
ERROR:pymule.validator.validator:Test 1: SoftLimitValidator(0)... failed fatally,
INFO:pymule.validator.validator:Test 2: SwapSymmetryValidator(2, 3)... passed
INFO:pymule.validator.validator:Test 3: SwapSymmetryValidator(4, 5)... passed
# Another graph appears!
ERROR:pymule.validator.validator:Test 4: SoftSwapValidator(0, 4, 5)... failed fatally,
...
ERROR:pymule.validator.validator:Overall pass rate: 7 / 9 with 2 fatal failures! 77.78%

The failure occurs because at the photon energy proportion \(\xi = 0\), the soft photon has \(0\) energy so only one photon exists so numphotons = 1. At any higher energy this photon is counted so numphotons = 2. This discrepancy means the limit above doesn’t hold, so the test fails. Notice that there is some additional functionality when a soft limit test (or soft-swap test) fails. The graph that is generated shows which of the observables caused the test to fail. In this case numpho is shown but not const. We define the \(Z = \|\frac{S(\xi) - S(0)}{tol}\|\) where \(tol\) is a relative tolerance to the value \(S(0)\) (usually \(S(0)\times10^{-3}\)) plus an absolute tolerance which is used for \(S(0) = 0\) (usually also \(10^{-3}\)).

../_images/SoftLimitGraph.svg

Figure 6 There is a line at \(Z=1\) indicating that a value must (consistently) fall below that line in order to be considered smooth and going to the expected limit. Note that numpho is shown but not const.

A measurement function more inline with physical reality is one where the photon is not counted if it is below some threshold, e.g. the detection limit of the detector.

Listing 16 An example measurement function which is soft limit safe
58! This should pass soft limit !
59FUNCTION QUANT(q1,q2,q3,q4,q5,q6,q7)
60    ! m(q1) -> e(q2) n(q3) n(q4) g(q5) R/0(q6)
61    real (kind=prec), intent(in) :: q1(4),q2(4),q3(4),q4(4), q5(4),q6(4),q7(4)
62    real (kind=prec) :: quant(nr_q)
63    integer :: numphotons
64    !! ==== keep the line below in any case ==== !!
65    call fix_mu
66
67    numphotons = 0
68
69    if (q5(4) > 1e-3) then ! NOTE: threshold for counting the photon
70        numphotons = numphotons + 1
71    endif
72    if (q6(4) > 1e-3) then ! NOTE: threshold for counting the photon
73        numphotons = numphotons + 1
74    endif
75
76    names(1) = 'numpho'
77    quant(1) = numphotons
78
79    names(2) = 'const'
80    quant(2) = 5
81
82END FUNCTION QUANT
$ pymule validate -w m2enngR -l INFO PassesSoftLimit.so
...
INFO:pymule.validator.validator:Testing whichpiece: m2enngR with 9 tests
INFO:pymule.validator.validator:Test 1: SoftLimitValidator(0)... passed
INFO:pymule.validator.validator:Test 2: SwapSymmetryValidator(2, 3)... passed
INFO:pymule.validator.validator:Test 3: SwapSymmetryValidator(4, 5)... passed
INFO:pymule.validator.validator:Test 4: SoftSwapValidator(0, 4, 5)... passed
...

If we still wish to see the graph (for example if it shouldn’t have passed but did) we can send the DEBUG flag to the logging and force the graph to be displayed.

$ pymule validate -w m2enngR -l DEBUG PassesSoftLimit.so
...
DEBUG:matplotlib ...
...
# Graph appears!
INFO:pymule.validator.validator:Test 1: SoftLimitValidator(0)... passed
...

There is unfortunately quite a bit of matplotlib spam (In the works to be fixed!) but the following graph is produced.

../_images/SoftLimitGraph_Debug.svg

Figure 7 We can see that at around \(\xi = 10^{-5}\) the graph starts to dip down. Numerically numpho has gone to \(1\) but the Richardson extrapolation which is used reach the limit faster for other observables lags slightly. After the value has been below the \(Z = 1\) threshold for three iterations it is said to converge. Note also that the observable const is also in the legend of the graph.

Note

DEBUG will force all observables onto the graph regardless of if they pass or not. In this case all values of const are at \(Z = 0\) which is not displayed on the log-log plot.

In the m2enngR which_piece the energy to be passed around is \(E \approx 100\) MeV, meaning the cut off for the photon counting happens at \(E \xi = 10^{-3}\) MeV, \(\xi \approx 10^{-5}\) as seen in the graph.

Rotation validator and Flip validator

These two tests are bundled together as it is difficult to only fail one and not both. When calling the validator two rotation tests are performed. The first involves only rotating around the Z-axis. This will be passed if, for example, the observable only looks at the momentum perpendicular to the Z-axis. But the test would be failed if the observable looks for example only at the momentum in the X-axis. The second rotation test is a rotation of an arbitrary set of angles, designed to test complete rotational invariance.

The flip tests are designed to check the symmetry of flipping the system in the X, Y, and Z axes (individually).

All these tests are non-fatal, meaning that failing them is not necessarily bad as McMule will not be performing arbitrary transformations of the phase space. Failing them however is still to be discouraged, if only for completeness.

An example failing the complete rotation test, but not the perpendicular rotation test is

Listing 17 An example measurement function which is not rotation safe
58! This should fail total rotation but pass Z rotation !
59FUNCTION QUANT(q1,q2,q3,q4,q5,q6,q7)
60    ! m(q1) 2 e(q2) n(q3) n(q4) !
61    real (kind=prec), intent(in) :: q1(4),q2(4),q3(4),q4(4), q5(4),q6(4),q7(4)
62    real (kind=prec) :: qHigh(4), qLow(4)
63    real (kind=prec) :: quant(nr_q)
64    !! ==== keep the line below in any case ==== !!
65    call fix_mu
66
67    names(1) = 'ElPerp'
68    quant(1) = sqrt(q2(1)**2 + q2(2)**2)
69
70    names(2) = 'ElLong'
71    quant(2) = q2(3)
72
73END FUNCTION QUANT
74
75
76SUBROUTINE USEREVENT(X, NDIM)
77    integer :: ndim
78    real(kind=prec) :: x(ndim)
79    userweight = 1.
80END SUBROUTINE USEREVENT
$ pymule validate -w m2enngR -l INFO FailsRotation.so
...
INFO:pymule.validator.validator:Test 5: RotationValidator(1.7834175407965058, 0, 0)... passed
WARNING:pymule.validator.validator:Test 6: RotationValidator(4.772496460705239, 3.8977085904094153, 2.9542522534242135)... failed,
INFO:pymule.validator.validator:Test 7: FlipValidator(x)... passed
INFO:pymule.validator.validator:Test 8: FlipValidator(y)... passed
WARNING:pymule.validator.validator:Test 9: FlipValidator(z)... failed,
WARNING:pymule.validator.validator:Overall pass rate: 7 / 9, 77.78%

As we can see, the rotation around Z (Test 5) passed as expected and the arbitrary rotation (Test 6) failed as expected. The flip validator across the Z axis also fails. This is because when looking at the longitudinal momentum the Z momentum is used and will flip sign, thus changing the observable and failing the test. An easy solution to this is to replace line 71 with

quant(2) = abs(q2(3))

thus becoming invariant under flipping across Z.

Minimum and Maximum values, and bin-width testing

Currently the minimum (min_val) and maximum values (max_val) for the histogram as well as the bin widths (determined by number of bins nrbins) are decided by the user before knowing the shape of the data output in the process. This can lead to a time consuming trial-and-error of runs until the correct histogram shape and size is established.

It is intended in the validator to have an additional test which generates many random results and sees if the parameters for the histogram are suitable. This will still (potentially) be a trial-and-error process but will be much quicker than setting up and running the entirety of McMule.

Distance function tests in detail

A distance function may optionally be included in the userfile. This function is used to determine the distance between two points in the momentum phasespace. The simplest example of this would perhaps be a euclidean distance between the spatial components of all the momenta. This distance function must obey certain criteria however.

Triangle validator

A test of the triangle inequality such that:

\[a \lt b + c,\]

where \(a\) is the direct distance between two events, and \(b\) and \(c\) are distances between a third event to the first two.

Event Swap validator

The distance from event 1 (\(e_1\)) to event 2 (\(e_2\)) should be equal to the distance from event 2 to event 1. If D is our distance function then:

\[D(e_1, e_2) = D(e_2, e_1),\]

in a similar way to the particle swapping in Swap symmetry validator

Same event validator

If an event is the same as (indistinguishable to) another event then the distance between these two events should be zero.