diff --git a/epoch1/oneapi/gg_ttgg/Cards/param_card.dat b/epoch1/oneapi/gg_ttgg/Cards/param_card.dat new file mode 100644 index 0000000000..637294342d --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/Cards/param_card.dat @@ -0,0 +1,78 @@ +###################################################################### +## PARAM_CARD AUTOMATICALY GENERATED BY MG5 FOLLOWING UFO MODEL #### +###################################################################### +## ## +## Width set on Auto will be computed following the information ## +## present in the decay.py files of the model. ## +## See arXiv:1402.1178 for more details. ## +## ## +###################################################################### + +################################### +## INFORMATION FOR MASS +################################### +Block mass + 5 4.700000e+00 # MB + 6 1.730000e+02 # MT + 15 1.777000e+00 # MTA + 23 9.118800e+01 # MZ + 25 1.250000e+02 # MH +## Dependent parameters, given by model restrictions. +## Those values should be edited following the +## analytical expression. MG5 ignores those values +## but they are important for interfacing the output of MG5 +## to external program such as Pythia. + 1 0.000000e+00 # d : 0.0 + 2 0.000000e+00 # u : 0.0 + 3 0.000000e+00 # s : 0.0 + 4 0.000000e+00 # c : 0.0 + 11 0.000000e+00 # e- : 0.0 + 12 0.000000e+00 # ve : 0.0 + 13 0.000000e+00 # mu- : 0.0 + 14 0.000000e+00 # vm : 0.0 + 16 0.000000e+00 # vt : 0.0 + 21 0.000000e+00 # g : 0.0 + 22 0.000000e+00 # a : 0.0 + 24 8.041900e+01 # w+ : cmath.sqrt(MZ__exp__2/2. + cmath.sqrt(MZ__exp__4/4. - (aEW*cmath.pi*MZ__exp__2)/(Gf*sqrt__2))) + +################################### +## INFORMATION FOR SMINPUTS +################################### +Block sminputs + 1 1.325070e+02 # aEWM1 + 2 1.166390e-05 # Gf + 3 1.180000e-01 # aS + +################################### +## INFORMATION FOR YUKAWA +################################### +Block yukawa + 5 4.700000e+00 # ymb + 6 1.730000e+02 # ymt + 15 1.777000e+00 # ymtau + +################################### +## INFORMATION FOR DECAY +################################### +DECAY 6 1.491500e+00 # WT +DECAY 23 2.441404e+00 # WZ +DECAY 24 2.047600e+00 # WW +DECAY 25 6.382339e-03 # WH +## Dependent parameters, given by model restrictions. +## Those values should be edited following the +## analytical expression. MG5 ignores those values +## but they are important for interfacing the output of MG5 +## to external program such as Pythia. +DECAY 1 0.000000e+00 # d : 0.0 +DECAY 2 0.000000e+00 # u : 0.0 +DECAY 3 0.000000e+00 # s : 0.0 +DECAY 4 0.000000e+00 # c : 0.0 +DECAY 5 0.000000e+00 # b : 0.0 +DECAY 11 0.000000e+00 # e- : 0.0 +DECAY 12 0.000000e+00 # ve : 0.0 +DECAY 13 0.000000e+00 # mu- : 0.0 +DECAY 14 0.000000e+00 # vm : 0.0 +DECAY 15 0.000000e+00 # ta- : 0.0 +DECAY 16 0.000000e+00 # vt : 0.0 +DECAY 21 0.000000e+00 # g : 0.0 +DECAY 22 0.000000e+00 # a : 0.0 diff --git a/epoch1/oneapi/gg_ttgg/SubProcesses/Makefile b/epoch1/oneapi/gg_ttgg/SubProcesses/Makefile new file mode 100644 index 0000000000..413fdd3277 --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/SubProcesses/Makefile @@ -0,0 +1,61 @@ + +LIBDIR=../../lib +INCDIR=../../src +MODELLIB=model_sm +CXXFLAGS=-Wno-sycl-strict -I$(INCDIR) -I. +#CXXFLAGS= -O3 -I$(INCDIR) -I. +#CUARCHFLAGS= -arch=compute_$(CUARCHNUM) #--gpu-architecture=compute_35 --gpu-code=sm_35 +#CUFLAGS= $(CUARCHFLAGS) -use_fast_math -lineinfo +LIBFLAGS= -L$(LIBDIR) -l$(MODELLIB) +#NVCC=nvcc +CXX=dpcpp +MAKEDEBUG= + +main=check.exe +cxx_objects=CPPProcess.o check_sa.o + +all: check + +check: $(main) + +debug: CXXFLAGS:=$(filter-out -O3,$(CXXFLAGS)) +debug: CXXFLAGS += -g -O0 -DDEBUG2 +debug: MAKEDEBUG := debug +debug: check + +$(LIBDIR)/lib$(MODELLIB).a: + @cd ../../src && make $(MAKEDEBUG) + +#check_sa.o: check_sa.cpp +# $(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $< -o $@ + +%.o : %.cpp %.h + $(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $< -o $@ + +$(main): $(LIBDIR)/lib$(MODELLIB).a $(cxx_objects) $(cu_objects) + $(CXX) -o $@ $(cxx_objects) $(LIBFLAGS) + +.PHONY: clean + +clean: + cd ../../src && make clean + rm -f $(main) + rm -f $(cxx_objects) + +memcheck: check + /usr/local/cuda/bin/cuda-memcheck --check-api-memory-access yes --check-deprecated-instr yes --check-device-heap yes --demangle full --language c --leak-check full --racecheck-report all --report-api-errors all --show-backtrace yes --tool memcheck --track-unused-memory yes ./check.exe 5 5 10 + +perf: force + make clean && make + time ./check.exe -p 1024 16 384 && date + +test: force + ./check.exe -v 1 1 10 + +force: + +#Allowed values for this option: 'compute_30', 'compute_32', 'compute_35', 'compute_37', 'compute_50', 'compute_52', 'compute_53', 'compute_60', 'compute_61', 'compute_62', 'compute_70', 'compute_72', 'compute_75', 'sm_30', 'sm_32', 'sm_35', 'sm_37', 'sm_50', 'sm_52', 'sm_53', 'sm_60', 'sm_61', 'sm_62', 'sm_70', 'sm_72', 'sm_75'. + +# Max compute architectures +# cern batch (tesla v100): 70 +# jetson nano (maxwell): 35 diff --git a/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/CPPProcess.cpp b/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/CPPProcess.cpp new file mode 100644 index 0000000000..9e65ce9b73 --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/CPPProcess.cpp @@ -0,0 +1,21 @@ +//========================================================================== +// This file has been automatically generated for C++ Standalone by +// MadGraph5_aMC@NLO v. 2.7.3.py3, 2020-06-28 +// By the MadGraph5_aMC@NLO Development Team +// Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch +//========================================================================== + +#include +#include "HelAmps_sm.h" +#include +#include +#include +#include + + +//========================================================================== +// Private class member functions + +//-------------------------------------------------------------------------- + + diff --git a/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/CPPProcess.h b/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/CPPProcess.h new file mode 100644 index 0000000000..180598d271 --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/CPPProcess.h @@ -0,0 +1,1964 @@ +//========================================================================== +// This file has been automatically generated for C++ Standalone +// MadGraph5_aMC@NLO v. 2.7.3.py3, 2020-06-28 +// By the MadGraph5_aMC@NLO Development Team +// Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch +// GPU Template +//========================================================================== + +#ifndef HelAmps_sm_H +#define HelAmps_sm_H + +#include +#include +#include "Parameters_sm.h" + +int cHel[64][6]; +double cmME[6]; // value hardcoded now +int cPerm[4]; +// +double cIPC[6]; // coupling ? +double cIPD[2]; + +///// + +void ixxxxx(double pvec[3], double fmass, int nhel, int nsf, + std::complex fi[6]) +{ + std::complex chi[2]; + double sf[2], sfomega[2], omega[2], pp, pp3, sqp0p3, sqm[2]; + int ip, im, nh; + + double p[4] = {0, pvec[0], pvec[1], pvec[2]}; + p[0] = sycl::sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3] + fmass * fmass); + fi[0] = std::complex(-p[0] * nsf, -p[3] * nsf); + fi[1] = std::complex(-p[1] * nsf, -p[2] * nsf); + nh = nhel * nsf; + if (fmass != 0.0) + { + pp = sycl::min(p[0], sycl::sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3])); + if (pp == 0.0) + { + sqm[0] = sycl::sqrt(sycl::fabs(fmass)); + sqm[1] = (fmass < 0) ? -sycl::fabs(sqm[0]) : sycl::fabs(sqm[0]); + ip = (1 + nh)/2; + im = (1 - nh)/2; + fi[2] = ip * sqm[ip]; + fi[3] = im * nsf * sqm[ip]; + fi[4] = ip * nsf * sqm[im]; + fi[5] = im * sqm[im]; + } + else + { + sf[0] = (1 + nsf + (1 - nsf) * nh) * 0.5; + sf[1] = (1 + nsf - (1 - nsf) * nh) * 0.5; + omega[0] = sycl::sqrt(p[0] + pp); + omega[1] = fmass/omega[0]; + ip = (1 + nh)/2; + im = (1 - nh)/2; + sfomega[0] = sf[0] * omega[ip]; + sfomega[1] = sf[1] * omega[im]; + pp3 = sycl::max((double)(pp + p[3]), 0.0); + chi[0] = std::complex(sycl::sqrt(pp3 * 0.5 / pp), 0); + if (pp3 == 0.0) + { + chi[1] = std::complex(-nh, 0); + } + else + { + chi[1] = + std::complex(nh * p[1], p[2]) / sycl::sqrt(2.0 * pp * pp3); + } + fi[2] = sfomega[0] * chi[im]; + fi[3] = sfomega[0] * chi[ip]; + fi[4] = sfomega[1] * chi[im]; + fi[5] = sfomega[1] * chi[ip]; + } + } + else + { + if (p[1] == 0.0 and p[2] == 0.0 and p[3] < 0.0) + { + sqp0p3 = 0.0; + } + else + { + sqp0p3 = sycl::sqrt(sycl::max((double)(p[0] + p[3]), 0.0)) * nsf; + } + chi[0] = std::complex(sqp0p3, 0.0); + if (sqp0p3 == 0.0) + { + chi[1] = std::complex(-nhel * sycl::sqrt(2.0 * p[0]), 0.0); + } + else + { + chi[1] = std::complex(nh * p[1], p[2]) / sqp0p3; + } + if (nh == 1) + { + fi[2] = std::complex(0.0, 0.0); + fi[3] = std::complex(0.0, 0.0); + fi[4] = chi[0]; + fi[5] = chi[1]; + } + else + { + fi[2] = chi[1]; + fi[3] = chi[0]; + fi[4] = std::complex(0.0, 0.0); + fi[5] = std::complex(0.0, 0.0); + } + } + + return; +} + +void txxxxx(double pvec[3], double tmass, int nhel, int nst, + std::complex tc[18]) +{ + std::complex ft[6][4], ep[4], em[4], e0[4]; + double pt, pt2, pp, pzpt, emp, sqh, sqs; + int i, j; + + double p[4] = {0, pvec[0], pvec[1], pvec[2]}; + p[0] = sycl::sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3] + tmass * tmass); + sqh = sycl::sqrt(0.5); + sqs = sycl::sqrt(0.5 / 3); + + pt2 = p[1] * p[1] + p[2] * p[2]; + pp = sycl::min(p[0], sycl::sqrt(pt2 + p[3] * p[3])); + pt = sycl::min(pp, sycl::sqrt(pt2)); + + ft[4][0] = std::complex(p[0] * nst, p[3] * nst); + ft[5][0] = std::complex(p[1] * nst, p[2] * nst); + + // construct eps+ + if (nhel >= 0) + { + if (pp == 0) + { + ep[0] = std::complex(0, 0); + ep[1] = std::complex(-sqh, 0); + ep[2] = std::complex(0, nst * sqh); + ep[3] = std::complex(0, 0); + } + else + { + ep[0] = std::complex(0, 0); + ep[3] = std::complex(pt / pp * sqh, 0); + + if (pt != 0) + { + pzpt = p[3]/(pp * pt) * sqh; + ep[1] = std::complex(-p[1] * pzpt, -nst * p[2] / pt * sqh); + ep[2] = std::complex(-p[2] * pzpt, nst * p[1] / pt * sqh); + } + else + { + ep[1] = std::complex(-sqh, 0); + ep[2] = std::complex(0, (nst * (p[3] < 0)) ? -sycl::fabs(sqh) + : sycl::fabs(sqh)); + } + } + } + + // construct eps- + if (nhel <= 0) + { + if (pp == 0) + { + em[0] = std::complex(0, 0); + em[1] = std::complex(sqh, 0); + em[2] = std::complex(0, nst * sqh); + em[3] = std::complex(0, 0); + } + else + { + em[0] = std::complex(0, 0); + em[3] = std::complex(-pt / pp * sqh, 0); + + if (pt != 0) + { + pzpt = -p[3]/(pp * pt) * sqh; + em[1] = std::complex(-p[1] * pzpt, -nst * p[2] / pt * sqh); + em[2] = std::complex(-p[2] * pzpt, nst * p[1] / pt * sqh); + } + else + { + em[1] = std::complex(sqh, 0); + em[2] = std::complex(0, (nst * (p[3] < 0)) ? -sycl::fabs(sqh) + : sycl::fabs(sqh)); + } + } + } + + // construct eps0 + if (std::labs(nhel) <= 1) + { + if (pp == 0) + { + e0[0] = std::complex(0, 0); + e0[1] = std::complex(0, 0); + e0[2] = std::complex(0, 0); + e0[3] = std::complex(1, 0); + } + else + { + emp = p[0]/(tmass * pp); + e0[0] = std::complex(pp / tmass, 0); + e0[3] = std::complex(p[3] * emp, 0); + + if (pt != 0) + { + e0[1] = std::complex(p[1] * emp, 0); + e0[2] = std::complex(p[2] * emp, 0); + } + else + { + e0[1] = std::complex(0, 0); + e0[2] = std::complex(0, 0); + } + } + } + + if (nhel == 2) + { + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + ft[i][j] = ep[i] * ep[j]; + } + } + else if (nhel == -2) + { + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + ft[i][j] = em[i] * em[j]; + } + } + else if (tmass == 0) + { + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + ft[i][j] = 0; + } + } + else if (tmass != 0) + { + if (nhel == 1) + { + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + ft[i][j] = sqh * (ep[i] * e0[j] + e0[i] * ep[j]); + } + } + else if (nhel == 0) + { + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + ft[i][j] = + sqs * (ep[i] * em[j] + em[i] * ep[j] + 2.0 * e0[i] * e0[j]); + } + } + else if (nhel == -1) + { + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + ft[i][j] = sqh * (em[i] * e0[j] + e0[i] * em[j]); + } + } + else + { + // sr fixme // std::cerr << "Invalid helicity in txxxxx.\n"; + // sr fixme // std::exit(1); + } + } + + tc[0] = ft[4][0]; + tc[1] = ft[5][0]; + + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + tc[j * 4 + i + 2] = ft[j][i]; + } +} + +void vxxxxx(double pvec[3], double vmass, int nhel, int nsv, + std::complex vc[6]) +{ + double hel, hel0, pt, pt2, pp, pzpt, emp, sqh; + int nsvahl; + + double p[4] = {0, pvec[0], pvec[1], pvec[2]}; + p[0] = sycl::sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3] + vmass * vmass); + + sqh = sycl::sqrt(0.5); + hel = double(nhel); + nsvahl = nsv * sycl::fabs(hel); + pt2 = (p[1] * p[1]) + (p[2] * p[2]); + pp = sycl::min(p[0], sycl::sqrt(pt2 + (p[3] * p[3]))); + pt = sycl::min(pp, sycl::sqrt(pt2)); + vc[0] = std::complex(p[0] * nsv, p[3] * nsv); + vc[1] = std::complex(p[1] * nsv, p[2] * nsv); + if (vmass != 0.0) + { + hel0 = 1.0 - sycl::fabs(hel); + if (pp == 0.0) + { + vc[2] = std::complex(0.0, 0.0); + vc[3] = std::complex(-hel * sqh, 0.0); + vc[4] = std::complex(0.0, nsvahl * sqh); + vc[5] = std::complex(hel0, 0.0); + } + else + { + emp = p[0]/(vmass * pp); + vc[2] = std::complex(hel0 * pp / vmass, 0.0); + vc[5] = + std::complex(hel0 * p[3] * emp + hel * pt / pp * sqh, 0.0); + if (pt != 0.0) + { + pzpt = p[3]/(pp * pt) * sqh * hel; + vc[3] = std::complex(hel0 * p[1] * emp - p[1] * pzpt, + -nsvahl * p[2] / pt * sqh); + vc[4] = std::complex(hel0 * p[2] * emp - p[2] * pzpt, + nsvahl * p[1] / pt * sqh); + } + else + { + vc[3] = std::complex(-hel * sqh, 0.0); + vc[4] = std::complex( + 0.0, (nsvahl * (p[3] < 0)) ? -sycl::fabs(sqh) : sycl::fabs(sqh)); + } + } + } + else + { + pp = p[0]; + pt = sycl::sqrt((p[1] * p[1]) + (p[2] * p[2])); + vc[2] = std::complex(0.0, 0.0); + vc[5] = std::complex(hel * pt / pp * sqh, 0.0); + if (pt != 0.0) + { + pzpt = p[3]/(pp * pt) * sqh * hel; + vc[3] = std::complex(-p[1] * pzpt, -nsv * p[2] / pt * sqh); + vc[4] = std::complex(-p[2] * pzpt, nsv * p[1] / pt * sqh); + } + else + { + vc[3] = std::complex(-hel * sqh, 0.0); + vc[4] = std::complex(0.0, (nsv * (p[3] < 0)) ? -sycl::fabs(sqh) + : sycl::fabs(sqh)); + } + } + return; +} + +void sxxxxx(double pvec[3], int nss, std::complex sc[3]) +{ + // double p[4] = {0, pvec[0], pvec[1], pvec[2]}; + // p[0] = sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3]+fmass*fmass); + double p[4] = {0, 0, 0, 0}; + sc[2] = std::complex(1.00, 0.00); + sc[0] = std::complex(p[0] * nss, p[3] * nss); + sc[1] = std::complex(p[1] * nss, p[2] * nss); + return; +} + +void oxxxxx(double pvec[3], double fmass, int nhel, int nsf, + std::complex fo[6]) +{ + std::complex chi[2]; + double sf[2], sfomeg[2], omega[2], pp, pp3, sqp0p3, sqm[2]; + int nh, ip, im; + + double p[4] = {0, pvec[0], pvec[1], pvec[2]}; + p[0] = sycl::sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3] + fmass * fmass); + + fo[0] = std::complex(p[0] * nsf, p[3] * nsf); + fo[1] = std::complex(p[1] * nsf, p[2] * nsf); + nh = nhel * nsf; + if (fmass != 0.000) + { + pp = sycl::min(p[0], + sycl::sqrt((p[1] * p[1]) + (p[2] * p[2]) + (p[3] * p[3]))); + if (pp == 0.000) + { + sqm[0] = sycl::sqrt(sycl::fabs(fmass)); + sqm[1] = (fmass < 0) ? -sycl::fabs(sqm[0]) : sycl::fabs(sqm[0]); + ip = -((1 - nh)/2) * nhel; + im = (1 + nh)/2 * nhel; + fo[2] = im * sqm[sycl::abs(ip)]; + fo[3] = ip * nsf * sqm[sycl::abs(ip)]; + fo[4] = im * nsf * sqm[sycl::abs(im)]; + fo[5] = ip * sqm[sycl::abs(im)]; + } + else + { + pp = sycl::min(p[0], + sycl::sqrt((p[1] * p[1]) + (p[2] * p[2]) + (p[3] * p[3]))); + sf[0] = double(1 + nsf + (1 - nsf) * nh) * 0.5; + sf[1] = double(1 + nsf - (1 - nsf) * nh) * 0.5; + omega[0] = sycl::sqrt(p[0] + pp); + omega[1] = fmass/omega[0]; + ip = (1 + nh)/2; + im = (1 - nh)/2; + sfomeg[0] = sf[0] * omega[ip]; + sfomeg[1] = sf[1] * omega[im]; + pp3 = sycl::max((double)(pp + p[3]), 0.00); + chi[0] = std::complex(sycl::sqrt(pp3 * 0.5 / pp), 0.00); + if (pp3 == 0.00) + { + chi[1] = std::complex(-nh, 0.00); + } + else + { + chi[1] = + std::complex(nh * p[1], -p[2]) / sycl::sqrt(2.0 * pp * pp3); + } + fo[2] = sfomeg[1] * chi[im]; + fo[3] = sfomeg[1] * chi[ip]; + fo[4] = sfomeg[0] * chi[im]; + fo[5] = sfomeg[0] * chi[ip]; + } + } + else + { + if ((p[1] == 0.00) and (p[2] == 0.00) and (p[3] < 0.00)) + { + sqp0p3 = 0.00; + } + else + { + sqp0p3 = sycl::sqrt(sycl::max((double)(p[0] + p[3]), 0.00)) * nsf; + } + chi[0] = std::complex(sqp0p3, 0.00); + if (sqp0p3 == 0.000) + { + chi[1] = std::complex(-nhel, 0.00) * sycl::sqrt(2.0 * p[0]); + } + else + { + chi[1] = std::complex(nh * p[1], -p[2]) / sqp0p3; + } + if (nh == 1) + { + fo[2] = chi[0]; + fo[3] = chi[1]; + fo[4] = std::complex(0.00, 0.00); + fo[5] = std::complex(0.00, 0.00); + } + else + { + fo[2] = std::complex(0.00, 0.00); + fo[3] = std::complex(0.00, 0.00); + fo[4] = chi[1]; + fo[5] = chi[0]; + } + } + return; +} +void VVVV3_0(std::complex V1[], const std::complex V2[], + const std::complex V3[], + const std::complex V4[], + const std::complex COUP, std::complex *vertex) +{ + std::complex cI = std::complex(0., 1.); + std::complex TMP0; + std::complex TMP1; + std::complex TMP2; + std::complex TMP3; + TMP3 = (V1[2] * V2[2] - V1[3] * V2[3] - V1[4] * V2[4] - V1[5] * V2[5]); + TMP1 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); + TMP2 = (V4[2] * V3[2] - V4[3] * V3[3] - V4[4] * V3[4] - V4[5] * V3[5]); + TMP0 = (V4[2] * V1[2] - V4[3] * V1[3] - V4[4] * V1[4] - V4[5] * V1[5]); + (*vertex) = COUP * (-cI * (TMP0 * TMP1) + cI * (TMP2 * TMP3)); +} + +void VVVV3P0_1(std::complex V2[], const std::complex V3[], + const std::complex V4[], + const std::complex COUP, const double M1, + const double W1, std::complex V1[]) +{ + std::complex cI = std::complex(0., 1.); + double P1[4]; + std::complex TMP1; + std::complex TMP2; + std::complex denom; + V1[0] = +V2[0] + V3[0] + V4[0]; + V1[1] = +V2[1] + V3[1] + V4[1]; + P1[0] = -V1[0].real(); + P1[1] = -V1[1].real(); + P1[2] = -V1[1].imag(); + P1[3] = -V1[0].imag(); + TMP1 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); + TMP2 = (V4[2] * V3[2] - V4[3] * V3[3] - V4[4] * V3[4] - V4[5] * V3[5]); + denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * + P1[3]) - M1 * (M1 - cI * W1)); + V1[2] = denom * (-cI * (V4[2] * TMP1) + cI * (V2[2] * TMP2)); + V1[3] = denom * (-cI * (V4[3] * TMP1) + cI * (V2[3] * TMP2)); + V1[4] = denom * (-cI * (V4[4] * TMP1) + cI * (V2[4] * TMP2)); + V1[5] = denom * (-cI * (V4[5] * TMP1) + cI * (V2[5] * TMP2)); +} + +void VVVV1_0(std::complex V1[], const std::complex V2[], + const std::complex V3[], + const std::complex V4[], + const std::complex COUP, std::complex *vertex) +{ + std::complex cI = std::complex(0., 1.); + std::complex TMP0; + std::complex TMP1; + std::complex TMP4; + std::complex TMP5; + TMP5 = (V1[2] * V3[2] - V1[3] * V3[3] - V1[4] * V3[4] - V1[5] * V3[5]); + TMP1 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); + TMP0 = (V4[2] * V1[2] - V4[3] * V1[3] - V4[4] * V1[4] - V4[5] * V1[5]); + TMP4 = (V4[2] * V2[2] - V4[3] * V2[3] - V4[4] * V2[4] - V4[5] * V2[5]); + (*vertex) = COUP * (-cI * (TMP0 * TMP1) + cI * (TMP4 * TMP5)); +} + +void VVVV1P0_1(std::complex V2[], const std::complex V3[], + const std::complex V4[], + const std::complex COUP, const double M1, + const double W1, std::complex V1[]) +{ + std::complex cI = std::complex(0., 1.); + double P1[4]; + std::complex TMP1; + std::complex TMP4; + std::complex denom; + V1[0] = +V2[0] + V3[0] + V4[0]; + V1[1] = +V2[1] + V3[1] + V4[1]; + P1[0] = -V1[0].real(); + P1[1] = -V1[1].real(); + P1[2] = -V1[1].imag(); + P1[3] = -V1[0].imag(); + TMP1 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); + TMP4 = (V4[2] * V2[2] - V4[3] * V2[3] - V4[4] * V2[4] - V4[5] * V2[5]); + denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * + P1[3]) - M1 * (M1 - cI * W1)); + V1[2] = denom * (-cI * (V4[2] * TMP1) + cI * (V3[2] * TMP4)); + V1[3] = denom * (-cI * (V4[3] * TMP1) + cI * (V3[3] * TMP4)); + V1[4] = denom * (-cI * (V4[4] * TMP1) + cI * (V3[4] * TMP4)); + V1[5] = denom * (-cI * (V4[5] * TMP1) + cI * (V3[5] * TMP4)); +} + +void FFV1_0(std::complex F1[], const std::complex F2[], + const std::complex V3[], + const std::complex COUP, std::complex *vertex) +{ + std::complex cI = std::complex(0., 1.); + std::complex TMP6; + TMP6 = (F1[2] * (F2[4] * (V3[2] + V3[5]) + F2[5] * (V3[3] + cI * (V3[4]))) + + (F1[3] * (F2[4] * (V3[3] - cI * (V3[4])) + F2[5] * (V3[2] - V3[5])) + + (F1[4] * (F2[2] * (V3[2] - V3[5]) - F2[3] * (V3[3] + cI * (V3[4]))) + + F1[5] * (F2[2] * (-V3[3] + cI * (V3[4])) + F2[3] * (V3[2] + V3[5]))))); + (*vertex) = COUP * - cI * TMP6; +} + +void FFV1_1(std::complex F2[], const std::complex V3[], + const std::complex COUP, const double M1, + const double W1, std::complex F1[]) +{ + std::complex cI = std::complex(0., 1.); + double P1[4]; + std::complex denom; + F1[0] = +F2[0] + V3[0]; + F1[1] = +F2[1] + V3[1]; + P1[0] = -F1[0].real(); + P1[1] = -F1[1].real(); + P1[2] = -F1[1].imag(); + P1[3] = -F1[0].imag(); + denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * + P1[3]) - M1 * (M1 - cI * W1)); + F1[2] = denom * cI * (F2[2] * (P1[0] * (-V3[2] + V3[5]) + (P1[1] * (V3[3] - + cI * (V3[4])) + (P1[2] * (+cI * (V3[3]) + V3[4]) + P1[3] * (-V3[2] + + V3[5])))) + (F2[3] * (P1[0] * (V3[3] + cI * (V3[4])) + (P1[1] * (-1.) * + (V3[2] + V3[5]) + (P1[2] * (-1.) * (+cI * (V3[2] + V3[5])) + P1[3] * + (V3[3] + cI * (V3[4]))))) + M1 * (F2[4] * (V3[2] + V3[5]) + F2[5] * + (V3[3] + cI * (V3[4]))))); + F1[3] = denom * (-cI) * (F2[2] * (P1[0] * (-V3[3] + cI * (V3[4])) + (P1[1] * + (V3[2] - V3[5]) + (P1[2] * (-cI * (V3[2]) + cI * (V3[5])) + P1[3] * + (V3[3] - cI * (V3[4]))))) + (F2[3] * (P1[0] * (V3[2] + V3[5]) + (P1[1] * + (-1.) * (V3[3] + cI * (V3[4])) + (P1[2] * (+cI * (V3[3]) - V3[4]) - P1[3] + * (V3[2] + V3[5])))) + M1 * (F2[4] * (-V3[3] + cI * (V3[4])) + F2[5] * + (-V3[2] + V3[5])))); + F1[4] = denom * (-cI) * (F2[4] * (P1[0] * (V3[2] + V3[5]) + (P1[1] * (-V3[3] + + cI * (V3[4])) + (P1[2] * (-1.) * (+cI * (V3[3]) + V3[4]) - P1[3] * + (V3[2] + V3[5])))) + (F2[5] * (P1[0] * (V3[3] + cI * (V3[4])) + (P1[1] * + (-V3[2] + V3[5]) + (P1[2] * (-cI * (V3[2]) + cI * (V3[5])) - P1[3] * + (V3[3] + cI * (V3[4]))))) + M1 * (F2[2] * (-V3[2] + V3[5]) + F2[3] * + (V3[3] + cI * (V3[4]))))); + F1[5] = denom * cI * (F2[4] * (P1[0] * (-V3[3] + cI * (V3[4])) + (P1[1] * + (V3[2] + V3[5]) + (P1[2] * (-1.) * (+cI * (V3[2] + V3[5])) + P1[3] * + (-V3[3] + cI * (V3[4]))))) + (F2[5] * (P1[0] * (-V3[2] + V3[5]) + (P1[1] + * (V3[3] + cI * (V3[4])) + (P1[2] * (-cI * (V3[3]) + V3[4]) + P1[3] * + (-V3[2] + V3[5])))) + M1 * (F2[2] * (-V3[3] + cI * (V3[4])) + F2[3] * + (V3[2] + V3[5])))); +} + +void FFV1_2(std::complex F1[], const std::complex V3[], + const std::complex COUP, const double M2, + const double W2, std::complex F2[]) +{ + std::complex cI = std::complex(0., 1.); + double P2[4]; + std::complex denom; + F2[0] = +F1[0] + V3[0]; + F2[1] = +F1[1] + V3[1]; + P2[0] = -F2[0].real(); + P2[1] = -F2[1].real(); + P2[2] = -F2[1].imag(); + P2[3] = -F2[0].imag(); + denom = COUP/((P2[0] * P2[0]) - (P2[1] * P2[1]) - (P2[2] * P2[2]) - (P2[3] * + P2[3]) - M2 * (M2 - cI * W2)); + F2[2] = denom * cI * (F1[2] * (P2[0] * (V3[2] + V3[5]) + (P2[1] * (-1.) * + (V3[3] + cI * (V3[4])) + (P2[2] * (+cI * (V3[3]) - V3[4]) - P2[3] * + (V3[2] + V3[5])))) + (F1[3] * (P2[0] * (V3[3] - cI * (V3[4])) + (P2[1] * + (-V3[2] + V3[5]) + (P2[2] * (+cI * (V3[2]) - cI * (V3[5])) + P2[3] * + (-V3[3] + cI * (V3[4]))))) + M2 * (F1[4] * (V3[2] - V3[5]) + F1[5] * + (-V3[3] + cI * (V3[4]))))); + F2[3] = denom * (-cI) * (F1[2] * (P2[0] * (-1.) * (V3[3] + cI * (V3[4])) + + (P2[1] * (V3[2] + V3[5]) + (P2[2] * (+cI * (V3[2] + V3[5])) - P2[3] * + (V3[3] + cI * (V3[4]))))) + (F1[3] * (P2[0] * (-V3[2] + V3[5]) + (P2[1] * + (V3[3] - cI * (V3[4])) + (P2[2] * (+cI * (V3[3]) + V3[4]) + P2[3] * + (-V3[2] + V3[5])))) + M2 * (F1[4] * (V3[3] + cI * (V3[4])) - F1[5] * + (V3[2] + V3[5])))); + F2[4] = denom * (-cI) * (F1[4] * (P2[0] * (-V3[2] + V3[5]) + (P2[1] * (V3[3] + + cI * (V3[4])) + (P2[2] * (-cI * (V3[3]) + V3[4]) + P2[3] * (-V3[2] + + V3[5])))) + (F1[5] * (P2[0] * (V3[3] - cI * (V3[4])) + (P2[1] * (-1.) * + (V3[2] + V3[5]) + (P2[2] * (+cI * (V3[2] + V3[5])) + P2[3] * (V3[3] - cI + * (V3[4]))))) + M2 * (F1[2] * (-1.) * (V3[2] + V3[5]) + F1[3] * (-V3[3] + + cI * (V3[4]))))); + F2[5] = denom * cI * (F1[4] * (P2[0] * (-1.) * (V3[3] + cI * (V3[4])) + + (P2[1] * (V3[2] - V3[5]) + (P2[2] * (+cI * (V3[2]) - cI * (V3[5])) + + P2[3] * (V3[3] + cI * (V3[4]))))) + (F1[5] * (P2[0] * (V3[2] + V3[5]) + + (P2[1] * (-V3[3] + cI * (V3[4])) + (P2[2] * (-1.) * (+cI * (V3[3]) + + V3[4]) - P2[3] * (V3[2] + V3[5])))) + M2 * (F1[2] * (V3[3] + cI * + (V3[4])) + F1[3] * (V3[2] - V3[5])))); +} + +void FFV1P0_3(std::complex F1[], const std::complex F2[], + const std::complex COUP, const double M3, + const double W3, std::complex V3[]) +{ + std::complex cI = std::complex(0., 1.); + double P3[4]; + std::complex denom; + V3[0] = +F1[0] + F2[0]; + V3[1] = +F1[1] + F2[1]; + P3[0] = -V3[0].real(); + P3[1] = -V3[1].real(); + P3[2] = -V3[1].imag(); + P3[3] = -V3[0].imag(); + denom = COUP/((P3[0] * P3[0]) - (P3[1] * P3[1]) - (P3[2] * P3[2]) - (P3[3] * + P3[3]) - M3 * (M3 - cI * W3)); + V3[2] = denom * (-cI) * (F1[2] * F2[4] + F1[3] * F2[5] + F1[4] * F2[2] + + F1[5] * F2[3]); + V3[3] = denom * (-cI) * (-F1[2] * F2[5] - F1[3] * F2[4] + F1[4] * F2[3] + + F1[5] * F2[2]); + V3[4] = denom * (-cI) * (-cI * (F1[2] * F2[5] + F1[5] * F2[2]) + cI * (F1[3] + * F2[4] + F1[4] * F2[3])); + V3[5] = denom * (-cI) * (-F1[2] * F2[4] - F1[5] * F2[3] + F1[3] * F2[5] + + F1[4] * F2[2]); +} + +void VVVV4_0(std::complex V1[], const std::complex V2[], + const std::complex V3[], + const std::complex V4[], + const std::complex COUP, std::complex *vertex) +{ + std::complex cI = std::complex(0., 1.); + std::complex TMP2; + std::complex TMP3; + std::complex TMP4; + std::complex TMP5; + TMP3 = (V1[2] * V2[2] - V1[3] * V2[3] - V1[4] * V2[4] - V1[5] * V2[5]); + TMP5 = (V1[2] * V3[2] - V1[3] * V3[3] - V1[4] * V3[4] - V1[5] * V3[5]); + TMP2 = (V4[2] * V3[2] - V4[3] * V3[3] - V4[4] * V3[4] - V4[5] * V3[5]); + TMP4 = (V4[2] * V2[2] - V4[3] * V2[3] - V4[4] * V2[4] - V4[5] * V2[5]); + (*vertex) = COUP * (-cI * (TMP4 * TMP5) + cI * (TMP2 * TMP3)); +} + +void VVVV4P0_1(std::complex V2[], const std::complex V3[], + const std::complex V4[], + const std::complex COUP, const double M1, + const double W1, std::complex V1[]) +{ + std::complex cI = std::complex(0., 1.); + double P1[4]; + std::complex TMP2; + std::complex TMP4; + std::complex denom; + V1[0] = +V2[0] + V3[0] + V4[0]; + V1[1] = +V2[1] + V3[1] + V4[1]; + P1[0] = -V1[0].real(); + P1[1] = -V1[1].real(); + P1[2] = -V1[1].imag(); + P1[3] = -V1[0].imag(); + TMP2 = (V4[2] * V3[2] - V4[3] * V3[3] - V4[4] * V3[4] - V4[5] * V3[5]); + TMP4 = (V4[2] * V2[2] - V4[3] * V2[3] - V4[4] * V2[4] - V4[5] * V2[5]); + denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * + P1[3]) - M1 * (M1 - cI * W1)); + V1[2] = denom * (-cI * (V3[2] * TMP4) + cI * (V2[2] * TMP2)); + V1[3] = denom * (-cI * (V3[3] * TMP4) + cI * (V2[3] * TMP2)); + V1[4] = denom * (-cI * (V3[4] * TMP4) + cI * (V2[4] * TMP2)); + V1[5] = denom * (-cI * (V3[5] * TMP4) + cI * (V2[5] * TMP2)); +} + +void VVV1_0(std::complex V1[], const std::complex V2[], + const std::complex V3[], + const std::complex COUP, std::complex *vertex) +{ + std::complex cI = std::complex(0., 1.); + double P1[4]; + double P2[4]; + double P3[4]; + std::complex TMP1; + std::complex TMP10; + std::complex TMP11; + std::complex TMP12; + std::complex TMP3; + std::complex TMP5; + std::complex TMP7; + std::complex TMP8; + std::complex TMP9; + P1[0] = V1[0].real(); + P1[1] = V1[1].real(); + P1[2] = V1[1].imag(); + P1[3] = V1[0].imag(); + P2[0] = V2[0].real(); + P2[1] = V2[1].real(); + P2[2] = V2[1].imag(); + P2[3] = V2[0].imag(); + P3[0] = V3[0].real(); + P3[1] = V3[1].real(); + P3[2] = V3[1].imag(); + P3[3] = V3[0].imag(); + TMP9 = (V2[2] * P1[0] - V2[3] * P1[1] - V2[4] * P1[2] - V2[5] * P1[3]); + TMP8 = (V3[2] * P2[0] - V3[3] * P2[1] - V3[4] * P2[2] - V3[5] * P2[3]); + TMP3 = (V1[2] * V2[2] - V1[3] * V2[3] - V1[4] * V2[4] - V1[5] * V2[5]); + TMP1 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); + TMP7 = (V3[2] * P1[0] - V3[3] * P1[1] - V3[4] * P1[2] - V3[5] * P1[3]); + TMP5 = (V1[2] * V3[2] - V1[3] * V3[3] - V1[4] * V3[4] - V1[5] * V3[5]); + TMP10 = (V2[2] * P3[0] - V2[3] * P3[1] - V2[4] * P3[2] - V2[5] * P3[3]); + TMP11 = (V1[2] * P2[0] - V1[3] * P2[1] - V1[4] * P2[2] - V1[5] * P2[3]); + TMP12 = (V1[2] * P3[0] - V1[3] * P3[1] - V1[4] * P3[2] - V1[5] * P3[3]); + (*vertex) = COUP * (TMP1 * (-cI * (TMP11) + cI * (TMP12)) + (TMP3 * (-cI * + (TMP7) + cI * (TMP8)) + TMP5 * (+cI * (TMP9) - cI * (TMP10)))); +} + +void VVV1P0_1(std::complex V2[], const std::complex V3[], + const std::complex COUP, const double M1, + const double W1, std::complex V1[]) +{ + std::complex cI = std::complex(0., 1.); + double P1[4]; + double P2[4]; + double P3[4]; + std::complex TMP1; + std::complex TMP10; + std::complex TMP7; + std::complex TMP8; + std::complex TMP9; + std::complex denom; + P2[0] = V2[0].real(); + P2[1] = V2[1].real(); + P2[2] = V2[1].imag(); + P2[3] = V2[0].imag(); + P3[0] = V3[0].real(); + P3[1] = V3[1].real(); + P3[2] = V3[1].imag(); + P3[3] = V3[0].imag(); + V1[0] = +V2[0] + V3[0]; + V1[1] = +V2[1] + V3[1]; + P1[0] = -V1[0].real(); + P1[1] = -V1[1].real(); + P1[2] = -V1[1].imag(); + P1[3] = -V1[0].imag(); + TMP9 = (V2[2] * P1[0] - V2[3] * P1[1] - V2[4] * P1[2] - V2[5] * P1[3]); + TMP8 = (V3[2] * P2[0] - V3[3] * P2[1] - V3[4] * P2[2] - V3[5] * P2[3]); + TMP1 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); + TMP7 = (V3[2] * P1[0] - V3[3] * P1[1] - V3[4] * P1[2] - V3[5] * P1[3]); + TMP10 = (V2[2] * P3[0] - V2[3] * P3[1] - V2[4] * P3[2] - V2[5] * P3[3]); + denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * + P1[3]) - M1 * (M1 - cI * W1)); + V1[2] = denom * (TMP1 * (-cI * (P2[0]) + cI * (P3[0])) + (V2[2] * (-cI * + (TMP7) + cI * (TMP8)) + V3[2] * (+cI * (TMP9) - cI * (TMP10)))); + V1[3] = denom * (TMP1 * (-cI * (P2[1]) + cI * (P3[1])) + (V2[3] * (-cI * + (TMP7) + cI * (TMP8)) + V3[3] * (+cI * (TMP9) - cI * (TMP10)))); + V1[4] = denom * (TMP1 * (-cI * (P2[2]) + cI * (P3[2])) + (V2[4] * (-cI * + (TMP7) + cI * (TMP8)) + V3[4] * (+cI * (TMP9) - cI * (TMP10)))); + V1[5] = denom * (TMP1 * (-cI * (P2[3]) + cI * (P3[3])) + (V2[5] * (-cI * + (TMP7) + cI * (TMP8)) + V3[5] * (+cI * (TMP9) - cI * (TMP10)))); +} + +void calculate_wavefunctions(int ihel, + double local_mom[6][3], + double &matrix, + int **cHel, + double *cIPC, + double *cIPD) +{ + std::complex amp[159]; + // Calculate wavefunctions for all processes + std::complex w[26][6]; + vxxxxx(local_mom[0], 0., cHel[ihel][0], -1, w[0]); + vxxxxx(local_mom[1], 0., cHel[ihel][1], -1, w[1]); + oxxxxx(local_mom[2], cIPD[0], cHel[ihel][2], +1, w[2]); + ixxxxx(local_mom[3], cIPD[0], cHel[ihel][3], -1, w[3]); + vxxxxx(local_mom[4], 0., cHel[ihel][4], +1, w[4]); + vxxxxx(local_mom[5], 0., cHel[ihel][5], +1, w[5]); + VVV1P0_1(w[0], w[1], std::complex(cIPC[0], cIPC[1]), 0., 0., w[6]); + FFV1P0_3(w[3], w[2], std::complex(cIPC[2], cIPC[3]), 0., 0., w[7]); + // Amplitude(s) for diagram number 1 + VVVV1_0(w[6], w[7], w[4], w[5], std::complex(cIPC[4], cIPC[5]), + &[0]); + VVVV3_0(w[6], w[7], w[4], w[5], std::complex(cIPC[4], cIPC[5]), + &[1]); + VVVV4_0(w[6], w[7], w[4], w[5], std::complex(cIPC[4], cIPC[5]), + &[2]); + VVV1P0_1(w[6], w[4], std::complex(cIPC[0], cIPC[1]), 0., 0., w[8]); + // Amplitude(s) for diagram number 2 + VVV1_0(w[7], w[5], w[8], std::complex(cIPC[0], cIPC[1]), &[3]); + VVV1P0_1(w[6], w[5], std::complex(cIPC[0], cIPC[1]), 0., 0., w[9]); + // Amplitude(s) for diagram number 3 + VVV1_0(w[7], w[4], w[9], std::complex(cIPC[0], cIPC[1]), &[4]); + VVV1P0_1(w[4], w[5], std::complex(cIPC[0], cIPC[1]), 0., 0., w[10]); + // Amplitude(s) for diagram number 4 + VVV1_0(w[6], w[7], w[10], std::complex(cIPC[0], cIPC[1]), &[5]); + FFV1_1(w[2], w[4], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[11]); + FFV1_2(w[3], w[6], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[12]); + // Amplitude(s) for diagram number 5 + FFV1_0(w[12], w[11], w[5], std::complex(cIPC[2], cIPC[3]), &[6]); + // Amplitude(s) for diagram number 6 + FFV1_0(w[3], w[11], w[9], std::complex(cIPC[2], cIPC[3]), &[7]); + FFV1_2(w[3], w[5], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[13]); + // Amplitude(s) for diagram number 7 + FFV1_0(w[13], w[11], w[6], std::complex(cIPC[2], cIPC[3]), &[8]); + FFV1_1(w[2], w[5], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[14]); + // Amplitude(s) for diagram number 8 + FFV1_0(w[12], w[14], w[4], std::complex(cIPC[2], cIPC[3]), &[9]); + // Amplitude(s) for diagram number 9 + FFV1_0(w[3], w[14], w[8], std::complex(cIPC[2], cIPC[3]), &[10]); + FFV1_2(w[3], w[4], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[15]); + // Amplitude(s) for diagram number 10 + FFV1_0(w[15], w[14], w[6], std::complex(cIPC[2], cIPC[3]), &[11]); + FFV1_1(w[2], w[6], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[16]); + // Amplitude(s) for diagram number 11 + FFV1_0(w[15], w[16], w[5], std::complex(cIPC[2], cIPC[3]), &[12]); + // Amplitude(s) for diagram number 12 + FFV1_0(w[15], w[2], w[9], std::complex(cIPC[2], cIPC[3]), &[13]); + // Amplitude(s) for diagram number 13 + FFV1_0(w[13], w[16], w[4], std::complex(cIPC[2], cIPC[3]), &[14]); + // Amplitude(s) for diagram number 14 + FFV1_0(w[13], w[2], w[8], std::complex(cIPC[2], cIPC[3]), &[15]); + // Amplitude(s) for diagram number 15 + FFV1_0(w[3], w[16], w[10], std::complex(cIPC[2], cIPC[3]), &[16]); + // Amplitude(s) for diagram number 16 + FFV1_0(w[12], w[2], w[10], std::complex(cIPC[2], cIPC[3]), &[17]); + FFV1_1(w[2], w[0], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[12]); + FFV1_2(w[3], w[1], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[16]); + FFV1_1(w[12], w[4], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[8]); + // Amplitude(s) for diagram number 17 + FFV1_0(w[16], w[8], w[5], std::complex(cIPC[2], cIPC[3]), &[18]); + FFV1_1(w[12], w[5], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[9]); + // Amplitude(s) for diagram number 18 + FFV1_0(w[16], w[9], w[4], std::complex(cIPC[2], cIPC[3]), &[19]); + // Amplitude(s) for diagram number 19 + FFV1_0(w[16], w[12], w[10], std::complex(cIPC[2], cIPC[3]), &[20]); + VVV1P0_1(w[1], w[4], std::complex(cIPC[0], cIPC[1]), 0., 0., w[6]); + FFV1P0_3(w[3], w[12], std::complex(cIPC[2], cIPC[3]), 0., 0., w[17]); + // Amplitude(s) for diagram number 20 + VVV1_0(w[6], w[5], w[17], std::complex(cIPC[0], cIPC[1]), &[21]); + // Amplitude(s) for diagram number 21 + FFV1_0(w[3], w[9], w[6], std::complex(cIPC[2], cIPC[3]), &[22]); + // Amplitude(s) for diagram number 22 + FFV1_0(w[13], w[12], w[6], std::complex(cIPC[2], cIPC[3]), &[23]); + VVV1P0_1(w[1], w[5], std::complex(cIPC[0], cIPC[1]), 0., 0., w[18]); + // Amplitude(s) for diagram number 23 + VVV1_0(w[18], w[4], w[17], std::complex(cIPC[0], cIPC[1]), &[24]); + // Amplitude(s) for diagram number 24 + FFV1_0(w[3], w[8], w[18], std::complex(cIPC[2], cIPC[3]), &[25]); + // Amplitude(s) for diagram number 25 + FFV1_0(w[15], w[12], w[18], std::complex(cIPC[2], cIPC[3]), &[26]); + FFV1_1(w[12], w[1], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[19]); + // Amplitude(s) for diagram number 26 + FFV1_0(w[15], w[19], w[5], std::complex(cIPC[2], cIPC[3]), &[27]); + // Amplitude(s) for diagram number 27 + FFV1_0(w[15], w[9], w[1], std::complex(cIPC[2], cIPC[3]), &[28]); + // Amplitude(s) for diagram number 28 + FFV1_0(w[13], w[19], w[4], std::complex(cIPC[2], cIPC[3]), &[29]); + // Amplitude(s) for diagram number 29 + FFV1_0(w[13], w[8], w[1], std::complex(cIPC[2], cIPC[3]), &[30]); + // Amplitude(s) for diagram number 30 + FFV1_0(w[3], w[19], w[10], std::complex(cIPC[2], cIPC[3]), &[31]); + // Amplitude(s) for diagram number 31 + VVV1_0(w[1], w[10], w[17], std::complex(cIPC[0], cIPC[1]), &[32]); + VVVV1P0_1(w[1], w[4], w[5], std::complex(cIPC[4], cIPC[5]), 0., 0., + w[17]); + VVVV3P0_1(w[1], w[4], w[5], std::complex(cIPC[4], cIPC[5]), 0., 0., + w[19]); + VVVV4P0_1(w[1], w[4], w[5], std::complex(cIPC[4], cIPC[5]), 0., 0., + w[8]); + // Amplitude(s) for diagram number 32 + FFV1_0(w[3], w[12], w[17], std::complex(cIPC[2], cIPC[3]), &[33]); + FFV1_0(w[3], w[12], w[19], std::complex(cIPC[2], cIPC[3]), &[34]); + FFV1_0(w[3], w[12], w[8], std::complex(cIPC[2], cIPC[3]), &[35]); + FFV1_2(w[3], w[0], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[12]); + FFV1_1(w[2], w[1], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[9]); + FFV1_2(w[12], w[4], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[20]); + // Amplitude(s) for diagram number 33 + FFV1_0(w[20], w[9], w[5], std::complex(cIPC[2], cIPC[3]), &[36]); + FFV1_2(w[12], w[5], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[21]); + // Amplitude(s) for diagram number 34 + FFV1_0(w[21], w[9], w[4], std::complex(cIPC[2], cIPC[3]), &[37]); + // Amplitude(s) for diagram number 35 + FFV1_0(w[12], w[9], w[10], std::complex(cIPC[2], cIPC[3]), &[38]); + FFV1P0_3(w[12], w[2], std::complex(cIPC[2], cIPC[3]), 0., 0., w[22]); + // Amplitude(s) for diagram number 36 + VVV1_0(w[6], w[5], w[22], std::complex(cIPC[0], cIPC[1]), &[39]); + // Amplitude(s) for diagram number 37 + FFV1_0(w[21], w[2], w[6], std::complex(cIPC[2], cIPC[3]), &[40]); + // Amplitude(s) for diagram number 38 + FFV1_0(w[12], w[14], w[6], std::complex(cIPC[2], cIPC[3]), &[41]); + // Amplitude(s) for diagram number 39 + VVV1_0(w[18], w[4], w[22], std::complex(cIPC[0], cIPC[1]), &[42]); + // Amplitude(s) for diagram number 40 + FFV1_0(w[20], w[2], w[18], std::complex(cIPC[2], cIPC[3]), &[43]); + // Amplitude(s) for diagram number 41 + FFV1_0(w[12], w[11], w[18], std::complex(cIPC[2], cIPC[3]), &[44]); + FFV1_2(w[12], w[1], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[23]); + // Amplitude(s) for diagram number 42 + FFV1_0(w[23], w[11], w[5], std::complex(cIPC[2], cIPC[3]), &[45]); + // Amplitude(s) for diagram number 43 + FFV1_0(w[21], w[11], w[1], std::complex(cIPC[2], cIPC[3]), &[46]); + // Amplitude(s) for diagram number 44 + FFV1_0(w[23], w[14], w[4], std::complex(cIPC[2], cIPC[3]), &[47]); + // Amplitude(s) for diagram number 45 + FFV1_0(w[20], w[14], w[1], std::complex(cIPC[2], cIPC[3]), &[48]); + // Amplitude(s) for diagram number 46 + FFV1_0(w[23], w[2], w[10], std::complex(cIPC[2], cIPC[3]), &[49]); + // Amplitude(s) for diagram number 47 + VVV1_0(w[1], w[10], w[22], std::complex(cIPC[0], cIPC[1]), &[50]); + // Amplitude(s) for diagram number 48 + FFV1_0(w[12], w[2], w[17], std::complex(cIPC[2], cIPC[3]), &[51]); + FFV1_0(w[12], w[2], w[19], std::complex(cIPC[2], cIPC[3]), &[52]); + FFV1_0(w[12], w[2], w[8], std::complex(cIPC[2], cIPC[3]), &[53]); + VVV1P0_1(w[0], w[4], std::complex(cIPC[0], cIPC[1]), 0., 0., w[12]); + FFV1_2(w[3], w[12], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[22]); + // Amplitude(s) for diagram number 49 + FFV1_0(w[22], w[9], w[5], std::complex(cIPC[2], cIPC[3]), &[54]); + VVV1P0_1(w[12], w[5], std::complex(cIPC[0], cIPC[1]), 0., 0., w[23]); + // Amplitude(s) for diagram number 50 + FFV1_0(w[3], w[9], w[23], std::complex(cIPC[2], cIPC[3]), &[55]); + // Amplitude(s) for diagram number 51 + FFV1_0(w[13], w[9], w[12], std::complex(cIPC[2], cIPC[3]), &[56]); + FFV1_1(w[2], w[12], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[20]); + // Amplitude(s) for diagram number 52 + FFV1_0(w[16], w[20], w[5], std::complex(cIPC[2], cIPC[3]), &[57]); + // Amplitude(s) for diagram number 53 + FFV1_0(w[16], w[2], w[23], std::complex(cIPC[2], cIPC[3]), &[58]); + // Amplitude(s) for diagram number 54 + FFV1_0(w[16], w[14], w[12], std::complex(cIPC[2], cIPC[3]), &[59]); + // Amplitude(s) for diagram number 55 + FFV1_0(w[3], w[20], w[18], std::complex(cIPC[2], cIPC[3]), &[60]); + // Amplitude(s) for diagram number 56 + FFV1_0(w[22], w[2], w[18], std::complex(cIPC[2], cIPC[3]), &[61]); + // Amplitude(s) for diagram number 57 + VVV1_0(w[12], w[18], w[7], std::complex(cIPC[0], cIPC[1]), &[62]); + // Amplitude(s) for diagram number 58 + VVVV1_0(w[12], w[1], w[7], w[5], std::complex(cIPC[4], cIPC[5]), + &[63]); + VVVV3_0(w[12], w[1], w[7], w[5], std::complex(cIPC[4], cIPC[5]), + &[64]); + VVVV4_0(w[12], w[1], w[7], w[5], std::complex(cIPC[4], cIPC[5]), + &[65]); + VVV1P0_1(w[12], w[1], std::complex(cIPC[0], cIPC[1]), 0., 0., w[21]); + // Amplitude(s) for diagram number 59 + VVV1_0(w[7], w[5], w[21], std::complex(cIPC[0], cIPC[1]), &[66]); + // Amplitude(s) for diagram number 60 + VVV1_0(w[1], w[7], w[23], std::complex(cIPC[0], cIPC[1]), &[67]); + // Amplitude(s) for diagram number 61 + FFV1_0(w[3], w[14], w[21], std::complex(cIPC[2], cIPC[3]), &[68]); + // Amplitude(s) for diagram number 62 + FFV1_0(w[22], w[14], w[1], std::complex(cIPC[2], cIPC[3]), &[69]); + // Amplitude(s) for diagram number 63 + FFV1_0(w[13], w[2], w[21], std::complex(cIPC[2], cIPC[3]), &[70]); + // Amplitude(s) for diagram number 64 + FFV1_0(w[13], w[20], w[1], std::complex(cIPC[2], cIPC[3]), &[71]); + VVV1P0_1(w[0], w[5], std::complex(cIPC[0], cIPC[1]), 0., 0., w[20]); + FFV1_2(w[3], w[20], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[21]); + // Amplitude(s) for diagram number 65 + FFV1_0(w[21], w[9], w[4], std::complex(cIPC[2], cIPC[3]), &[72]); + VVV1P0_1(w[20], w[4], std::complex(cIPC[0], cIPC[1]), 0., 0., w[22]); + // Amplitude(s) for diagram number 66 + FFV1_0(w[3], w[9], w[22], std::complex(cIPC[2], cIPC[3]), &[73]); + // Amplitude(s) for diagram number 67 + FFV1_0(w[15], w[9], w[20], std::complex(cIPC[2], cIPC[3]), &[74]); + FFV1_1(w[2], w[20], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[23]); + // Amplitude(s) for diagram number 68 + FFV1_0(w[16], w[23], w[4], std::complex(cIPC[2], cIPC[3]), &[75]); + // Amplitude(s) for diagram number 69 + FFV1_0(w[16], w[2], w[22], std::complex(cIPC[2], cIPC[3]), &[76]); + // Amplitude(s) for diagram number 70 + FFV1_0(w[16], w[11], w[20], std::complex(cIPC[2], cIPC[3]), &[77]); + // Amplitude(s) for diagram number 71 + FFV1_0(w[3], w[23], w[6], std::complex(cIPC[2], cIPC[3]), &[78]); + // Amplitude(s) for diagram number 72 + FFV1_0(w[21], w[2], w[6], std::complex(cIPC[2], cIPC[3]), &[79]); + // Amplitude(s) for diagram number 73 + VVV1_0(w[20], w[6], w[7], std::complex(cIPC[0], cIPC[1]), &[80]); + // Amplitude(s) for diagram number 74 + VVVV1_0(w[20], w[1], w[7], w[4], std::complex(cIPC[4], cIPC[5]), + &[81]); + VVVV3_0(w[20], w[1], w[7], w[4], std::complex(cIPC[4], cIPC[5]), + &[82]); + VVVV4_0(w[20], w[1], w[7], w[4], std::complex(cIPC[4], cIPC[5]), + &[83]); + VVV1P0_1(w[20], w[1], std::complex(cIPC[0], cIPC[1]), 0., 0., w[12]); + // Amplitude(s) for diagram number 75 + VVV1_0(w[7], w[4], w[12], std::complex(cIPC[0], cIPC[1]), &[84]); + // Amplitude(s) for diagram number 76 + VVV1_0(w[1], w[7], w[22], std::complex(cIPC[0], cIPC[1]), &[85]); + // Amplitude(s) for diagram number 77 + FFV1_0(w[3], w[11], w[12], std::complex(cIPC[2], cIPC[3]), &[86]); + // Amplitude(s) for diagram number 78 + FFV1_0(w[21], w[11], w[1], std::complex(cIPC[2], cIPC[3]), &[87]); + // Amplitude(s) for diagram number 79 + FFV1_0(w[15], w[2], w[12], std::complex(cIPC[2], cIPC[3]), &[88]); + // Amplitude(s) for diagram number 80 + FFV1_0(w[15], w[23], w[1], std::complex(cIPC[2], cIPC[3]), &[89]); + FFV1_1(w[9], w[0], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[23]); + // Amplitude(s) for diagram number 81 + FFV1_0(w[15], w[23], w[5], std::complex(cIPC[2], cIPC[3]), &[90]); + FFV1_2(w[15], w[0], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[12]); + // Amplitude(s) for diagram number 82 + FFV1_0(w[12], w[9], w[5], std::complex(cIPC[2], cIPC[3]), &[91]); + // Amplitude(s) for diagram number 83 + FFV1_0(w[13], w[23], w[4], std::complex(cIPC[2], cIPC[3]), &[92]); + FFV1_2(w[13], w[0], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[21]); + // Amplitude(s) for diagram number 84 + FFV1_0(w[21], w[9], w[4], std::complex(cIPC[2], cIPC[3]), &[93]); + // Amplitude(s) for diagram number 85 + FFV1_0(w[3], w[23], w[10], std::complex(cIPC[2], cIPC[3]), &[94]); + VVV1P0_1(w[0], w[10], std::complex(cIPC[0], cIPC[1]), 0., 0., w[23]); + // Amplitude(s) for diagram number 86 + FFV1_0(w[3], w[9], w[23], std::complex(cIPC[2], cIPC[3]), &[95]); + FFV1_2(w[16], w[0], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[22]); + // Amplitude(s) for diagram number 87 + FFV1_0(w[22], w[11], w[5], std::complex(cIPC[2], cIPC[3]), &[96]); + FFV1_1(w[11], w[0], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[20]); + // Amplitude(s) for diagram number 88 + FFV1_0(w[16], w[20], w[5], std::complex(cIPC[2], cIPC[3]), &[97]); + // Amplitude(s) for diagram number 89 + FFV1_0(w[22], w[14], w[4], std::complex(cIPC[2], cIPC[3]), &[98]); + FFV1_1(w[14], w[0], std::complex(cIPC[2], cIPC[3]), cIPD[0], cIPD[1], + w[24]); + // Amplitude(s) for diagram number 90 + FFV1_0(w[16], w[24], w[4], std::complex(cIPC[2], cIPC[3]), &[99]); + // Amplitude(s) for diagram number 91 + FFV1_0(w[22], w[2], w[10], std::complex(cIPC[2], cIPC[3]), &[100]); + // Amplitude(s) for diagram number 92 + FFV1_0(w[16], w[2], w[23], std::complex(cIPC[2], cIPC[3]), &[101]); + // Amplitude(s) for diagram number 93 + VVVV1_0(w[0], w[6], w[7], w[5], std::complex(cIPC[4], cIPC[5]), + &[102]); + VVVV3_0(w[0], w[6], w[7], w[5], std::complex(cIPC[4], cIPC[5]), + &[103]); + VVVV4_0(w[0], w[6], w[7], w[5], std::complex(cIPC[4], cIPC[5]), + &[104]); + VVV1P0_1(w[0], w[6], std::complex(cIPC[0], cIPC[1]), 0., 0., w[22]); + // Amplitude(s) for diagram number 94 + VVV1_0(w[7], w[5], w[22], std::complex(cIPC[0], cIPC[1]), &[105]); + VVV1P0_1(w[0], w[7], std::complex(cIPC[0], cIPC[1]), 0., 0., w[25]); + // Amplitude(s) for diagram number 95 + VVV1_0(w[6], w[5], w[25], std::complex(cIPC[0], cIPC[1]), &[106]); + // Amplitude(s) for diagram number 96 + FFV1_0(w[3], w[14], w[22], std::complex(cIPC[2], cIPC[3]), &[107]); + // Amplitude(s) for diagram number 97 + FFV1_0(w[3], w[24], w[6], std::complex(cIPC[2], cIPC[3]), &[108]); + // Amplitude(s) for diagram number 98 + FFV1_0(w[13], w[2], w[22], std::complex(cIPC[2], cIPC[3]), &[109]); + // Amplitude(s) for diagram number 99 + FFV1_0(w[21], w[2], w[6], std::complex(cIPC[2], cIPC[3]), &[110]); + // Amplitude(s) for diagram number 100 + VVVV1_0(w[0], w[18], w[7], w[4], std::complex(cIPC[4], cIPC[5]), + &[111]); + VVVV3_0(w[0], w[18], w[7], w[4], std::complex(cIPC[4], cIPC[5]), + &[112]); + VVVV4_0(w[0], w[18], w[7], w[4], std::complex(cIPC[4], cIPC[5]), + &[113]); + VVV1P0_1(w[0], w[18], std::complex(cIPC[0], cIPC[1]), 0., 0., w[6]); + // Amplitude(s) for diagram number 101 + VVV1_0(w[7], w[4], w[6], std::complex(cIPC[0], cIPC[1]), &[114]); + // Amplitude(s) for diagram number 102 + VVV1_0(w[18], w[4], w[25], std::complex(cIPC[0], cIPC[1]), &[115]); + // Amplitude(s) for diagram number 103 + FFV1_0(w[3], w[11], w[6], std::complex(cIPC[2], cIPC[3]), &[116]); + // Amplitude(s) for diagram number 104 + FFV1_0(w[3], w[20], w[18], std::complex(cIPC[2], cIPC[3]), &[117]); + // Amplitude(s) for diagram number 105 + FFV1_0(w[15], w[2], w[6], std::complex(cIPC[2], cIPC[3]), &[118]); + // Amplitude(s) for diagram number 106 + FFV1_0(w[12], w[2], w[18], std::complex(cIPC[2], cIPC[3]), &[119]); + // Amplitude(s) for diagram number 107 + VVVV1_0(w[0], w[1], w[7], w[10], std::complex(cIPC[4], cIPC[5]), + &[120]); + VVVV3_0(w[0], w[1], w[7], w[10], std::complex(cIPC[4], cIPC[5]), + &[121]); + VVVV4_0(w[0], w[1], w[7], w[10], std::complex(cIPC[4], cIPC[5]), + &[122]); + // Amplitude(s) for diagram number 108 + VVV1_0(w[1], w[10], w[25], std::complex(cIPC[0], cIPC[1]), &[123]); + // Amplitude(s) for diagram number 109 + VVV1_0(w[1], w[7], w[23], std::complex(cIPC[0], cIPC[1]), &[124]); + // Amplitude(s) for diagram number 110 + FFV1_0(w[13], w[20], w[1], std::complex(cIPC[2], cIPC[3]), &[125]); + // Amplitude(s) for diagram number 111 + FFV1_0(w[21], w[11], w[1], std::complex(cIPC[2], cIPC[3]), &[126]); + // Amplitude(s) for diagram number 112 + FFV1_0(w[15], w[24], w[1], std::complex(cIPC[2], cIPC[3]), &[127]); + // Amplitude(s) for diagram number 113 + FFV1_0(w[12], w[14], w[1], std::complex(cIPC[2], cIPC[3]), &[128]); + VVVV1P0_1(w[0], w[1], w[4], std::complex(cIPC[4], cIPC[5]), 0., 0., + w[12]); + VVVV3P0_1(w[0], w[1], w[4], std::complex(cIPC[4], cIPC[5]), 0., 0., + w[24]); + VVVV4P0_1(w[0], w[1], w[4], std::complex(cIPC[4], cIPC[5]), 0., 0., + w[21]); + // Amplitude(s) for diagram number 114 + VVV1_0(w[12], w[7], w[5], std::complex(cIPC[0], cIPC[1]), &[129]); + VVV1_0(w[24], w[7], w[5], std::complex(cIPC[0], cIPC[1]), &[130]); + VVV1_0(w[21], w[7], w[5], std::complex(cIPC[0], cIPC[1]), &[131]); + // Amplitude(s) for diagram number 115 + FFV1_0(w[3], w[14], w[12], std::complex(cIPC[2], cIPC[3]), &[132]); + FFV1_0(w[3], w[14], w[24], std::complex(cIPC[2], cIPC[3]), &[133]); + FFV1_0(w[3], w[14], w[21], std::complex(cIPC[2], cIPC[3]), &[134]); + // Amplitude(s) for diagram number 116 + FFV1_0(w[13], w[2], w[12], std::complex(cIPC[2], cIPC[3]), &[135]); + FFV1_0(w[13], w[2], w[24], std::complex(cIPC[2], cIPC[3]), &[136]); + FFV1_0(w[13], w[2], w[21], std::complex(cIPC[2], cIPC[3]), &[137]); + VVVV1P0_1(w[0], w[1], w[5], std::complex(cIPC[4], cIPC[5]), 0., 0., + w[21]); + VVVV3P0_1(w[0], w[1], w[5], std::complex(cIPC[4], cIPC[5]), 0., 0., + w[13]); + VVVV4P0_1(w[0], w[1], w[5], std::complex(cIPC[4], cIPC[5]), 0., 0., + w[24]); + // Amplitude(s) for diagram number 117 + VVV1_0(w[21], w[7], w[4], std::complex(cIPC[0], cIPC[1]), &[138]); + VVV1_0(w[13], w[7], w[4], std::complex(cIPC[0], cIPC[1]), &[139]); + VVV1_0(w[24], w[7], w[4], std::complex(cIPC[0], cIPC[1]), &[140]); + // Amplitude(s) for diagram number 118 + FFV1_0(w[3], w[11], w[21], std::complex(cIPC[2], cIPC[3]), &[141]); + FFV1_0(w[3], w[11], w[13], std::complex(cIPC[2], cIPC[3]), &[142]); + FFV1_0(w[3], w[11], w[24], std::complex(cIPC[2], cIPC[3]), &[143]); + // Amplitude(s) for diagram number 119 + FFV1_0(w[15], w[2], w[21], std::complex(cIPC[2], cIPC[3]), &[144]); + FFV1_0(w[15], w[2], w[13], std::complex(cIPC[2], cIPC[3]), &[145]); + FFV1_0(w[15], w[2], w[24], std::complex(cIPC[2], cIPC[3]), &[146]); + VVVV1P0_1(w[0], w[4], w[5], std::complex(cIPC[4], cIPC[5]), 0., 0., + w[24]); + VVVV3P0_1(w[0], w[4], w[5], std::complex(cIPC[4], cIPC[5]), 0., 0., + w[15]); + VVVV4P0_1(w[0], w[4], w[5], std::complex(cIPC[4], cIPC[5]), 0., 0., + w[13]); + // Amplitude(s) for diagram number 120 + FFV1_0(w[3], w[9], w[24], std::complex(cIPC[2], cIPC[3]), &[147]); + FFV1_0(w[3], w[9], w[15], std::complex(cIPC[2], cIPC[3]), &[148]); + FFV1_0(w[3], w[9], w[13], std::complex(cIPC[2], cIPC[3]), &[149]); + // Amplitude(s) for diagram number 121 + FFV1_0(w[16], w[2], w[24], std::complex(cIPC[2], cIPC[3]), &[150]); + FFV1_0(w[16], w[2], w[15], std::complex(cIPC[2], cIPC[3]), &[151]); + FFV1_0(w[16], w[2], w[13], std::complex(cIPC[2], cIPC[3]), &[152]); + // Amplitude(s) for diagram number 122 + VVV1_0(w[24], w[1], w[7], std::complex(cIPC[0], cIPC[1]), &[153]); + VVV1_0(w[15], w[1], w[7], std::complex(cIPC[0], cIPC[1]), &[154]); + VVV1_0(w[13], w[1], w[7], std::complex(cIPC[0], cIPC[1]), &[155]); + // Amplitude(s) for diagram number 123 + VVV1_0(w[0], w[17], w[7], std::complex(cIPC[0], cIPC[1]), &[156]); + VVV1_0(w[0], w[19], w[7], std::complex(cIPC[0], cIPC[1]), &[157]); + VVV1_0(w[0], w[8], w[7], std::complex(cIPC[0], cIPC[1]), &[158]); + // double CPPProcess::matrix_1_gg_ttxgg() { + int i, j; + // Local variables + + // const int ngraphs = 2; + const int ncolor = 24; + std::complex ztemp; + std::complex jamp[ncolor]; + // The color matrix; + static const double denom[ncolor] = {54, 54, 54, 54, 54, 54, 54, 54, 54, 54, + 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54}; + static const double cf[ncolor][ncolor] = {{512, -64, -64, 8, 8, 80, -64, 8, + 8, -1, -1, -10, 8, -1, 80, -10, 71, 62, -1, -10, -10, 62, 62, -28}, {-64, + 512, 8, 80, -64, 8, 8, -64, -1, -10, 8, -1, -1, -10, -10, 62, 62, -28, 8, + -1, 80, -10, 71, 62}, {-64, 8, 512, -64, 80, 8, 8, -1, 80, -10, 71, 62, + -64, 8, 8, -1, -1, -10, -10, -1, 62, -28, -10, 62}, {8, 80, -64, 512, 8, + -64, -1, -10, -10, 62, 62, -28, 8, -64, -1, -10, 8, -1, -1, 8, 71, 62, + 80, -10}, {8, -64, 80, 8, 512, -64, -1, 8, 71, 62, 80, -10, -10, -1, 62, + -28, -10, 62, -64, 8, 8, -1, -1, -10}, {80, 8, 8, -64, -64, 512, -10, -1, + 62, -28, -10, 62, -1, 8, 71, 62, 80, -10, 8, -64, -1, -10, 8, -1}, {-64, + 8, 8, -1, -1, -10, 512, -64, -64, 8, 8, 80, 80, -10, 8, -1, 62, 71, -10, + 62, -1, -10, -28, 62}, {8, -64, -1, -10, 8, -1, -64, 512, 8, 80, -64, 8, + -10, 62, -1, -10, -28, 62, 80, -10, 8, -1, 62, 71}, {8, -1, 80, -10, 71, + 62, -64, 8, 512, -64, 80, 8, 8, -1, -64, 8, -10, -1, 62, -28, -10, -1, + 62, -10}, {-1, -10, -10, 62, 62, -28, 8, 80, -64, 512, 8, -64, -1, -10, + 8, -64, -1, 8, 71, 62, -1, 8, -10, 80}, {-1, 8, 71, 62, 80, -10, 8, -64, + 80, 8, 512, -64, 62, -28, -10, -1, 62, -10, 8, -1, -64, 8, -10, -1}, + {-10, -1, 62, -28, -10, 62, 80, 8, 8, -64, -64, 512, 71, 62, -1, 8, -10, + 80, -1, -10, 8, -64, -1, 8}, {8, -1, -64, 8, -10, -1, 80, -10, 8, -1, 62, + 71, 512, -64, -64, 8, 8, 80, 62, -10, -28, 62, -1, -10}, {-1, -10, 8, + -64, -1, 8, -10, 62, -1, -10, -28, 62, -64, 512, 8, 80, -64, 8, -10, 80, + 62, 71, 8, -1}, {80, -10, 8, -1, 62, 71, 8, -1, -64, 8, -10, -1, -64, 8, + 512, -64, 80, 8, -28, 62, 62, -10, -10, -1}, {-10, 62, -1, -10, -28, 62, + -1, -10, 8, -64, -1, 8, 8, 80, -64, 512, 8, -64, 62, 71, -10, 80, -1, 8}, + {71, 62, -1, 8, -10, 80, 62, -28, -10, -1, 62, -10, 8, -64, 80, 8, 512, + -64, -1, 8, -10, -1, -64, 8}, {62, -28, -10, -1, 62, -10, 71, 62, -1, 8, + -10, 80, 80, 8, 8, -64, -64, 512, -10, -1, -1, 8, 8, -64}, {-1, 8, -10, + -1, -64, 8, -10, 80, 62, 71, 8, -1, 62, -10, -28, 62, -1, -10, 512, -64, + -64, 8, 8, 80}, {-10, -1, -1, 8, 8, -64, 62, -10, -28, 62, -1, -10, -10, + 80, 62, 71, 8, -1, -64, 512, 8, 80, -64, 8}, {-10, 80, 62, 71, 8, -1, -1, + 8, -10, -1, -64, 8, -28, 62, 62, -10, -10, -1, -64, 8, 512, -64, 80, 8}, + {62, -10, -28, 62, -1, -10, -10, -1, -1, 8, 8, -64, 62, 71, -10, 80, -1, + 8, 8, 80, -64, 512, 8, -64}, {62, 71, -10, 80, -1, 8, -28, 62, 62, -10, + -10, -1, -1, 8, -10, -1, -64, 8, 8, -64, 80, 8, 512, -64}, {-28, 62, 62, + -10, -10, -1, 62, 71, -10, 80, -1, 8, -10, -1, -1, 8, 8, -64, 80, 8, 8, + -64, -64, 512}}; + + // Calculate color flows + jamp[0] = +std::complex(0, 1) * amp[0] + + std::complex(0, 1) * amp[1] + + std::complex(0, 1) * amp[3] + + std::complex(0, 1) * amp[5] + + std::complex(0, 1) * amp[14] + amp[15] + amp[16] + amp[21] + + std::complex(0, 1) * amp[23] - amp[29] + + std::complex(0, 1) * amp[31] + amp[32] + amp[33] - amp[35] + + std::complex(0, 1) * amp[102] + + std::complex(0, 1) * amp[103] + + std::complex(0, 1) * amp[105] + + std::complex(0, 1) * amp[106] + amp[109] + + std::complex(0, 1) * amp[120] + + std::complex(0, 1) * amp[121] + + std::complex(0, 1) * amp[123] + + std::complex(0, 1) * amp[129] - + std::complex(0, 1) * amp[131] + amp[135] - amp[137] - + std::complex(0, 1) * amp[156] + + std::complex(0, 1) * amp[158]; + jamp[1] = -std::complex(0, 1) * amp[0] + + std::complex(0, 1) * amp[2] + + std::complex(0, 1) * amp[4] - + std::complex(0, 1) * amp[5] + + std::complex(0, 1) * amp[12] + amp[13] - amp[16] + amp[24] + + std::complex(0, 1) * amp[26] - amp[27] - + std::complex(0, 1) * amp[31] - amp[32] - amp[33] - amp[34] + + std::complex(0, 1) * amp[111] + + std::complex(0, 1) * amp[112] + + std::complex(0, 1) * amp[114] + + std::complex(0, 1) * amp[115] + amp[118] - + std::complex(0, 1) * amp[120] - + std::complex(0, 1) * amp[121] - + std::complex(0, 1) * amp[123] + + std::complex(0, 1) * amp[138] - + std::complex(0, 1) * amp[140] + amp[144] - amp[146] + + std::complex(0, 1) * amp[156] + + std::complex(0, 1) * amp[157]; + jamp[2] = -amp[21] - std::complex(0, 1) * amp[23] - amp[24] + + std::complex(0, 1) * amp[25] - amp[30] + amp[34] + amp[35] + + amp[60] - std::complex(0, 1) * amp[62] + + std::complex(0, 1) * amp[63] + + std::complex(0, 1) * amp[64] + + std::complex(0, 1) * amp[66] + amp[70] + + std::complex(0, 1) * amp[71] - + std::complex(0, 1) * amp[102] - + std::complex(0, 1) * amp[103] - + std::complex(0, 1) * amp[105] - + std::complex(0, 1) * amp[106] - amp[109] - + std::complex(0, 1) * amp[112] - + std::complex(0, 1) * amp[113] - + std::complex(0, 1) * amp[115] - + std::complex(0, 1) * amp[129] - + std::complex(0, 1) * amp[130] - amp[135] - amp[136] - + std::complex(0, 1) * amp[157] - + std::complex(0, 1) * amp[158]; + jamp[3] = -amp[18] + std::complex(0, 1) * amp[20] + amp[24] - + std::complex(0, 1) * amp[25] - amp[32] - amp[33] - amp[34] + + std::complex(0, 1) * amp[57] + amp[58] - amp[60] + + std::complex(0, 1) * amp[62] - + std::complex(0, 1) * amp[64] - + std::complex(0, 1) * amp[65] - + std::complex(0, 1) * amp[67] + amp[101] + + std::complex(0, 1) * amp[112] + + std::complex(0, 1) * amp[113] + + std::complex(0, 1) * amp[115] - + std::complex(0, 1) * amp[121] - + std::complex(0, 1) * amp[122] - + std::complex(0, 1) * amp[123] - + std::complex(0, 1) * amp[124] + amp[150] - amp[152] - + std::complex(0, 1) * amp[153] + + std::complex(0, 1) * amp[155] + + std::complex(0, 1) * amp[156] + + std::complex(0, 1) * amp[157]; + jamp[4] = -amp[21] + std::complex(0, 1) * amp[22] - amp[24] - + std::complex(0, 1) * amp[26] - amp[28] + amp[34] + amp[35] + + amp[78] - std::complex(0, 1) * amp[80] + + std::complex(0, 1) * amp[81] + + std::complex(0, 1) * amp[82] + + std::complex(0, 1) * amp[84] + amp[88] + + std::complex(0, 1) * amp[89] - + std::complex(0, 1) * amp[103] - + std::complex(0, 1) * amp[104] - + std::complex(0, 1) * amp[106] - + std::complex(0, 1) * amp[111] - + std::complex(0, 1) * amp[112] - + std::complex(0, 1) * amp[114] - + std::complex(0, 1) * amp[115] - amp[118] - + std::complex(0, 1) * amp[138] - + std::complex(0, 1) * amp[139] - amp[144] - amp[145] - + std::complex(0, 1) * amp[157] - + std::complex(0, 1) * amp[158]; + jamp[5] = -amp[19] - std::complex(0, 1) * amp[20] + amp[21] - + std::complex(0, 1) * amp[22] + amp[32] + amp[33] - amp[35] + + std::complex(0, 1) * amp[75] + amp[76] - amp[78] + + std::complex(0, 1) * amp[80] - + std::complex(0, 1) * amp[82] - + std::complex(0, 1) * amp[83] - + std::complex(0, 1) * amp[85] - amp[101] + + std::complex(0, 1) * amp[103] + + std::complex(0, 1) * amp[104] + + std::complex(0, 1) * amp[106] + + std::complex(0, 1) * amp[121] + + std::complex(0, 1) * amp[122] + + std::complex(0, 1) * amp[123] + + std::complex(0, 1) * amp[124] - amp[150] - amp[151] + + std::complex(0, 1) * amp[153] + + std::complex(0, 1) * amp[154] - + std::complex(0, 1) * amp[156] + + std::complex(0, 1) * amp[158]; + jamp[6] = -std::complex(0, 1) * amp[0] - + std::complex(0, 1) * amp[1] - + std::complex(0, 1) * amp[3] - + std::complex(0, 1) * amp[5] - + std::complex(0, 1) * amp[14] - amp[15] - amp[16] + amp[55] + + std::complex(0, 1) * amp[56] - + std::complex(0, 1) * amp[63] + + std::complex(0, 1) * amp[65] - + std::complex(0, 1) * amp[66] + + std::complex(0, 1) * amp[67] - amp[70] - amp[92] + + std::complex(0, 1) * amp[94] + amp[95] - + std::complex(0, 1) * amp[120] + + std::complex(0, 1) * amp[122] + + std::complex(0, 1) * amp[124] + + std::complex(0, 1) * amp[130] + + std::complex(0, 1) * amp[131] + amp[136] + amp[137] + + amp[147] - amp[149] + std::complex(0, 1) * amp[153] - + std::complex(0, 1) * amp[155]; + jamp[7] = +std::complex(0, 1) * amp[0] - + std::complex(0, 1) * amp[2] - + std::complex(0, 1) * amp[4] + + std::complex(0, 1) * amp[5] - + std::complex(0, 1) * amp[12] - amp[13] + amp[16] + amp[73] + + std::complex(0, 1) * amp[74] - + std::complex(0, 1) * amp[81] + + std::complex(0, 1) * amp[83] - + std::complex(0, 1) * amp[84] + + std::complex(0, 1) * amp[85] - amp[88] - amp[90] - + std::complex(0, 1) * amp[94] - amp[95] + + std::complex(0, 1) * amp[120] - + std::complex(0, 1) * amp[122] - + std::complex(0, 1) * amp[124] + + std::complex(0, 1) * amp[139] + + std::complex(0, 1) * amp[140] + amp[145] + amp[146] - + amp[147] - amp[148] - std::complex(0, 1) * amp[153] - + std::complex(0, 1) * amp[154]; + jamp[8] = -amp[55] - std::complex(0, 1) * amp[56] + + std::complex(0, 1) * amp[63] - + std::complex(0, 1) * amp[65] + + std::complex(0, 1) * amp[66] - + std::complex(0, 1) * amp[67] + amp[70] + + std::complex(0, 1) * amp[72] - amp[73] + amp[79] + + std::complex(0, 1) * amp[80] - + std::complex(0, 1) * amp[82] - + std::complex(0, 1) * amp[83] - + std::complex(0, 1) * amp[85] - amp[93] - + std::complex(0, 1) * amp[102] + + std::complex(0, 1) * amp[104] - + std::complex(0, 1) * amp[105] - amp[109] + + std::complex(0, 1) * amp[110] - + std::complex(0, 1) * amp[129] - + std::complex(0, 1) * amp[130] - amp[135] - amp[136] + + amp[148] + amp[149] + std::complex(0, 1) * amp[154] + + std::complex(0, 1) * amp[155]; + jamp[9] = -amp[37] + std::complex(0, 1) * amp[38] + amp[39] + + std::complex(0, 1) * amp[40] + amp[50] + amp[51] - amp[53] - + std::complex(0, 1) * amp[72] + amp[73] - amp[79] - + std::complex(0, 1) * amp[80] + + std::complex(0, 1) * amp[82] + + std::complex(0, 1) * amp[83] + + std::complex(0, 1) * amp[85] - amp[95] - + std::complex(0, 1) * amp[103] - + std::complex(0, 1) * amp[104] - + std::complex(0, 1) * amp[106] - + std::complex(0, 1) * amp[121] - + std::complex(0, 1) * amp[122] - + std::complex(0, 1) * amp[123] - + std::complex(0, 1) * amp[124] - amp[147] - amp[148] - + std::complex(0, 1) * amp[153] - + std::complex(0, 1) * amp[154] + + std::complex(0, 1) * amp[156] - + std::complex(0, 1) * amp[158]; + jamp[10] = +std::complex(0, 1) * amp[54] - amp[55] + amp[61] + + std::complex(0, 1) * amp[62] - + std::complex(0, 1) * amp[64] - + std::complex(0, 1) * amp[65] - + std::complex(0, 1) * amp[67] - amp[73] - + std::complex(0, 1) * amp[74] + + std::complex(0, 1) * amp[81] - + std::complex(0, 1) * amp[83] + + std::complex(0, 1) * amp[84] - + std::complex(0, 1) * amp[85] + amp[88] - amp[91] - + std::complex(0, 1) * amp[111] + + std::complex(0, 1) * amp[113] - + std::complex(0, 1) * amp[114] - amp[118] + + std::complex(0, 1) * amp[119] - + std::complex(0, 1) * amp[138] - + std::complex(0, 1) * amp[139] - amp[144] - amp[145] + + amp[148] + amp[149] + std::complex(0, 1) * amp[154] + + std::complex(0, 1) * amp[155]; + jamp[11] = -amp[36] - std::complex(0, 1) * amp[38] + amp[42] + + std::complex(0, 1) * amp[43] - amp[50] - amp[51] - + amp[52] - std::complex(0, 1) * amp[54] + amp[55] - + amp[61] - std::complex(0, 1) * amp[62] + + std::complex(0, 1) * amp[64] + + std::complex(0, 1) * amp[65] + + std::complex(0, 1) * amp[67] + amp[95] - + std::complex(0, 1) * amp[112] - + std::complex(0, 1) * amp[113] - + std::complex(0, 1) * amp[115] + + std::complex(0, 1) * amp[121] + + std::complex(0, 1) * amp[122] + + std::complex(0, 1) * amp[123] + + std::complex(0, 1) * amp[124] + amp[147] - amp[149] + + std::complex(0, 1) * amp[153] - + std::complex(0, 1) * amp[155] - + std::complex(0, 1) * amp[156] - + std::complex(0, 1) * amp[157]; + jamp[12] = -std::complex(0, 1) * amp[1] - + std::complex(0, 1) * amp[2] - + std::complex(0, 1) * amp[3] - + std::complex(0, 1) * amp[4] + amp[7] + + std::complex(0, 1) * amp[8] - amp[15] - amp[60] + + std::complex(0, 1) * amp[62] - + std::complex(0, 1) * amp[63] - + std::complex(0, 1) * amp[64] - + std::complex(0, 1) * amp[66] - amp[70] - + std::complex(0, 1) * amp[71] - + std::complex(0, 1) * amp[111] + + std::complex(0, 1) * amp[113] - + std::complex(0, 1) * amp[114] + amp[116] + + std::complex(0, 1) * amp[117] - amp[125] + + std::complex(0, 1) * amp[130] + + std::complex(0, 1) * amp[131] + amp[136] + amp[137] - + std::complex(0, 1) * amp[138] + + std::complex(0, 1) * amp[140] + amp[141] - amp[143]; + jamp[13] = -std::complex(0, 1) * amp[57] - amp[58] + amp[60] - + std::complex(0, 1) * amp[62] + + std::complex(0, 1) * amp[64] + + std::complex(0, 1) * amp[65] + + std::complex(0, 1) * amp[67] - amp[76] + + std::complex(0, 1) * amp[77] - + std::complex(0, 1) * amp[81] + + std::complex(0, 1) * amp[83] - + std::complex(0, 1) * amp[84] + + std::complex(0, 1) * amp[85] + amp[86] - amp[97] + + std::complex(0, 1) * amp[111] - + std::complex(0, 1) * amp[113] + + std::complex(0, 1) * amp[114] - amp[116] - + std::complex(0, 1) * amp[117] + + std::complex(0, 1) * amp[138] + + std::complex(0, 1) * amp[139] - amp[141] - amp[142] + + amp[151] + amp[152] - std::complex(0, 1) * amp[154] - + std::complex(0, 1) * amp[155]; + jamp[14] = +std::complex(0, 1) * amp[1] + + std::complex(0, 1) * amp[2] + + std::complex(0, 1) * amp[3] + + std::complex(0, 1) * amp[4] - amp[7] - + std::complex(0, 1) * amp[8] + amp[15] - amp[79] - + std::complex(0, 1) * amp[80] + + std::complex(0, 1) * amp[81] + + std::complex(0, 1) * amp[82] + + std::complex(0, 1) * amp[84] - amp[86] + + std::complex(0, 1) * amp[87] + + std::complex(0, 1) * amp[102] - + std::complex(0, 1) * amp[104] + + std::complex(0, 1) * amp[105] + amp[109] - + std::complex(0, 1) * amp[110] - amp[126] + + std::complex(0, 1) * amp[129] - + std::complex(0, 1) * amp[131] + amp[135] - amp[137] - + std::complex(0, 1) * amp[139] - + std::complex(0, 1) * amp[140] + amp[142] + amp[143]; + jamp[15] = -amp[39] - std::complex(0, 1) * amp[40] - amp[42] + + std::complex(0, 1) * amp[44] - amp[46] + amp[52] + + amp[53] + amp[79] + std::complex(0, 1) * amp[80] - + std::complex(0, 1) * amp[81] - + std::complex(0, 1) * amp[82] - + std::complex(0, 1) * amp[84] + amp[86] - + std::complex(0, 1) * amp[87] + + std::complex(0, 1) * amp[103] + + std::complex(0, 1) * amp[104] + + std::complex(0, 1) * amp[106] + + std::complex(0, 1) * amp[111] + + std::complex(0, 1) * amp[112] + + std::complex(0, 1) * amp[114] + + std::complex(0, 1) * amp[115] - amp[116] + + std::complex(0, 1) * amp[138] + + std::complex(0, 1) * amp[139] - amp[141] - amp[142] + + std::complex(0, 1) * amp[157] + + std::complex(0, 1) * amp[158]; + jamp[16] = -std::complex(0, 1) * amp[0] + + std::complex(0, 1) * amp[2] + + std::complex(0, 1) * amp[4] - + std::complex(0, 1) * amp[5] + + std::complex(0, 1) * amp[6] - amp[7] + amp[17] + amp[76] - + std::complex(0, 1) * amp[77] + + std::complex(0, 1) * amp[81] - + std::complex(0, 1) * amp[83] + + std::complex(0, 1) * amp[84] - + std::complex(0, 1) * amp[85] - amp[86] - amp[96] + + std::complex(0, 1) * amp[100] - amp[101] - + std::complex(0, 1) * amp[120] + + std::complex(0, 1) * amp[122] + + std::complex(0, 1) * amp[124] - + std::complex(0, 1) * amp[139] - + std::complex(0, 1) * amp[140] + amp[142] + amp[143] - + amp[150] - amp[151] + std::complex(0, 1) * amp[153] + + std::complex(0, 1) * amp[154]; + jamp[17] = +std::complex(0, 1) * amp[0] - + std::complex(0, 1) * amp[2] - + std::complex(0, 1) * amp[4] + + std::complex(0, 1) * amp[5] - + std::complex(0, 1) * amp[6] + amp[7] - amp[17] + amp[42] - + std::complex(0, 1) * amp[44] - amp[45] + + std::complex(0, 1) * amp[49] - amp[50] - amp[51] - + amp[52] - std::complex(0, 1) * amp[111] - + std::complex(0, 1) * amp[112] - + std::complex(0, 1) * amp[114] - + std::complex(0, 1) * amp[115] + amp[116] + + std::complex(0, 1) * amp[120] + + std::complex(0, 1) * amp[121] + + std::complex(0, 1) * amp[123] - + std::complex(0, 1) * amp[138] + + std::complex(0, 1) * amp[140] + amp[141] - amp[143] - + std::complex(0, 1) * amp[156] - + std::complex(0, 1) * amp[157]; + jamp[18] = -std::complex(0, 1) * amp[1] - + std::complex(0, 1) * amp[2] - + std::complex(0, 1) * amp[3] - + std::complex(0, 1) * amp[4] + amp[10] + + std::complex(0, 1) * amp[11] - amp[13] - amp[78] + + std::complex(0, 1) * amp[80] - + std::complex(0, 1) * amp[81] - + std::complex(0, 1) * amp[82] - + std::complex(0, 1) * amp[84] - amp[88] - + std::complex(0, 1) * amp[89] - + std::complex(0, 1) * amp[102] + + std::complex(0, 1) * amp[104] - + std::complex(0, 1) * amp[105] + amp[107] + + std::complex(0, 1) * amp[108] - amp[127] - + std::complex(0, 1) * amp[129] + + std::complex(0, 1) * amp[131] + amp[132] - amp[134] + + std::complex(0, 1) * amp[139] + + std::complex(0, 1) * amp[140] + amp[145] + amp[146]; + jamp[19] = -amp[58] + std::complex(0, 1) * amp[59] - + std::complex(0, 1) * amp[63] + + std::complex(0, 1) * amp[65] - + std::complex(0, 1) * amp[66] + + std::complex(0, 1) * amp[67] + amp[68] - + std::complex(0, 1) * amp[75] - amp[76] + amp[78] - + std::complex(0, 1) * amp[80] + + std::complex(0, 1) * amp[82] + + std::complex(0, 1) * amp[83] + + std::complex(0, 1) * amp[85] - amp[99] + + std::complex(0, 1) * amp[102] - + std::complex(0, 1) * amp[104] + + std::complex(0, 1) * amp[105] - amp[107] - + std::complex(0, 1) * amp[108] + + std::complex(0, 1) * amp[129] + + std::complex(0, 1) * amp[130] - amp[132] - amp[133] + + amp[151] + amp[152] - std::complex(0, 1) * amp[154] - + std::complex(0, 1) * amp[155]; + jamp[20] = +std::complex(0, 1) * amp[1] + + std::complex(0, 1) * amp[2] + + std::complex(0, 1) * amp[3] + + std::complex(0, 1) * amp[4] - amp[10] - + std::complex(0, 1) * amp[11] + amp[13] - amp[61] - + std::complex(0, 1) * amp[62] + + std::complex(0, 1) * amp[63] + + std::complex(0, 1) * amp[64] + + std::complex(0, 1) * amp[66] - amp[68] + + std::complex(0, 1) * amp[69] + + std::complex(0, 1) * amp[111] - + std::complex(0, 1) * amp[113] + + std::complex(0, 1) * amp[114] + amp[118] - + std::complex(0, 1) * amp[119] - amp[128] - + std::complex(0, 1) * amp[130] - + std::complex(0, 1) * amp[131] + amp[133] + amp[134] + + std::complex(0, 1) * amp[138] - + std::complex(0, 1) * amp[140] + amp[144] - amp[146]; + jamp[21] = -amp[39] + std::complex(0, 1) * amp[41] - amp[42] - + std::complex(0, 1) * amp[43] - amp[48] + amp[52] + + amp[53] + amp[61] + std::complex(0, 1) * amp[62] - + std::complex(0, 1) * amp[63] - + std::complex(0, 1) * amp[64] - + std::complex(0, 1) * amp[66] + amp[68] - + std::complex(0, 1) * amp[69] + + std::complex(0, 1) * amp[102] + + std::complex(0, 1) * amp[103] + + std::complex(0, 1) * amp[105] + + std::complex(0, 1) * amp[106] - amp[107] + + std::complex(0, 1) * amp[112] + + std::complex(0, 1) * amp[113] + + std::complex(0, 1) * amp[115] + + std::complex(0, 1) * amp[129] + + std::complex(0, 1) * amp[130] - amp[132] - amp[133] + + std::complex(0, 1) * amp[157] + + std::complex(0, 1) * amp[158]; + jamp[22] = +std::complex(0, 1) * amp[0] + + std::complex(0, 1) * amp[1] + + std::complex(0, 1) * amp[3] + + std::complex(0, 1) * amp[5] + + std::complex(0, 1) * amp[9] - amp[10] - amp[17] + amp[58] - + std::complex(0, 1) * amp[59] + + std::complex(0, 1) * amp[63] - + std::complex(0, 1) * amp[65] + + std::complex(0, 1) * amp[66] - + std::complex(0, 1) * amp[67] - amp[68] - amp[98] - + std::complex(0, 1) * amp[100] + amp[101] + + std::complex(0, 1) * amp[120] - + std::complex(0, 1) * amp[122] - + std::complex(0, 1) * amp[124] - + std::complex(0, 1) * amp[130] - + std::complex(0, 1) * amp[131] + amp[133] + amp[134] + + amp[150] - amp[152] - std::complex(0, 1) * amp[153] + + std::complex(0, 1) * amp[155]; + jamp[23] = -std::complex(0, 1) * amp[0] - + std::complex(0, 1) * amp[1] - + std::complex(0, 1) * amp[3] - + std::complex(0, 1) * amp[5] - + std::complex(0, 1) * amp[9] + amp[10] + amp[17] + amp[39] - + std::complex(0, 1) * amp[41] - amp[47] - + std::complex(0, 1) * amp[49] + amp[50] + amp[51] - + amp[53] - std::complex(0, 1) * amp[102] - + std::complex(0, 1) * amp[103] - + std::complex(0, 1) * amp[105] - + std::complex(0, 1) * amp[106] + amp[107] - + std::complex(0, 1) * amp[120] - + std::complex(0, 1) * amp[121] - + std::complex(0, 1) * amp[123] - + std::complex(0, 1) * amp[129] + + std::complex(0, 1) * amp[131] + amp[132] - amp[134] + + std::complex(0, 1) * amp[156] - + std::complex(0, 1) * amp[158]; + + // Sum and square the color flows to get the matrix element + for(i = 0; i < ncolor; i++ ) + { + ztemp = 0.; + for(j = 0; j < ncolor; j++ ) + ztemp = ztemp + cf[i][j] * jamp[j]; + /* + DPCT1007:3: Migration of this CUDA API is not supported by the Intel(R) + DPC++ Compatibility Tool. + */ + matrix = matrix + (ztemp * conj(jamp[i])).real() / denom[i]; + } + + // Store the leading color flows for choice of color + // for(i=0;i < ncolor; i++) + // jamp2[0][i] += real(jamp[i]*conj(jamp[i])); + +} + + +//-------------------------------------------------------------------------- +// Evaluate |M|^2, part independent of incoming flavour. + +class sigmaKin { + + public: + // constructor + sigmaKin(double *allmomenta, + double *output, + int **cHel, + double *cIPC, + double *cIPD + ) + : m_allmomenta(allmomenta), + m_output(output), + m_cHel(cHel), + m_cIPC(cIPC), + m_cIPD(cIPD) + {} + + void operator()(sycl::nd_item<3> item) const { + + // Set the parameters which change event by event + // Need to discuss this with Stefan + // pars->setDependentParameters(); + // pars->setDependentCouplings(); + + // Reset color flows + + // for (int xx = 0; xx < 384; ++xx) { + const int nprocesses = 1; + int tid = item.get_group(2) * item.get_local_range().get(2) + item.get_local_id(2); + + // char *devPtr = (char *)tp.ptr; + // size_t dpt = tp.pitch; + // size_t slicePitch = dpt * 6; + + // char *dps = devPtr + dim * slicePitch; + double matrix_element[nprocesses]; + + std::complex amp[159]; + + double local_m[6][3]; + + int DIM = item.get_local_range().get(2) * item.get_group_range().get(2); + // for (int i=0; i<20;i++){ + // printf(" %f ", allmomenta[i]); + // } + // printf("\n"); + // printf("DIM is %i/%i\n", tid, DIM); + for (int i = 0; i < 6; i++ ) + { + for (int j = 0; j < 3; j++ ) + { + local_m[i][j] = m_allmomenta[i * 3 * DIM + j * DIM + tid]; + // printf(" %f ", local_m[i][j]); + } + // printf("\n"); + } + + + // Local variables and constants + const int ncomb = 64; + // static bool goodhel[ncomb] = {ncomb * false}; + // static int ntry = 0, sum_hel = 0, ngood = 0; + // static int igood[ncomb]; + // static int jhel; + // std::complex **wfs; + // double t[1]; + // Helicities for the process + // static const int helicities[ncomb][nexternal] = + // {{-1,-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1,1},{-1,-1,-1,-1,1,-1},{-1,-1,-1,-1,1,1 + // },{-1,-1,-1,1,-1,-1},{-1,-1,-1,1,-1,1},{-1,-1,-1,1,1,-1},{-1,-1,-1,1,1,1},{ + // -1,-1,1,-1,-1,-1},{-1,-1,1,-1,-1,1},{-1,-1,1,-1,1,-1},{-1,-1,1,-1,1,1},{-1, + // -1,1,1,-1,-1},{-1,-1,1,1,-1,1},{-1,-1,1,1,1,-1},{-1,-1,1,1,1,1},{-1,1,-1,-1 + // ,-1,-1},{-1,1,-1,-1,-1,1},{-1,1,-1,-1,1,-1},{-1,1,-1,-1,1,1},{-1,1,-1,1,-1, + // -1},{-1,1,-1,1,-1,1},{-1,1,-1,1,1,-1},{-1,1,-1,1,1,1},{-1,1,1,-1,-1,-1},{-1 + // ,1,1,-1,-1,1},{-1,1,1,-1,1,-1},{-1,1,1,-1,1,1},{-1,1,1,1,-1,-1},{-1,1,1,1,- + // 1,1},{-1,1,1,1,1,-1},{-1,1,1,1,1,1},{1,-1,-1,-1,-1,-1},{1,-1,-1,-1,-1,1},{1 + // ,-1,-1,-1,1,-1},{1,-1,-1,-1,1,1},{1,-1,-1,1,-1,-1},{1,-1,-1,1,-1,1},{1,-1,- + // 1,1,1,-1},{1,-1,-1,1,1,1},{1,-1,1,-1,-1,-1},{1,-1,1,-1,-1,1},{1,-1,1,-1,1,- + // 1},{1,-1,1,-1,1,1},{1,-1,1,1,-1,-1},{1,-1,1,1,-1,1},{1,-1,1,1,1,-1},{1,-1,1 + // ,1,1,1},{1,1,-1,-1,-1,-1},{1,1,-1,-1,-1,1},{1,1,-1,-1,1,-1},{1,1,-1,-1,1,1} + // ,{1,1,-1,1,-1,-1},{1,1,-1,1,-1,1},{1,1,-1,1,1,-1},{1,1,-1,1,1,1},{1,1,1,-1, + // -1,-1},{1,1,1,-1,-1,1},{1,1,1,-1,1,-1},{1,1,1,-1,1,1},{1,1,1,1,-1,-1},{1,1, + // 1,1,-1,1},{1,1,1,1,1,-1},{1,1,1,1,1,1}}; + // Denominators: spins, colors and identical particles + const int denominators[1] = {512}; + + + // Reset the matrix elements + for(int i = 0; i < nprocesses; i++ ) + { + matrix_element[i] = 0.; + } + // Define permutation + // int perm[nexternal]; + // for(int i = 0; i < nexternal; i++){ + // perm[i]=i; + // } + + + for (int ihel = 0; ihel < ncomb; ihel++ ) + { + calculate_wavefunctions(ihel, local_m, matrix_element[0], m_cHel, m_cIPC, m_cIPD); + } + + for (int i = 0; i < nprocesses; ++ i) + { + matrix_element[i] /= denominators[i]; + } + for (int i = 0; i < nprocesses; ++ i) + { + m_output[i * nprocesses + tid] = matrix_element[i]; + // printf("output %i %i %i %f", tid, i, i*nprocesses+tid, + // output[i*nprocesses+tid]); + + } + } + + private: + double *m_allmomenta; + double *m_output; + int **m_cHel; + double *m_cIPC; + double *m_cIPD; +}; + +class CPPProcess +{ + public: + + CPPProcess(int numiterations, int gpublocks, int gputhreads, + bool verbose, bool debug): + m_numiterations(numiterations), + gpu_nblocks(gpublocks), + gpu_nthreads(gputhreads), + dim(gpu_nblocks * gpu_nthreads) + { + + + // Helicities for the process - nodim + static const int tHel[ncomb][nexternal] = {{-1, -1, -1, -1, -1, -1}, {-1, -1, + -1, -1, -1, 1}, {-1, -1, -1, -1, 1, -1}, {-1, -1, -1, -1, 1, 1}, {-1, -1, + -1, 1, -1, -1}, {-1, -1, -1, 1, -1, 1}, {-1, -1, -1, 1, 1, -1}, {-1, -1, + -1, 1, 1, 1}, {-1, -1, 1, -1, -1, -1}, {-1, -1, 1, -1, -1, 1}, {-1, -1, + 1, -1, 1, -1}, {-1, -1, 1, -1, 1, 1}, {-1, -1, 1, 1, -1, -1}, {-1, -1, 1, + 1, -1, 1}, {-1, -1, 1, 1, 1, -1}, {-1, -1, 1, 1, 1, 1}, {-1, 1, -1, -1, + -1, -1}, {-1, 1, -1, -1, -1, 1}, {-1, 1, -1, -1, 1, -1}, {-1, 1, -1, -1, + 1, 1}, {-1, 1, -1, 1, -1, -1}, {-1, 1, -1, 1, -1, 1}, {-1, 1, -1, 1, 1, + -1}, {-1, 1, -1, 1, 1, 1}, {-1, 1, 1, -1, -1, -1}, {-1, 1, 1, -1, -1, 1}, + {-1, 1, 1, -1, 1, -1}, {-1, 1, 1, -1, 1, 1}, {-1, 1, 1, 1, -1, -1}, {-1, + 1, 1, 1, -1, 1}, {-1, 1, 1, 1, 1, -1}, {-1, 1, 1, 1, 1, 1}, {1, -1, -1, + -1, -1, -1}, {1, -1, -1, -1, -1, 1}, {1, -1, -1, -1, 1, -1}, {1, -1, -1, + -1, 1, 1}, {1, -1, -1, 1, -1, -1}, {1, -1, -1, 1, -1, 1}, {1, -1, -1, 1, + 1, -1}, {1, -1, -1, 1, 1, 1}, {1, -1, 1, -1, -1, -1}, {1, -1, 1, -1, -1, + 1}, {1, -1, 1, -1, 1, -1}, {1, -1, 1, -1, 1, 1}, {1, -1, 1, 1, -1, -1}, + {1, -1, 1, 1, -1, 1}, {1, -1, 1, 1, 1, -1}, {1, -1, 1, 1, 1, 1}, {1, 1, + -1, -1, -1, -1}, {1, 1, -1, -1, -1, 1}, {1, 1, -1, -1, 1, -1}, {1, 1, -1, + -1, 1, 1}, {1, 1, -1, 1, -1, -1}, {1, 1, -1, 1, -1, 1}, {1, 1, -1, 1, 1, + -1}, {1, 1, -1, 1, 1, 1}, {1, 1, 1, -1, -1, -1}, {1, 1, 1, -1, -1, 1}, + {1, 1, 1, -1, 1, -1}, {1, 1, 1, -1, 1, 1}, {1, 1, 1, 1, -1, -1}, {1, 1, + 1, 1, -1, 1}, {1, 1, 1, 1, 1, -1}, {1, 1, 1, 1, 1, 1}}; + // perm - nodim + // static int perm[nexternal] = {0, 1, 2, 3}; + } + + ~CPPProcess() = default; + + const std::vector &getMasses() const {return mME;} + + //-------------------------------------------------------------------------- + // Initialize process. + void initProc(string param_card_name) + { + sycl::device dev_ct1 = sycl::device(); + sycl::queue q_ct1 = sycl::queue(); + // Instantiate the model class and set parameters that stay fixed during run + pars = Parameters_sm::getInstance(); + SLHAReader slha(param_card_name); + pars->setIndependentParameters(slha); + pars->setIndependentCouplings(); + pars->printIndependentParameters(); + pars->printIndependentCouplings(); + pars->setDependentParameters(); + pars->setDependentCouplings(); + // Set external particle masses for this matrix element + mME.push_back(pars->ZERO); + mME.push_back(pars->ZERO); + mME.push_back(pars->mdl_MT); + mME.push_back(pars->mdl_MT); + mME.push_back(pars->ZERO); + mME.push_back(pars->ZERO); + static std::complex tIPC[3] = {pars->GC_10, pars->GC_11, pars->GC_12}; + static double tIPD[2] = {pars->mdl_MT, pars->mdl_WT}; + + } + + /* void setInitial(int inid1, int inid2) */ + /* { */ + /* id1 = inid1; */ + /* id2 = inid2; */ + /* } */ + + /* int getDim() const {return dim;} */ + + /* int getNIOParticles() const {return nexternal;} */ + + + // Constants for array limits + static const int ninitial = 2; + static const int nexternal = 6; + static const int nprocesses = 1; + + + + private: + int m_numiterations; + // gpu variables + int gpu_nblocks; + int gpu_nthreads; + int dim; // gpu_nblocks * gpu_nthreads; + + // print verbose info + bool m_verbose; + + // print debug info + bool m_debug; + + static const int nwavefuncs = 6; + static const int namplitudes = 159; + static const int ncomb = 64; + static const int wrows = 63; + // static const int nioparticles = 6; + + std::complex **amp; + + // Pointer to the model parameters + Parameters_sm * pars; + + // vector with external particle masses + std::vector mME; + + // Initial particle ids + int id1, id2; + +}; + + +#endif // MG5_Sigma_sm_gg_ttxgg_H diff --git a/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/Makefile b/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/Makefile new file mode 100644 index 0000000000..2b4f760488 --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/Makefile @@ -0,0 +1,60 @@ +#CUARCHNUM=70 +LIBDIR=../../lib +INCDIR=../../src +MODELLIB=model_sm +CXXFLAGS= -O3 -I$(INCDIR) -I. +#CUARCHFLAGS= -arch=compute_$(CUARCHNUM) #--gpu-architecture=compute_35 --gpu-code=sm_35 +#CUFLAGS= $(CUARCHFLAGS) -use_fast_math -lineinfo +LIBFLAGS= -L$(LIBDIR) -l$(MODELLIB) +CXX=dpcpp +MAKEDEBUG= + +main=check.exe +#cxx_objects=check_sa.o +cxx_objects=CPPProcess.o check_sa.o + +all: check + +check: $(main) + +debug: CXXFLAGS:=$(filter-out -O3,$(CXXFLAGS)) +debug: CXXFLAGS += -g -O0 -DDEBUG2 +debug: MAKEDEBUG := debug +debug: check + +$(LIBDIR)/lib$(MODELLIB).a: + @cd ../../src && make $(MAKEDEBUG) + +check_sa.o: check_sa.cpp + $(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $< -o $@ + +%.o : %.cpp %.h + $(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $< -o $@ + +$(main): $(LIBDIR)/lib$(MODELLIB).a $(cxx_objects) $(cu_objects) + $(CXX) -o $@ $(cxx_objects) $(LIBFLAGS) -lsycl + +.PHONY: clean + +clean: + cd ../../src && make clean + rm -f $(main) + rm -f $(cxx_objects) + +memcheck: check + /usr/local/cuda/bin/cuda-memcheck --check-api-memory-access yes --check-deprecated-instr yes --check-device-heap yes --demangle full --language c --leak-check full --racecheck-report all --report-api-errors all --show-backtrace yes --tool memcheck --track-unused-memory yes ./check.exe 5 5 10 + +perf: force + make clean && make + time ./check.exe -p 1024 16 384 && date + +test: force + ./check.exe -v 1 1 10 + +force: + +#Allowed values for this option: 'compute_30', 'compute_32', 'compute_35', 'compute_37', 'compute_50', 'compute_52', 'compute_53', 'compute_60', 'compute_61', 'compute_62', 'compute_70', 'compute_72', 'compute_75', 'sm_30', 'sm_32', 'sm_35', 'sm_37', 'sm_50', 'sm_52', 'sm_53', 'sm_60', 'sm_61', 'sm_62', 'sm_70', 'sm_72', 'sm_75'. + +# Max compute architectures +# cern batch (tesla v100): 70 +# jetson nano (maxwell): 35 diff --git a/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/check_sa.cpp b/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/check_sa.cpp new file mode 100644 index 0000000000..1627e0be55 --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/check_sa.cpp @@ -0,0 +1,286 @@ +#include +#include // perf stats +#include +#include +#include +#include // perf stats +#include +#include + +#include "CPPProcess.h" +#include "HelAmps_sm.h" + +#include "rambo.h" +#include "timer.h" +#include + +#define gpuErrchk3(ans) \ + { gpuAssert3((ans), __FILE__, __LINE__); } + +class sycl_kernel; + +inline void gpuAssert3(int code, const char *file, int line, + bool abort = true) { +} + +#define TIMERTYPE std::chrono::high_resolution_clock + +bool is_number(const char *s) { + const char *t = s; + while (*t != '\0' && isdigit(*t)) + ++t; + return strlen(s) == t - s; +} + +int usage(char* argv0, int ret = 1) { + std::cout << "Usage: " << argv0 + << " [--verbose|-v] [--debug|-d] [--performance|-p]" + << " [#gpuBlocksPerGrid #gpuThreadsPerBlock] #iterations" << std::endl; + return ret; +} + +int main(int argc, char **argv) { + sycl::cpu_selector selector; + sycl::queue q_ct1(selector); + sycl::device dev_ct1 = q_ct1.get_device();; + std::cout << "Device: " << dev_ct1.get_info() << "\n"; + bool verbose = false, debug = false, perf = false; + int numiter = 0, gpublocks = 1, gputhreads = 1; + std::vector numvec; + Timer timer; + std::vector wavetimes; + + + for (int argn = 1; argn < argc; ++argn) { + if (strcmp(argv[argn], "--verbose") == 0 || strcmp(argv[argn], "-v") == 0) + verbose = true; + else if (strcmp(argv[argn], "--debug") == 0 || + strcmp(argv[argn], "-d") == 0) + debug = true; + else if (strcmp(argv[argn], "--performance") == 0 || + strcmp(argv[argn], "-p") == 0) + perf = true; + else if (is_number(argv[argn])) + numvec.push_back(atoi(argv[argn])); + else + return usage(argv[0]); + } + int veclen = numvec.size(); + if (veclen == 3) { + gpublocks = numvec[0]; + gputhreads = numvec[1]; + numiter = numvec[2]; + } else if (veclen == 1) { + numiter = numvec[0]; + } else { + return usage(argv[0]); + } + + if (numiter == 0) + return usage(argv[0]); + + if (verbose) + std::cout << "# iterations: " << numiter << std::endl; + + // Create a process object + CPPProcess process(numiter, gpublocks, gputhreads, verbose, debug); + + // Read param_card and set parameters + process.initProc("../../Cards/param_card.dat"); + + double energy = 1500; + double weight; + + int meGeVexponent = -(2 * process.nexternal - 8); + + int dim = gpublocks * gputhreads; + + // Local Memory + //typedef double arr_t[6][4]; + double* lp = new double[6*3*dim]; + + std::vector matrixelementvector; + + + + for (int x = 0; x < numiter; ++x) { + // Get phase space point + std::vector > p = get_momenta(process.ninitial, energy, process.getMasses(), weight, dim); + + // Set momenta for this event + for (int d = 0; d < dim; ++d) { + for (int i = 0; i < 6; ++i) { + for (int j = 0; j < 3; ++j) { + lp[i*dim*3+j*dim+d] = p[d][i][1+j]; + } + } + } + //new + double * allmomenta = sycl::malloc_shared(3*6*dim*sizeof(double) ,q_ct1); + double * me = sycl::malloc_shared(dim*sizeof(double), q_ct1); + + if (perf) { + timer.Start(); + } + + // Evaluate matrix element + // later process.sigmaKin(ncomb, goodhel, ntry, sum_hel, ngood, igood, + // jhel); + + /* Setup Kernel Variables */ + sycl::nd_range nd_range = sycl::nd_range<3>(sycl::range<3>(1, 1, gpublocks) * + sycl::range<3>(1, 1, gputhreads), + sycl::range<3>(1, 1, gputhreads)); + + // USM Allocations, copy host -> device + const int cHEL_ROW = 64; + const int cHEL_COL = 6; + int **cHel_local = sycl::malloc_device(cHEL_ROW *sizeof(int*), q_ct1); + // allocate elements using host pointers, copy + for (int i=0; i< cHEL_ROW ; i++){ + cHel_local[i] = sycl::malloc_device(cHEL_COL * sizeof(int), q_ct1); + q_ct1.memcpy(cHel_local[i], &cHel[i][0], cHEL_COL * sizeof(int)).wait_and_throw(); + } + // copy host-side pointers + q_ct1.memcpy(cHel_local, &cHel[0], cHEL_ROW * cHEL_COL * sizeof(int)).wait_and_throw(); + + double *cIPC_local = sycl::malloc_device(sizeof(cIPC), q_ct1); + q_ct1.memcpy(cIPC_local, cIPC, sizeof(cIPC)).wait_and_throw(); + + double *cIPD_local = sycl::malloc_device(sizeof(cIPD), q_ct1); + q_ct1.memcpy(cIPD_local, cIPD, sizeof(cIPD)).wait_and_throw(); + + try{ + q_ct1.submit([&](sycl::handler &cgh) { + // instantiate kernel + auto k = sigmaKin(allmomenta, + me, + cHel_local, + cIPC_local, + cIPD_local + ); + + cgh.parallel_for(nd_range, k); + + }).wait_and_throw(); //, debug, verbose); + } catch(const sycl::exception& e){ + std::cout << "error: " << e.what(); + } + + // Shared matrix element means no need to memcpy to host + + if (verbose) + std::cout << "***********************************" << std::endl + << "Iteration #" << x+1 << " of " << numiter << std::endl; + + if (perf) { + float gputime = timer.GetDuration(); + wavetimes.push_back(gputime); + if (verbose) + std::cout << "Wave function time: " << gputime << std::endl; + } + + if (verbose || perf) { + + for (int d = 0; d < dim; ++d) { + + if (verbose) { + std::cout << "Momenta:" << std::endl; + for (int i = 0; i < process.nexternal; i++) + std::cout << std::setw(4) << i + 1 + << setiosflags(std::ios::scientific) << std::setw(14) + << p[d][i][0] << setiosflags(std::ios::scientific) + << std::setw(14) << p[d][i][1] + << setiosflags(std::ios::scientific) << std::setw(14) + << p[d][i][2] << setiosflags(std::ios::scientific) + << std::setw(14) << p[d][i][3] << std::endl; + std::cout << std::string(80, '-') << std::endl; + } + + // Display matrix elements + for (int i = 0; i < process.nprocesses; i++) { + if (verbose) + std::cout << " Matrix element = " + // << setiosflags(ios::fixed) << setprecision(17) + << me[i*1 + d] << " GeV^" << meGeVexponent << std::endl; + if (perf) + matrixelementvector.push_back(me[i*1 + d]); + } + + if (verbose) + std::cout << std::string(80, '-') << std::endl; + } + } else if (!debug) { + std::cout << "."; + } + + for (std::vector >::iterator it = p.begin(); + it != p.end(); ++it) { + for (std::vector::iterator jt = it->begin(); jt != it->end(); + ++jt) { + delete[] & (**jt); + } + } + } + + if (!(verbose || debug || perf)) { + std::cout << std::endl; + } + + if (perf) { + float sum = std::accumulate(wavetimes.begin(), wavetimes.end(), 0.0); + int num_wts = wavetimes.size(); + float mean = sum / num_wts; + float sq_sum = std::inner_product(wavetimes.begin(), wavetimes.end(), + wavetimes.begin(), 0.0); + float stdev = std::sqrt(sq_sum / num_wts - mean * mean); + std::vector::iterator mintime = + std::min_element(wavetimes.begin(), wavetimes.end()); + std::vector::iterator maxtime = + std::max_element(wavetimes.begin(), wavetimes.end()); + + int num_mes = matrixelementvector.size(); + float sumelem = std::accumulate(matrixelementvector.begin(), matrixelementvector.end(), 0.0); + float meanelem = sumelem / num_mes; + float sqselem = std::inner_product(matrixelementvector.begin(), matrixelementvector.end(), + matrixelementvector.begin(), 0.0); + float stdelem = std::sqrt(sqselem / num_mes - meanelem * meanelem); + std::vector::iterator maxelem = std::max_element( + matrixelementvector.begin(), matrixelementvector.end()); + std::vector::iterator minelem = std::min_element( + matrixelementvector.begin(), matrixelementvector.end()); + + std::cout << "***********************************" << std::endl + << "NumIterations = " << numiter << std::endl + << "NumThreadsPerBlock = " << gputhreads << std::endl + << "NumBlocksPerGrid = " << gpublocks << std::endl + << "-----------------------------------" << std::endl + << "NumberOfEntries = " << num_wts << std::endl + << std::scientific + << "TotalTimeInWaveFuncs = " << sum << " sec" << std::endl + << "MeanTimeInWaveFuncs = " << mean << " sec" << std::endl + << "StdDevTimeInWaveFuncs = " << stdev << " sec" << std::endl + << "MinTimeInWaveFuncs = " << *mintime << " sec" << std::endl + << "MaxTimeInWaveFuncs = " << *maxtime << " sec" << std::endl + << "-----------------------------------" << std::endl + << "ProcessID: = " << getpid() << std::endl + << "NProcesses = " << process.nprocesses << std::endl + << "NumMatrixElements = " << num_mes << std::endl + << "MatrixElementsPerSec = " << num_mes/sum << " sec^-1" << std::endl; + + std::cout << "***********************************" << std::endl + << "NumMatrixElements = " << num_mes << std::endl + << std::scientific << "MeanMatrixElemValue = " << meanelem + << " GeV^" << meGeVexponent << std::endl + << "StdErrMatrixElemValue = " << stdelem / sqrt(num_mes) + << " GeV^" << meGeVexponent << std::endl + << "StdDevMatrixElemValue = " << stdelem << " GeV^" + << meGeVexponent << std::endl + << "MinMatrixElemValue = " << *minelem << " GeV^" + << meGeVexponent << std::endl + << "MaxMatrixElemValue = " << *maxelem << " GeV^" + << meGeVexponent << std::endl; + } + delete[] lp; + +} diff --git a/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/timer.h b/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/timer.h new file mode 100644 index 0000000000..8cb692bd0a --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/SubProcesses/P1_Sigma_sm_gg_ttxgg/timer.h @@ -0,0 +1,52 @@ +#include +#include + +/* +high_resolution_clock +steady_clock +system_clock + +from https://www.modernescpp.com/index.php/the-three-clocks +and https://codereview.stackexchange.com/questions/196245/extremely-simple-timer-class-in-c +*/ + +template class Timer { + + public: + void Start(); + float GetDuration(); + void Info(); + + private: + typedef typename T::time_point TTP; + TTP m_StartTime; + +}; + +template void Timer::Start() +{ + m_StartTime = T::now(); +} + +template float Timer::GetDuration() +{ + std::chrono::duration duration = T::now() - m_StartTime; + return duration.count(); +} + + +template void Timer::Info() { + typedef typename T::period TPER; + typedef typename std::ratio_multiply MilliSec; + typedef typename std::ratio_multiply MicroSec; + + std::cout << std::boolalpha << std::endl; + std::cout << "clock info: " << std::endl; + std::cout << " is steady: " << T::is_steady << std::endl; + std::cout << " precision: " << TPER::num << "/" << TPER::den << " second " << std::endl; + std::cout << std::fixed; + std::cout << " " << static_cast(MilliSec::num)/MilliSec::den << " milliseconds " << std::endl; + std::cout << " " << static_cast(MicroSec::num)/MicroSec::den << " microseconds " << std::endl; + std::cout << std::endl; +} + diff --git a/epoch1/oneapi/gg_ttgg/SubProcesses/timer.h b/epoch1/oneapi/gg_ttgg/SubProcesses/timer.h new file mode 100644 index 0000000000..8cb692bd0a --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/SubProcesses/timer.h @@ -0,0 +1,52 @@ +#include +#include + +/* +high_resolution_clock +steady_clock +system_clock + +from https://www.modernescpp.com/index.php/the-three-clocks +and https://codereview.stackexchange.com/questions/196245/extremely-simple-timer-class-in-c +*/ + +template class Timer { + + public: + void Start(); + float GetDuration(); + void Info(); + + private: + typedef typename T::time_point TTP; + TTP m_StartTime; + +}; + +template void Timer::Start() +{ + m_StartTime = T::now(); +} + +template float Timer::GetDuration() +{ + std::chrono::duration duration = T::now() - m_StartTime; + return duration.count(); +} + + +template void Timer::Info() { + typedef typename T::period TPER; + typedef typename std::ratio_multiply MilliSec; + typedef typename std::ratio_multiply MicroSec; + + std::cout << std::boolalpha << std::endl; + std::cout << "clock info: " << std::endl; + std::cout << " is steady: " << T::is_steady << std::endl; + std::cout << " precision: " << TPER::num << "/" << TPER::den << " second " << std::endl; + std::cout << std::fixed; + std::cout << " " << static_cast(MilliSec::num)/MilliSec::den << " milliseconds " << std::endl; + std::cout << " " << static_cast(MicroSec::num)/MicroSec::den << " microseconds " << std::endl; + std::cout << std::endl; +} + diff --git a/epoch1/oneapi/gg_ttgg/known_issues.md b/epoch1/oneapi/gg_ttgg/known_issues.md new file mode 100644 index 0000000000..bfaaa932ed --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/known_issues.md @@ -0,0 +1,3 @@ +The function `ranmar` within rambo is not thread safe. Variables jranmr and iranmr are shared in a thread and the equal to zero check can fail due to a race condition where the value becomes less than zero. This causes negative values within square root functions when calculating the ME, leading to a segmentation fault. This is present in the original CUDA, so left as-is to be true to the original code. The ranmar function has been removed from epoch 2 code. + +This issue tracked in git issue [#79](https://github.com/madgraph5/madgraph4gpu/issues/79). \ No newline at end of file diff --git a/epoch1/oneapi/gg_ttgg/src/HelAmps_sm.cpp b/epoch1/oneapi/gg_ttgg/src/HelAmps_sm.cpp new file mode 100644 index 0000000000..68abec3fc8 --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/src/HelAmps_sm.cpp @@ -0,0 +1,806 @@ +//========================================================================== +// This file has been automatically generated for C++ Standalone by +// MadGraph5_aMC@NLO v. 2.7.3.py3, 2020-06-28 +// By the MadGraph5_aMC@NLO Development Team +// Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch +//========================================================================== + +#include +#include +#include "HelAmps_sm.h" +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +namespace MG5_sm +{ + +void ixxxxx(const double *pvec, double fmass, int nhel, int nsf, + std::complex fi[6]) +{ + std::complex chi[2]; + double sf[2], sfomega[2], omega[2], pp, pp3, sqp0p3, sqm[2]; + int ip, im, nh; + + double p[4] = {0, pvec[0], pvec[1], pvec[2]}; + p[0] = sycl::sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3] + fmass * fmass); + fi[0] = std::complex(-p[0] * nsf, -p[3] * nsf); + fi[1] = std::complex(-p[1] * nsf, -p[2] * nsf); + nh = nhel * nsf; + if (fmass != 0.0) + { + pp = sycl::min(p[0], sycl::sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3])); + if (pp == 0.0) + { + sqm[0] = sycl::sqrt(sycl::fabs(fmass)); + sqm[1] = (fmass < 0) ? -sycl::fabs(sqm[0]) : sycl::fabs(sqm[0]); + ip = (1 + nh)/2; + im = (1 - nh)/2; + fi[2] = ip * sqm[ip]; + fi[3] = im * nsf * sqm[ip]; + fi[4] = ip * nsf * sqm[im]; + fi[5] = im * sqm[im]; + } + else + { + sf[0] = (1 + nsf + (1 - nsf) * nh) * 0.5; + sf[1] = (1 + nsf - (1 - nsf) * nh) * 0.5; + omega[0] = sycl::sqrt(p[0] + pp); + omega[1] = fmass/omega[0]; + ip = (1 + nh)/2; + im = (1 - nh)/2; + sfomega[0] = sf[0] * omega[ip]; + sfomega[1] = sf[1] * omega[im]; + pp3 = sycl::max((double)(pp + p[3]), 0.0); + chi[0] = std::complex(sycl::sqrt(pp3 * 0.5 / pp), 0); + if (pp3 == 0.0) + { + chi[1] = std::complex(-nh, 0); + } + else + { + chi[1] = + std::complex(nh * p[1], p[2]) / sycl::sqrt(2.0 * pp * pp3); + } + fi[2] = sfomega[0] * chi[im]; + fi[3] = sfomega[0] * chi[ip]; + fi[4] = sfomega[1] * chi[im]; + fi[5] = sfomega[1] * chi[ip]; + } + } + else + { + if (p[1] == 0.0 and p[2] == 0.0 and p[3] < 0.0) + { + sqp0p3 = 0.0; + } + else + { + sqp0p3 = sycl::sqrt(sycl::max((double)(p[0] + p[3]), 0.0)) * nsf; + } + chi[0] = std::complex(sqp0p3, 0.0); + if (sqp0p3 == 0.0) + { + chi[1] = std::complex(-nhel * sycl::sqrt(2.0 * p[0]), 0.0); + } + else + { + chi[1] = std::complex(nh * p[1], p[2]) / sqp0p3; + } + if (nh == 1) + { + fi[2] = std::complex(0.0, 0.0); + fi[3] = std::complex(0.0, 0.0); + fi[4] = chi[0]; + fi[5] = chi[1]; + } + else + { + fi[2] = chi[1]; + fi[3] = chi[0]; + fi[4] = std::complex(0.0, 0.0); + fi[5] = std::complex(0.0, 0.0); + } + } + + return; +} + +void txxxxx(double pvec[3], double tmass, int nhel, int nst, + std::complex tc[18]) +{ + std::complex ft[6][4], ep[4], em[4], e0[4]; + double pt, pt2, pp, pzpt, emp, sqh, sqs; + int i, j; + + double p[4] = {0, pvec[0], pvec[1], pvec[2]}; + p[0] = sycl::sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3] + tmass * tmass); + sqh = sycl::sqrt(0.5); + sqs = sycl::sqrt(0.5 / 3); + + pt2 = p[1] * p[1] + p[2] * p[2]; + pp = sycl::min(p[0], sycl::sqrt(pt2 + p[3] * p[3])); + pt = sycl::min(pp, sycl::sqrt(pt2)); + + ft[4][0] = std::complex(p[0] * nst, p[3] * nst); + ft[5][0] = std::complex(p[1] * nst, p[2] * nst); + + // construct eps+ + if (nhel >= 0) + { + if (pp == 0) + { + ep[0] = std::complex(0, 0); + ep[1] = std::complex(-sqh, 0); + ep[2] = std::complex(0, nst * sqh); + ep[3] = std::complex(0, 0); + } + else + { + ep[0] = std::complex(0, 0); + ep[3] = std::complex(pt / pp * sqh, 0); + + if (pt != 0) + { + pzpt = p[3]/(pp * pt) * sqh; + ep[1] = std::complex(-p[1] * pzpt, -nst * p[2] / pt * sqh); + ep[2] = std::complex(-p[2] * pzpt, nst * p[1] / pt * sqh); + } + else + { + ep[1] = std::complex(-sqh, 0); + ep[2] = std::complex(0, nst * (p[3] < 0) ? -sycl::fabs(sqh) + : sycl::fabs(sqh)); + } + } + } + + // construct eps- + if (nhel <= 0) + { + if (pp == 0) + { + em[0] = std::complex(0, 0); + em[1] = std::complex(sqh, 0); + em[2] = std::complex(0, nst * sqh); + em[3] = std::complex(0, 0); + } + else + { + em[0] = std::complex(0, 0); + em[3] = std::complex(-pt / pp * sqh, 0); + + if (pt != 0) + { + pzpt = -p[3]/(pp * pt) * sqh; + em[1] = std::complex(-p[1] * pzpt, -nst * p[2] / pt * sqh); + em[2] = std::complex(-p[2] * pzpt, nst * p[1] / pt * sqh); + } + else + { + em[1] = std::complex(sqh, 0); + em[2] = std::complex(0, nst * (p[3] < 0) ? -sycl::fabs(sqh) + : sycl::fabs(sqh)); + } + } + } + + // construct eps0 + if (std::labs(nhel) <= 1) + { + if (pp == 0) + { + e0[0] = std::complex(0, 0); + e0[1] = std::complex(0, 0); + e0[2] = std::complex(0, 0); + e0[3] = std::complex(1, 0); + } + else + { + emp = p[0]/(tmass * pp); + e0[0] = std::complex(pp / tmass, 0); + e0[3] = std::complex(p[3] * emp, 0); + + if (pt != 0) + { + e0[1] = std::complex(p[1] * emp, 0); + e0[2] = std::complex(p[2] * emp, 0); + } + else + { + e0[1] = std::complex(0, 0); + e0[2] = std::complex(0, 0); + } + } + } + + if (nhel == 2) + { + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + ft[i][j] = ep[i] * ep[j]; + } + } + else if (nhel == -2) + { + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + ft[i][j] = em[i] * em[j]; + } + } + else if (tmass == 0) + { + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + ft[i][j] = 0; + } + } + else if (tmass != 0) + { + if (nhel == 1) + { + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + ft[i][j] = sqh * (ep[i] * e0[j] + e0[i] * ep[j]); + } + } + else if (nhel == 0) + { + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + ft[i][j] = + sqs * (ep[i] * em[j] + em[i] * ep[j] + 2.0 * e0[i] * e0[j]); + } + } + else if (nhel == -1) + { + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + ft[i][j] = sqh * (em[i] * e0[j] + e0[i] * em[j]); + } + } + else + { + // sr fixme // std::cerr << "Invalid helicity in txxxxx.\n"; + // sr fixme // std::exit(1); + } + } + + tc[0] = ft[4][0]; + tc[1] = ft[5][0]; + + for (j = 0; j < 4; j++ ) + { + for (i = 0; i < 4; i++ ) + tc[j * 4 + i + 2] = ft[j][i]; + } +} + +void vxxxxx(double pvec[3], double vmass, int nhel, int nsv, + std::complex vc[6]) +{ + double hel, hel0, pt, pt2, pp, pzpt, emp, sqh; + int nsvahl; + + double p[4] = {0, pvec[0], pvec[1], pvec[2]}; + p[0] = sycl::sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3] + vmass * vmass); + + sqh = sycl::sqrt(0.5); + hel = double(nhel); + nsvahl = nsv * sycl::fabs(hel); + pt2 = (p[1] * p[1]) + (p[2] * p[2]); + pp = sycl::min(p[0], sycl::sqrt(pt2 + (p[3] * p[3]))); + pt = sycl::min(pp, sycl::sqrt(pt2)); + vc[0] = std::complex(p[0] * nsv, p[3] * nsv); + vc[1] = std::complex(p[1] * nsv, p[2] * nsv); + if (vmass != 0.0) + { + hel0 = 1.0 - sycl::fabs(hel); + if (pp == 0.0) + { + vc[2] = std::complex(0.0, 0.0); + vc[3] = std::complex(-hel * sqh, 0.0); + vc[4] = std::complex(0.0, nsvahl * sqh); + vc[5] = std::complex(hel0, 0.0); + } + else + { + emp = p[0]/(vmass * pp); + vc[2] = std::complex(hel0 * pp / vmass, 0.0); + vc[5] = + std::complex(hel0 * p[3] * emp + hel * pt / pp * sqh, 0.0); + if (pt != 0.0) + { + pzpt = p[3]/(pp * pt) * sqh * hel; + vc[3] = std::complex(hel0 * p[1] * emp - p[1] * pzpt, + -nsvahl * p[2] / pt * sqh); + vc[4] = std::complex(hel0 * p[2] * emp - p[2] * pzpt, + nsvahl * p[1] / pt * sqh); + } + else + { + vc[3] = std::complex(-hel * sqh, 0.0); + vc[4] = std::complex( + 0.0, nsvahl * (p[3] < 0) ? -sycl::fabs(sqh) : sycl::fabs(sqh)); + } + } + } + else + { + pp = p[0]; + pt = sycl::sqrt((p[1] * p[1]) + (p[2] * p[2])); + vc[2] = std::complex(0.0, 0.0); + vc[5] = std::complex(hel * pt / pp * sqh, 0.0); + if (pt != 0.0) + { + pzpt = p[3]/(pp * pt) * sqh * hel; + vc[3] = std::complex(-p[1] * pzpt, -nsv * p[2] / pt * sqh); + vc[4] = std::complex(-p[2] * pzpt, nsv * p[1] / pt * sqh); + } + else + { + vc[3] = std::complex(-hel * sqh, 0.0); + vc[4] = std::complex(0.0, nsv * (p[3] < 0) ? -sycl::fabs(sqh) + : sycl::fabs(sqh)); + } + } + return; +} + +void sxxxxx(double pvec[3], int nss, std::complex sc[3], + sycl::stream stream_ct1) +{ + // double p[4] = {0, pvec[0], pvec[1], pvec[2]}; + // p[0] = sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3]+fmass*fmass); + double p[4] = {0, 0, 0, 0}; + stream_ct1 << "scalar not supported so far. to do: fix mass issue"; + sc[2] = std::complex(1.00, 0.00); + sc[0] = std::complex(p[0] * nss, p[3] * nss); + sc[1] = std::complex(p[1] * nss, p[2] * nss); + return; +} + +void oxxxxx(double pvec[3], double fmass, int nhel, int nsf, + std::complex fo[6]) +{ + std::complex chi[2]; + double sf[2], sfomeg[2], omega[2], pp, pp3, sqp0p3, sqm[2]; + int nh, ip, im; + + double p[4] = {0, pvec[0], pvec[1], pvec[2]}; + p[0] = sycl::sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3] + fmass * fmass); + + fo[0] = std::complex(p[0] * nsf, p[3] * nsf); + fo[1] = std::complex(p[1] * nsf, p[2] * nsf); + nh = nhel * nsf; + if (fmass != 0.000) + { + pp = sycl::min(p[0], + sycl::sqrt((p[1] * p[1]) + (p[2] * p[2]) + (p[3] * p[3]))); + if (pp == 0.000) + { + sqm[0] = sycl::sqrt(sycl::fabs(fmass)); + sqm[1] = (fmass < 0) ? -sycl::fabs(sqm[0]) : sycl::fabs(sqm[0]); + ip = -((1 - nh)/2) * nhel; + im = (1 + nh)/2 * nhel; + fo[2] = im * sqm[sycl::abs(ip)]; + fo[3] = ip * nsf * sqm[sycl::abs(ip)]; + fo[4] = im * nsf * sqm[sycl::abs(im)]; + fo[5] = ip * sqm[sycl::abs(im)]; + } + else + { + pp = sycl::min(p[0], + sycl::sqrt((p[1] * p[1]) + (p[2] * p[2]) + (p[3] * p[3]))); + sf[0] = double(1 + nsf + (1 - nsf) * nh) * 0.5; + sf[1] = double(1 + nsf - (1 - nsf) * nh) * 0.5; + omega[0] = sycl::sqrt(p[0] + pp); + omega[1] = fmass/omega[0]; + ip = (1 + nh)/2; + im = (1 - nh)/2; + sfomeg[0] = sf[0] * omega[ip]; + sfomeg[1] = sf[1] * omega[im]; + pp3 = sycl::max((double)(pp + p[3]), 0.00); + chi[0] = std::complex(sycl::sqrt(pp3 * 0.5 / pp), 0.00); + if (pp3 == 0.00) + { + chi[1] = std::complex(-nh, 0.00); + } + else + { + chi[1] = + std::complex(nh * p[1], -p[2]) / sycl::sqrt(2.0 * pp * pp3); + } + fo[2] = sfomeg[1] * chi[im]; + fo[3] = sfomeg[1] * chi[ip]; + fo[4] = sfomeg[0] * chi[im]; + fo[5] = sfomeg[0] * chi[ip]; + } + } + else + { + if ((p[1] == 0.00) and (p[2] == 0.00) and (p[3] < 0.00)) + { + sqp0p3 = 0.00; + } + else + { + sqp0p3 = sycl::sqrt(sycl::max((double)(p[0] + p[3]), 0.00)) * nsf; + } + chi[0] = std::complex(sqp0p3, 0.00); + if (sqp0p3 == 0.000) + { + chi[1] = std::complex(-nhel, 0.00) * sycl::sqrt(2.0 * p[0]); + } + else + { + chi[1] = std::complex(nh * p[1], -p[2]) / sqp0p3; + } + if (nh == 1) + { + fo[2] = chi[0]; + fo[3] = chi[1]; + fo[4] = std::complex(0.00, 0.00); + fo[5] = std::complex(0.00, 0.00); + } + else + { + fo[2] = std::complex(0.00, 0.00); + fo[3] = std::complex(0.00, 0.00); + fo[4] = chi[1]; + fo[5] = chi[0]; + } + } + return; +} +void FFV1_0(std::complex F1[], const thrust::complex F2[], + const thrust::complex V3[], + const thrust::complex COUP, std::complex *vertex) +{ + std::complex cI = std::complex(0., 1.); + std::complex TMP0; + TMP0 = (F1[2] * (F2[4] * (V3[2] + V3[5]) + F2[5] * (V3[3] + cI * (V3[4]))) + + (F1[3] * (F2[4] * (V3[3] - cI * (V3[4])) + F2[5] * (V3[2] - V3[5])) + + (F1[4] * (F2[2] * (V3[2] - V3[5]) - F2[3] * (V3[3] + cI * (V3[4]))) + + F1[5] * (F2[2] * (-V3[3] + cI * (V3[4])) + F2[3] * (V3[2] + V3[5]))))); + (*vertex) = COUP * - cI * TMP0; +} + +void FFV1_1(std::complex F2[], const thrust::complex V3[], + const thrust::complex COUP, const double M1, + const double W1, std::complex F1[]) +{ + std::complex cI = std::complex(0., 1.); + double P1[4]; + std::complex denom; + F1[0] = +F2[0] + V3[0]; + F1[1] = +F2[1] + V3[1]; + P1[0] = -F1[0].real(); + P1[1] = -F1[1].real(); + P1[2] = -F1[1].imag(); + P1[3] = -F1[0].imag(); + denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * + P1[3]) - M1 * (M1 - cI * W1)); + F1[2] = denom * cI * (F2[2] * (P1[0] * (-V3[2] + V3[5]) + (P1[1] * (V3[3] - + cI * (V3[4])) + (P1[2] * (+cI * (V3[3]) + V3[4]) + P1[3] * (-V3[2] + + V3[5])))) + (F2[3] * (P1[0] * (V3[3] + cI * (V3[4])) + (P1[1] * (-1.) * + (V3[2] + V3[5]) + (P1[2] * (-1.) * (+cI * (V3[2] + V3[5])) + P1[3] * + (V3[3] + cI * (V3[4]))))) + M1 * (F2[4] * (V3[2] + V3[5]) + F2[5] * + (V3[3] + cI * (V3[4]))))); + F1[3] = denom * (-cI) * (F2[2] * (P1[0] * (-V3[3] + cI * (V3[4])) + (P1[1] * + (V3[2] - V3[5]) + (P1[2] * (-cI * (V3[2]) + cI * (V3[5])) + P1[3] * + (V3[3] - cI * (V3[4]))))) + (F2[3] * (P1[0] * (V3[2] + V3[5]) + (P1[1] * + (-1.) * (V3[3] + cI * (V3[4])) + (P1[2] * (+cI * (V3[3]) - V3[4]) - P1[3] + * (V3[2] + V3[5])))) + M1 * (F2[4] * (-V3[3] + cI * (V3[4])) + F2[5] * + (-V3[2] + V3[5])))); + F1[4] = denom * (-cI) * (F2[4] * (P1[0] * (V3[2] + V3[5]) + (P1[1] * (-V3[3] + + cI * (V3[4])) + (P1[2] * (-1.) * (+cI * (V3[3]) + V3[4]) - P1[3] * + (V3[2] + V3[5])))) + (F2[5] * (P1[0] * (V3[3] + cI * (V3[4])) + (P1[1] * + (-V3[2] + V3[5]) + (P1[2] * (-cI * (V3[2]) + cI * (V3[5])) - P1[3] * + (V3[3] + cI * (V3[4]))))) + M1 * (F2[2] * (-V3[2] + V3[5]) + F2[3] * + (V3[3] + cI * (V3[4]))))); + F1[5] = denom * cI * (F2[4] * (P1[0] * (-V3[3] + cI * (V3[4])) + (P1[1] * + (V3[2] + V3[5]) + (P1[2] * (-1.) * (+cI * (V3[2] + V3[5])) + P1[3] * + (-V3[3] + cI * (V3[4]))))) + (F2[5] * (P1[0] * (-V3[2] + V3[5]) + (P1[1] + * (V3[3] + cI * (V3[4])) + (P1[2] * (-cI * (V3[3]) + V3[4]) + P1[3] * + (-V3[2] + V3[5])))) + M1 * (F2[2] * (-V3[3] + cI * (V3[4])) + F2[3] * + (V3[2] + V3[5])))); +} + +void FFV1_2(std::complex F1[], const thrust::complex V3[], + const thrust::complex COUP, const double M2, + const double W2, std::complex F2[]) +{ + std::complex cI = std::complex(0., 1.); + double P2[4]; + std::complex denom; + F2[0] = +F1[0] + V3[0]; + F2[1] = +F1[1] + V3[1]; + P2[0] = -F2[0].real(); + P2[1] = -F2[1].real(); + P2[2] = -F2[1].imag(); + P2[3] = -F2[0].imag(); + denom = COUP/((P2[0] * P2[0]) - (P2[1] * P2[1]) - (P2[2] * P2[2]) - (P2[3] * + P2[3]) - M2 * (M2 - cI * W2)); + F2[2] = denom * cI * (F1[2] * (P2[0] * (V3[2] + V3[5]) + (P2[1] * (-1.) * + (V3[3] + cI * (V3[4])) + (P2[2] * (+cI * (V3[3]) - V3[4]) - P2[3] * + (V3[2] + V3[5])))) + (F1[3] * (P2[0] * (V3[3] - cI * (V3[4])) + (P2[1] * + (-V3[2] + V3[5]) + (P2[2] * (+cI * (V3[2]) - cI * (V3[5])) + P2[3] * + (-V3[3] + cI * (V3[4]))))) + M2 * (F1[4] * (V3[2] - V3[5]) + F1[5] * + (-V3[3] + cI * (V3[4]))))); + F2[3] = denom * (-cI) * (F1[2] * (P2[0] * (-1.) * (V3[3] + cI * (V3[4])) + + (P2[1] * (V3[2] + V3[5]) + (P2[2] * (+cI * (V3[2] + V3[5])) - P2[3] * + (V3[3] + cI * (V3[4]))))) + (F1[3] * (P2[0] * (-V3[2] + V3[5]) + (P2[1] * + (V3[3] - cI * (V3[4])) + (P2[2] * (+cI * (V3[3]) + V3[4]) + P2[3] * + (-V3[2] + V3[5])))) + M2 * (F1[4] * (V3[3] + cI * (V3[4])) - F1[5] * + (V3[2] + V3[5])))); + F2[4] = denom * (-cI) * (F1[4] * (P2[0] * (-V3[2] + V3[5]) + (P2[1] * (V3[3] + + cI * (V3[4])) + (P2[2] * (-cI * (V3[3]) + V3[4]) + P2[3] * (-V3[2] + + V3[5])))) + (F1[5] * (P2[0] * (V3[3] - cI * (V3[4])) + (P2[1] * (-1.) * + (V3[2] + V3[5]) + (P2[2] * (+cI * (V3[2] + V3[5])) + P2[3] * (V3[3] - cI + * (V3[4]))))) + M2 * (F1[2] * (-1.) * (V3[2] + V3[5]) + F1[3] * (-V3[3] + + cI * (V3[4]))))); + F2[5] = denom * cI * (F1[4] * (P2[0] * (-1.) * (V3[3] + cI * (V3[4])) + + (P2[1] * (V3[2] - V3[5]) + (P2[2] * (+cI * (V3[2]) - cI * (V3[5])) + + P2[3] * (V3[3] + cI * (V3[4]))))) + (F1[5] * (P2[0] * (V3[2] + V3[5]) + + (P2[1] * (-V3[3] + cI * (V3[4])) + (P2[2] * (-1.) * (+cI * (V3[3]) + + V3[4]) - P2[3] * (V3[2] + V3[5])))) + M2 * (F1[2] * (V3[3] + cI * + (V3[4])) + F1[3] * (V3[2] - V3[5])))); +} + +void FFV1P0_3(std::complex F1[], const thrust::complex F2[], + const thrust::complex COUP, const double M3, + const double W3, std::complex V3[]) +{ + std::complex cI = std::complex(0., 1.); + double P3[4]; + std::complex denom; + V3[0] = +F1[0] + F2[0]; + V3[1] = +F1[1] + F2[1]; + P3[0] = -V3[0].real(); + P3[1] = -V3[1].real(); + P3[2] = -V3[1].imag(); + P3[3] = -V3[0].imag(); + denom = COUP/((P3[0] * P3[0]) - (P3[1] * P3[1]) - (P3[2] * P3[2]) - (P3[3] * + P3[3]) - M3 * (M3 - cI * W3)); + V3[2] = denom * (-cI) * (F1[2] * F2[4] + F1[3] * F2[5] + F1[4] * F2[2] + + F1[5] * F2[3]); + V3[3] = denom * (-cI) * (-F1[2] * F2[5] - F1[3] * F2[4] + F1[4] * F2[3] + + F1[5] * F2[2]); + V3[4] = denom * (-cI) * (-cI * (F1[2] * F2[5] + F1[5] * F2[2]) + cI * (F1[3] + * F2[4] + F1[4] * F2[3])); + V3[5] = denom * (-cI) * (-F1[2] * F2[4] - F1[5] * F2[3] + F1[3] * F2[5] + + F1[4] * F2[2]); +} + +void VVVV3_0(std::complex V1[], const thrust::complex V2[], + const thrust::complex V3[], + const thrust::complex V4[], + const thrust::complex COUP, std::complex *vertex) +{ + std::complex cI = std::complex(0., 1.); + std::complex TMP1; + std::complex TMP2; + std::complex TMP3; + std::complex TMP4; + TMP4 = (V1[2] * V2[2] - V1[3] * V2[3] - V1[4] * V2[4] - V1[5] * V2[5]); + TMP1 = (V4[2] * V1[2] - V4[3] * V1[3] - V4[4] * V1[4] - V4[5] * V1[5]); + TMP2 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); + TMP3 = (V3[2] * V4[2] - V3[3] * V4[3] - V3[4] * V4[4] - V3[5] * V4[5]); + (*vertex) = COUP * (-cI * (TMP1 * TMP2) + cI * (TMP3 * TMP4)); +} + +void VVVV3P0_1(std::complex V2[], const thrust::complex V3[], + const thrust::complex V4[], + const thrust::complex COUP, const double M1, + const double W1, std::complex V1[]) +{ + std::complex cI = std::complex(0., 1.); + double P1[4]; + std::complex TMP2; + std::complex TMP3; + std::complex denom; + V1[0] = +V2[0] + V3[0] + V4[0]; + V1[1] = +V2[1] + V3[1] + V4[1]; + P1[0] = -V1[0].real(); + P1[1] = -V1[1].real(); + P1[2] = -V1[1].imag(); + P1[3] = -V1[0].imag(); + TMP2 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); + TMP3 = (V3[2] * V4[2] - V3[3] * V4[3] - V3[4] * V4[4] - V3[5] * V4[5]); + denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * + P1[3]) - M1 * (M1 - cI * W1)); + V1[2] = denom * (-cI * (V4[2] * TMP2) + cI * (V2[2] * TMP3)); + V1[3] = denom * (-cI * (V4[3] * TMP2) + cI * (V2[3] * TMP3)); + V1[4] = denom * (-cI * (V4[4] * TMP2) + cI * (V2[4] * TMP3)); + V1[5] = denom * (-cI * (V4[5] * TMP2) + cI * (V2[5] * TMP3)); +} + +void VVVV1_0(std::complex V1[], const thrust::complex V2[], + const thrust::complex V3[], + const thrust::complex V4[], + const thrust::complex COUP, std::complex *vertex) +{ + std::complex cI = std::complex(0., 1.); + std::complex TMP1; + std::complex TMP2; + std::complex TMP5; + std::complex TMP6; + TMP6 = (V3[2] * V1[2] - V3[3] * V1[3] - V3[4] * V1[4] - V3[5] * V1[5]); + TMP5 = (V4[2] * V2[2] - V4[3] * V2[3] - V4[4] * V2[4] - V4[5] * V2[5]); + TMP1 = (V4[2] * V1[2] - V4[3] * V1[3] - V4[4] * V1[4] - V4[5] * V1[5]); + TMP2 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); + (*vertex) = COUP * (-cI * (TMP1 * TMP2) + cI * (TMP5 * TMP6)); +} + +void VVVV1P0_1(std::complex V2[], const thrust::complex V3[], + const thrust::complex V4[], + const thrust::complex COUP, const double M1, + const double W1, std::complex V1[]) +{ + std::complex cI = std::complex(0., 1.); + double P1[4]; + std::complex TMP2; + std::complex TMP5; + std::complex denom; + V1[0] = +V2[0] + V3[0] + V4[0]; + V1[1] = +V2[1] + V3[1] + V4[1]; + P1[0] = -V1[0].real(); + P1[1] = -V1[1].real(); + P1[2] = -V1[1].imag(); + P1[3] = -V1[0].imag(); + TMP5 = (V4[2] * V2[2] - V4[3] * V2[3] - V4[4] * V2[4] - V4[5] * V2[5]); + TMP2 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); + denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * + P1[3]) - M1 * (M1 - cI * W1)); + V1[2] = denom * (-cI * (V4[2] * TMP2) + cI * (V3[2] * TMP5)); + V1[3] = denom * (-cI * (V4[3] * TMP2) + cI * (V3[3] * TMP5)); + V1[4] = denom * (-cI * (V4[4] * TMP2) + cI * (V3[4] * TMP5)); + V1[5] = denom * (-cI * (V4[5] * TMP2) + cI * (V3[5] * TMP5)); +} + +void VVVV4_0(std::complex V1[], const thrust::complex V2[], + const thrust::complex V3[], + const thrust::complex V4[], + const thrust::complex COUP, std::complex *vertex) +{ + std::complex cI = std::complex(0., 1.); + std::complex TMP3; + std::complex TMP4; + std::complex TMP5; + std::complex TMP6; + TMP4 = (V1[2] * V2[2] - V1[3] * V2[3] - V1[4] * V2[4] - V1[5] * V2[5]); + TMP5 = (V4[2] * V2[2] - V4[3] * V2[3] - V4[4] * V2[4] - V4[5] * V2[5]); + TMP6 = (V3[2] * V1[2] - V3[3] * V1[3] - V3[4] * V1[4] - V3[5] * V1[5]); + TMP3 = (V3[2] * V4[2] - V3[3] * V4[3] - V3[4] * V4[4] - V3[5] * V4[5]); + (*vertex) = COUP * (-cI * (TMP5 * TMP6) + cI * (TMP3 * TMP4)); +} + +void VVVV4P0_1(std::complex V2[], const thrust::complex V3[], + const thrust::complex V4[], + const thrust::complex COUP, const double M1, + const double W1, std::complex V1[]) +{ + std::complex cI = std::complex(0., 1.); + double P1[4]; + std::complex TMP3; + std::complex TMP5; + std::complex denom; + V1[0] = +V2[0] + V3[0] + V4[0]; + V1[1] = +V2[1] + V3[1] + V4[1]; + P1[0] = -V1[0].real(); + P1[1] = -V1[1].real(); + P1[2] = -V1[1].imag(); + P1[3] = -V1[0].imag(); + TMP5 = (V4[2] * V2[2] - V4[3] * V2[3] - V4[4] * V2[4] - V4[5] * V2[5]); + TMP3 = (V3[2] * V4[2] - V3[3] * V4[3] - V3[4] * V4[4] - V3[5] * V4[5]); + denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * + P1[3]) - M1 * (M1 - cI * W1)); + V1[2] = denom * (-cI * (V3[2] * TMP5) + cI * (V2[2] * TMP3)); + V1[3] = denom * (-cI * (V3[3] * TMP5) + cI * (V2[3] * TMP3)); + V1[4] = denom * (-cI * (V3[4] * TMP5) + cI * (V2[4] * TMP3)); + V1[5] = denom * (-cI * (V3[5] * TMP5) + cI * (V2[5] * TMP3)); +} + +void VVV1_0(std::complex V1[], const thrust::complex V2[], + const thrust::complex V3[], + const thrust::complex COUP, std::complex *vertex) +{ + std::complex cI = std::complex(0., 1.); + double P1[4]; + double P2[4]; + double P3[4]; + std::complex TMP10; + std::complex TMP11; + std::complex TMP12; + std::complex TMP2; + std::complex TMP4; + std::complex TMP6; + std::complex TMP7; + std::complex TMP8; + std::complex TMP9; + P1[0] = V1[0].real(); + P1[1] = V1[1].real(); + P1[2] = V1[1].imag(); + P1[3] = V1[0].imag(); + P2[0] = V2[0].real(); + P2[1] = V2[1].real(); + P2[2] = V2[1].imag(); + P2[3] = V2[0].imag(); + P3[0] = V3[0].real(); + P3[1] = V3[1].real(); + P3[2] = V3[1].imag(); + P3[3] = V3[0].imag(); + TMP8 = (V3[2] * P2[0] - V3[3] * P2[1] - V3[4] * P2[2] - V3[5] * P2[3]); + TMP10 = (V2[2] * P3[0] - V2[3] * P3[1] - V2[4] * P3[2] - V2[5] * P3[3]); + TMP9 = (P1[0] * V2[2] - P1[1] * V2[3] - P1[2] * V2[4] - P1[3] * V2[5]); + TMP2 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); + TMP4 = (V1[2] * V2[2] - V1[3] * V2[3] - V1[4] * V2[4] - V1[5] * V2[5]); + TMP6 = (V3[2] * V1[2] - V3[3] * V1[3] - V3[4] * V1[4] - V3[5] * V1[5]); + TMP12 = (V1[2] * P3[0] - V1[3] * P3[1] - V1[4] * P3[2] - V1[5] * P3[3]); + TMP7 = (V3[2] * P1[0] - V3[3] * P1[1] - V3[4] * P1[2] - V3[5] * P1[3]); + TMP11 = (P2[0] * V1[2] - P2[1] * V1[3] - P2[2] * V1[4] - P2[3] * V1[5]); + (*vertex) = COUP * (TMP2 * (-cI * (TMP11) + cI * (TMP12)) + (TMP4 * (-cI * + (TMP7) + cI * (TMP8)) + TMP6 * (+cI * (TMP9) - cI * (TMP10)))); +} + +void VVV1P0_1(std::complex V2[], const thrust::complex V3[], + const thrust::complex COUP, const double M1, + const double W1, std::complex V1[]) +{ + std::complex cI = std::complex(0., 1.); + double P1[4]; + double P2[4]; + double P3[4]; + std::complex TMP10; + std::complex TMP2; + std::complex TMP7; + std::complex TMP8; + std::complex TMP9; + std::complex denom; + P2[0] = V2[0].real(); + P2[1] = V2[1].real(); + P2[2] = V2[1].imag(); + P2[3] = V2[0].imag(); + P3[0] = V3[0].real(); + P3[1] = V3[1].real(); + P3[2] = V3[1].imag(); + P3[3] = V3[0].imag(); + V1[0] = +V2[0] + V3[0]; + V1[1] = +V2[1] + V3[1]; + P1[0] = -V1[0].real(); + P1[1] = -V1[1].real(); + P1[2] = -V1[1].imag(); + P1[3] = -V1[0].imag(); + TMP8 = (V3[2] * P2[0] - V3[3] * P2[1] - V3[4] * P2[2] - V3[5] * P2[3]); + TMP10 = (V2[2] * P3[0] - V2[3] * P3[1] - V2[4] * P3[2] - V2[5] * P3[3]); + TMP9 = (P1[0] * V2[2] - P1[1] * V2[3] - P1[2] * V2[4] - P1[3] * V2[5]); + TMP2 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); + TMP7 = (V3[2] * P1[0] - V3[3] * P1[1] - V3[4] * P1[2] - V3[5] * P1[3]); + denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * + P1[3]) - M1 * (M1 - cI * W1)); + V1[2] = denom * (TMP2 * (-cI * (P2[0]) + cI * (P3[0])) + (V2[2] * (-cI * + (TMP7) + cI * (TMP8)) + V3[2] * (+cI * (TMP9) - cI * (TMP10)))); + V1[3] = denom * (TMP2 * (-cI * (P2[1]) + cI * (P3[1])) + (V2[3] * (-cI * + (TMP7) + cI * (TMP8)) + V3[3] * (+cI * (TMP9) - cI * (TMP10)))); + V1[4] = denom * (TMP2 * (-cI * (P2[2]) + cI * (P3[2])) + (V2[4] * (-cI * + (TMP7) + cI * (TMP8)) + V3[4] * (+cI * (TMP9) - cI * (TMP10)))); + V1[5] = denom * (TMP2 * (-cI * (P2[3]) + cI * (P3[3])) + (V2[5] * (-cI * + (TMP7) + cI * (TMP8)) + V3[5] * (+cI * (TMP9) - cI * (TMP10)))); +} + + +} // end namespace $(namespace)s_sm + diff --git a/epoch1/oneapi/gg_ttgg/src/HelAmps_sm.h b/epoch1/oneapi/gg_ttgg/src/HelAmps_sm.h new file mode 100644 index 0000000000..c92e8deacd --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/src/HelAmps_sm.h @@ -0,0 +1,88 @@ +//========================================================================== +// This file has been automatically generated for C++ Standalone +// MadGraph5_aMC@NLO v. 2.7.3.py3, 2020-06-28 +// By the MadGraph5_aMC@NLO Development Team +// Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch +// GPU Template +//========================================================================== + +#ifndef HelAmps_sm_H +#define HelAmps_sm_H + +#include +#include +#include + +namespace MG5_sm +{ +void oxxxxx(double p[4], double fmass, int nhel, int nsf, +std::complex fo[6]); + +void sxxxxx(double p[4], int nss, std::complex sc[3]); + +void ixxxxx(double p[4], double fmass, int nhel, int nsf, +std::complex fi[6]); + +void txxxxx(double p[4], double tmass, int nhel, int nst, +std::complex fi[18]); + +void vxxxxx(double p[4], double vmass, int nhel, int nsv, +std::complex v[6]); + +void FFV1_0(std::complex F1[], const + std::complex F2[], const std::complex V3[], const + std::complex COUP, std::complex * vertex); + +void FFV1_1(std::complex F2[], const + std::complex V3[], const std::complex COUP, const + double M1, const double W1, std::complex F1[]); + +void FFV1_2(std::complex F1[], const + std::complex V3[], const std::complex COUP, const + double M2, const double W2, std::complex F2[]); + +void FFV1P0_3(std::complex F1[], const + std::complex F2[], const std::complex COUP, const + double M3, const double W3, std::complex V3[]); + +void VVVV3_0(std::complex V1[], const + std::complex V2[], const std::complex V3[], const + std::complex V4[], const std::complex COUP, + std::complex * vertex); + +void VVVV3P0_1(std::complex V2[], const + std::complex V3[], const std::complex V4[], const + std::complex COUP, const double M1, const double W1, + std::complex V1[]); + +void VVVV1_0(std::complex V1[], const + std::complex V2[], const std::complex V3[], const + std::complex V4[], const std::complex COUP, + std::complex * vertex); + +void VVVV1P0_1(std::complex V2[], const + std::complex V3[], const std::complex V4[], const + std::complex COUP, const double M1, const double W1, + std::complex V1[]); + +void VVVV4_0(std::complex V1[], const + std::complex V2[], const std::complex V3[], const + std::complex V4[], const std::complex COUP, + std::complex * vertex); + +void VVVV4P0_1(std::complex V2[], const + std::complex V3[], const std::complex V4[], const + std::complex COUP, const double M1, const double W1, + std::complex V1[]); + +void VVV1_0(std::complex V1[], const + std::complex V2[], const std::complex V3[], const + std::complex COUP, std::complex * vertex); + +void VVV1P0_1(std::complex V2[], const + std::complex V3[], const std::complex COUP, const + double M1, const double W1, std::complex V1[]); + +} // end namespace MG5_sm + +#endif // HelAmps_sm_H diff --git a/epoch1/oneapi/gg_ttgg/src/Makefile b/epoch1/oneapi/gg_ttgg/src/Makefile new file mode 100644 index 0000000000..d87385f442 --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/src/Makefile @@ -0,0 +1,38 @@ +#CUARCHNUM=70 +LIBDIR=../lib +#CXXFLAGS= -O3 -I. +CXXFLAGS=-Wno-sycl-strict -I. +#CUARCHFLAGS= -arch=compute_$(CUARCHNUM) +#CUFLAGS= $(CUARCHFLAGS) -use_fast_math -lineinfo +#CUINC=/usr/local/cuda/targets/x86_64-linux/include +CXX=dpcpp + + +target=$(LIBDIR)/libmodel_sm.a +cxx_objects=Parameters_sm.o read_slha.o rambo.o + + +all: $(target) + +debug: CXXFLAGS:=$(filter-out -O3,$(CXXFLAGS)) +debug: CXXFLAGS += -g -O0 -DDEBUG2 +debug: $(target) + +# sr fixme # +# not sure including the cuda includes here is a good idea... +# ... needed for thrust/complex which in principal is a C++ class +%.o : %.cc %.h + $(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $< -o $@ + + + +$(target): $(cxx_objects) + if [ ! -d $(LIBDIR) ]; then mkdir $(LIBDIR); fi + $(AR) cru $@ $(cxx_objects) + ranlib $(target) + +.PHONY: clean + +clean: + rm -f $(target) + rm -f $(cxx_objects) diff --git a/epoch1/oneapi/gg_ttgg/src/Parameters_sm.cc b/epoch1/oneapi/gg_ttgg/src/Parameters_sm.cc new file mode 100644 index 0000000000..4567def136 --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/src/Parameters_sm.cc @@ -0,0 +1,225 @@ +//========================================================================== +// This file has been automatically generated for C++ by +// MadGraph5_aMC@NLO v. 2.7.3.py3, 2020-06-28 +// By the MadGraph5_aMC@NLO Development Team +// Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch +//========================================================================== + +#include +#include +#include "Parameters_sm.h" + +// Initialize static instance +Parameters_sm * Parameters_sm::instance = 0; + +// Function to get static instance - only one instance per program +Parameters_sm * Parameters_sm::getInstance() +{ + if (instance == 0) + instance = new Parameters_sm(); + + return instance; +} + +void Parameters_sm::setIndependentParameters(SLHAReader& slha) +{ + // Define "zero" + zero = 0; + ZERO = 0; + // Prepare a vector for indices + vector indices(2, 0); + mdl_WH = slha.get_block_entry("decay", 25, 6.382339e-03); + mdl_WW = slha.get_block_entry("decay", 24, 2.047600e+00); + mdl_WZ = slha.get_block_entry("decay", 23, 2.441404e+00); + mdl_WT = slha.get_block_entry("decay", 6, 1.491500e+00); + mdl_ymtau = slha.get_block_entry("yukawa", 15, 1.777000e+00); + mdl_ymt = slha.get_block_entry("yukawa", 6, 1.730000e+02); + mdl_ymb = slha.get_block_entry("yukawa", 5, 4.700000e+00); + aS = slha.get_block_entry("sminputs", 3, 1.180000e-01); + mdl_Gf = slha.get_block_entry("sminputs", 2, 1.166390e-05); + aEWM1 = slha.get_block_entry("sminputs", 1, 1.325070e+02); + mdl_MH = slha.get_block_entry("mass", 25, 1.250000e+02); + mdl_MZ = slha.get_block_entry("mass", 23, 9.118800e+01); + mdl_MTA = slha.get_block_entry("mass", 15, 1.777000e+00); + mdl_MT = slha.get_block_entry("mass", 6, 1.730000e+02); + mdl_MB = slha.get_block_entry("mass", 5, 4.700000e+00); + mdl_conjg__CKM1x1 = 1.; + mdl_conjg__CKM3x3 = 1.; + mdl_CKM3x3 = 1.; + mdl_complexi = std::complex (0., 1.); + mdl_MZ__exp__2 = ((mdl_MZ) * (mdl_MZ)); + mdl_MZ__exp__4 = ((mdl_MZ) * (mdl_MZ) * (mdl_MZ) * (mdl_MZ)); + mdl_sqrt__2 = sqrt(2.); + mdl_MH__exp__2 = ((mdl_MH) * (mdl_MH)); + mdl_aEW = 1./aEWM1; + mdl_MW = sqrt(mdl_MZ__exp__2/2. + sqrt(mdl_MZ__exp__4/4. - (mdl_aEW * M_PI * + mdl_MZ__exp__2)/(mdl_Gf * mdl_sqrt__2))); + mdl_sqrt__aEW = sqrt(mdl_aEW); + mdl_ee = 2. * mdl_sqrt__aEW * sqrt(M_PI); + mdl_MW__exp__2 = ((mdl_MW) * (mdl_MW)); + mdl_sw2 = 1. - mdl_MW__exp__2/mdl_MZ__exp__2; + mdl_cw = sqrt(1. - mdl_sw2); + mdl_sqrt__sw2 = sqrt(mdl_sw2); + mdl_sw = mdl_sqrt__sw2; + mdl_g1 = mdl_ee/mdl_cw; + mdl_gw = mdl_ee/mdl_sw; + mdl_vev = (2. * mdl_MW * mdl_sw)/mdl_ee; + mdl_vev__exp__2 = ((mdl_vev) * (mdl_vev)); + mdl_lam = mdl_MH__exp__2/(2. * mdl_vev__exp__2); + mdl_yb = (mdl_ymb * mdl_sqrt__2)/mdl_vev; + mdl_yt = (mdl_ymt * mdl_sqrt__2)/mdl_vev; + mdl_ytau = (mdl_ymtau * mdl_sqrt__2)/mdl_vev; + mdl_muH = sqrt(mdl_lam * mdl_vev__exp__2); + mdl_I1x33 = mdl_yb * mdl_conjg__CKM3x3; + mdl_I2x33 = mdl_yt * mdl_conjg__CKM3x3; + mdl_I3x33 = mdl_CKM3x3 * mdl_yt; + mdl_I4x33 = mdl_CKM3x3 * mdl_yb; + mdl_ee__exp__2 = ((mdl_ee) * (mdl_ee)); + mdl_sw__exp__2 = ((mdl_sw) * (mdl_sw)); + mdl_cw__exp__2 = ((mdl_cw) * (mdl_cw)); +} +void Parameters_sm::setIndependentCouplings() +{ + +} +void Parameters_sm::setDependentParameters() +{ + mdl_sqrt__aS = sqrt(aS); + G = 2. * mdl_sqrt__aS * sqrt(M_PI); + mdl_G__exp__2 = ((G) * (G)); +} +void Parameters_sm::setDependentCouplings() +{ + GC_10 = -G; + GC_11 = mdl_complexi * G; + GC_12 = mdl_complexi * mdl_G__exp__2; +} + +// Routines for printing out parameters +void Parameters_sm::printIndependentParameters() +{ + cout << "sm model parameters independent of event kinematics:" << endl; + cout << setw(20) << "mdl_WH " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_WH << endl; + cout << setw(20) << "mdl_WW " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_WW << endl; + cout << setw(20) << "mdl_WZ " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_WZ << endl; + cout << setw(20) << "mdl_WT " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_WT << endl; + cout << setw(20) << "mdl_ymtau " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_ymtau << endl; + cout << setw(20) << "mdl_ymt " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_ymt << endl; + cout << setw(20) << "mdl_ymb " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_ymb << endl; + cout << setw(20) << "aS " << "= " << setiosflags(ios::scientific) << + setw(10) << aS << endl; + cout << setw(20) << "mdl_Gf " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_Gf << endl; + cout << setw(20) << "aEWM1 " << "= " << setiosflags(ios::scientific) << + setw(10) << aEWM1 << endl; + cout << setw(20) << "mdl_MH " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_MH << endl; + cout << setw(20) << "mdl_MZ " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_MZ << endl; + cout << setw(20) << "mdl_MTA " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_MTA << endl; + cout << setw(20) << "mdl_MT " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_MT << endl; + cout << setw(20) << "mdl_MB " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_MB << endl; + cout << setw(20) << "mdl_conjg__CKM1x1 " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_conjg__CKM1x1 << endl; + cout << setw(20) << "mdl_conjg__CKM3x3 " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_conjg__CKM3x3 << endl; + cout << setw(20) << "mdl_CKM3x3 " << "= " << setiosflags(ios::scientific) + << setw(10) << mdl_CKM3x3 << endl; + cout << setw(20) << "mdl_complexi " << "= " << setiosflags(ios::scientific) + << setw(10) << mdl_complexi << endl; + cout << setw(20) << "mdl_MZ__exp__2 " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_MZ__exp__2 << endl; + cout << setw(20) << "mdl_MZ__exp__4 " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_MZ__exp__4 << endl; + cout << setw(20) << "mdl_sqrt__2 " << "= " << setiosflags(ios::scientific) + << setw(10) << mdl_sqrt__2 << endl; + cout << setw(20) << "mdl_MH__exp__2 " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_MH__exp__2 << endl; + cout << setw(20) << "mdl_aEW " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_aEW << endl; + cout << setw(20) << "mdl_MW " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_MW << endl; + cout << setw(20) << "mdl_sqrt__aEW " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_sqrt__aEW << endl; + cout << setw(20) << "mdl_ee " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_ee << endl; + cout << setw(20) << "mdl_MW__exp__2 " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_MW__exp__2 << endl; + cout << setw(20) << "mdl_sw2 " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_sw2 << endl; + cout << setw(20) << "mdl_cw " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_cw << endl; + cout << setw(20) << "mdl_sqrt__sw2 " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_sqrt__sw2 << endl; + cout << setw(20) << "mdl_sw " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_sw << endl; + cout << setw(20) << "mdl_g1 " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_g1 << endl; + cout << setw(20) << "mdl_gw " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_gw << endl; + cout << setw(20) << "mdl_vev " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_vev << endl; + cout << setw(20) << "mdl_vev__exp__2 " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_vev__exp__2 << endl; + cout << setw(20) << "mdl_lam " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_lam << endl; + cout << setw(20) << "mdl_yb " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_yb << endl; + cout << setw(20) << "mdl_yt " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_yt << endl; + cout << setw(20) << "mdl_ytau " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_ytau << endl; + cout << setw(20) << "mdl_muH " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_muH << endl; + cout << setw(20) << "mdl_I1x33 " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_I1x33 << endl; + cout << setw(20) << "mdl_I2x33 " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_I2x33 << endl; + cout << setw(20) << "mdl_I3x33 " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_I3x33 << endl; + cout << setw(20) << "mdl_I4x33 " << "= " << setiosflags(ios::scientific) << + setw(10) << mdl_I4x33 << endl; + cout << setw(20) << "mdl_ee__exp__2 " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_ee__exp__2 << endl; + cout << setw(20) << "mdl_sw__exp__2 " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_sw__exp__2 << endl; + cout << setw(20) << "mdl_cw__exp__2 " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_cw__exp__2 << endl; +} +void Parameters_sm::printIndependentCouplings() +{ + cout << "sm model couplings independent of event kinematics:" << endl; + +} +void Parameters_sm::printDependentParameters() +{ + cout << "sm model parameters dependent on event kinematics:" << endl; + cout << setw(20) << "mdl_sqrt__aS " << "= " << setiosflags(ios::scientific) + << setw(10) << mdl_sqrt__aS << endl; + cout << setw(20) << "G " << "= " << setiosflags(ios::scientific) << + setw(10) << G << endl; + cout << setw(20) << "mdl_G__exp__2 " << "= " << + setiosflags(ios::scientific) << setw(10) << mdl_G__exp__2 << endl; +} +void Parameters_sm::printDependentCouplings() +{ + cout << "sm model couplings dependent on event kinematics:" << endl; + cout << setw(20) << "GC_10 " << "= " << setiosflags(ios::scientific) << + setw(10) << GC_10 << endl; + cout << setw(20) << "GC_11 " << "= " << setiosflags(ios::scientific) << + setw(10) << GC_11 << endl; + cout << setw(20) << "GC_12 " << "= " << setiosflags(ios::scientific) << + setw(10) << GC_12 << endl; +} + + diff --git a/epoch1/oneapi/gg_ttgg/src/Parameters_sm.h b/epoch1/oneapi/gg_ttgg/src/Parameters_sm.h new file mode 100644 index 0000000000..607595677c --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/src/Parameters_sm.h @@ -0,0 +1,66 @@ +//========================================================================== +// This file has been automatically generated for C++ +// MadGraph5_aMC@NLO v. 2.7.3.py3, 2020-06-28 +// By the MadGraph5_aMC@NLO Development Team +// Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch +//========================================================================== + +#ifndef Parameters_sm_H +#define Parameters_sm_H + +#include + +#include "read_slha.h" +using namespace std; + +class Parameters_sm +{ + public: + + static Parameters_sm * getInstance(); + + // Define "zero" + double zero, ZERO; + // Model parameters independent of aS + double mdl_WH, mdl_WW, mdl_WZ, mdl_WT, mdl_ymtau, mdl_ymt, mdl_ymb, aS, + mdl_Gf, aEWM1, mdl_MH, mdl_MZ, mdl_MTA, mdl_MT, mdl_MB, + mdl_conjg__CKM1x1, mdl_conjg__CKM3x3, mdl_CKM3x3, mdl_MZ__exp__2, + mdl_MZ__exp__4, mdl_sqrt__2, mdl_MH__exp__2, mdl_aEW, mdl_MW, + mdl_sqrt__aEW, mdl_ee, mdl_MW__exp__2, mdl_sw2, mdl_cw, mdl_sqrt__sw2, + mdl_sw, mdl_g1, mdl_gw, mdl_vev, mdl_vev__exp__2, mdl_lam, mdl_yb, + mdl_yt, mdl_ytau, mdl_muH, mdl_ee__exp__2, mdl_sw__exp__2, + mdl_cw__exp__2; + std::complex mdl_complexi, mdl_I1x33, mdl_I2x33, mdl_I3x33, + mdl_I4x33; + // Model parameters dependent on aS + double mdl_sqrt__aS, G, mdl_G__exp__2; + // Model couplings independent of aS + + // Model couplings dependent on aS + std::complex GC_10, GC_11, GC_12; + + // Set parameters that are unchanged during the run + void setIndependentParameters(SLHAReader& slha); + // Set couplings that are unchanged during the run + void setIndependentCouplings(); + // Set parameters that are changed event by event + void setDependentParameters(); + // Set couplings that are changed event by event + void setDependentCouplings(); + + // Print parameters that are unchanged during the run + void printIndependentParameters(); + // Print couplings that are unchanged during the run + void printIndependentCouplings(); + // Print parameters that are changed event by event + void printDependentParameters(); + // Print couplings that are changed event by event + void printDependentCouplings(); + + + private: + static Parameters_sm * instance; +}; + +#endif // Parameters_sm_H + diff --git a/epoch1/oneapi/gg_ttgg/src/rambo.cc b/epoch1/oneapi/gg_ttgg/src/rambo.cc new file mode 100644 index 0000000000..de74a018f8 --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/src/rambo.cc @@ -0,0 +1,382 @@ +#include +#include +#include +#include +#include +#include + +#include "rambo.h" + +double Random::ranmar() { + /* ----------------- + * universal random number generator proposed by marsaglia and zaman + * in report fsu-scri-87-50 + * in this version rvec is a double precision variable. */ + double uni = ranu[iranmr] - ranu[jranmr]; + if (uni < 0) + uni = uni + 1; + ranu[iranmr] = uni; + iranmr = iranmr - 1; + jranmr = jranmr - 1; + if (iranmr < 0) + iranmr = 97; + if (jranmr < 0) + jranmr = 97; + ranc = ranc - rancd; + if (ranc < 0) + ranc = ranc + rancm; + uni = uni - ranc; + if (uni < 0) + uni = uni + 1; + return uni; +} + +void Random::rmarin(int ij, int kl) { + /* ----------------- + * initializing routine for ranmar, must be called before generating + * any pseudorandom numbers with ranmar. the input values should be in + * the ranges 0<=ij<=31328 ; 0<=kl<=30081 */ + /* this shows correspondence between the simplified input seeds ij, kl + * and the original marsaglia-zaman seeds i,j,k,l. + * to get the standard values in the marsaglia-zaman paper (i=12,j=34 + * k=56,l=78) put ij=1802, kl=9373 */ + int i = ij / 177 % 177 + 2; + int j = ij % 177 + 2; + int k = (kl / 169) % 178 + 1; + int l = kl % 169; + for (int ii = 1; ii < 98; ii++) { + double s = 0; + double t = .5; + for (int jj = 1; jj < 25; jj++) { + int m = ((i * j % 179) * k) % 179; + i = j; + j = k; + k = m; + l = (53 * l + 1) % 169; + if ((l * m) % 64 >= 32) + s = s + t; + t = .5 * t; + } + ranu[ii] = s; + } + ranc = 362436. / 16777216.; + rancd = 7654321. / 16777216.; + rancm = 16777213. / 16777216.; + iranmr = 97; + jranmr = 33; +} + +double rn(int idummy) { + static Random rand; + double ran; + static int init = 1; + // Prevent unused variable warning + if (false) + idummy = idummy; + if (init == 1) { + init = 0; + rand.rmarin(1802, 9373); + } + + while (true) { + ran = rand.ranmar(); + if (ran > 1e-16) + break; + } + return ran; +} + + +bool pass_cuts(const std::vector& pout){ + // return true if the events pass the cuts. (fails if not) + // two cut are implemented + // - pt cut on gluon (particle 3 and 4): 20 GeV + // - delta_R between the gluon (0.1) + + // pt cut on first gluon + double pt2 =pout[2][1]*pout[2][1] + pout[2][2]*pout[2][2]; + if (pt2 < 400){ + return false; + } + + // pt cut on second gluon + pt2 =pout[3][1]*pout[3][1] + pout[3][2]*pout[3][2]; + if (pt2 < 400){ + return false; + } + + //delta_R between the gluon + // delta_R**2 = delta_phi**2 + delta_eta**2 + // computing delta_phi first + double denom = sqrt(pout[2][1]*pout[2][1] + pout[2][2]*pout[2][2]) * sqrt(pout[3][1]*pout[3][1] + pout[3][2]*pout[3][2]); + double delta_phi = acos((pout[2][1]*pout[3][1] + pout[2][2]*pout[3][2]) / denom); + + // computing delta_eta (note rap1/2 are not lorentz invariant but the difference is. + // This evaluates rapidity in the center of mass frame (not the same as lab frame) + double rap1 = 0.5 * log( (pout[2][0] + pout[2][3])/(pout[2][0] - pout[2][3])); + double rap2 = 0.5 * log( (pout[3][0] + pout[3][3])/(pout[3][0] - pout[3][3])); + + //computing deltaR and check bounds + double delta_rap2 = delta_phi*delta_phi + (rap1 -rap2) * (rap1 - rap2); + if (delta_rap2 < 0.01){ + return false; + } + + + return true; +} + +std::vector> get_momenta(int ninitial, double energy, + std::vector masses, + double &wgt, int dim) { + //---- auxiliary function to change convention between MadGraph5_aMC@NLO and + // rambo + //---- four momenta. + int nexternal = masses.size(); + int nfinal = nexternal - ninitial; + double e2 = pow(energy, 2); + double m1 = masses[0]; + std::vector> p2; + + if (ninitial == 1) { + int done = 0; + while (done < dim){ + // for (int d = 0; d < dim; ++d) { + // Momenta for the incoming particle + std::vector p(1, new double[4]); + p[0][0] = m1; + p[0][1] = 0.; + p[0][2] = 0.; + p[0][3] = 0.; + + std::vector finalmasses(++masses.begin(), masses.end()); + std::vector p_rambo = rambo(m1, finalmasses, wgt); + if (pass_cuts(p_rambo)){ + p.insert(++p.begin(), p_rambo.begin(), p_rambo.end()); + p2.push_back(p); + done++; + } + } + return p2; + } + + else if (ninitial != 2) { + std::cout << "Rambo needs 1 or 2 incoming particles" << std::endl; + exit(-1); + } + + if (nfinal == 1) + energy = m1; + + double m2 = masses[1]; + + double mom = + sqrt((pow(e2, 2) - 2 * e2 * pow(m1, 2) + pow(m1, 4) - + 2 * e2 * pow(m2, 2) - 2 * pow(m1, 2) * pow(m2, 2) + pow(m2, 4)) / + (4 * e2)); + double energy1 = sqrt(pow(mom, 2) + pow(m1, 2)); + double energy2 = sqrt(pow(mom, 2) + pow(m2, 2)); + // Set momenta for incoming particles + for (int d = 0; d < dim; ++d) { + + std::vector p(1, new double[4]); + p[0][0] = energy1; + p[0][1] = 0; + p[0][2] = 0; + p[0][3] = mom; + p.push_back(new double[4]); + p[1][0] = energy2; + p[1][1] = 0; + p[1][2] = 0; + p[1][3] = -mom; + + if (nfinal == 1) { + p.push_back(new double[4]); + p[2][0] = energy; + wgt = 1; + p2.push_back(p); + } + std::vector finalmasses(++(++masses.begin()), masses.end()); + std::vector p_rambo = rambo(energy, finalmasses, wgt); + p.insert(++(++p.begin()), p_rambo.begin(), p_rambo.end()); + p2.push_back(p); + } + return p2; +} + +std::vector rambo(double et, std::vector &xm, double &wt) { + /********************************************************************** + * rambo * + * ra(ndom) m(omenta) b(eautifully) o(rganized) * + * * + * a democratic multi-particle phase space generator * + * authors: s.d. ellis, r. kleiss, w.j. stirling * + * this is version 1.0 - written by r. kleiss * + * -- adjusted by hans kuijf, weights are logarithmic (20-08-90) * + * * + * n = number of particles * + * et = total centre-of-mass energy * + * xm = particle masses ( dim=nexternal-nincoming ) * + * p = particle momenta ( dim=(4,nexternal-nincoming) ) * + * wt = weight of the event * + ***********************************************************************/ + int n = xm.size(); + std::vector q, p; + std::vector z(n), r(4), b(3), p2(n), xm2(n), e(n), v(n); + static std::vector iwarn(5, 0); + static double acc = 1e-14; + static int itmax = 6, ibegin = 0; + static double twopi = 8. * atan(1.); + static double po2log = log(twopi / 4.); + + for (int i = 0; i < n; i++) { + q.push_back(new double[4]); + p.push_back(new double[4]); + } + // initialization step: factorials for the phase space weight + if (ibegin == 0) { + ibegin = 1; + z[1] = po2log; + for (int k = 2; k < n; k++) + z[k] = z[k - 1] + po2log - 2. * log(double(k - 1)); + for (int k = 2; k < n; k++) + z[k] = (z[k] - log(double(k))); + } + // check on the number of particles + if (n < 1 || n > 101) { + std::cout << "Too few or many particles: " << n << std::endl; + exit(-1); + } + // check whether total energy is sufficient; count nonzero masses + double xmt = 0.; + int nm = 0; + for (int i = 0; i < n; i++) { + if (xm[i] != 0.) + nm = nm + 1; + xmt = xmt + sycl::abs(xm[i]); + } + if (xmt > et) { + std::cout << "Too low energy: " << et << " needed " << xmt << std::endl; + exit(-1); + } + // the parameter values are now accepted + + // generate n massless momenta in infinite phase space + for (int i = 0; i < n; i++) { + double r1 = rn(1); + double c = 2. * r1 - 1.; + double s = sqrt(1. - c * c); + double f = twopi * rn(2); + r1 = rn(3); + double r2 = rn(4); + q[i][0] = -log(r1 * r2); + q[i][3] = q[i][0] * c; + q[i][2] = q[i][0] * s * cos(f); + q[i][1] = q[i][0] * s * sin(f); + } + // calculate the parameters of the conformal transformation + for (int i = 0; i < 4; i++) + r[i] = 0.; + for (int i = 0; i < n; i++) { + for (int k = 0; k < 4; k++) + r[k] = r[k] + q[i][k]; + } + double rmas = sqrt(pow(r[0], 2) - pow(r[3], 2) - pow(r[2], 2) - pow(r[1], 2)); + for (int k = 1; k < 4; k++) + b[k - 1] = -r[k] / rmas; + double g = r[0] / rmas; + double a = 1. / (1. + g); + double x = et / rmas; + + // transform the q's conformally into the p's + for (int i = 0; i < n; i++) { + double bq = b[0] * q[i][1] + b[1] * q[i][2] + b[2] * q[i][3]; + for (int k = 1; k < 4; k++) + p[i][k] = x * (q[i][k] + b[k - 1] * (q[i][0] + a * bq)); + p[i][0] = x * (g * q[i][0] + bq); + } + + for (int i = 0; i < n; ++i) { + delete[] q[i]; + } + + // calculate weight and possible warnings + wt = po2log; + if (n != 2) + wt = (2. * n - 4.) * log(et) + z[n - 1]; + if (wt < -180.) { + if (iwarn[0] <= 5) + std::cout << "Too small wt, risk for underflow: " << wt << std::endl; + iwarn[0] = iwarn[0] + 1; + } + if (wt > 174.) { + if (iwarn[1] <= 5) + std::cout << "Too large wt, risk for overflow: " << wt << std::endl; + iwarn[1] = iwarn[1] + 1; + } + + // return for weighted massless momenta + if (nm == 0) { + // return log of weight + return p; + } + + // massive particles: rescale the momenta by a factor x + double xmax = sqrt(1. - pow(xmt / et, 2)); + for (int i = 0; i < n; i++) { + xm2[i] = pow(xm[i], 2); + p2[i] = pow(p[i][0], 2); + } + int iter = 0; + x = xmax; + double accu = et * acc; + while (true) { + double f0 = -et; + double g0 = 0.; + double x2 = x * x; + for (int i = 0; i < n; i++) { + e[i] = sqrt(xm2[i] + x2 * p2[i]); + f0 = f0 + e[i]; + g0 = g0 + p2[i] / e[i]; + } + if (sycl::abs(f0) <= accu) + break; + iter = iter + 1; + if (iter > itmax) { + std::cout << "Too many iterations without desired accuracy: " << itmax + << std::endl; + break; + } + x = x - f0 / (x * g0); + } + for (int i = 0; i < n; i++) { + v[i] = x * p[i][0]; + for (int k = 1; k < 4; k++) + p[i][k] = x * p[i][k]; + p[i][0] = e[i]; + } + + // calculate the mass-effect weight factor + double wt2 = 1.; + double wt3 = 0.; + for (int i = 0; i < n; i++) { + wt2 = wt2 * v[i] / e[i]; + wt3 = wt3 + pow(v[i], 2) / e[i]; + } + double wtm = (2. * n - 3.) * log(x) + log(wt2 / wt3 * et); + + // return for weighted massive momenta + wt = wt + wtm; + if (wt < -180.) { + if (iwarn[2] <= 5) + std::cout << "Too small wt, risk for underflow: " << wt << std::endl; + iwarn[2] = iwarn[2] + 1; + } + if (wt > 174.) { + if (iwarn[3] <= 5) + std::cout << "Too large wt, risk for overflow: " << wt << std::endl; + iwarn[3] = iwarn[3] + 1; + } + // return log of weight + return p; +} diff --git a/epoch1/oneapi/gg_ttgg/src/rambo.h b/epoch1/oneapi/gg_ttgg/src/rambo.h new file mode 100644 index 0000000000..6c675e80d7 --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/src/rambo.h @@ -0,0 +1,20 @@ +#include + +class Random { +public: + double ranmar(); + void rmarin(int ij, int kl); + +private: + double ranu[98]; + double ranc, rancd, rancm; + int iranmr, jranmr; +}; + +double rn(int idummy); + +std::vector> get_momenta(int ninitial, double energy, + std::vector masses, + double &wgt, int dim); + +std::vector rambo(double et, std::vector &xm, double &wt); diff --git a/epoch1/oneapi/gg_ttgg/src/read_slha.cc b/epoch1/oneapi/gg_ttgg/src/read_slha.cc new file mode 100644 index 0000000000..f83b2b6056 --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/src/read_slha.cc @@ -0,0 +1,144 @@ +#include +#include +#include +#include "read_slha.h" + +void SLHABlock::set_entry(std::vector indices, double value) +{ + if (_entries.size() == 0) + _indices = indices.size(); + else if(indices.size() != _indices) + throw "Wrong number of indices in set_entry"; + + _entries[indices] = value; +} + +double SLHABlock::get_entry(std::vector indices, double def_val) +{ + if (_entries.find(indices) == _entries.end()){ + std::cout << "Warning: No such entry in " << _name << ", using default value " + << def_val << std::endl; + return def_val; + } + return _entries[indices]; +} + +void SLHAReader::read_slha_file(std::string file_name) +{ + std::ifstream param_card; + param_card.open(file_name.c_str(), std::ifstream::in); + if(!param_card.good()) + throw "Error while opening param card"; + std::cout << "Opened slha file " << file_name << " for reading" << std::endl; + char buf[200]; + std::string line; + std::string block(""); + + while(param_card.good()){ + param_card.getline(buf, 200); + line = buf; + // Change to lowercase + transform(line.begin(), line.end(), line.begin(), (int(*)(int)) tolower); + if(line != "" && line[0] != '#'){ + if(block != ""){ + // Look for double index blocks + double dindex1, dindex2; + double value; + std::stringstream linestr2(line); + if (linestr2 >> dindex1 >> dindex2 >> value && + dindex1 == int(dindex1) and dindex2 == int(dindex2)) + { + std::vector indices; + indices.push_back(int(dindex1)); + indices.push_back(int(dindex2)); + set_block_entry(block, indices, value); + // Done with this line, read next + continue; + } + std::stringstream linestr1(line); + // Look for single index blocks + if(linestr1 >> dindex1 >> value && dindex1 == int(dindex1)) + { + std::vector indices; + indices.push_back(int(dindex1)); + set_block_entry(block, indices, value); + // Done with this line, read next + continue; + } + } + // Look for block + if(line.find("block ") != line.npos){ + line = line.substr(6); + // Get rid of spaces between block and block name + while (line[0] == ' ') + line = line.substr(1); + // Now find end of block name + int space_pos = line.find(' '); + if(space_pos != ((int) line.npos)) + line = line.substr(0, space_pos); + block = line; + continue; + } + // Look for decay + if(line.find("decay ") == 0){ + line = line.substr(6); + block = ""; + std::stringstream linestr(line); + int pdg_code; + double value; + if(linestr >> pdg_code >> value) + set_block_entry("decay", pdg_code, value); + else + std::cout << "Warning: Wrong format for decay block " << line << std::endl; + continue; + } + } + } + + if (_blocks.size() == 0) + throw "No information read from SLHA card"; + + param_card.close(); +} + +double SLHAReader::get_block_entry(std::string block_name, std::vector indices, + double def_val) +{ + if (_blocks.find(block_name) == _blocks.end()){ + std::cout << "No such block " << block_name << ", using default value " + << def_val << std::endl; + return def_val; + } + return _blocks[block_name].get_entry(indices); +} + +double SLHAReader::get_block_entry(std::string block_name, int index, + double def_val) +{ + std::vector indices; + indices.push_back(index); + return get_block_entry(block_name, indices, def_val); +} + + +void SLHAReader::set_block_entry(std::string block_name, std::vector indices, + double value) +{ + if (_blocks.find(block_name) == _blocks.end()){ + SLHABlock block(block_name); + _blocks[block_name] = block; + } + _blocks[block_name].set_entry(indices, value); + /* cout << "Set block " << block_name << " entry "; + for (int i=0;i < indices.size();i++) + cout << indices[i] << " "; + cout << "to " << _blocks[block_name].get_entry(indices) << endl;*/ +} + +void SLHAReader::set_block_entry(std::string block_name, int index, + double value) +{ + std::vector indices; + indices.push_back(index); + set_block_entry(block_name, indices, value); +} diff --git a/epoch1/oneapi/gg_ttgg/src/read_slha.h b/epoch1/oneapi/gg_ttgg/src/read_slha.h new file mode 100644 index 0000000000..658a51049d --- /dev/null +++ b/epoch1/oneapi/gg_ttgg/src/read_slha.h @@ -0,0 +1,46 @@ +#ifndef READ_SLHA_H +#define READ_SLHA_H + +#include +#include +#include +#include + +class SLHABlock +{ + public: + SLHABlock(std::string name = ""){_name = name;} + ~SLHABlock(){} + + void set_entry(std::vector indices, double value); + double get_entry(std::vector indices, double def_val = 0); + void set_name(std::string name) {_name = name;} + std::string get_name(){return _name;} + int get_indices() { return _indices;} + + private: + std::string _name; + std::map, double> _entries; + unsigned int _indices; +}; + +class SLHAReader +{ + public: + SLHAReader(std::string file_name = "") + {if(file_name != "") read_slha_file(file_name);} + + void read_slha_file(std::string file_name); + double get_block_entry(std::string block_name, std::vector indices, + double def_val = 0); + double get_block_entry(std::string block_name, int index, + double def_val = 0); + void set_block_entry(std::string block_name, std::vector indices, + double value); + void set_block_entry(std::string block_name, int index, + double value); + private: + std::map _blocks; +}; + +#endif diff --git a/epoch1/oneapi/remove.me b/epoch1/oneapi/remove.me deleted file mode 100644 index f05e81d02a..0000000000 --- a/epoch1/oneapi/remove.me +++ /dev/null @@ -1 +0,0 @@ -placeholder file, git rm it once the directory is filled with code