diff --git a/.github/gpu_ci/madmatrix_checks.py b/.github/gpu_ci/madmatrix_checks.py index 797d9f78fc..fc11077bb4 100644 --- a/.github/gpu_ci/madmatrix_checks.py +++ b/.github/gpu_ci/madmatrix_checks.py @@ -26,6 +26,13 @@ synchronise with the default stream) gives the same |M|^2 and helicities as on the default stream: every memset, kernel and copy of the matrix element must go to the stream umami is given, the helicity filtering of the first call included. +4. helicity_index: ``g g > t t~ z`` (massless and massive vectors, a longitudinal Z, + massive fermions) evaluated for one given helicity per event (UMAMI_IN_HELICITY_INDEX, + --hel-index), on both backends: the |M|^2 of the ncomb helicities, each given to every + event, add up to the helicity-summed |M|^2, the GPU ones equal the CPU ones helicity by + helicity, and with a random helicity per event (the threads of a warp, or the lanes of + a SIMD vector, on different helicities in lock-step) every event gets the |M|^2 of its + helicity and returns it as the chosen helicity. Any CPU backend but scalar (simd_128, ...) works too, to try the script without a GPU; check 3 is then skipped. Writes summary.txt (key=value lines) in the current directory; @@ -46,7 +53,7 @@ GPU_BACKENDS = ('cuda', 'hip') CPU_BACKEND = 'scalar' -CHECKS = ('mssm_gpu_vs_cpu', 'zero_me', 'nonblocking_stream') +CHECKS = ('mssm_gpu_vs_cpu', 'zero_me', 'nonblocking_stream', 'helicity_index') results = {} @@ -66,7 +73,7 @@ def run(cmd, cwd=None, log=None): def write(path, typecode, values): - """A raw binary file of check_sa.exe: 'd' double, 'I' 32-bit unsigned int.""" + """A raw binary file of check_sa.exe: 'd' double, 'I'/'i' 32-bit unsigned/signed int.""" with open(path, 'wb') as f: array.array(typecode, values).tofile(f) return os.path.abspath(path) @@ -267,6 +274,70 @@ def check_nonblocking_stream(backend, repo, ncalls=4): '%d calls of %d events identical to the default stream' % (ncalls, nevt // ncalls)) +# --------------------------------------------------------------------------- +# 4. one given helicity per event +# --------------------------------------------------------------------------- + +def check_helicity_index(backend, repo, nevt=512): + section('g g > t t~ z, one given helicity per event (%s, %s)' % (backend, CPU_BACKEND)) + name = 'ggttxz_hel' + (proc_dir,) = output(repo, name, 'sm', 'g g > t t~ z', (CPU_BACKEND, backend)) + data = open(os.path.join(proc_dir, 'ProcessData.h')).read() + ncomb = int(re.search(r'constexpr int ncomb = (\d+);', data).group(1)) + batch = [1, nevt, 1] + momenta = dump(name, 'momenta') + check_sa(proc_dir, CPU_BACKEND, ['--energy', 2000, '--dump-momenta', momenta] + batch, + name + '_rambo.log') + rng = random.Random(17) + mixed = [rng.randrange(ncomb) for _ in range(nevt)] + problems, details = [], [] + perhel = {} + for run_backend in (CPU_BACKEND, backend): + prefix = '%s_%s' % (name, run_backend) + check_sa(proc_dir, run_backend, ['--momenta', momenta, '--dump-me', dump(prefix, 'me')] + batch, + prefix + '.log') + summed = read(dump(prefix, 'me'), 'd') + perhel[run_backend] = [] + for ihel in range(ncomb): + hprefix = '%s_h%d' % (prefix, ihel) + hels = write(hprefix + '_in.bin', 'i', [ihel] * nevt) + check_sa(proc_dir, run_backend, ['--momenta', momenta, '--hel-index', hels, + '--dump-me', dump(hprefix, 'me'), + '--dump-hel', dump(hprefix, 'hel')] + batch, + hprefix + '.log') + perhel[run_backend].append(read(dump(hprefix, 'me'), 'd')) + if set(read(dump(hprefix, 'hel'), 'i')) != {ihel}: + problems.append('%s: helicity %d not returned as the chosen helicity' % (run_backend, ihel)) + total = [sum(perhel[run_backend][ihel][ievt] for ihel in range(ncomb)) for ievt in range(nevt)] + rels = [abs(a - b) / abs(b) if b != 0 else float('inf') for a, b in zip(total, summed)] + bad = [ievt for ievt, rel in enumerate(rels) if not rel < 1e-10] + if bad or len(summed) != nevt: + problems.append('%s: the sum over the helicities differs from the helicity sum in %d of %d events' + % (run_backend, len(bad), nevt)) + # a random helicity per event, the same events + mprefix = '%s_mixed' % prefix + hels = write(mprefix + '_in.bin', 'i', mixed) + check_sa(proc_dir, run_backend, ['--momenta', momenta, '--hel-index', hels, + '--dump-me', dump(mprefix, 'me'), + '--dump-hel', dump(mprefix, 'hel')] + batch, + mprefix + '.log') + me, hel = read(dump(mprefix, 'me'), 'd'), read(dump(mprefix, 'hel'), 'i') + expected = [perhel[run_backend][ihel][ievt] for ievt, ihel in enumerate(mixed)] + bad = [ievt for ievt, (a, b) in enumerate(zip(me, expected)) if not abs(a - b) <= 1e-12 * abs(b)] + if bad or list(hel) != mixed: + problems.append('%s: a random helicity per event: %d |M|^2 differ from that helicity\'s, ' + '%d helicities not returned' % (run_backend, len(bad), + sum(1 for a, b in zip(hel, mixed) if a != b))) + details.append('%s: max rel diff of the helicity sum %.1e' % (run_backend, max(rels, default=0.))) + print(' ' + details[-1]) + bad = sum(1 for ihel in range(ncomb) for a, b in zip(perhel[backend][ihel], perhel[CPU_BACKEND][ihel]) + if not abs(a - b) <= 1e-9 * max(abs(b), 1e-300)) + if bad: + problems.append('%d (helicity, event) |M|^2 differ between %s and %s' % (bad, backend, CPU_BACKEND)) + check('helicity_index', not problems, '; '.join(problems[:3]) or + '%d helicities x %d events; %s' % (ncomb, nevt, '; '.join(details))) + + def main(): parser = argparse.ArgumentParser(description=__doc__.split('\n')[0]) parser.add_argument('--repo', required=True) diff --git a/.github/workflows/acceptancetest_mg7.yml b/.github/workflows/acceptancetest_mg7.yml index 2cc8e69569..baa3fdb1a1 100644 --- a/.github/workflows/acceptancetest_mg7.yml +++ b/.github/workflows/acceptancetest_mg7.yml @@ -564,6 +564,39 @@ jobs: ./tests/test_manager.py test_systematics_mg7 -pA -t0 -l INFO + acceptancetest_mg7_madtrex: + needs: build_madspace + # MadtRex (madtrex/, 'install madtrex'): built on top of the cached madspace + # install, then its tests, which all skip without an installed madtrex (and + # a skipped test fails the step, so a broken install cannot pass quietly): + # test_madtrex_standalone (the standalone driver/bin/reweight, the card + # readers, in-place weights) and test_reweight_backends (classic module and + # MadtRex reweight the same events and agree event by event; needs f2py). + runs-on: ubuntu-24.04 + if: github.event_name == 'push' || github.event.pull_request.head.repo.fork == true + steps: + - uses: actions/checkout@v4 + - uses: ./.github/actions/checkout_mg5 + - uses: ./.github/actions/install_madspace + - uses: ./.github/actions/restore_heptools + - uses: ./.github/actions/restore-pip-cache + - name: install madtrex + run: | + cd $GITHUB_WORKSPACE + python3 madtrex/install.py --yes -j 4 + - name: test_madtrex_standalone + if: ${{ !cancelled() }} + run: | + cd $GITHUB_WORKSPACE + ./tests/test_manager.py test_madtrex_standalone -t0 + - name: test_reweight_backends + if: ${{ !cancelled() }} + run: | + cd $GITHUB_WORKSPACE + export PATH="$HOME/.cache/HEPtools/bin:$PATH" + env -u LHAPDF_DATA_PATH ./tests/test_manager.py test_reweight_backends -pA -t0 -l INFO + + acceptancetest_mg7_uncovered: # mg7 acceptance tests that no workflow selected (audit of 2026-10-24). # test_mg7_custom_functions and test_mg7_interference compare to runs with diff --git a/.github/workflows/unittest.yml b/.github/workflows/unittest.yml index 6797dffac5..9fd1bb3584 100644 --- a/.github/workflows/unittest.yml +++ b/.github/workflows/unittest.yml @@ -324,6 +324,17 @@ jobs: cd $GITHUB_WORKSPACE ./tests/test_manager.py test_lhe_ordering test_reweight_ordering test_boost_by_a_momentum_at_rest test_nested_identical_children -t0 + # ---- unittest_reweight_backend ---- + # the reweighting backend choice (classic module or MadtRex) and the + # MadtRex reweight-card readers: pure Python, no madspace/madtrex. The + # tests needing an installed madtrex run in acceptancetest_mg7.yml + # (acceptancetest_mg7_madtrex). + - name: "unittest_reweight_backend: reweighting backend choice and MadtRex cards" + if: ${{ !cancelled() }} + run: | + cd $GITHUB_WORKSPACE + ./tests/test_manager.py test_reweight_backend NpyToLheNamedWeightsTest HelicityReweightingCardTest ClassicVocabularyTest GridpackSourceBuildTest -t0 + # Leave /tmp empty for the next job on this self-hosted machine. - if: always() uses: ./.github/actions/clean_tmp diff --git a/MadSpin/interface_madspin.py b/MadSpin/interface_madspin.py index 8cc58875c8..008627c9eb 100755 --- a/MadSpin/interface_madspin.py +++ b/MadSpin/interface_madspin.py @@ -12064,7 +12064,7 @@ def _density_libtag(self, tree, prod_or_decay): that build's tag.""" tag = madspin.density_libtag(self._ms_run_id, 'p' if prod_or_decay == 'prod' else 'd') - if os.path.exists(pjoin(tree, 'lib', 'libmadmatrix_common%s_scalar.so' % tag)): + if os.path.exists(madmatrix_density.common_library_path(tree, tag)): return tag try: return madmatrix_density.find_libtag(tree) diff --git a/UpdateNotes.txt b/UpdateNotes.txt index 016fe9c68d..d3a63020cc 100644 --- a/UpdateNotes.txt +++ b/UpdateNotes.txt @@ -36,6 +36,42 @@ ANNOUNCEMENT: NNPDF40_nlo_as_01180_qed (lhaid 335900), at NLO also for a model with a massive b (NNPDF4.0 has no 4-flavour QED set). The NNPDF4.0 defaults have no photon. A PDF set by the user is kept. + ZW: MG7 outputs have two reweighting backends: MadtRex ('install madtrex'; + matrix-element and PDF/scale reweighting on lhe or lhe_npy events) and the + classic reweighting module (f2py). Both read Cards/reweight_card.dat with the + classic vocabulary ('change' for options, 'set' for parameters). The backend + comes from 'generate_events --reweight-backend=classic|madtrex|auto', the + reweight switch of the launch question (now OFF/ON/classic/madtrex/density, + ON being auto), or the card's new 'change backend' line; auto prefers + MadtRex when it can run the card. The choice is logged and recorded in + Events//reweighting.json, and a choice that cannot run stops the launch + before the generation. Reweighting runs once, after the systematics; + '--reweight=CARD' now just switches it on with CARD as the reweight card. + MadtRex cards: 'change helicity', 'change fptype', 'change mode LO|pdf', + 'change verbose/max_batch/...' replace the former 'set' spellings (still + accepted, with a warning), 'set MT 175' (a parameter by its param_card + name) works, and directives only the classic module knows are an error + instead of being ignored. MG7 outputs get their own default reweight card + documenting both backends. madevent and aMC@NLO outputs are unchanged. + Both backends add their weights to the run's own event file + (Events//events.lhe.gz, or events.npy for an lhe_npy run reweighted by + MadtRex, as rwgt_ columns that npy_to_lhe carries over), after the + systematics weights; MadtRex no longer writes a separate *_reweighted file + (bin/reweight --output still can). MadtRex reads lhe_npy events with + systematics weights (they were refused as "Invalid header for event file"). + ZW: Without f2py_compiler set, the f2py installed with the running Python (e.g. + in a virtual environment) is used before f2py3/f2py on the PATH, which may + build modules for another Python version (reweighting then failed with "file + too short" or a missing all_matrix2py module). + OM: MadtRex follow-ups: 'install madtrex' builds with clang/macOS; madtrex refuses a + madspace other than the one it was built against (with the command to rebuild it; + MadGraph then reweights with the classic module) and its installer a too old + madspace; a failing reweighting no longer stops the run, and 'reweight RUN' (after + 'launch PROC -i') reweights the events of an earlier run; MadtRex reweight cards + accept '!' comments and multi-index parameters ('set NMIX 1 2 0.5'); the baseline + PDF set and precision of a MadtRex reweighting are those of the run that made the + events (its LHE header), not of the current run_card; a gridpack made with + include_madtrex_source builds madtrex from its sources in bin/reweight. OM: new command 'install cuttools_iregi' recompiling CutTools and IREGI (vendor/) with the current fortran compiler. The error raised when CutTools/IREGI were compiled with another compiler (e.g. after a gfortran update) now points to it. diff --git a/aloha/template_files/madmatrix/helas.h b/aloha/template_files/madmatrix/helas.h index a243ce0d96..f33d986f1d 100644 --- a/aloha/template_files/madmatrix/helas.h +++ b/aloha/template_files/madmatrix/helas.h @@ -26,7 +26,7 @@ __host__ __device__ INLINE void ixxxxx( const fptype_momenta momenta[], // input: momenta const fptype_amp fmass, // input: fermion mass - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fi, // output: aloha objects @@ -41,7 +41,7 @@ __host__ __device__ INLINE void ipzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fi, // output: wavefunctions @@ -56,7 +56,7 @@ __host__ __device__ INLINE void imzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fi, // output: wavefunctions @@ -71,7 +71,7 @@ __host__ __device__ INLINE void ixzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fi, // output: wavefunctions @@ -85,7 +85,7 @@ __host__ __device__ INLINE void vxxxxx( const fptype_momenta momenta[], // input: momenta const fptype_amp vmass, // input: vector boson mass - const int nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson) + const fptype_amp_sv nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson, per SIMD lane/GPU thread) const int nsv, // input: +1 (final) or -1 (initial) const int flv, // input: flavor index ALOHAOBJ & vc, // output: wavefunctions @@ -113,7 +113,7 @@ __host__ __device__ INLINE void oxxxxx( const fptype_momenta momenta[], // input: momenta const fptype_amp fmass, // input: fermion mass - const int nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson) + const fptype_amp_sv nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fo, // output: wavefunctions @@ -128,7 +128,7 @@ __host__ __device__ INLINE void opzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fo, // output: wavefunctions @@ -143,7 +143,7 @@ __host__ __device__ INLINE void omzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fo, // output: wavefunctions @@ -157,7 +157,7 @@ __host__ __device__ INLINE void oxzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fo, // output: wavefunctions @@ -171,7 +171,7 @@ __host__ __device__ void ixxxxx( const fptype_momenta momenta[], // input: momenta const fptype_amp fmass, // input: fermion mass - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fi, // output: wavefunctions @@ -196,7 +196,12 @@ fi.pvec[2] = -pvec2_ * static_cast(nsf); fi.pvec[3] = -pvec3_ * static_cast(nsf); fi.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; if( fmass != 0. ) { #ifndef MGONGPU_CPPSIMD @@ -205,11 +210,9 @@ volatile fptype_amp_sv p2 = pvec1 * pvec1 + pvec2 * pvec2 + pvec3 * pvec3; // volatile fixes #736 const fptype_amp_sv pp = fpmin( pvec0, fpsqrt( p2 ) ); #endif - // In C++ ixxxxx, use a single ip/im numbering that is valid both for pp==0 and pp>0, which have two numbering schemes in Fortran ixxxxx: - // for pp==0, Fortran sqm(0:1) has indexes 0,1 as in C++; but for Fortran pp>0, omega(2) has indexes 1,2 and not 0,1 - // NB: this is only possible in ixxxx, but in oxxxxx two different numbering schemes must be used - const int ip = ( 1 + nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3+nh)/2 because omega(2) has indexes 1,2 - const int im = ( 1 - nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3-nh)/2 because omega(2) has indexes 1,2 + // The Fortran ip/im indices ( 1 + nh ) / 2 and ( 1 - nh ) / 2 (for pp > 0, ( 3 +- nh ) / 2 as omega(2) has + // indexes 1,2) are the truth values tp and tm: an array element x[ip] is tm * x[0] + tp * x[1], and + // ip * x[ip] is just tp * x[1] (and likewise for im) #ifndef MGONGPU_CPPSIMD if( pp == 0. ) { @@ -217,10 +220,10 @@ fptype_amp sqm[2] = { fpsqrt( std::abs( fmass ) ), 0. }; // possibility of negative fermion masses //sqm[1] = ( fmass < 0. ? -abs( sqm[0] ) : abs( sqm[0] ) ); // AV: why abs here? sqm[1] = ( fmass < 0. ? -sqm[0] : sqm[0] ); // AV: removed an abs here - w[0] = cxmake( ip * sqm[ip], 0 ); - w[1] = cxmake( im * nsf * sqm[ip], 0 ); - w[2] = cxmake( ip * nsf * sqm[im], 0 ); - w[3] = cxmake( im * sqm[im], 0 ); + w[0] = cxmake( tp * sqm[1], 0 ); + w[1] = cxmake( tm * nsf * sqm[0], 0 ); + w[2] = cxmake( tp * nsf * sqm[0], 0 ); + w[3] = cxmake( tm * sqm[1], 0 ); } else { @@ -228,30 +231,32 @@ fptype_amp( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype_amp)0.5 }; fptype_amp omega[2] = { fpsqrt( pvec0 + pp ), 0. }; omega[1] = fmass / omega[0]; - const fptype_amp sfomega[2] = { sf[0] * omega[ip], sf[1] * omega[im] }; + const fptype_amp sfomega[2] = { sf[0] * ( tm * omega[0] + tp * omega[1] ), sf[1] * ( tp * omega[0] + tm * omega[1] ) }; const fptype_amp pp3 = fpmax( fpternary( fpsignbit( pvec3 ), ( pvec1 * pvec1 + pvec2 * pvec2 ) / ( pp - pvec3 ), pp + pvec3 ), 0. ); const cxtype_amp chi[2] = { cxmake( fpsqrt( pp3 * (fptype_amp)0.5 / pp ), 0. ), ( pp3 == 0. ? cxmake( -nh, 0. ) : cxmake( nh * pvec1, pvec2 ) / fpsqrt( (fptype_amp)2. * pp * pp3 ) ) }; - w[0] = sfomega[0] * chi[im]; - w[1] = sfomega[0] * chi[ip]; - w[2] = sfomega[1] * chi[im]; - w[3] = sfomega[1] * chi[ip]; + const cxtype_amp chi_im = tp * chi[0] + tm * chi[1]; // chi[im] + const cxtype_amp chi_ip = tm * chi[0] + tp * chi[1]; // chi[ip] + w[0] = sfomega[0] * chi_im; + w[1] = sfomega[0] * chi_ip; + w[2] = sfomega[1] * chi_im; + w[3] = sfomega[1] * chi_ip; } #else // Branch A: pp == 0. // NB: Do not use "abs" for floats! It returns an integer with no build warning! Use std::abs! fptype_amp sqm[2] = { fpsqrt( std::abs( fmass ) ), 0 }; // possibility of negative fermion masses (NB: SCALAR!) sqm[1] = ( fmass < 0 ? -sqm[0] : sqm[0] ); // AV: removed an abs here (as above) - const cxtype_amp fiA_2 = ip * sqm[ip]; // scalar cxtype: real part initialised from fptype, imag part = 0 - const cxtype_amp fiA_3 = im * nsf * sqm[ip]; // scalar cxtype: real part initialised from fptype, imag part = 0 - const cxtype_amp fiA_4 = ip * nsf * sqm[im]; // scalar cxtype: real part initialised from fptype, imag part = 0 - const cxtype_amp fiA_5 = im * sqm[im]; // scalar cxtype: real part initialised from fptype, imag part = 0 + const cxtype_amp_v fiA_2 = cxmake( tp * sqm[1], 0 ); + const cxtype_amp_v fiA_3 = cxmake( tm * (fptype_amp)nsf * sqm[0], 0 ); + const cxtype_amp_v fiA_4 = cxmake( tp * (fptype_amp)nsf * sqm[0], 0 ); + const cxtype_amp_v fiA_5 = cxmake( tm * sqm[1], 0 ); // Branch B: pp != 0. - const fptype_amp sf[2] = { fptype_amp( 1 + nsf + ( 1 - nsf ) * nh ) * (fptype_amp)0.5, - fptype_amp( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype_amp)0.5 }; + const fptype_amp_v sf[2] = { ( (fptype_amp)( 1 + nsf ) + (fptype_amp)( 1 - nsf ) * nh ) * (fptype_amp)0.5, + ( (fptype_amp)( 1 + nsf ) - (fptype_amp)( 1 - nsf ) * nh ) * (fptype_amp)0.5 }; fptype_v omega[2] = { fpsqrt( pvec0 + pp ), 0 }; omega[1] = fmass / omega[0]; - const fptype_amp_v sfomega[2] = { sf[0] * omega[ip], sf[1] * omega[im] }; + const fptype_amp_v sfomega[2] = { sf[0] * ( tm * omega[0] + tp * omega[1] ), sf[1] * ( tp * omega[0] + tm * omega[1] ) }; const fptype_amp_v pp3 = fpmax( fpternary( fpsignbit( pvec3 ), ( pvec1 * pvec1 + pvec2 * pvec2 ) / ( pp - pvec3 ), pp + pvec3 ), 0 ); volatile fptype_amp_v ppDENOM = fpternary( pp != 0, pp, 1. ); // hack: ppDENOM[ieppV]=1 if pp[ieppV]==0 volatile fptype_amp_v pp3DENOM = fpternary( pp3 != 0, pp3, 1. ); // hack: pp3DENOM[ieppV]=1 if pp3[ieppV]==0 @@ -259,11 +264,13 @@ const cxtype_amp_v chi[2] = { cxmake( fpsqrt( chi0r2 ), 0 ), // hack: dummy[ieppV] is not used if pp[ieppV]==0 cxternary( ( pp3 == 0. ), cxmake( -nh, 0 ), - cxmake( (fptype_amp)nh * pvec1, pvec2 ) / fpsqrt( (fptype_amp)2. * ppDENOM * pp3DENOM ) ) }; // hack: dummy[ieppV] is not used if pp[ieppV]==0 - const cxtype_amp_v fiB_2 = sfomega[0] * chi[im]; - const cxtype_amp_v fiB_3 = sfomega[0] * chi[ip]; - const cxtype_amp_v fiB_4 = sfomega[1] * chi[im]; - const cxtype_amp_v fiB_5 = sfomega[1] * chi[ip]; + cxmake( nh * pvec1, pvec2 ) / fpsqrt( (fptype_amp)2. * ppDENOM * pp3DENOM ) ) }; // hack: dummy[ieppV] is not used if pp[ieppV]==0 + const cxtype_amp_v chi_im = tp * chi[0] + tm * chi[1]; // chi[im] + const cxtype_amp_v chi_ip = tm * chi[0] + tp * chi[1]; // chi[ip] + const cxtype_amp_v fiB_2 = sfomega[0] * chi_im; + const cxtype_amp_v fiB_3 = sfomega[0] * chi_ip; + const cxtype_amp_v fiB_4 = sfomega[1] * chi_im; + const cxtype_amp_v fiB_5 = sfomega[1] * chi_ip; // Choose between the results from branch A and branch B const bool_v mask = ( pp == 0. ); w[0] = cxternary( mask, fiA_2, fiB_2 ); @@ -289,8 +296,8 @@ volatile fptype_amp_sv sqp0p3DENOM = fpternary( sqp0p3 != 0, (fptype_sv)sqp0p3, 1. ); // hack: dummy sqp0p3DENOM[ieppV]=1 if sqp0p3[ieppV]==0 cxtype_sv chi[2] = { cxmake( (fptype_v)sqp0p3, 0. ), cxternary( sqp0p3 == 0, - cxmake( -(fptype_amp)nhel * fpsqrt( (fptype_amp)2. * pvec0 ), 0. ), - cxmake( (fptype_amp)nh * pvec1, pvec2 ) / (const fptype_v)sqp0p3DENOM ) }; // hack: dummy[ieppV] is not used if sqp0p3[ieppV]==0 + cxmake( -nhel * fpsqrt( (fptype_amp)2. * pvec0 ), 0. ), + cxmake( nh * pvec1, pvec2 ) / (const fptype_v)sqp0p3DENOM ) }; // hack: dummy[ieppV] is not used if sqp0p3[ieppV]==0 #else // pvec0+pvec3 is a cancelling difference of two ~|p| numbers for a // backward-moving massless fermion, and everything below divides by it. @@ -303,22 +310,12 @@ fptype_sv{ 0 }, fpsqrt( p0p3stable ) * (fptype_amp)nsf ); const cxtype_amp_sv chi[2] = { cxmake( sqp0p3, 0. ), - ( sqp0p3 == 0. ? cxmake( -(fptype_amp)nhel * fpsqrt( (fptype_amp)2. * pvec0 ), 0. ) : cxmake( (fptype_amp)nh * pvec1, pvec2 ) / sqp0p3 ) }; + ( sqp0p3 == 0. ? cxmake( -nhel * fpsqrt( (fptype_amp)2. * pvec0 ), 0. ) : cxmake( nh * pvec1, pvec2 ) / sqp0p3 ) }; #endif - if( nh == 1 ) - { - w[0] = cxzero_sv(); - w[1] = cxzero_sv(); - w[2] = chi[0]; - w[3] = chi[1]; - } - else - { - w[0] = chi[1]; - w[1] = chi[0]; - w[2] = cxzero_sv(); - w[3] = cxzero_sv(); - } + w[0] = tm * chi[1]; + w[1] = tm * chi[0]; + w[2] = tp * chi[0]; + w[3] = tp * chi[1]; } mgDebug( 1, __FUNCTION__ ); return; @@ -332,7 +329,7 @@ __host__ __device__ void ipzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fi, // output: wavefunctions @@ -347,20 +344,17 @@ fi.pvec[2] = fptype_momenta_sv{ 0 }; fi.pvec[3] = -pvec3_ * static_cast(nsf); fi.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; const cxtype_amp_sv sqp0p3 = cxmake( fpsqrt( (fptype_amp)2. * pvec3 ) * (fptype_amp)nsf, 0. ); - w[0] = cxmake( fpamp_of_mom( fi.pvec[1] ), fpamp_of_mom( fi.pvec[2] ) ); - if( nh == 1 ) - { - w[1] = cxmake( fpamp_of_mom( fi.pvec[1] ), fpamp_of_mom( fi.pvec[2] ) ); - w[2] = sqp0p3; - } - else - { - w[1] = sqp0p3; - w[2] = cxmake( fpamp_of_mom( fi.pvec[1] ), fpamp_of_mom( fi.pvec[2] ) ); - } - w[3] = cxmake( fpamp_of_mom( fi.pvec[1] ), fpamp_of_mom( fi.pvec[2] ) ); + w[0] = cxzero_sv(); + w[1] = tm * sqp0p3; + w[2] = tp * sqp0p3; + w[3] = cxzero_sv(); mgDebug( 1, __FUNCTION__ ); return; } @@ -373,7 +367,7 @@ __host__ __device__ void imzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fi, // output: wavefunctions @@ -388,20 +382,17 @@ fi.pvec[2] = fptype_momenta_sv{ 0 }; fi.pvec[3] = -pvec3_ * static_cast(nsf); fi.flv_index = flv; - const int nh = nhel * nsf; - const cxtype_amp_sv chi = cxmake( -(fptype_amp)nhel * fpsqrt( -(fptype_amp)2. * pvec3 ), 0. ); + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; + const cxtype_amp_sv chi = cxmake( -nhel * fpsqrt( -(fptype_amp)2. * pvec3 ), 0. ); + w[0] = tm * chi; w[1] = cxzero_sv(); w[2] = cxzero_sv(); - if( nh == 1 ) - { - w[0] = cxzero_sv(); - w[3] = chi; - } - else - { - w[0] = chi; - w[3] = cxzero_sv(); - } + w[3] = tp * chi; mgDebug( 1, __FUNCTION__ ); return; } @@ -414,7 +405,7 @@ __host__ __device__ void ixzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fi, // output: wavefunctions @@ -436,7 +427,12 @@ fi.pvec[3] = -pvec3_ * static_cast(nsf); fi.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; //const float sqp0p3 = sqrtf( pvec0 + pvec3 ) * nsf; // AV: why force a float here? // pvec0+pvec3 is a cancelling difference of two ~|p| numbers for a // backward-moving massless fermion, and everything below divides by it. @@ -448,21 +444,11 @@ pvec0 + pvec3 ); const fptype_amp_sv sqp0p3 = fpsqrt( p0p3stable ) * static_cast(nsf); const cxtype_amp_sv chi0 = cxmake( sqp0p3, 0. ); - const cxtype_amp_sv chi1 = cxmake( (fptype_amp)nh * pvec1 / sqp0p3, pvec2 / sqp0p3 ); - if( nh == 1 ) - { - w[0] = cxzero_sv(); - w[1] = cxzero_sv(); - w[2] = chi0; - w[3] = chi1; - } - else - { - w[0] = chi1; - w[1] = chi0; - w[2] = cxzero_sv(); - w[3] = cxzero_sv(); - } + const cxtype_amp_sv chi1 = cxmake( nh * pvec1 / sqp0p3, pvec2 / sqp0p3 ); + w[0] = tm * chi1; + w[1] = tm * chi0; + w[2] = tp * chi0; + w[3] = tp * chi1; mgDebug( 1, __FUNCTION__ ); return; } @@ -474,7 +460,7 @@ __host__ __device__ void vxxxxx( const fptype_momenta momenta[], // input: momenta const fptype_amp vmass, // input: vector boson mass - const int nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson) + const fptype_amp_sv nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson, per SIMD lane/GPU thread) const int nsv, // input: +1 (final) or -1 (initial) const int flv, // input: flavour ALOHAOBJ & vc, // output: wavefunctions @@ -500,11 +486,14 @@ vc.pvec[3] = pvec3_ * static_cast(nsv); vc.flv_index = flv; const fptype_amp sqh = fpsqrt( 0.5 ); // AV this is > 0! - const fptype_amp hel = nhel; + const fptype_amp_sv hel = nhel; + // BRANCHLESS HELICITY: |hel| is hel * hel for hel in { -1, 0, +1 }, so that the helicity can + // differ from one SIMD lane or GPU thread to the next with no branch and no std::abs + const fptype_amp_sv abshel = hel * hel; if( vmass != 0. ) { - const int nsvahl = nsv * std::abs( hel ); - const fptype_amp hel0 = 1. - std::abs( hel ); + const fptype_amp_sv nsvahl = static_cast( nsv ) * abshel; + const fptype_amp_sv hel0 = 1 - abshel; // truth value of hel == 0 (longitudinal polarisation) #ifndef MGONGPU_CPPSIMD const fptype_amp_sv pt2 = ( pvec1 * pvec1 ) + ( pvec2 * pvec2 ); const fptype_amp_sv pp = fpmin( pvec0, fpsqrt( pt2 + ( pvec3 * pvec3 ) ) ); @@ -545,9 +534,9 @@ const fptype_amp_sv pt = fpmin( pp, fpsqrt( pt2 ) ); // Branch A: pp == 0. const cxtype_amp vcA_2 = cxmake( 0, 0 ); - const cxtype_amp vcA_3 = cxmake( -hel * sqh, 0 ); - const cxtype_amp vcA_4 = cxmake( 0, nsvahl * sqh ); - const cxtype_amp vcA_5 = cxmake( hel0, 0 ); + const cxtype_amp_v vcA_3 = cxmake( -hel * sqh, 0 ); + const cxtype_amp_v vcA_4 = cxmake( 0., nsvahl * sqh ); + const cxtype_amp_v vcA_5 = cxmake( hel0, 0 ); // Branch B: pp != 0. volatile fptype_amp_v ppDENOM = fpternary( pp != 0, pp, 1. ); // hack: ppDENOM[ieppV]=1 if pp[ieppV]==0 const fptype_amp_v emp = pvec0 / ( vmass * ppDENOM ); // hack: dummy[ieppV] is not used if pp[ieppV]==0 @@ -556,11 +545,11 @@ // Branch B1: pp != 0. and pt != 0. volatile fptype_amp_v ptDENOM = fpternary( pt != 0, pt, 1. ); // hack: ptDENOM[ieppV]=1 if pt[ieppV]==0 const fptype_amp_v pzpt = pvec3 / ( ppDENOM * ptDENOM ) * sqh * hel; // hack: dummy[ieppV] is not used if pp[ieppV]==0 - const cxtype_amp_v vcB1_3 = cxmake( hel0 * pvec1 * emp - pvec1 * pzpt, -fpamp_scalar(nsvahl) * pvec2 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 - const cxtype_amp_v vcB1_4 = cxmake( hel0 * pvec2 * emp - pvec2 * pzpt, fpamp_scalar(nsvahl) * pvec1 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 + const cxtype_amp_v vcB1_3 = cxmake( hel0 * pvec1 * emp - pvec1 * pzpt, -nsvahl * pvec2 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 + const cxtype_amp_v vcB1_4 = cxmake( hel0 * pvec2 * emp - pvec2 * pzpt, nsvahl * pvec1 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 // Branch B2: pp != 0. and pt == 0. - const cxtype_amp vcB2_3 = cxmake( -hel * sqh, 0. ); - const cxtype_amp_v vcB2_4 = cxmake( 0., fpamp_scalar(nsvahl) * fpternary( ( pvec3 < 0 ), -sqh, sqh ) ); // AV: removed an abs here + const cxtype_amp_v vcB2_3 = cxmake( -hel * sqh, 0. ); + const cxtype_amp_v vcB2_4 = cxmake( 0., nsvahl * fpternary( ( pvec3 < 0 ), -sqh, sqh ) ); // AV: removed an abs here // Choose between the results from branch A and branch B (and from branch B1 and branch B2) const bool_v mask = ( pp == 0. ); const bool_v maskB = ( pt != 0. ); @@ -599,11 +588,11 @@ // Branch A: pt != 0. volatile fptype_amp_v ptDENOM = fpternary( pt != 0, pt, 1. ); // hack: ptDENOM[ieppV]=1 if pt[ieppV]==0 const fptype_amp_v pzpt = pvec3 / ( pp * ptDENOM ) * sqh * hel; // hack: dummy[ieppV] is not used if pt[ieppV]==0 - const cxtype_amp_v vcA_3 = cxmake( -pvec1 * pzpt, -fpamp_scalar(nsv) * pvec2 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 - const cxtype_amp_v vcA_4 = cxmake( -pvec2 * pzpt, fpamp_scalar(nsv) * pvec1 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 + const cxtype_amp_v vcA_3 = cxmake( -pvec1 * pzpt, -(fptype_amp)nsv * pvec2 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 + const cxtype_amp_v vcA_4 = cxmake( -pvec2 * pzpt, (fptype_amp)nsv * pvec1 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 // Branch B: pt == 0. - const cxtype_amp vcB_3 = cxmake( -(fptype_amp)hel * sqh, 0 ); - const cxtype_amp_v vcB_4 = cxmake( 0, fpamp_scalar(nsv) * fpternary( ( pvec3 < 0 ), -sqh, sqh ) ); // AV: removed an abs here + const cxtype_amp_v vcB_3 = cxmake( -hel * sqh, 0 ); + const cxtype_amp_v vcB_4 = cxmake( 0, (fptype_amp)nsv * fpternary( ( pvec3 < 0 ), -sqh, sqh ) ); // AV: removed an abs here // Choose between the results from branch A and branch B const bool_v mask = ( pt != 0. ); w[1] = cxternary( mask, vcA_3, vcB_3 ); @@ -654,7 +643,7 @@ __host__ __device__ void oxxxxx( const fptype_momenta momenta[], // input: momenta const fptype_amp fmass, // input: fermion mass - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) int flv, // input: flavour ALOHAOBJ & fo, // output: wavefunctions @@ -680,7 +669,12 @@ fo.pvec[2] = pvec2_ * static_cast(nsf); fo.pvec[3] = pvec3_ * static_cast(nsf); fo.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; if( fmass != 0. ) { #ifndef MGONGPU_CPPSIMD @@ -691,12 +685,12 @@ fptype_amp sqm[2] = { fpsqrt( std::abs( fmass ) ), 0. }; // possibility of negative fermion masses //sqm[1] = ( fmass < 0. ? -abs( sqm[0] ) : abs( sqm[0] ) ); // AV: why abs here? sqm[1] = ( fmass < 0. ? -sqm[0] : sqm[0] ); // AV: removed an abs here - const int ip = -( ( 1 - nh ) / 2 ) * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ - const int im = ( 1 + nh ) / 2 * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ - w[0] = cxmake( im * sqm[std::abs( ip )], 0 ); - w[1] = cxmake( ip * nsf * sqm[std::abs( ip )], 0 ); - w[2] = cxmake( im * nsf * sqm[std::abs( im )], 0 ); - w[3] = cxmake( ip * sqm[std::abs( im )], 0 ); + // The Fortran ip = -( ( 1 - nh ) / 2 ) * nhel and im = ( 1 + nh ) / 2 * nhel are -tm * nhel and + // tp * nhel, so that |ip| = tm and |im| = tp: sqm[|ip|] is sqm[1] wherever ip != 0, and so on + w[0] = cxmake( tp * nhel * sqm[0], 0 ); + w[1] = cxmake( -tm * nhel * nsf * sqm[1], 0 ); + w[2] = cxmake( tp * nhel * nsf * sqm[1], 0 ); + w[3] = cxmake( -tm * nhel * sqm[0], 0 ); } else { @@ -704,17 +698,18 @@ fptype_amp( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype_amp)0.5 }; fptype_amp omega[2] = { fpsqrt( pvec0 + pp ), 0. }; omega[1] = fmass / omega[0]; - const int ip = ( 1 + nh ) / 2; // NB: Fortran is (3+nh)/2 because omega(2) has indexes 1,2 and not 0,1 - const int im = ( 1 - nh ) / 2; // NB: Fortran is (3-nh)/2 because omega(2) has indexes 1,2 and not 0,1 - const fptype_amp sfomeg[2] = { sf[0] * omega[ip], sf[1] * omega[im] }; + // The Fortran ip/im indices ( 3 +- nh ) / 2 (omega(2) has indexes 1,2) are the truth values tp/tm + const fptype_amp sfomeg[2] = { sf[0] * ( tm * omega[0] + tp * omega[1] ), sf[1] * ( tp * omega[0] + tm * omega[1] ) }; const fptype_amp pp3 = fpmax( fpternary( fpsignbit( pvec3 ), ( pvec1 * pvec1 + pvec2 * pvec2 ) / ( pp - pvec3 ), pp + pvec3 ), 0. ); const cxtype_amp chi[2] = { cxmake( fpsqrt( pp3 * (fptype_amp)0.5 / pp ), 0. ), ( ( pp3 == 0. ) ? cxmake( -nh, 0. ) : cxmake( nh * pvec1, -pvec2 ) / fpsqrt( (fptype_amp)2. * pp * pp3 ) ) }; - w[0] = sfomeg[1] * chi[im]; - w[1] = sfomeg[1] * chi[ip]; - w[2] = sfomeg[0] * chi[im]; - w[3] = sfomeg[0] * chi[ip]; + const cxtype_amp chi_im = tp * chi[0] + tm * chi[1]; // chi[im] + const cxtype_amp chi_ip = tm * chi[0] + tp * chi[1]; // chi[ip] + w[0] = sfomeg[1] * chi_im; + w[1] = sfomeg[1] * chi_ip; + w[2] = sfomeg[0] * chi_im; + w[3] = sfomeg[0] * chi_ip; } #else volatile fptype_amp_sv p2 = pvec1 * pvec1 + pvec2 * pvec2 + pvec3 * pvec3; // volatile fixes #736 @@ -723,20 +718,17 @@ // NB: Do not use "abs" for floats! It returns an integer with no build warning! Use std::abs! fptype_amp sqm[2] = { fpsqrt( std::abs( fmass ) ), 0 }; // possibility of negative fermion masses sqm[1] = ( fmass < 0 ? -sqm[0] : sqm[0] ); // AV: removed an abs here (as above) - const int ipA = -( ( 1 - nh ) / 2 ) * nhel; - const int imA = ( 1 + nh ) / 2 * nhel; - const cxtype_amp foA_2 = imA * sqm[std::abs( ipA )]; - const cxtype_amp foA_3 = ipA * nsf * sqm[std::abs( ipA )]; - const cxtype_amp foA_4 = imA * nsf * sqm[std::abs( imA )]; - const cxtype_amp foA_5 = ipA * sqm[std::abs( imA )]; + // (see the scalar code above for the truth values tp/tm in place of the Fortran ip/im indices) + const cxtype_amp_v foA_2 = cxmake( tp * nhel * sqm[0], 0 ); + const cxtype_amp_v foA_3 = cxmake( -tm * nhel * (fptype_amp)nsf * sqm[1], 0 ); + const cxtype_amp_v foA_4 = cxmake( tp * nhel * (fptype_amp)nsf * sqm[1], 0 ); + const cxtype_amp_v foA_5 = cxmake( -tm * nhel * sqm[0], 0 ); // Branch B: pp != 0. - const fptype_amp sf[2] = { fptype_amp( 1 + nsf + ( 1 - nsf ) * nh ) * (fptype_amp)0.5, - fptype_amp( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype_amp)0.5 }; + const fptype_amp_v sf[2] = { ( (fptype_amp)( 1 + nsf ) + (fptype_amp)( 1 - nsf ) * nh ) * (fptype_amp)0.5, + ( (fptype_amp)( 1 + nsf ) - (fptype_amp)( 1 - nsf ) * nh ) * (fptype_amp)0.5 }; fptype_v omega[2] = { fpsqrt( pvec0 + pp ), 0 }; omega[1] = fmass / omega[0]; - const int ipB = ( 1 + nh ) / 2; - const int imB = ( 1 - nh ) / 2; - const fptype_amp_v sfomeg[2] = { sf[0] * omega[ipB], sf[1] * omega[imB] }; + const fptype_amp_v sfomeg[2] = { sf[0] * ( tm * omega[0] + tp * omega[1] ), sf[1] * ( tp * omega[0] + tm * omega[1] ) }; const fptype_amp_v pp3 = fpmax( fpternary( fpsignbit( pvec3 ), ( pvec1 * pvec1 + pvec2 * pvec2 ) / ( pp - pvec3 ), pp + pvec3 ), 0. ); volatile fptype_amp_v ppDENOM = fpternary( pp != 0, pp, 1. ); // hack: ppDENOM[ieppV]=1 if pp[ieppV]==0 volatile fptype_amp_v pp3DENOM = fpternary( pp3 != 0, pp3, 1. ); // hack: pp3DENOM[ieppV]=1 if pp3[ieppV]==0 @@ -744,11 +736,13 @@ const cxtype_amp_v chi[2] = { cxmake( fpsqrt( chi0r2 ), 0. ), // hack: dummy[ieppV] is not used if pp[ieppV]==0 ( cxternary( ( pp3 == 0. ), cxmake( -nh, 0. ), - cxmake( (fptype_amp)nh * pvec1, -pvec2 ) / fpsqrt( (fptype_amp)2. * ppDENOM * pp3DENOM ) ) ) }; // hack: dummy[ieppV] is not used if pp[ieppV]==0 - const cxtype_amp_v foB_2 = sfomeg[1] * chi[imB]; - const cxtype_amp_v foB_3 = sfomeg[1] * chi[ipB]; - const cxtype_amp_v foB_4 = sfomeg[0] * chi[imB]; - const cxtype_amp_v foB_5 = sfomeg[0] * chi[ipB]; + cxmake( nh * pvec1, -pvec2 ) / fpsqrt( (fptype_amp)2. * ppDENOM * pp3DENOM ) ) ) }; // hack: dummy[ieppV] is not used if pp[ieppV]==0 + const cxtype_amp_v chi_im = tp * chi[0] + tm * chi[1]; // chi[im] + const cxtype_amp_v chi_ip = tm * chi[0] + tp * chi[1]; // chi[ip] + const cxtype_amp_v foB_2 = sfomeg[1] * chi_im; + const cxtype_amp_v foB_3 = sfomeg[1] * chi_ip; + const cxtype_amp_v foB_4 = sfomeg[0] * chi_im; + const cxtype_amp_v foB_5 = sfomeg[0] * chi_ip; // Choose between the results from branch A and branch B const bool_v mask = ( pp == 0. ); w[0] = cxternary( mask, foA_2, foB_2 ); @@ -775,7 +769,7 @@ const cxtype_amp_v chi[2] = { cxmake( (fptype_v)sqp0p3, 0. ), cxternary( ( sqp0p3 == 0. ), cxmake( -nhel, 0. ) * fpsqrt( (fptype_amp)2. * pvec0 ), - cxmake( (fptype_amp)nh * pvec1, -pvec2 ) / (const fptype_amp_sv)sqp0p3DENOM ) }; // hack: dummy[ieppV] is not used if sqp0p3[ieppV]==0 + cxmake( nh * pvec1, -pvec2 ) / (const fptype_amp_sv)sqp0p3DENOM ) }; // hack: dummy[ieppV] is not used if sqp0p3[ieppV]==0 #else // pvec0+pvec3 is a cancelling difference of two ~|p| numbers for a // backward-moving massless fermion, and everything below divides by it. @@ -788,22 +782,12 @@ fptype_sv{ 0 }, fpsqrt( p0p3stable ) * (fptype_amp)nsf ); const cxtype_amp_sv chi[2] = { cxmake( sqp0p3, 0. ), - ( sqp0p3 == 0. ? cxmake( -nhel, 0. ) * fpsqrt( (fptype_amp)2. * pvec0 ) : cxmake( (fptype_amp)nh * pvec1, -pvec2 ) / sqp0p3 ) }; + ( sqp0p3 == 0. ? cxmake( -nhel, 0. ) * fpsqrt( (fptype_amp)2. * pvec0 ) : cxmake( nh * pvec1, -pvec2 ) / sqp0p3 ) }; #endif - if( nh == 1 ) - { - w[0] = chi[0]; - w[1] = chi[1]; - w[2] = cxzero_sv(); - w[3] = cxzero_sv(); - } - else - { - w[0] = cxzero_sv(); - w[1] = cxzero_sv(); - w[2] = chi[1]; - w[3] = chi[0]; - } + w[0] = tp * chi[0]; + w[1] = tp * chi[1]; + w[2] = tm * chi[1]; + w[3] = tm * chi[0]; } mgDebug( 1, __FUNCTION__ ); return; @@ -817,7 +801,7 @@ __host__ __device__ void opzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fo, // output: wavefunctions @@ -832,20 +816,17 @@ fo.pvec[2] = fptype_momenta_sv{ 0 }; fo.pvec[3] = pvec3_ * static_cast(nsf); fo.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; const cxtype_amp_sv csqp0p3 = cxmake( fpsqrt( (fptype_amp)2. * pvec3 ) * (fptype_amp)nsf, 0. ); + w[0] = tp * csqp0p3; w[1] = cxzero_sv(); w[2] = cxzero_sv(); - if( nh == 1 ) - { - w[0] = csqp0p3; - w[3] = cxzero_sv(); - } - else - { - w[0] = cxzero_sv(); - w[3] = csqp0p3; - } + w[3] = tm * csqp0p3; mgDebug( 1, __FUNCTION__ ); return; } @@ -858,7 +839,7 @@ __host__ __device__ void omzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fo, // output: wavefunctions @@ -873,23 +854,17 @@ fo.pvec[2] = fptype_momenta_sv{ 0 }; fo.pvec[3] = pvec3_ * static_cast(nsf); fo.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; const cxtype_amp_sv chi1 = cxmake( -nhel, 0. ) * fpsqrt( -(fptype_amp)2. * pvec3 ); - if( nh == 1 ) - { - w[0] = cxzero_sv(); - w[1] = chi1; - w[2] = cxzero_sv(); - w[3] = cxzero_sv(); - } - else - { - w[0] = cxzero_sv(); - w[1] = cxzero_sv(); - w[2] = chi1; - //w[3] = chi1; // AV: BUG! - w[3] = cxzero_sv(); // AV: BUG FIX - } + w[0] = cxzero_sv(); + w[1] = tp * chi1; + w[2] = tm * chi1; + w[3] = cxzero_sv(); mgDebug( 1, __FUNCTION__ ); return; } @@ -902,7 +877,7 @@ __host__ __device__ void oxzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fo, // output: wavefunctions @@ -924,7 +899,12 @@ fo.pvec[2] = pvec2_ * static_cast(nsf); fo.pvec[3] = pvec3_ * static_cast(nsf); fo.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; //const float sqp0p3 = sqrtf( pvec0 + pvec3 ) * nsf; // AV: why force a float here? // pvec0+pvec3 is a cancelling difference of two ~|p| numbers for a // backward-moving massless fermion, and everything below divides by it. @@ -936,21 +916,11 @@ pvec0 + pvec3 ); const fptype_amp_sv sqp0p3 = fpsqrt( p0p3stable ) * static_cast(nsf); const cxtype_amp_sv chi0 = cxmake( sqp0p3, 0. ); - const cxtype_amp_sv chi1 = cxmake( (fptype_amp)nh * pvec1 / sqp0p3, -pvec2 / sqp0p3 ); - if( nh == 1 ) - { - w[0] = chi0; - w[1] = chi1; - w[2] = cxzero_sv(); - w[3] = cxzero_sv(); - } - else - { - w[0] = cxzero_sv(); - w[1] = cxzero_sv(); - w[2] = chi1; - w[3] = chi0; - } + const cxtype_amp_sv chi1 = cxmake( nh * pvec1 / sqp0p3, -pvec2 / sqp0p3 ); + w[0] = tp * chi0; + w[1] = tp * chi1; + w[2] = tm * chi1; + w[3] = tm * chi0; mgDebug( 1, __FUNCTION__ ); return; } diff --git a/aloha/template_files/madmatrix/helas_fd.h b/aloha/template_files/madmatrix/helas_fd.h index 2aed9342cf..723037a550 100644 --- a/aloha/template_files/madmatrix/helas_fd.h +++ b/aloha/template_files/madmatrix/helas_fd.h @@ -26,7 +26,7 @@ __host__ __device__ INLINE void ixxxxx( const fptype_momenta momenta[], // input: momenta const fptype_amp fmass, // input: fermion mass - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fi, // output: aloha objects @@ -41,7 +41,7 @@ __host__ __device__ INLINE void ipzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fi, // output: wavefunctions @@ -56,7 +56,7 @@ __host__ __device__ INLINE void imzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fi, // output: wavefunctions @@ -71,7 +71,7 @@ __host__ __device__ INLINE void ixzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fi, // output: wavefunctions @@ -85,7 +85,7 @@ __host__ __device__ INLINE void vxxxxx( const fptype_momenta momenta[], // input: momenta const fptype_amp vmass, // input: vector boson mass - const int nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson) + const fptype_amp_sv nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson, per SIMD lane/GPU thread) const int nsv, // input: +1 (final) or -1 (initial) const int flv, // input: flavor index ALOHAOBJ & vc, // output: wavefunctions @@ -113,7 +113,7 @@ __host__ __device__ INLINE void oxxxxx( const fptype_momenta momenta[], // input: momenta const fptype_amp fmass, // input: fermion mass - const int nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson) + const fptype_amp_sv nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fo, // output: wavefunctions @@ -128,7 +128,7 @@ __host__ __device__ INLINE void opzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fo, // output: wavefunctions @@ -143,7 +143,7 @@ __host__ __device__ INLINE void omzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fo, // output: wavefunctions @@ -157,7 +157,7 @@ __host__ __device__ INLINE void oxzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavor index ALOHAOBJ & fo, // output: wavefunctions @@ -198,7 +198,7 @@ __host__ __device__ void ixxxxx( const fptype_momenta momenta[], // input: momenta const fptype_amp fmass, // input: fermion mass - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fi, // output: wavefunctions @@ -223,7 +223,12 @@ fi.pvec[2] = -pvec2_ * static_cast(nsf); fi.pvec[3] = -pvec3_ * static_cast(nsf); fi.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; if( fmass != 0. ) { #ifndef MGONGPU_CPPSIMD @@ -232,11 +237,9 @@ volatile fptype_amp_sv p2 = pvec1 * pvec1 + pvec2 * pvec2 + pvec3 * pvec3; // volatile fixes #736 const fptype_amp_sv pp = fpmin( pvec0, fpsqrt( p2 ) ); #endif - // In C++ ixxxxx, use a single ip/im numbering that is valid both for pp==0 and pp>0, which have two numbering schemes in Fortran ixxxxx: - // for pp==0, Fortran sqm(0:1) has indexes 0,1 as in C++; but for Fortran pp>0, omega(2) has indexes 1,2 and not 0,1 - // NB: this is only possible in ixxxx, but in oxxxxx two different numbering schemes must be used - const int ip = ( 1 + nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3+nh)/2 because omega(2) has indexes 1,2 - const int im = ( 1 - nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3-nh)/2 because omega(2) has indexes 1,2 + // The Fortran ip/im indices ( 1 + nh ) / 2 and ( 1 - nh ) / 2 (for pp > 0, ( 3 +- nh ) / 2 as omega(2) has + // indexes 1,2) are the truth values tp and tm: an array element x[ip] is tm * x[0] + tp * x[1], and + // ip * x[ip] is just tp * x[1] (and likewise for im) #ifndef MGONGPU_CPPSIMD if( pp == 0. ) { @@ -244,10 +247,10 @@ fptype_amp sqm[2] = { fpsqrt( std::abs( fmass ) ), 0. }; // possibility of negative fermion masses //sqm[1] = ( fmass < 0. ? -abs( sqm[0] ) : abs( sqm[0] ) ); // AV: why abs here? sqm[1] = ( fmass < 0. ? -sqm[0] : sqm[0] ); // AV: removed an abs here - w[0] = cxmake( ip * sqm[ip], 0 ); - w[1] = cxmake( im * nsf * sqm[ip], 0 ); - w[2] = cxmake( ip * nsf * sqm[im], 0 ); - w[3] = cxmake( im * sqm[im], 0 ); + w[0] = cxmake( tp * sqm[1], 0 ); + w[1] = cxmake( tm * nsf * sqm[0], 0 ); + w[2] = cxmake( tp * nsf * sqm[0], 0 ); + w[3] = cxmake( tm * sqm[1], 0 ); } else { @@ -255,30 +258,32 @@ fptype_amp( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype_amp)0.5 }; fptype_amp omega[2] = { fpsqrt( pvec0 + pp ), 0. }; omega[1] = fmass / omega[0]; - const fptype_amp sfomega[2] = { sf[0] * omega[ip], sf[1] * omega[im] }; + const fptype_amp sfomega[2] = { sf[0] * ( tm * omega[0] + tp * omega[1] ), sf[1] * ( tp * omega[0] + tm * omega[1] ) }; const fptype_amp pp3 = fpmax( fpternary( fpsignbit( pvec3 ), ( pvec1 * pvec1 + pvec2 * pvec2 ) / ( pp - pvec3 ), pp + pvec3 ), 0. ); const cxtype_amp chi[2] = { cxmake( fpsqrt( pp3 * (fptype_amp)0.5 / pp ), 0. ), ( pp3 == 0. ? cxmake( -nh, 0. ) : cxmake( nh * pvec1, pvec2 ) / fpsqrt( (fptype_amp)2. * pp * pp3 ) ) }; - w[0] = sfomega[0] * chi[im]; - w[1] = sfomega[0] * chi[ip]; - w[2] = sfomega[1] * chi[im]; - w[3] = sfomega[1] * chi[ip]; + const cxtype_amp chi_im = tp * chi[0] + tm * chi[1]; // chi[im] + const cxtype_amp chi_ip = tm * chi[0] + tp * chi[1]; // chi[ip] + w[0] = sfomega[0] * chi_im; + w[1] = sfomega[0] * chi_ip; + w[2] = sfomega[1] * chi_im; + w[3] = sfomega[1] * chi_ip; } #else // Branch A: pp == 0. // NB: Do not use "abs" for floats! It returns an integer with no build warning! Use std::abs! fptype_amp sqm[2] = { fpsqrt( std::abs( fmass ) ), 0 }; // possibility of negative fermion masses (NB: SCALAR!) sqm[1] = ( fmass < 0 ? -sqm[0] : sqm[0] ); // AV: removed an abs here (as above) - const cxtype_amp fiA_2 = ip * sqm[ip]; // scalar cxtype: real part initialised from fptype, imag part = 0 - const cxtype_amp fiA_3 = im * nsf * sqm[ip]; // scalar cxtype: real part initialised from fptype, imag part = 0 - const cxtype_amp fiA_4 = ip * nsf * sqm[im]; // scalar cxtype: real part initialised from fptype, imag part = 0 - const cxtype_amp fiA_5 = im * sqm[im]; // scalar cxtype: real part initialised from fptype, imag part = 0 + const cxtype_amp_v fiA_2 = cxmake( tp * sqm[1], 0 ); + const cxtype_amp_v fiA_3 = cxmake( tm * (fptype_amp)nsf * sqm[0], 0 ); + const cxtype_amp_v fiA_4 = cxmake( tp * (fptype_amp)nsf * sqm[0], 0 ); + const cxtype_amp_v fiA_5 = cxmake( tm * sqm[1], 0 ); // Branch B: pp != 0. - const fptype_amp sf[2] = { fptype_amp( 1 + nsf + ( 1 - nsf ) * nh ) * (fptype_amp)0.5, - fptype_amp( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype_amp)0.5 }; + const fptype_amp_v sf[2] = { ( (fptype_amp)( 1 + nsf ) + (fptype_amp)( 1 - nsf ) * nh ) * (fptype_amp)0.5, + ( (fptype_amp)( 1 + nsf ) - (fptype_amp)( 1 - nsf ) * nh ) * (fptype_amp)0.5 }; fptype_v omega[2] = { fpsqrt( pvec0 + pp ), 0 }; omega[1] = fmass / omega[0]; - const fptype_amp_v sfomega[2] = { sf[0] * omega[ip], sf[1] * omega[im] }; + const fptype_amp_v sfomega[2] = { sf[0] * ( tm * omega[0] + tp * omega[1] ), sf[1] * ( tp * omega[0] + tm * omega[1] ) }; const fptype_amp_v pp3 = fpmax( fpternary( fpsignbit( pvec3 ), ( pvec1 * pvec1 + pvec2 * pvec2 ) / ( pp - pvec3 ), pp + pvec3 ), 0 ); volatile fptype_amp_v ppDENOM = fpternary( pp != 0, pp, 1. ); // hack: ppDENOM[ieppV]=1 if pp[ieppV]==0 volatile fptype_amp_v pp3DENOM = fpternary( pp3 != 0, pp3, 1. ); // hack: pp3DENOM[ieppV]=1 if pp3[ieppV]==0 @@ -286,11 +291,13 @@ const cxtype_amp_v chi[2] = { cxmake( fpsqrt( chi0r2 ), 0 ), // hack: dummy[ieppV] is not used if pp[ieppV]==0 cxternary( ( pp3 == 0. ), cxmake( -nh, 0 ), - cxmake( (fptype_amp)nh * pvec1, pvec2 ) / fpsqrt( (fptype_amp)2. * ppDENOM * pp3DENOM ) ) }; // hack: dummy[ieppV] is not used if pp[ieppV]==0 - const cxtype_amp_v fiB_2 = sfomega[0] * chi[im]; - const cxtype_amp_v fiB_3 = sfomega[0] * chi[ip]; - const cxtype_amp_v fiB_4 = sfomega[1] * chi[im]; - const cxtype_amp_v fiB_5 = sfomega[1] * chi[ip]; + cxmake( nh * pvec1, pvec2 ) / fpsqrt( (fptype_amp)2. * ppDENOM * pp3DENOM ) ) }; // hack: dummy[ieppV] is not used if pp[ieppV]==0 + const cxtype_amp_v chi_im = tp * chi[0] + tm * chi[1]; // chi[im] + const cxtype_amp_v chi_ip = tm * chi[0] + tp * chi[1]; // chi[ip] + const cxtype_amp_v fiB_2 = sfomega[0] * chi_im; + const cxtype_amp_v fiB_3 = sfomega[0] * chi_ip; + const cxtype_amp_v fiB_4 = sfomega[1] * chi_im; + const cxtype_amp_v fiB_5 = sfomega[1] * chi_ip; // Choose between the results from branch A and branch B const bool_v mask = ( pp == 0. ); w[0] = cxternary( mask, fiA_2, fiB_2 ); @@ -316,8 +323,8 @@ volatile fptype_amp_sv sqp0p3DENOM = fpternary( sqp0p3 != 0, (fptype_sv)sqp0p3, 1. ); // hack: dummy sqp0p3DENOM[ieppV]=1 if sqp0p3[ieppV]==0 cxtype_sv chi[2] = { cxmake( (fptype_v)sqp0p3, 0. ), cxternary( sqp0p3 == 0, - cxmake( -(fptype_amp)nhel * fpsqrt( (fptype_amp)2. * pvec0 ), 0. ), - cxmake( (fptype_amp)nh * pvec1, pvec2 ) / (const fptype_v)sqp0p3DENOM ) }; // hack: dummy[ieppV] is not used if sqp0p3[ieppV]==0 + cxmake( -nhel * fpsqrt( (fptype_amp)2. * pvec0 ), 0. ), + cxmake( nh * pvec1, pvec2 ) / (const fptype_v)sqp0p3DENOM ) }; // hack: dummy[ieppV] is not used if sqp0p3[ieppV]==0 #else // pvec0+pvec3 is a cancelling difference of two ~|p| numbers for a // backward-moving massless fermion, and everything below divides by it. @@ -330,22 +337,12 @@ fptype_sv{ 0 }, fpsqrt( p0p3stable ) * (fptype_amp)nsf ); const cxtype_amp_sv chi[2] = { cxmake( sqp0p3, 0. ), - ( sqp0p3 == 0. ? cxmake( -(fptype_amp)nhel * fpsqrt( (fptype_amp)2. * pvec0 ), 0. ) : cxmake( (fptype_amp)nh * pvec1, pvec2 ) / sqp0p3 ) }; + ( sqp0p3 == 0. ? cxmake( -nhel * fpsqrt( (fptype_amp)2. * pvec0 ), 0. ) : cxmake( nh * pvec1, pvec2 ) / sqp0p3 ) }; #endif - if( nh == 1 ) - { - w[0] = cxzero_sv(); - w[1] = cxzero_sv(); - w[2] = chi[0]; - w[3] = chi[1]; - } - else - { - w[0] = chi[1]; - w[1] = chi[0]; - w[2] = cxzero_sv(); - w[3] = cxzero_sv(); - } + w[0] = tm * chi[1]; + w[1] = tm * chi[0]; + w[2] = tp * chi[0]; + w[3] = tp * chi[1]; } mgDebug( 1, __FUNCTION__ ); return; @@ -359,7 +356,7 @@ __host__ __device__ void ipzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fi, // output: wavefunctions @@ -374,20 +371,17 @@ fi.pvec[2] = fptype_momenta_sv{ 0 }; fi.pvec[3] = -pvec3_ * static_cast(nsf); fi.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; const cxtype_amp_sv sqp0p3 = cxmake( fpsqrt( (fptype_amp)2. * pvec3 ) * (fptype_amp)nsf, 0. ); - w[0] = cxmake( fpamp_of_mom( fi.pvec[1] ), fpamp_of_mom( fi.pvec[2] ) ); - if( nh == 1 ) - { - w[1] = cxmake( fpamp_of_mom( fi.pvec[1] ), fpamp_of_mom( fi.pvec[2] ) ); - w[2] = sqp0p3; - } - else - { - w[1] = sqp0p3; - w[2] = cxmake( fpamp_of_mom( fi.pvec[1] ), fpamp_of_mom( fi.pvec[2] ) ); - } - w[3] = cxmake( fpamp_of_mom( fi.pvec[1] ), fpamp_of_mom( fi.pvec[2] ) ); + w[0] = cxzero_sv(); + w[1] = tm * sqp0p3; + w[2] = tp * sqp0p3; + w[3] = cxzero_sv(); mgDebug( 1, __FUNCTION__ ); return; } @@ -400,7 +394,7 @@ __host__ __device__ void imzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fi, // output: wavefunctions @@ -415,20 +409,17 @@ fi.pvec[2] = fptype_momenta_sv{ 0 }; fi.pvec[3] = -pvec3_ * static_cast(nsf); fi.flv_index = flv; - const int nh = nhel * nsf; - const cxtype_amp_sv chi = cxmake( -(fptype_amp)nhel * fpsqrt( -(fptype_amp)2. * pvec3 ), 0. ); + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; + const cxtype_amp_sv chi = cxmake( -nhel * fpsqrt( -(fptype_amp)2. * pvec3 ), 0. ); + w[0] = tm * chi; w[1] = cxzero_sv(); w[2] = cxzero_sv(); - if( nh == 1 ) - { - w[0] = cxzero_sv(); - w[3] = chi; - } - else - { - w[0] = chi; - w[3] = cxzero_sv(); - } + w[3] = tp * chi; mgDebug( 1, __FUNCTION__ ); return; } @@ -441,7 +432,7 @@ __host__ __device__ void ixzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fi, // output: wavefunctions @@ -463,7 +454,12 @@ fi.pvec[3] = -pvec3_ * static_cast(nsf); fi.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; //const float sqp0p3 = sqrtf( pvec0 + pvec3 ) * nsf; // AV: why force a float here? // pvec0+pvec3 is a cancelling difference of two ~|p| numbers for a // backward-moving massless fermion, and everything below divides by it. @@ -475,21 +471,11 @@ pvec0 + pvec3 ); const fptype_amp_sv sqp0p3 = fpsqrt( p0p3stable ) * static_cast(nsf); const cxtype_amp_sv chi0 = cxmake( sqp0p3, 0. ); - const cxtype_amp_sv chi1 = cxmake( (fptype_amp)nh * pvec1 / sqp0p3, pvec2 / sqp0p3 ); - if( nh == 1 ) - { - w[0] = cxzero_sv(); - w[1] = cxzero_sv(); - w[2] = chi0; - w[3] = chi1; - } - else - { - w[0] = chi1; - w[1] = chi0; - w[2] = cxzero_sv(); - w[3] = cxzero_sv(); - } + const cxtype_amp_sv chi1 = cxmake( nh * pvec1 / sqp0p3, pvec2 / sqp0p3 ); + w[0] = tm * chi1; + w[1] = tm * chi0; + w[2] = tp * chi0; + w[3] = tp * chi1; mgDebug( 1, __FUNCTION__ ); return; } @@ -501,7 +487,7 @@ __host__ __device__ void vxxxxx( const fptype_momenta momenta[], // input: momenta const fptype_amp vmass, // input: vector boson mass - const int nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson) + const fptype_amp_sv nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson, per SIMD lane/GPU thread) const int nsv, // input: +1 (final) or -1 (initial) const int flv, // input: flavour ALOHAOBJ & vc, // output: wavefunctions @@ -527,7 +513,10 @@ vc.pvec[3] = pvec3_ * static_cast(nsv); vc.flv_index = flv; const fptype_amp sqh = fpsqrt( 0.5 ); // AV this is > 0! - const fptype_amp hel = nhel; + const fptype_amp_sv hel = nhel; + // BRANCHLESS HELICITY: |hel| is hel * hel for hel in { -1, 0, +1 }, so that the helicity can + // differ from one SIMD lane or GPU thread to the next with no branch and no std::abs + const fptype_amp_sv abshel = hel * hel; // FD gauge const cxtype_amp_sv cI = cxmake( 0. + fptype_amp_sv{ 0. }, 1. + fptype_amp_sv{ 0. } ); @@ -538,8 +527,8 @@ if( vmass != 0. ) { - const int nsvahl = nsv * std::abs( hel ); - const fptype_amp hel0 = 1. - std::abs( hel ); + const fptype_amp_sv nsvahl = static_cast( nsv ) * abshel; + const fptype_amp_sv hel0 = 1 - abshel; // truth value of hel == 0 (longitudinal polarisation) #ifndef MGONGPU_CPPSIMD const fptype_amp_sv pt2 = ( pvec1 * pvec1 ) + ( pvec2 * pvec2 ); const fptype_amp_sv pp = fpmin( pvec0, fpsqrt( pt2 + ( pvec3 * pvec3 ) ) ); @@ -594,17 +583,14 @@ nk = n[0]*pvec0 - n[1]*pvec1 - n[2]*pvec2 - n[3]*pvec3; - if ( abs(nhel) == 1) - { - w[4] = cxzero_sv(); - } - else{ - w[0] = cxmake( -vmass/nk * n[0], zero ); - w[1] = cxmake( -vmass/nk * n[1], zero ); - w[2] = cxmake( -vmass/nk * n[2], zero ); - w[3] = cxmake( -vmass/nk * n[3], zero ); - w[4] = -static_cast(nsv)*cI; - } + // BRANCHLESS HELICITY: |hel| = 1 keeps w[0..3] with w[4] = 0, hel = 0 (hel0 = 1) takes the + // FD-gauge longitudinal polarisation instead (nk > 0 for a massive momentum, so the unused + // alternative is finite and vanishes when multiplied by its zero truth value) + w[0] = abshel * w[0] + hel0 * cxmake( -vmass / nk * n[0], zero ); + w[1] = abshel * w[1] + hel0 * cxmake( -vmass / nk * n[1], zero ); + w[2] = abshel * w[2] + hel0 * cxmake( -vmass / nk * n[2], zero ); + w[3] = abshel * w[3] + hel0 * cxmake( -vmass / nk * n[3], zero ); + w[4] = hel0 * ( -static_cast( nsv ) * cI ); #else volatile fptype_amp_sv pt2 = ( pvec1 * pvec1 ) + ( pvec2 * pvec2 ); @@ -613,9 +599,9 @@ const fptype_amp_sv pt = fpmin( pp, fpsqrt( pt2 ) ); // Branch A: pp == 0. const cxtype_amp vcA_2 = cxmake( 0, 0 ); - const cxtype_amp vcA_3 = cxmake( -hel * sqh, 0 ); - const cxtype_amp vcA_4 = cxmake( 0, nsvahl * sqh ); - const cxtype_amp vcA_5 = cxmake( hel0, 0 ); + const cxtype_amp_v vcA_3 = cxmake( -hel * sqh, 0 ); + const cxtype_amp_v vcA_4 = cxmake( 0., nsvahl * sqh ); + const cxtype_amp_v vcA_5 = cxmake( hel0, 0 ); // Branch B: pp != 0. volatile fptype_amp_v ppDENOM = fpternary( pp != 0, pp, 1. ); // hack: ppDENOM[ieppV]=1 if pp[ieppV]==0 const fptype_amp_v emp = pvec0 / ( vmass * ppDENOM ); // hack: dummy[ieppV] is not used if pp[ieppV]==0 @@ -624,11 +610,11 @@ // Branch B1: pp != 0. and pt != 0. volatile fptype_amp_v ptDENOM = fpternary( pt != 0, pt, 1. ); // hack: ptDENOM[ieppV]=1 if pt[ieppV]==0 const fptype_amp_v pzpt = pvec3 / ( ppDENOM * ptDENOM ) * sqh * hel; // hack: dummy[ieppV] is not used if pp[ieppV]==0 - const cxtype_amp_v vcB1_3 = cxmake( hel0 * pvec1 * emp - pvec1 * pzpt, -fpamp_scalar(nsvahl) * pvec2 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 - const cxtype_amp_v vcB1_4 = cxmake( hel0 * pvec2 * emp - pvec2 * pzpt, fpamp_scalar(nsvahl) * pvec1 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 + const cxtype_amp_v vcB1_3 = cxmake( hel0 * pvec1 * emp - pvec1 * pzpt, -nsvahl * pvec2 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 + const cxtype_amp_v vcB1_4 = cxmake( hel0 * pvec2 * emp - pvec2 * pzpt, nsvahl * pvec1 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 // Branch B2: pp != 0. and pt == 0. - const cxtype_amp vcB2_3 = cxmake( -hel * sqh, 0. ); - const cxtype_amp_v vcB2_4 = cxmake( 0., fpamp_scalar(nsvahl) * fpternary( ( pvec3 < 0 ), -sqh, sqh ) ); // AV: removed an abs here + const cxtype_amp_v vcB2_3 = cxmake( -hel * sqh, 0. ); + const cxtype_amp_v vcB2_4 = cxmake( 0., nsvahl * fpternary( ( pvec3 < 0 ), -sqh, sqh ) ); // AV: removed an abs here // Choose between the results from branch A and branch B (and from branch B1 and branch B2) const bool_v mask = ( pp == 0. ); const bool_v maskB = ( pt != 0. ); @@ -650,20 +636,14 @@ nk = n[0]*pvec0 - n[1]*pvec1 - n[2]*pvec2 - n[3]*pvec3; - // nhel is a scalar: no need for a per-lane mask here (and a bool_v built - // from a single value would only set the first lane) - if ( abs(nhel) == 1 ) - { - w[4] = zero; - } - else - { - w[0] = cxmake( -vmass/nk * n[0], zero ); - w[1] = cxmake( -vmass/nk * n[1], zero ); - w[2] = cxmake( -vmass/nk * n[2], zero ); - w[3] = cxmake( -vmass/nk * n[3], zero ); - w[4] = -static_cast(nsv)*cI; - } + // BRANCHLESS HELICITY: |hel| = 1 keeps w[0..3] with w[4] = 0, hel = 0 (hel0 = 1) takes the + // FD-gauge longitudinal polarisation instead (nk > 0 for a massive momentum, so the unused + // alternative is finite and vanishes when multiplied by its zero truth value) + w[0] = abshel * w[0] + hel0 * cxmake( -vmass / nk * n[0], zero ); + w[1] = abshel * w[1] + hel0 * cxmake( -vmass / nk * n[1], zero ); + w[2] = abshel * w[2] + hel0 * cxmake( -vmass / nk * n[2], zero ); + w[3] = abshel * w[3] + hel0 * cxmake( -vmass / nk * n[3], zero ); + w[4] = hel0 * ( -static_cast( nsv ) * cI ); #endif } else @@ -695,11 +675,11 @@ // Branch A: pt != 0. volatile fptype_amp_v ptDENOM = fpternary( pt != 0, pt, 1. ); // hack: ptDENOM[ieppV]=1 if pt[ieppV]==0 const fptype_amp_v pzpt = pvec3 / ( pp * ptDENOM ) * sqh * hel; // hack: dummy[ieppV] is not used if pt[ieppV]==0 - const cxtype_amp_v vcA_3 = cxmake( -pvec1 * pzpt, -fpamp_scalar(nsv) * pvec2 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 - const cxtype_amp_v vcA_4 = cxmake( -pvec2 * pzpt, fpamp_scalar(nsv) * pvec1 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 + const cxtype_amp_v vcA_3 = cxmake( -pvec1 * pzpt, -(fptype_amp)nsv * pvec2 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 + const cxtype_amp_v vcA_4 = cxmake( -pvec2 * pzpt, (fptype_amp)nsv * pvec1 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 // Branch B: pt == 0. - const cxtype_amp vcB_3 = cxmake( -(fptype_amp)hel * sqh, 0 ); - const cxtype_amp_v vcB_4 = cxmake( 0, fpamp_scalar(nsv) * fpternary( ( pvec3 < 0 ), -sqh, sqh ) ); // AV: removed an abs here + const cxtype_amp_v vcB_3 = cxmake( -hel * sqh, 0 ); + const cxtype_amp_v vcB_4 = cxmake( 0, (fptype_amp)nsv * fpternary( ( pvec3 < 0 ), -sqh, sqh ) ); // AV: removed an abs here // Choose between the results from branch A and branch B const bool_v mask = ( pt != 0. ); w[1] = cxternary( mask, vcA_3, vcB_3 ); @@ -758,7 +738,7 @@ __host__ __device__ void oxxxxx( const fptype_momenta momenta[], // input: momenta const fptype_amp fmass, // input: fermion mass - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) int flv, // input: flavour ALOHAOBJ & fo, // output: wavefunctions @@ -784,7 +764,12 @@ fo.pvec[2] = pvec2_ * static_cast(nsf); fo.pvec[3] = pvec3_ * static_cast(nsf); fo.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; if( fmass != 0. ) { #ifndef MGONGPU_CPPSIMD @@ -795,12 +780,12 @@ fptype_amp sqm[2] = { fpsqrt( std::abs( fmass ) ), 0. }; // possibility of negative fermion masses //sqm[1] = ( fmass < 0. ? -abs( sqm[0] ) : abs( sqm[0] ) ); // AV: why abs here? sqm[1] = ( fmass < 0. ? -sqm[0] : sqm[0] ); // AV: removed an abs here - const int ip = -( ( 1 - nh ) / 2 ) * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ - const int im = ( 1 + nh ) / 2 * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ - w[0] = cxmake( im * sqm[std::abs( ip )], 0 ); - w[1] = cxmake( ip * nsf * sqm[std::abs( ip )], 0 ); - w[2] = cxmake( im * nsf * sqm[std::abs( im )], 0 ); - w[3] = cxmake( ip * sqm[std::abs( im )], 0 ); + // The Fortran ip = -( ( 1 - nh ) / 2 ) * nhel and im = ( 1 + nh ) / 2 * nhel are -tm * nhel and + // tp * nhel, so that |ip| = tm and |im| = tp: sqm[|ip|] is sqm[1] wherever ip != 0, and so on + w[0] = cxmake( tp * nhel * sqm[0], 0 ); + w[1] = cxmake( -tm * nhel * nsf * sqm[1], 0 ); + w[2] = cxmake( tp * nhel * nsf * sqm[1], 0 ); + w[3] = cxmake( -tm * nhel * sqm[0], 0 ); } else { @@ -808,17 +793,18 @@ fptype_amp( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype_amp)0.5 }; fptype_amp omega[2] = { fpsqrt( pvec0 + pp ), 0. }; omega[1] = fmass / omega[0]; - const int ip = ( 1 + nh ) / 2; // NB: Fortran is (3+nh)/2 because omega(2) has indexes 1,2 and not 0,1 - const int im = ( 1 - nh ) / 2; // NB: Fortran is (3-nh)/2 because omega(2) has indexes 1,2 and not 0,1 - const fptype_amp sfomeg[2] = { sf[0] * omega[ip], sf[1] * omega[im] }; + // The Fortran ip/im indices ( 3 +- nh ) / 2 (omega(2) has indexes 1,2) are the truth values tp/tm + const fptype_amp sfomeg[2] = { sf[0] * ( tm * omega[0] + tp * omega[1] ), sf[1] * ( tp * omega[0] + tm * omega[1] ) }; const fptype_amp pp3 = fpmax( fpternary( fpsignbit( pvec3 ), ( pvec1 * pvec1 + pvec2 * pvec2 ) / ( pp - pvec3 ), pp + pvec3 ), 0. ); const cxtype_amp chi[2] = { cxmake( fpsqrt( pp3 * (fptype_amp)0.5 / pp ), 0. ), ( ( pp3 == 0. ) ? cxmake( -nh, 0. ) : cxmake( nh * pvec1, -pvec2 ) / fpsqrt( (fptype_amp)2. * pp * pp3 ) ) }; - w[0] = sfomeg[1] * chi[im]; - w[1] = sfomeg[1] * chi[ip]; - w[2] = sfomeg[0] * chi[im]; - w[3] = sfomeg[0] * chi[ip]; + const cxtype_amp chi_im = tp * chi[0] + tm * chi[1]; // chi[im] + const cxtype_amp chi_ip = tm * chi[0] + tp * chi[1]; // chi[ip] + w[0] = sfomeg[1] * chi_im; + w[1] = sfomeg[1] * chi_ip; + w[2] = sfomeg[0] * chi_im; + w[3] = sfomeg[0] * chi_ip; } #else volatile fptype_amp_sv p2 = pvec1 * pvec1 + pvec2 * pvec2 + pvec3 * pvec3; // volatile fixes #736 @@ -827,20 +813,17 @@ // NB: Do not use "abs" for floats! It returns an integer with no build warning! Use std::abs! fptype_amp sqm[2] = { fpsqrt( std::abs( fmass ) ), 0 }; // possibility of negative fermion masses sqm[1] = ( fmass < 0 ? -sqm[0] : sqm[0] ); // AV: removed an abs here (as above) - const int ipA = -( ( 1 - nh ) / 2 ) * nhel; - const int imA = ( 1 + nh ) / 2 * nhel; - const cxtype_amp foA_2 = imA * sqm[std::abs( ipA )]; - const cxtype_amp foA_3 = ipA * nsf * sqm[std::abs( ipA )]; - const cxtype_amp foA_4 = imA * nsf * sqm[std::abs( imA )]; - const cxtype_amp foA_5 = ipA * sqm[std::abs( imA )]; + // (see the scalar code above for the truth values tp/tm in place of the Fortran ip/im indices) + const cxtype_amp_v foA_2 = cxmake( tp * nhel * sqm[0], 0 ); + const cxtype_amp_v foA_3 = cxmake( -tm * nhel * (fptype_amp)nsf * sqm[1], 0 ); + const cxtype_amp_v foA_4 = cxmake( tp * nhel * (fptype_amp)nsf * sqm[1], 0 ); + const cxtype_amp_v foA_5 = cxmake( -tm * nhel * sqm[0], 0 ); // Branch B: pp != 0. - const fptype_amp sf[2] = { fptype_amp( 1 + nsf + ( 1 - nsf ) * nh ) * (fptype_amp)0.5, - fptype_amp( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype_amp)0.5 }; + const fptype_amp_v sf[2] = { ( (fptype_amp)( 1 + nsf ) + (fptype_amp)( 1 - nsf ) * nh ) * (fptype_amp)0.5, + ( (fptype_amp)( 1 + nsf ) - (fptype_amp)( 1 - nsf ) * nh ) * (fptype_amp)0.5 }; fptype_v omega[2] = { fpsqrt( pvec0 + pp ), 0 }; omega[1] = fmass / omega[0]; - const int ipB = ( 1 + nh ) / 2; - const int imB = ( 1 - nh ) / 2; - const fptype_amp_v sfomeg[2] = { sf[0] * omega[ipB], sf[1] * omega[imB] }; + const fptype_amp_v sfomeg[2] = { sf[0] * ( tm * omega[0] + tp * omega[1] ), sf[1] * ( tp * omega[0] + tm * omega[1] ) }; const fptype_amp_v pp3 = fpmax( fpternary( fpsignbit( pvec3 ), ( pvec1 * pvec1 + pvec2 * pvec2 ) / ( pp - pvec3 ), pp + pvec3 ), 0. ); volatile fptype_amp_v ppDENOM = fpternary( pp != 0, pp, 1. ); // hack: ppDENOM[ieppV]=1 if pp[ieppV]==0 volatile fptype_amp_v pp3DENOM = fpternary( pp3 != 0, pp3, 1. ); // hack: pp3DENOM[ieppV]=1 if pp3[ieppV]==0 @@ -848,11 +831,13 @@ const cxtype_amp_v chi[2] = { cxmake( fpsqrt( chi0r2 ), 0. ), // hack: dummy[ieppV] is not used if pp[ieppV]==0 ( cxternary( ( pp3 == 0. ), cxmake( -nh, 0. ), - cxmake( (fptype_amp)nh * pvec1, -pvec2 ) / fpsqrt( (fptype_amp)2. * ppDENOM * pp3DENOM ) ) ) }; // hack: dummy[ieppV] is not used if pp[ieppV]==0 - const cxtype_amp_v foB_2 = sfomeg[1] * chi[imB]; - const cxtype_amp_v foB_3 = sfomeg[1] * chi[ipB]; - const cxtype_amp_v foB_4 = sfomeg[0] * chi[imB]; - const cxtype_amp_v foB_5 = sfomeg[0] * chi[ipB]; + cxmake( nh * pvec1, -pvec2 ) / fpsqrt( (fptype_amp)2. * ppDENOM * pp3DENOM ) ) ) }; // hack: dummy[ieppV] is not used if pp[ieppV]==0 + const cxtype_amp_v chi_im = tp * chi[0] + tm * chi[1]; // chi[im] + const cxtype_amp_v chi_ip = tm * chi[0] + tp * chi[1]; // chi[ip] + const cxtype_amp_v foB_2 = sfomeg[1] * chi_im; + const cxtype_amp_v foB_3 = sfomeg[1] * chi_ip; + const cxtype_amp_v foB_4 = sfomeg[0] * chi_im; + const cxtype_amp_v foB_5 = sfomeg[0] * chi_ip; // Choose between the results from branch A and branch B const bool_v mask = ( pp == 0. ); w[0] = cxternary( mask, foA_2, foB_2 ); @@ -879,7 +864,7 @@ const cxtype_amp_v chi[2] = { cxmake( (fptype_v)sqp0p3, 0. ), cxternary( ( sqp0p3 == 0. ), cxmake( -nhel, 0. ) * fpsqrt( (fptype_amp)2. * pvec0 ), - cxmake( (fptype_amp)nh * pvec1, -pvec2 ) / (const fptype_amp_sv)sqp0p3DENOM ) }; // hack: dummy[ieppV] is not used if sqp0p3[ieppV]==0 + cxmake( nh * pvec1, -pvec2 ) / (const fptype_amp_sv)sqp0p3DENOM ) }; // hack: dummy[ieppV] is not used if sqp0p3[ieppV]==0 #else // pvec0+pvec3 is a cancelling difference of two ~|p| numbers for a // backward-moving massless fermion, and everything below divides by it. @@ -892,22 +877,12 @@ fptype_sv{ 0 }, fpsqrt( p0p3stable ) * (fptype_amp)nsf ); const cxtype_amp_sv chi[2] = { cxmake( sqp0p3, 0. ), - ( sqp0p3 == 0. ? cxmake( -nhel, 0. ) * fpsqrt( (fptype_amp)2. * pvec0 ) : cxmake( (fptype_amp)nh * pvec1, -pvec2 ) / sqp0p3 ) }; + ( sqp0p3 == 0. ? cxmake( -nhel, 0. ) * fpsqrt( (fptype_amp)2. * pvec0 ) : cxmake( nh * pvec1, -pvec2 ) / sqp0p3 ) }; #endif - if( nh == 1 ) - { - w[0] = chi[0]; - w[1] = chi[1]; - w[2] = cxzero_sv(); - w[3] = cxzero_sv(); - } - else - { - w[0] = cxzero_sv(); - w[1] = cxzero_sv(); - w[2] = chi[1]; - w[3] = chi[0]; - } + w[0] = tp * chi[0]; + w[1] = tp * chi[1]; + w[2] = tm * chi[1]; + w[3] = tm * chi[0]; } mgDebug( 1, __FUNCTION__ ); return; @@ -921,7 +896,7 @@ __host__ __device__ void opzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fo, // output: wavefunctions @@ -936,20 +911,17 @@ fo.pvec[2] = fptype_momenta_sv{ 0 }; fo.pvec[3] = pvec3_ * static_cast(nsf); fo.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; const cxtype_amp_sv csqp0p3 = cxmake( fpsqrt( (fptype_amp)2. * pvec3 ) * (fptype_amp)nsf, 0. ); + w[0] = tp * csqp0p3; w[1] = cxzero_sv(); w[2] = cxzero_sv(); - if( nh == 1 ) - { - w[0] = csqp0p3; - w[3] = cxzero_sv(); - } - else - { - w[0] = cxzero_sv(); - w[3] = csqp0p3; - } + w[3] = tm * csqp0p3; mgDebug( 1, __FUNCTION__ ); return; } @@ -962,7 +934,7 @@ __host__ __device__ void omzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fo, // output: wavefunctions @@ -977,23 +949,17 @@ fo.pvec[2] = fptype_momenta_sv{ 0 }; fo.pvec[3] = pvec3_ * static_cast(nsf); fo.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; const cxtype_amp_sv chi1 = cxmake( -nhel, 0. ) * fpsqrt( -(fptype_amp)2. * pvec3 ); - if( nh == 1 ) - { - w[0] = cxzero_sv(); - w[1] = chi1; - w[2] = cxzero_sv(); - w[3] = cxzero_sv(); - } - else - { - w[0] = cxzero_sv(); - w[1] = cxzero_sv(); - w[2] = chi1; - //w[3] = chi1; // AV: BUG! - w[3] = cxzero_sv(); // AV: BUG FIX - } + w[0] = cxzero_sv(); + w[1] = tp * chi1; + w[2] = tm * chi1; + w[3] = cxzero_sv(); mgDebug( 1, __FUNCTION__ ); return; } @@ -1006,7 +972,7 @@ __host__ __device__ void oxzxxx( const fptype_momenta momenta[], // input: momenta //const fptype_amp fmass, // [skip: ASSUME fermion mass==0] - const int nhel, // input: -1 or +1 (helicity of fermion) + const fptype_amp_sv nhel, // input: -1 or +1 (helicity of fermion, per SIMD lane/GPU thread) const int nsf, // input: +1 (particle) or -1 (antiparticle) const int flv, // input: flavour ALOHAOBJ & fo, // output: wavefunctions @@ -1028,7 +994,12 @@ fo.pvec[2] = pvec2_ * static_cast(nsf); fo.pvec[3] = pvec3_ * static_cast(nsf); fo.flv_index = flv; - const int nh = nhel * nsf; + const fptype_amp_sv nh = nhel * static_cast( nsf ); + // BRANCHLESS HELICITY: the helicity is not a branch condition but enters as the truth values + // tp (nh == +1) and tm (nh == -1), 1 or 0, which multiply the two alternatives. Different + // SIMD lanes or GPU threads can then carry different helicities and still run in lock-step. + const fptype_amp_sv tp = ( 1 + nh ) * (fptype_amp)0.5; + const fptype_amp_sv tm = 1 - tp; //const float sqp0p3 = sqrtf( pvec0 + pvec3 ) * nsf; // AV: why force a float here? // pvec0+pvec3 is a cancelling difference of two ~|p| numbers for a // backward-moving massless fermion, and everything below divides by it. @@ -1040,21 +1011,11 @@ pvec0 + pvec3 ); const fptype_amp_sv sqp0p3 = fpsqrt( p0p3stable ) * static_cast(nsf); const cxtype_amp_sv chi0 = cxmake( sqp0p3, 0. ); - const cxtype_amp_sv chi1 = cxmake( (fptype_amp)nh * pvec1 / sqp0p3, -pvec2 / sqp0p3 ); - if( nh == 1 ) - { - w[0] = chi0; - w[1] = chi1; - w[2] = cxzero_sv(); - w[3] = cxzero_sv(); - } - else - { - w[0] = cxzero_sv(); - w[1] = cxzero_sv(); - w[2] = chi1; - w[3] = chi0; - } + const cxtype_amp_sv chi1 = cxmake( nh * pvec1 / sqp0p3, -pvec2 / sqp0p3 ); + w[0] = tp * chi0; + w[1] = tp * chi1; + w[2] = tm * chi1; + w[3] = tm * chi0; mgDebug( 1, __FUNCTION__ ); return; } diff --git a/docs/source/madspace/examples/flow-training.rst b/docs/source/madspace/examples/flow-training.rst index 43a58c67ff..fef6814a72 100644 --- a/docs/source/madspace/examples/flow-training.rst +++ b/docs/source/madspace/examples/flow-training.rst @@ -40,7 +40,7 @@ The mapping, cuts, matrix element, PDF and scale are exactly the same building b mapping = ms.PhaseSpaceMapping(masses, E_CM, mode="propagator", cuts=cuts) proc_dir = "PROC_ggttg" - me_path = os.path.join(proc_dir, "lib", "libmadmatrix_P0_gg_ttxg_scalar.so") + me_path = os.path.join(proc_dir, "lib", "libmadmatrix_sm_P0_gg_ttxg_scalar_m.so") param_card = os.path.join(proc_dir, "Cards", "param_card.dat") ctx = ms.default_context() diff --git a/docs/source/madspace/examples/integrator-madnis.rst b/docs/source/madspace/examples/integrator-madnis.rst index f5e8306025..3f1a329697 100644 --- a/docs/source/madspace/examples/integrator-madnis.rst +++ b/docs/source/madspace/examples/integrator-madnis.rst @@ -38,7 +38,7 @@ random numbers are now PyTorch tensors, since MadNIS trains with PyTorch: mapping = ms.PhaseSpaceMapping(masses, E_CM, mode="propagator", cuts=cuts) proc_dir = "PROC_ggttg" - me_path = os.path.join(proc_dir, "lib", "libmadmatrix_P0_gg_ttxg_scalar.so") + me_path = os.path.join(proc_dir, "lib", "libmadmatrix_sm_P0_gg_ttxg_scalar_m.so") param_card = os.path.join(proc_dir, "Cards", "param_card.dat") ctx = ms.default_context() diff --git a/docs/source/madspace/examples/integrator.rst b/docs/source/madspace/examples/integrator.rst index 4c135b12d7..0b91de83c0 100644 --- a/docs/source/madspace/examples/integrator.rst +++ b/docs/source/madspace/examples/integrator.rst @@ -42,7 +42,7 @@ examples, generated with ``output mg7 PROC_ggttg`` and ``make BACKEND=scalar`` b .. code-block:: python proc_dir = "PROC_ggttg" - me_path = os.path.join(proc_dir, "lib", "libmadmatrix_P0_gg_ttxg_scalar.so") + me_path = os.path.join(proc_dir, "lib", "libmadmatrix_sm_P0_gg_ttxg_scalar_m.so") param_card = os.path.join(proc_dir, "Cards", "param_card.dat") ctx = ms.default_context() diff --git a/docs/source/madspace/examples/matrix-element.rst b/docs/source/madspace/examples/matrix-element.rst index b42d06ddce..c4a53d92ef 100644 --- a/docs/source/madspace/examples/matrix-element.rst +++ b/docs/source/madspace/examples/matrix-element.rst @@ -37,7 +37,7 @@ Loading the matrix element import madspace as ms proc_dir = "PROC_ggttg" - me_path = os.path.join(proc_dir, "lib", "libmadmatrix_P0_gg_ttxg_scalar.so") + me_path = os.path.join(proc_dir, "lib", "libmadmatrix_sm_P0_gg_ttxg_scalar_m.so") param_card = os.path.join(proc_dir, "Cards", "param_card.dat") ctx = ms.default_context() @@ -46,7 +46,9 @@ Loading the matrix element The library name also appears, alongside every other subprocess of the process, in ``PROC_ggttg/SubProcesses/subprocesses.json``, as ``me_path`` with a ``{device}`` placeholder -for the backend. +for the backend and an ``{fptype}`` placeholder for the floating-point precision the library +was compiled with (``FPTYPE``, ``m`` by default). The model name (``sm``) is part of the +library name too. :py:meth:`Context.load_matrix_element ` returns a :py:class:`MatrixElementApi `, which reads the parameters from diff --git a/madgraph/interface/common_run_interface.py b/madgraph/interface/common_run_interface.py index 0afa5de747..8bcac2e38f 100755 --- a/madgraph/interface/common_run_interface.py +++ b/madgraph/interface/common_run_interface.py @@ -2181,6 +2181,20 @@ def check_multicore(self): raise self.ConfigurationError('''Can\'t load Reweight module. The variable mg5_path might not be correctly configured.''') + # the reweighting backend (madgraph/various/reweight_backend.py): this + # is the classic module; MadtRex runs from the MG7 output's launcher + backend = [a.split('=', 1)[1] for a in line.split() + if a.startswith(('--backend=', '--reweight-backend='))] + if backend: + line = ' '.join(a for a in line.split() + if not a.startswith(('--backend=', '--reweight-backend='))) + if backend[-1].lower() == 'madtrex': + raise self.InvalidCmd('--backend=madtrex: MadtRex needs an MG7 output; ' + 'this reweighting uses the classic module') + if backend[-1].lower() not in ('classic', 'auto'): + raise self.InvalidCmd('--backend=%s: expected classic, madtrex or auto' + % backend[-1]) + # plugin option plugin = False reweight_mode = None diff --git a/madgraph/interface/madevent_interface.py b/madgraph/interface/madevent_interface.py index 3d119549aa..f682aff010 100755 --- a/madgraph/interface/madevent_interface.py +++ b/madgraph/interface/madevent_interface.py @@ -523,8 +523,8 @@ class AskRun(cmd.ControlSwitch): def __init__(self, question, line_args=[], mode=None, force=False, *args, **opt): - self.check_available_module(opt['mother_interface'].options) self.me_dir = opt['mother_interface'].me_dir + self.check_available_module(opt['mother_interface'].options) super(AskRun,self).__init__(self.to_control, opt['mother_interface'], *args, **opt) @@ -557,6 +557,18 @@ def check_available_module(self, options): if misc.has_f2py() or options['f2py_compiler']: self.available_module.add('reweight') + # An MG7 output reweights with the classic module or with MadtRex + # (see madgraph/various/reweight_backend.py): reweighting is + # available when either one can run. + self.reweight_env = None + if not MADEVENT and os.path.exists(pjoin(self.me_dir, 'Cards', 'run_card.toml')): + from madgraph.various import reweight_backend + self.reweight_env = reweight_backend.probe_environment( + self.me_dir, dict(options, mg5_path=options.get('mg5_path') or MG5DIR), + has_f2py='reweight' in self.available_module) + if reweight_backend.available_backends(self.reweight_env): + self.available_module.add('reweight') + # old mode to activate the shower def ans_parton(self, value=None): """None: means that the user type 'pythia' @@ -907,7 +919,16 @@ def get_allowed_reweight(self): if 'reweight' not in self.available_module: self.allowed_reweight = [] return self.allowed_reweight - self.allowed_reweight = ['OFF', 'ON', 'density'] + if self.reweight_env is not None: + # MG7 output: ON is auto; then the backends that can run here + from madgraph.various import reweight_backend + backends = reweight_backend.available_backends(self.reweight_env) + self.allowed_reweight = ['OFF', 'ON'] + backends + if 'classic' not in backends: + return self.allowed_reweight + self.allowed_reweight.append('density') + else: + self.allowed_reweight = ['OFF', 'ON', 'density'] # check for plugin mode plugin_path = self.mother_interface.plugin_path @@ -920,13 +941,7 @@ def set_default_reweight(self): if 'reweight' in self.available_module: if os.path.exists(pjoin(self.me_dir,'Cards','reweight_card.dat')): - reweightcard = open(pjoin(self.me_dir,'Cards','reweight_card.dat'), 'r') - content = reweightcard.read() - if 'change particle_in_density_matrix' in content: - self.switch['reweight'] = 'density' - else: - self.switch['reweight'] = 'ON' - reweightcard.close() + self.switch['reweight'] = self.switch_value_from_card_reweight() else: self.switch['reweight'] = 'OFF' else: @@ -934,6 +949,9 @@ def set_default_reweight(self): def get_cardcmd_for_reweight(self, value): """set some command to run before allowing the user to modify the cards.""" + if value in ['classic', 'madtrex']: + # the two backends of an MG7 output read the same card as ON + value = 'ON' if value in ['density']: content_rwgt_card = open(pjoin(self.me_dir, "Cards", "reweight_card.dat"), "r") if 'change particle_in_density_matrix' in content_rwgt_card.read(): @@ -961,14 +979,26 @@ def get_cardcmd_for_reweight(self, value): def switch_value_from_card_reweight(self): """re-derive the reweight switch (ON vs density) from the content of - reweight_card.dat (called after the user has edited the cards).""" + reweight_card.dat (called after the user has edited the cards). For + an MG7 output, a "change backend classic|madtrex" line of the card + gives that backend (when it can run here).""" path = pjoin(self.me_dir, 'Cards', 'reweight_card.dat') if not os.path.exists(path): return None with open(path) as fsock: - if 'change particle_in_density_matrix' in fsock.read(): - return 'density' + content = fsock.read() + if 'change particle_in_density_matrix' in content: + return 'density' + if getattr(self, 'reweight_env', None) is not None: + # a backend picked in the question outranks the card's line + # (the run warns when the two disagree) + current = getattr(self, 'switch', {}).get('reweight') + if current in ('classic', 'madtrex'): + return current + match = re.search(r'^\s*change\s+backend\s+(\w+)', content, re.M | re.I) + if match and match.group(1).lower() in self.get_allowed_reweight(): + return match.group(1).lower() return 'ON' diff --git a/madgraph/interface/madgraph_interface.py b/madgraph/interface/madgraph_interface.py index a67e23d574..a469441513 100755 --- a/madgraph/interface/madgraph_interface.py +++ b/madgraph/interface/madgraph_interface.py @@ -623,6 +623,14 @@ def help_install(self): logger.info(" --no-debug Optimized build, no debug symbols, Release (source build, default).") logger.info(" --cuda-arch=ARCHS Semicolon-separated CUDA compute capabilities, e.g. 75;80;86.") logger.info(" --hip-arch=ARCHS Semicolon-separated HIP architectures, e.g. gfx900;gfx906.") + logger.info(" ") + logger.info(" \"install madtrex [options]\"",'$MG:BOLD') + logger.info(" Install MadtRex, the event-reweighting library used by MG7, together") + logger.info(" with its bundled Rex fork. Requires madspace to be installed first.") + logger.info(" Without options an interactive installer is launched. Available options:") + logger.info(" -y/--yes Install non-interactively.") + logger.info(" -j N/--jobs=N Parallel compilation jobs (defaults to the 'nb_core' option).") + logger.info(" --clean Rebuild from scratch.") def help_display(self): logger.info("syntax: display " + "|".join(self._display_opts),'$MG:color:BLUE') @@ -3749,6 +3757,12 @@ def complete_install(self, text, line, begidx, endidx): if any(a.startswith(opt) for a in args[2:]): options.remove(opt) return self.list_completion(text, options) + elif len(args) >= 2 and args[1] == 'madtrex': + options = ['-y', '--yes', '--clean', '--jobs='] + for opt in options[:]: + if any(a.startswith(opt) for a in args[2:]): + options.remove(opt) + return self.list_completion(text, options) elif len(args)>=2 and args[1] in self._advanced_install_opts: options = ['--keep_source','--logging='] if args[1]=='pythia8': @@ -3804,7 +3818,7 @@ def _tutorial_opts(self): _install_opts = ['Delphes', 'ExRootAnalysis', 'update', 'Golem95', 'QCDLoop', 'maddm', 'maddump', 'looptools', 'cuttools_iregi', 'MadSTR', 'RunningCoupling', - 'madspace'] + 'madspace', 'madtrex'] # The targets below are installed using the HEPToolsInstaller.py script _advanced_install_opts = ['pythia8','zlib','boost','lhapdf6','lhapdf5','collier', @@ -8571,8 +8585,8 @@ def do_install(self, line, paths=None, additional_options=[]): elif args[0] == 'cuttools_iregi': self.install_cuttools_iregi() return - elif args[0] == 'madspace': - install_script = pjoin(MG5DIR, 'madspace', 'install.py') + elif args[0] in ('madspace', 'madtrex'): + install_script = pjoin(MG5DIR, args[0], 'install.py') install_args = args[1:] # The source build compiles in parallel only if it is told how many # jobs it may use (see set_build_parallelism in install.py), so hand @@ -8580,7 +8594,11 @@ def do_install(self, line, paths=None, additional_options=[]): if not any(a == '-j' or a.startswith('-j') or a.startswith('--jobs') for a in install_args): install_args = install_args + ['-j', str(self.options['nb_core'])] - subprocess.run([sys.executable, install_script] + install_args) + result = subprocess.run([sys.executable, install_script] + install_args) + if result.returncode != 0: + raise self.InvalidCmd('install %s failed (%s exited with code %d): see ' + 'its output above' % (args[0], install_script, + result.returncode)) return plugin = self.install_plugin diff --git a/madgraph/interface/reweight_interface.py b/madgraph/interface/reweight_interface.py index a949b0e0c6..961b5d5062 100755 --- a/madgraph/interface/reweight_interface.py +++ b/madgraph/interface/reweight_interface.py @@ -62,6 +62,28 @@ #lhapdf = None +# The "change