diff --git a/MELA/COLLIER/setup.sh b/MELA/COLLIER/setup.sh
index 575e9c3d..20db717b 100644
--- a/MELA/COLLIER/setup.sh
+++ b/MELA/COLLIER/setup.sh
@@ -2,6 +2,8 @@
 
 {
 
+set -euo pipefail
+
 scriptdir=$(dirname $0)
 curdir=$(pwd)
 
@@ -28,22 +30,24 @@ if [[ $# > 0 ]] && [[ "$1" == *"clean"* ]];then
     fi
   done
 
-  rm -f "../data/"$SCRAM_ARCH"/"$libname
+  rm -f ../data/*/$libname
 
 else
 
-  wget $tarweb
-  mkdir $tmpdir
-  tar -xvzf $tarname -C $tmpdir
-  rm $tarname
-  mv $tmpdir"/"$pkgdir"/src/"* ./
-  rm -rf $tmpdir
+  if ! [ -f "../data/$SCRAM_ARCH/$libname" ]; then
+    wget $tarweb
+    mkdir $tmpdir
+    tar -xvzf $tarname -C $tmpdir
+    rm $tarname
+    mv $tmpdir"/"$pkgdir"/src/"* ./
+    rm -rf $tmpdir
 
-  make
-  mv $libname "../data/"$SCRAM_ARCH"/"$libname
+    make
+    cp $libname "../data/$SCRAM_ARCH/$libname"
+  fi
 
 fi
 
 cd $curdir
 
-}
\ No newline at end of file
+}
diff --git a/MELA/data/slc6_amd64_gcc700/download.url b/MELA/data/slc6_amd64_gcc700/download.url
new file mode 100644
index 00000000..4fbf4729
--- /dev/null
+++ b/MELA/data/slc6_amd64_gcc700/download.url
@@ -0,0 +1 @@
+http://spin.pha.jhu.edu/Generator/MCFM-precompiled/v2/
diff --git a/MELA/data/slc7_amd64_gcc630/download.url b/MELA/data/slc7_amd64_gcc630/download.url
new file mode 100644
index 00000000..4fbf4729
--- /dev/null
+++ b/MELA/data/slc7_amd64_gcc630/download.url
@@ -0,0 +1 @@
+http://spin.pha.jhu.edu/Generator/MCFM-precompiled/v2/
diff --git a/MELA/data/slc7_amd64_gcc700/download.url b/MELA/data/slc7_amd64_gcc700/download.url
new file mode 100644
index 00000000..4fbf4729
--- /dev/null
+++ b/MELA/data/slc7_amd64_gcc700/download.url
@@ -0,0 +1 @@
+http://spin.pha.jhu.edu/Generator/MCFM-precompiled/v2/
diff --git a/MELA/fortran/mod_HashCollection.F90 b/MELA/fortran/mod_HashCollection.F90
index 6361a2ca..31643e80 100644
--- a/MELA/fortran/mod_HashCollection.F90
+++ b/MELA/fortran/mod_HashCollection.F90
@@ -5,7 +5,7 @@ MODULE ModHashCollection
 public
 
 ! MCFM qq_VVqq hash - JHU conventions
-integer, parameter :: Hash_MCFM_qqVVqq_Size = 204
+integer, parameter :: Hash_MCFM_qqVVqq_Size = 164
 integer, target :: Hash_MCFM_qqVVqq(1:Hash_MCFM_qqVVqq_Size,1:4)
 
 ! MCFM qq_VVqq hash (generation) - PDF/PDG conventions
@@ -71,191 +71,155 @@ subroutine Init_Hash_MCFM_qqVVqq()
    Hash_MCFM_qqVVqq(22,:) = (/ Str_ , Up_ , Up_ , Str_ /)
    Hash_MCFM_qqVVqq(23,:) = (/ Bot_ , Up_ , Up_ , Bot_ /)
    Hash_MCFM_qqVVqq(24,:) = (/ Bot_ , Chm_ , Chm_ , Bot_ /)
+
    Hash_MCFM_qqVVqq(25,:) = (/ Chm_ , Dn_ , Dn_ , Chm_ /)
    Hash_MCFM_qqVVqq(26,:) = (/ Up_ , Dn_ , Dn_ , Up_ /)
    Hash_MCFM_qqVVqq(27,:) = (/ Chm_ , Str_ , Str_ , Chm_ /)
    Hash_MCFM_qqVVqq(28,:) = (/ Chm_ , Dn_ , Up_ , Str_ /)
    Hash_MCFM_qqVVqq(29,:) = (/ Str_ , Up_ , Dn_ , Chm_ /)
-   Hash_MCFM_qqVVqq(30,:) = (/ Up_ , Up_ , Up_ , Up_ /)
-   Hash_MCFM_qqVVqq(31,:) = (/ Chm_ , Chm_ , Chm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(32,:) = (/ Dn_ , Dn_ , Dn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(33,:) = (/ Str_ , Str_ , Str_ , Str_ /)
-   Hash_MCFM_qqVVqq(34,:) = (/ Bot_ , Bot_ , Bot_ , Bot_ /)
-
-   Hash_MCFM_qqVVqq(35,:) = (/ AChm_ , AUp_ , AChm_ , AUp_ /)
-   Hash_MCFM_qqVVqq(36,:) = (/ AStr_ , ADn_ , AStr_ , ADn_ /)
-   Hash_MCFM_qqVVqq(37,:) = (/ ABot_ , ADn_ , ABot_ , ADn_ /)
-   Hash_MCFM_qqVVqq(38,:) = (/ ABot_ , AStr_ , ABot_ , AStr_ /)
-   Hash_MCFM_qqVVqq(39,:) = (/ AStr_ , AUp_ , AStr_ , AUp_ /)
-   Hash_MCFM_qqVVqq(40,:) = (/ ABot_ , AUp_ , ABot_ , AUp_ /)
-   Hash_MCFM_qqVVqq(41,:) = (/ ABot_ , AChm_ , ABot_ , AChm_ /)
-   Hash_MCFM_qqVVqq(42,:) = (/ AChm_ , ADn_ , AChm_ , ADn_ /)
-   Hash_MCFM_qqVVqq(43,:) = (/ AUp_ , ADn_ , AUp_ , ADn_ /)
-   Hash_MCFM_qqVVqq(44,:) = (/ AChm_ , AStr_ , AChm_ , AStr_ /)
-   Hash_MCFM_qqVVqq(45,:) = (/ AStr_ , AUp_ , AChm_ , ADn_ /)
-   Hash_MCFM_qqVVqq(46,:) = (/ AChm_ , ADn_ , AStr_ , AUp_ /)
-   Hash_MCFM_qqVVqq(47,:) = (/ AUp_ , AUp_ , AUp_ , AUp_ /)
-   Hash_MCFM_qqVVqq(48,:) = (/ AChm_ , AChm_ , AChm_ , AChm_ /)
-   Hash_MCFM_qqVVqq(49,:) = (/ ADn_ , ADn_ , ADn_ , ADn_ /)
-   Hash_MCFM_qqVVqq(50,:) = (/ AStr_ , AStr_ , AStr_ , AStr_ /)
-   Hash_MCFM_qqVVqq(51,:) = (/ ABot_ , ABot_ , ABot_ , ABot_ /)
-   Hash_MCFM_qqVVqq(52,:) = (/ AUp_ , AChm_ , AChm_ , AUp_ /)
-   Hash_MCFM_qqVVqq(53,:) = (/ ADn_ , AStr_ , AStr_ , ADn_ /)
-   Hash_MCFM_qqVVqq(54,:) = (/ ADn_ , ABot_ , ABot_ , ADn_ /)
-   Hash_MCFM_qqVVqq(55,:) = (/ AStr_ , ABot_ , ABot_ , AStr_ /)
-   Hash_MCFM_qqVVqq(56,:) = (/ AUp_ , AStr_ , AStr_ , AUp_ /)
-   Hash_MCFM_qqVVqq(57,:) = (/ AUp_ , ABot_ , ABot_ , AUp_ /)
-   Hash_MCFM_qqVVqq(58,:) = (/ AChm_ , ABot_ , ABot_ , AChm_ /)
-   Hash_MCFM_qqVVqq(59,:) = (/ ADn_ , AChm_ , AChm_ , ADn_ /)
-   Hash_MCFM_qqVVqq(60,:) = (/ ADn_ , AUp_ , AUp_ , ADn_ /)
-   Hash_MCFM_qqVVqq(61,:) = (/ AStr_ , AChm_ , AChm_ , AStr_ /)
-   Hash_MCFM_qqVVqq(62,:) = (/ AUp_ , AStr_ , AChm_ , ADn_ /)
-   Hash_MCFM_qqVVqq(63,:) = (/ ADn_ , AChm_ , AStr_ , AUp_ /)
-   Hash_MCFM_qqVVqq(64,:) = (/ AUp_ , AUp_ , AUp_ , AUp_ /)
-   Hash_MCFM_qqVVqq(65,:) = (/ AChm_ , AChm_ , AChm_ , AChm_ /)
-   Hash_MCFM_qqVVqq(66,:) = (/ ADn_ , ADn_ , ADn_ , ADn_ /)
-   Hash_MCFM_qqVVqq(67,:) = (/ AStr_ , AStr_ , AStr_ , AStr_ /)
-   Hash_MCFM_qqVVqq(68,:) = (/ ABot_ , ABot_ , ABot_ , ABot_ /)
-
-   Hash_MCFM_qqVVqq(69,:) = (/ AUp_ , Chm_ , AUp_ , Chm_ /)
-   Hash_MCFM_qqVVqq(70,:) = (/ ADn_ , Str_ , ADn_ , Str_ /)
-   Hash_MCFM_qqVVqq(71,:) = (/ ADn_ , Bot_ , ADn_ , Bot_ /)
-   Hash_MCFM_qqVVqq(72,:) = (/ AStr_ , Bot_ , AStr_ , Bot_ /)
-   Hash_MCFM_qqVVqq(73,:) = (/ AUp_ , Str_ , AUp_ , Str_ /)
-   Hash_MCFM_qqVVqq(74,:) = (/ AUp_ , Bot_ , AUp_ , Bot_ /)
-   Hash_MCFM_qqVVqq(75,:) = (/ AChm_ , Bot_ , AChm_ , Bot_ /)
-   Hash_MCFM_qqVVqq(76,:) = (/ ADn_ , Chm_ , ADn_ , Chm_ /)
-   Hash_MCFM_qqVVqq(77,:) = (/ ADn_ , Up_ , ADn_ , Up_ /)
-   Hash_MCFM_qqVVqq(78,:) = (/ AStr_ , Chm_ , AStr_ , Chm_ /)
-   Hash_MCFM_qqVVqq(79,:) = (/ AUp_ , Chm_ , ADn_ , Str_ /)
-   Hash_MCFM_qqVVqq(80,:) = (/ ADn_ , Str_ , AUp_ , Chm_ /)
-   Hash_MCFM_qqVVqq(81,:) = (/ AUp_ , Up_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(82,:) = (/ AChm_ , Chm_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(83,:) = (/ ADn_ , Dn_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(84,:) = (/ AStr_ , Str_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(85,:) = (/ ABot_ , Bot_ , ABot_ , Bot_ /)
-   Hash_MCFM_qqVVqq(86,:) = (/ AChm_ , Up_ , AChm_ , Up_ /)
-   Hash_MCFM_qqVVqq(87,:) = (/ AStr_ , Dn_ , AStr_ , Dn_ /)
-   Hash_MCFM_qqVVqq(88,:) = (/ ABot_ , Dn_ , ABot_ , Dn_ /)
-   Hash_MCFM_qqVVqq(89,:) = (/ ABot_ , Str_ , ABot_ , Str_ /)
-   Hash_MCFM_qqVVqq(90,:) = (/ AStr_ , Up_ , AStr_ , Up_ /)
-   Hash_MCFM_qqVVqq(91,:) = (/ ABot_ , Up_ , ABot_ , Up_ /)
-   Hash_MCFM_qqVVqq(92,:) = (/ ABot_ , Chm_ , ABot_ , Chm_ /)
-   Hash_MCFM_qqVVqq(93,:) = (/ AChm_ , Dn_ , AChm_ , Dn_ /)
-   Hash_MCFM_qqVVqq(94,:) = (/ AUp_ , Dn_ , AUp_ , Dn_ /)
-   Hash_MCFM_qqVVqq(95,:) = (/ AChm_ , Str_ , AChm_ , Str_ /)
-   Hash_MCFM_qqVVqq(96,:) = (/ AStr_ , Dn_ , AChm_ , Up_ /)
-   Hash_MCFM_qqVVqq(97,:) = (/ AChm_ , Up_ , AStr_ , Dn_ /)
-   Hash_MCFM_qqVVqq(98,:) = (/ AUp_ , Up_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(99,:) = (/ AChm_ , Chm_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(100,:) = (/ ADn_ , Dn_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(101,:) = (/ AStr_ , Str_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(102,:) = (/ ABot_ , Bot_ , ABot_ , Bot_ /)
-
-   Hash_MCFM_qqVVqq(103,:) = (/ Chm_ , AUp_ , AUp_ , Chm_ /)
-   Hash_MCFM_qqVVqq(104,:) = (/ Str_ , ADn_ , ADn_ , Str_ /)
-   Hash_MCFM_qqVVqq(105,:) = (/ Bot_ , ADn_ , ADn_ , Bot_ /)
-   Hash_MCFM_qqVVqq(106,:) = (/ Bot_ , AStr_ , AStr_ , Bot_ /)
-   Hash_MCFM_qqVVqq(107,:) = (/ Str_ , AUp_ , AUp_ , Str_ /)
-   Hash_MCFM_qqVVqq(108,:) = (/ Bot_ , AUp_ , AUp_ , Bot_ /)
-   Hash_MCFM_qqVVqq(109,:) = (/ Bot_ , AChm_ , AChm_ , Bot_ /)
-   Hash_MCFM_qqVVqq(110,:) = (/ Chm_ , ADn_ , ADn_ , Chm_ /)
-   Hash_MCFM_qqVVqq(111,:) = (/ Up_ , ADn_ , ADn_ , Up_ /)
-   Hash_MCFM_qqVVqq(112,:) = (/ Chm_ , AStr_ , AStr_ , Chm_ /)
-   Hash_MCFM_qqVVqq(113,:) = (/ Chm_ , AUp_ , ADn_ , Str_ /)
-   Hash_MCFM_qqVVqq(114,:) = (/ Str_ , ADn_ , AUp_ , Chm_ /)
-   Hash_MCFM_qqVVqq(115,:) = (/ Up_ , AUp_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(116,:) = (/ Chm_ , AChm_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(117,:) = (/ Dn_ , ADn_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(118,:) = (/ Str_ , AStr_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(119,:) = (/ Bot_ , ABot_ , ABot_ , Bot_ /)
-   Hash_MCFM_qqVVqq(120,:) = (/ Up_ , AChm_ , AChm_ , Up_ /)
-   Hash_MCFM_qqVVqq(121,:) = (/ Dn_ , AStr_ , AStr_ , Dn_ /)
-   Hash_MCFM_qqVVqq(122,:) = (/ Dn_ , ABot_ , ABot_ , Dn_ /)
-   Hash_MCFM_qqVVqq(123,:) = (/ Str_ , ABot_ , ABot_ , Str_ /)
-   Hash_MCFM_qqVVqq(124,:) = (/ Up_ , AStr_ , AStr_ , Up_ /)
-   Hash_MCFM_qqVVqq(125,:) = (/ Up_ , ABot_ , ABot_ , Up_ /)
-   Hash_MCFM_qqVVqq(126,:) = (/ Chm_ , ABot_ , ABot_ , Chm_ /)
-   Hash_MCFM_qqVVqq(127,:) = (/ Dn_ , AChm_ , AChm_ , Dn_ /)
-   Hash_MCFM_qqVVqq(128,:) = (/ Dn_ , AUp_ , AUp_ , Dn_ /)
-   Hash_MCFM_qqVVqq(129,:) = (/ Str_ , AChm_ , AChm_ , Str_ /)
-   Hash_MCFM_qqVVqq(130,:) = (/ Dn_ , AStr_ , AChm_ , Up_ /)
-   Hash_MCFM_qqVVqq(131,:) = (/ Up_ , AChm_ , AStr_ , Dn_ /)
-   Hash_MCFM_qqVVqq(132,:) = (/ Up_ , AUp_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(133,:) = (/ Chm_ , AChm_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(134,:) = (/ Dn_ , ADn_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(135,:) = (/ Str_ , AStr_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(136,:) = (/ Bot_ , ABot_ , ABot_ , Bot_ /)
-
-   Hash_MCFM_qqVVqq(137,:) = (/ Up_ , AUp_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(138,:) = (/ Dn_ , ADn_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(139,:) = (/ Dn_ , ADn_ , ABot_ , Bot_ /)
-   Hash_MCFM_qqVVqq(140,:) = (/ Str_ , AStr_ , ABot_ , Bot_ /)
-   Hash_MCFM_qqVVqq(141,:) = (/ Up_ , AUp_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(142,:) = (/ Up_ , AUp_ , ABot_ , Bot_ /)
-   Hash_MCFM_qqVVqq(143,:) = (/ Chm_ , AChm_ , ABot_ , Bot_ /)
-   Hash_MCFM_qqVVqq(144,:) = (/ Dn_ , ADn_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(145,:) = (/ Dn_ , ADn_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(146,:) = (/ Str_ , AStr_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(147,:) = (/ Dn_ , AUp_ , AChm_ , Str_ /)
-   Hash_MCFM_qqVVqq(148,:) = (/ Up_ , ADn_ , AStr_ , Chm_ /)
-   Hash_MCFM_qqVVqq(149,:) = (/ Up_ , AUp_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(150,:) = (/ Chm_ , AChm_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(151,:) = (/ Dn_ , ADn_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(152,:) = (/ Str_ , AStr_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(153,:) = (/ Bot_ , ABot_ , ABot_ , Bot_ /)
-   Hash_MCFM_qqVVqq(154,:) = (/ Chm_ , AChm_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(155,:) = (/ Str_ , AStr_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(156,:) = (/ Bot_ , ABot_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(157,:) = (/ Bot_ , ABot_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(158,:) = (/ Str_ , AStr_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(159,:) = (/ Bot_ , ABot_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(160,:) = (/ Bot_ , ABot_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(161,:) = (/ Chm_ , AChm_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(162,:) = (/ Up_ , AUp_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(163,:) = (/ Chm_ , AChm_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(164,:) = (/ Chm_ , AStr_ , ADn_ , Up_ /)
-   Hash_MCFM_qqVVqq(165,:) = (/ Str_ , AChm_ , AUp_ , Dn_ /)
-   Hash_MCFM_qqVVqq(166,:) = (/ Up_ , AUp_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(167,:) = (/ Chm_ , AChm_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(168,:) = (/ Dn_ , ADn_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(169,:) = (/ Str_ , AStr_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(170,:) = (/ Bot_ , ABot_ , ABot_ , Bot_ /)
-
-   Hash_MCFM_qqVVqq(171,:) = (/ AUp_ , Up_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(172,:) = (/ ADn_ , Dn_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(173,:) = (/ ADn_ , Dn_ , ABot_ , Bot_ /)
-   Hash_MCFM_qqVVqq(174,:) = (/ AStr_ , Str_ , ABot_ , Bot_ /)
-   Hash_MCFM_qqVVqq(175,:) = (/ AUp_ , Up_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(176,:) = (/ AUp_ , Up_ , ABot_ , Bot_ /)
-   Hash_MCFM_qqVVqq(177,:) = (/ AChm_ , Chm_ , ABot_ , Bot_ /)
-   Hash_MCFM_qqVVqq(178,:) = (/ ADn_ , Dn_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(179,:) = (/ ADn_ , Dn_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(180,:) = (/ AStr_ , Str_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(181,:) = (/ AUp_ , Dn_ , AChm_ , Str_ /)
-   Hash_MCFM_qqVVqq(182,:) = (/ ADn_ , Up_ , AStr_ , Chm_ /)
-   Hash_MCFM_qqVVqq(183,:) = (/ AUp_ , Up_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(184,:) = (/ AChm_ , Chm_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(185,:) = (/ ADn_ , Dn_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(186,:) = (/ AStr_ , Str_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(187,:) = (/ ABot_ , Bot_ , ABot_ , Bot_ /)
-   Hash_MCFM_qqVVqq(188,:) = (/ AChm_ , Chm_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(189,:) = (/ AStr_ , Str_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(190,:) = (/ ABot_ , Bot_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(191,:) = (/ ABot_ , Bot_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(192,:) = (/ AStr_ , Str_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(193,:) = (/ ABot_ , Bot_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(194,:) = (/ ABot_ , Bot_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(195,:) = (/ AChm_ , Chm_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(196,:) = (/ AUp_ , Up_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(197,:) = (/ AChm_ , Chm_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(198,:) = (/ AStr_ , Chm_ , ADn_ , Up_ /)
-   Hash_MCFM_qqVVqq(199,:) = (/ AChm_ , Str_ , AUp_ , Dn_ /)
-   Hash_MCFM_qqVVqq(200,:) = (/ AUp_ , Up_ , AUp_ , Up_ /)
-   Hash_MCFM_qqVVqq(201,:) = (/ AChm_ , Chm_ , AChm_ , Chm_ /)
-   Hash_MCFM_qqVVqq(202,:) = (/ ADn_ , Dn_ , ADn_ , Dn_ /)
-   Hash_MCFM_qqVVqq(203,:) = (/ AStr_ , Str_ , AStr_ , Str_ /)
-   Hash_MCFM_qqVVqq(204,:) = (/ ABot_ , Bot_ , ABot_ , Bot_ /)
+
+   Hash_MCFM_qqVVqq(30,:) = (/ AChm_ , AUp_ , AChm_ , AUp_ /)
+   Hash_MCFM_qqVVqq(31,:) = (/ AStr_ , ADn_ , AStr_ , ADn_ /)
+   Hash_MCFM_qqVVqq(32,:) = (/ ABot_ , ADn_ , ABot_ , ADn_ /)
+   Hash_MCFM_qqVVqq(33,:) = (/ ABot_ , AStr_ , ABot_ , AStr_ /)
+   Hash_MCFM_qqVVqq(34,:) = (/ AStr_ , AUp_ , AStr_ , AUp_ /)
+   Hash_MCFM_qqVVqq(35,:) = (/ ABot_ , AUp_ , ABot_ , AUp_ /)
+   Hash_MCFM_qqVVqq(36,:) = (/ ABot_ , AChm_ , ABot_ , AChm_ /)
+   Hash_MCFM_qqVVqq(37,:) = (/ AChm_ , ADn_ , AChm_ , ADn_ /)
+   Hash_MCFM_qqVVqq(38,:) = (/ AUp_ , ADn_ , AUp_ , ADn_ /)
+   Hash_MCFM_qqVVqq(39,:) = (/ AChm_ , AStr_ , AChm_ , AStr_ /)
+   Hash_MCFM_qqVVqq(40,:) = (/ AStr_ , AUp_ , AChm_ , ADn_ /)
+   Hash_MCFM_qqVVqq(41,:) = (/ AChm_ , ADn_ , AStr_ , AUp_ /)
+   Hash_MCFM_qqVVqq(42,:) = (/ AUp_ , AUp_ , AUp_ , AUp_ /)
+   Hash_MCFM_qqVVqq(43,:) = (/ AChm_ , AChm_ , AChm_ , AChm_ /)
+   Hash_MCFM_qqVVqq(44,:) = (/ ADn_ , ADn_ , ADn_ , ADn_ /)
+   Hash_MCFM_qqVVqq(45,:) = (/ AStr_ , AStr_ , AStr_ , AStr_ /)
+   Hash_MCFM_qqVVqq(46,:) = (/ ABot_ , ABot_ , ABot_ , ABot_ /)
+
+   Hash_MCFM_qqVVqq(47,:) = (/ AUp_ , AChm_ , AChm_ , AUp_ /)
+   Hash_MCFM_qqVVqq(48,:) = (/ ADn_ , AStr_ , AStr_ , ADn_ /)
+   Hash_MCFM_qqVVqq(49,:) = (/ ADn_ , ABot_ , ABot_ , ADn_ /)
+   Hash_MCFM_qqVVqq(50,:) = (/ AStr_ , ABot_ , ABot_ , AStr_ /)
+   Hash_MCFM_qqVVqq(51,:) = (/ AUp_ , AStr_ , AStr_ , AUp_ /)
+   Hash_MCFM_qqVVqq(52,:) = (/ AUp_ , ABot_ , ABot_ , AUp_ /)
+   Hash_MCFM_qqVVqq(53,:) = (/ AChm_ , ABot_ , ABot_ , AChm_ /)
+   Hash_MCFM_qqVVqq(54,:) = (/ ADn_ , AChm_ , AChm_ , ADn_ /)
+   Hash_MCFM_qqVVqq(55,:) = (/ ADn_ , AUp_ , AUp_ , ADn_ /)
+   Hash_MCFM_qqVVqq(56,:) = (/ AStr_ , AChm_ , AChm_ , AStr_ /)
+   Hash_MCFM_qqVVqq(57,:) = (/ AUp_ , AStr_ , AChm_ , ADn_ /)
+   Hash_MCFM_qqVVqq(58,:) = (/ ADn_ , AChm_ , AStr_ , AUp_ /)
+
+   Hash_MCFM_qqVVqq(59,:) = (/ AUp_ , Chm_ , AUp_ , Chm_ /)
+   Hash_MCFM_qqVVqq(60,:) = (/ ADn_ , Str_ , ADn_ , Str_ /)
+   Hash_MCFM_qqVVqq(61,:) = (/ ADn_ , Bot_ , ADn_ , Bot_ /)
+   Hash_MCFM_qqVVqq(62,:) = (/ AStr_ , Bot_ , AStr_ , Bot_ /)
+   Hash_MCFM_qqVVqq(63,:) = (/ AUp_ , Str_ , AUp_ , Str_ /)
+   Hash_MCFM_qqVVqq(64,:) = (/ AUp_ , Bot_ , AUp_ , Bot_ /)
+   Hash_MCFM_qqVVqq(65,:) = (/ AChm_ , Bot_ , AChm_ , Bot_ /)
+   Hash_MCFM_qqVVqq(66,:) = (/ ADn_ , Chm_ , ADn_ , Chm_ /)
+   Hash_MCFM_qqVVqq(67,:) = (/ ADn_ , Up_ , ADn_ , Up_ /)
+   Hash_MCFM_qqVVqq(68,:) = (/ AStr_ , Chm_ , AStr_ , Chm_ /)
+   Hash_MCFM_qqVVqq(69,:) = (/ AUp_ , Chm_ , ADn_ , Str_ /)
+   Hash_MCFM_qqVVqq(70,:) = (/ ADn_ , Str_ , AUp_ , Chm_ /)
+   Hash_MCFM_qqVVqq(71,:) = (/ AUp_ , Up_ , AUp_ , Up_ /)
+   Hash_MCFM_qqVVqq(72,:) = (/ AChm_ , Chm_ , AChm_ , Chm_ /)
+   Hash_MCFM_qqVVqq(73,:) = (/ ADn_ , Dn_ , ADn_ , Dn_ /)
+   Hash_MCFM_qqVVqq(74,:) = (/ AStr_ , Str_ , AStr_ , Str_ /)
+   Hash_MCFM_qqVVqq(75,:) = (/ ABot_ , Bot_ , ABot_ , Bot_ /)
+
+   Hash_MCFM_qqVVqq(76,:) = (/ AChm_ , Up_ , AChm_ , Up_ /)
+   Hash_MCFM_qqVVqq(77,:) = (/ AStr_ , Dn_ , AStr_ , Dn_ /)
+   Hash_MCFM_qqVVqq(78,:) = (/ ABot_ , Dn_ , ABot_ , Dn_ /)
+   Hash_MCFM_qqVVqq(79,:) = (/ ABot_ , Str_ , ABot_ , Str_ /)
+   Hash_MCFM_qqVVqq(80,:) = (/ AStr_ , Up_ , AStr_ , Up_ /)
+   Hash_MCFM_qqVVqq(81,:) = (/ ABot_ , Up_ , ABot_ , Up_ /)
+   Hash_MCFM_qqVVqq(82,:) = (/ ABot_ , Chm_ , ABot_ , Chm_ /)
+   Hash_MCFM_qqVVqq(83,:) = (/ AChm_ , Dn_ , AChm_ , Dn_ /)
+   Hash_MCFM_qqVVqq(84,:) = (/ AUp_ , Dn_ , AUp_ , Dn_ /)
+   Hash_MCFM_qqVVqq(85,:) = (/ AChm_ , Str_ , AChm_ , Str_ /)
+   Hash_MCFM_qqVVqq(86,:) = (/ AStr_ , Dn_ , AChm_ , Up_ /)
+   Hash_MCFM_qqVVqq(87,:) = (/ AChm_ , Up_ , AStr_ , Dn_ /)
+
+   Hash_MCFM_qqVVqq(88,:) = (/ Chm_ , AUp_ , AUp_ , Chm_ /)
+   Hash_MCFM_qqVVqq(89,:) = (/ Str_ , ADn_ , ADn_ , Str_ /)
+   Hash_MCFM_qqVVqq(90,:) = (/ Bot_ , ADn_ , ADn_ , Bot_ /)
+   Hash_MCFM_qqVVqq(91,:) = (/ Bot_ , AStr_ , AStr_ , Bot_ /)
+   Hash_MCFM_qqVVqq(92,:) = (/ Str_ , AUp_ , AUp_ , Str_ /)
+   Hash_MCFM_qqVVqq(93,:) = (/ Bot_ , AUp_ , AUp_ , Bot_ /)
+   Hash_MCFM_qqVVqq(94,:) = (/ Bot_ , AChm_ , AChm_ , Bot_ /)
+   Hash_MCFM_qqVVqq(95,:) = (/ Chm_ , ADn_ , ADn_ , Chm_ /)
+   Hash_MCFM_qqVVqq(96,:) = (/ Up_ , ADn_ , ADn_ , Up_ /)
+   Hash_MCFM_qqVVqq(97,:) = (/ Chm_ , AStr_ , AStr_ , Chm_ /)
+   Hash_MCFM_qqVVqq(98,:) = (/ Chm_ , AUp_ , ADn_ , Str_ /)
+   Hash_MCFM_qqVVqq(99,:) = (/ Str_ , ADn_ , AUp_ , Chm_ /)
+   Hash_MCFM_qqVVqq(100,:) = (/ Up_ , AUp_ , AUp_ , Up_ /)
+   Hash_MCFM_qqVVqq(101,:) = (/ Chm_ , AChm_ , AChm_ , Chm_ /)
+   Hash_MCFM_qqVVqq(102,:) = (/ Dn_ , ADn_ , ADn_ , Dn_ /)
+   Hash_MCFM_qqVVqq(103,:) = (/ Str_ , AStr_ , AStr_ , Str_ /)
+   Hash_MCFM_qqVVqq(104,:) = (/ Bot_ , ABot_ , ABot_ , Bot_ /)
+
+   Hash_MCFM_qqVVqq(105,:) = (/ Up_ , AChm_ , AChm_ , Up_ /)
+   Hash_MCFM_qqVVqq(106,:) = (/ Dn_ , AStr_ , AStr_ , Dn_ /)
+   Hash_MCFM_qqVVqq(107,:) = (/ Dn_ , ABot_ , ABot_ , Dn_ /)
+   Hash_MCFM_qqVVqq(108,:) = (/ Str_ , ABot_ , ABot_ , Str_ /)
+   Hash_MCFM_qqVVqq(109,:) = (/ Up_ , AStr_ , AStr_ , Up_ /)
+   Hash_MCFM_qqVVqq(110,:) = (/ Up_ , ABot_ , ABot_ , Up_ /)
+   Hash_MCFM_qqVVqq(111,:) = (/ Chm_ , ABot_ , ABot_ , Chm_ /)
+   Hash_MCFM_qqVVqq(112,:) = (/ Dn_ , AChm_ , AChm_ , Dn_ /)
+   Hash_MCFM_qqVVqq(113,:) = (/ Dn_ , AUp_ , AUp_ , Dn_ /)
+   Hash_MCFM_qqVVqq(114,:) = (/ Str_ , AChm_ , AChm_ , Str_ /)
+   Hash_MCFM_qqVVqq(115,:) = (/ Dn_ , AStr_ , AChm_ , Up_ /)
+   Hash_MCFM_qqVVqq(116,:) = (/ Up_ , AChm_ , AStr_ , Dn_ /)
+   Hash_MCFM_qqVVqq(117,:) = (/ Up_ , AUp_ , AChm_ , Chm_ /)
+   Hash_MCFM_qqVVqq(118,:) = (/ Dn_ , ADn_ , AStr_ , Str_ /)
+   Hash_MCFM_qqVVqq(119,:) = (/ Dn_ , ADn_ , ABot_ , Bot_ /)
+   Hash_MCFM_qqVVqq(120,:) = (/ Str_ , AStr_ , ABot_ , Bot_ /)
+   Hash_MCFM_qqVVqq(121,:) = (/ Up_ , AUp_ , AStr_ , Str_ /)
+   Hash_MCFM_qqVVqq(122,:) = (/ Up_ , AUp_ , ABot_ , Bot_ /)
+   Hash_MCFM_qqVVqq(123,:) = (/ Chm_ , AChm_ , ABot_ , Bot_ /)
+
+   Hash_MCFM_qqVVqq(124,:) = (/ Dn_ , ADn_ , AChm_ , Chm_ /)
+   Hash_MCFM_qqVVqq(125,:) = (/ Dn_ , ADn_ , AUp_ , Up_ /)
+   Hash_MCFM_qqVVqq(126,:) = (/ Str_ , AStr_ , AChm_ , Chm_ /)
+   Hash_MCFM_qqVVqq(127,:) = (/ Dn_ , AUp_ , AChm_ , Str_ /)
+   Hash_MCFM_qqVVqq(128,:) = (/ Up_ , ADn_ , AStr_ , Chm_ /)
+   Hash_MCFM_qqVVqq(129,:) = (/ Chm_ , AChm_ , AUp_ , Up_ /)
+   Hash_MCFM_qqVVqq(130,:) = (/ Str_ , AStr_ , ADn_ , Dn_ /)
+   Hash_MCFM_qqVVqq(131,:) = (/ Bot_ , ABot_ , ADn_ , Dn_ /)
+   Hash_MCFM_qqVVqq(132,:) = (/ Bot_ , ABot_ , AStr_ , Str_ /)
+   Hash_MCFM_qqVVqq(133,:) = (/ Str_ , AStr_ , AUp_ , Up_ /)
+   Hash_MCFM_qqVVqq(134,:) = (/ Bot_ , ABot_ , AUp_ , Up_ /)
+   Hash_MCFM_qqVVqq(135,:) = (/ Bot_ , ABot_ , AChm_ , Chm_ /)
+   Hash_MCFM_qqVVqq(136,:) = (/ Chm_ , AChm_ , ADn_ , Dn_ /)
+   Hash_MCFM_qqVVqq(137,:) = (/ Up_ , AUp_ , ADn_ , Dn_ /)
+   Hash_MCFM_qqVVqq(138,:) = (/ Chm_ , AChm_ , AStr_ , Str_ /)
+   Hash_MCFM_qqVVqq(139,:) = (/ Chm_ , AStr_ , ADn_ , Up_ /)
+   Hash_MCFM_qqVVqq(140,:) = (/ Str_ , AChm_ , AUp_ , Dn_ /)
+
+   Hash_MCFM_qqVVqq(141,:) = (/ AUp_ , Up_ , AChm_ , Chm_ /)
+   Hash_MCFM_qqVVqq(142,:) = (/ ADn_ , Dn_ , AStr_ , Str_ /)
+   Hash_MCFM_qqVVqq(143,:) = (/ ADn_ , Dn_ , ABot_ , Bot_ /)
+   Hash_MCFM_qqVVqq(144,:) = (/ AStr_ , Str_ , ABot_ , Bot_ /)
+   Hash_MCFM_qqVVqq(145,:) = (/ AUp_ , Up_ , AStr_ , Str_ /)
+   Hash_MCFM_qqVVqq(146,:) = (/ AUp_ , Up_ , ABot_ , Bot_ /)
+   Hash_MCFM_qqVVqq(147,:) = (/ AChm_ , Chm_ , ABot_ , Bot_ /)
+   Hash_MCFM_qqVVqq(148,:) = (/ ADn_ , Dn_ , AChm_ , Chm_ /)
+   Hash_MCFM_qqVVqq(149,:) = (/ ADn_ , Dn_ , AUp_ , Up_ /)
+   Hash_MCFM_qqVVqq(150,:) = (/ AStr_ , Str_ , AChm_ , Chm_ /)
+   Hash_MCFM_qqVVqq(151,:) = (/ AUp_ , Dn_ , AChm_ , Str_ /)
+   Hash_MCFM_qqVVqq(152,:) = (/ ADn_ , Up_ , AStr_ , Chm_ /)
+   Hash_MCFM_qqVVqq(153,:) = (/ AChm_ , Chm_ , AUp_ , Up_ /)
+   Hash_MCFM_qqVVqq(154,:) = (/ AStr_ , Str_ , ADn_ , Dn_ /)
+   Hash_MCFM_qqVVqq(155,:) = (/ ABot_ , Bot_ , ADn_ , Dn_ /)
+   Hash_MCFM_qqVVqq(156,:) = (/ ABot_ , Bot_ , AStr_ , Str_ /)
+   Hash_MCFM_qqVVqq(157,:) = (/ AStr_ , Str_ , AUp_ , Up_ /)
+   Hash_MCFM_qqVVqq(158,:) = (/ ABot_ , Bot_ , AUp_ , Up_ /)
+   Hash_MCFM_qqVVqq(159,:) = (/ ABot_ , Bot_ , AChm_ , Chm_ /)
+   Hash_MCFM_qqVVqq(160,:) = (/ AChm_ , Chm_ , ADn_ , Dn_ /)
+   Hash_MCFM_qqVVqq(161,:) = (/ AUp_ , Up_ , ADn_ , Dn_ /)
+   Hash_MCFM_qqVVqq(162,:) = (/ AChm_ , Chm_ , AStr_ , Str_ /)
+   Hash_MCFM_qqVVqq(163,:) = (/ AStr_ , Chm_ , ADn_ , Up_ /)
+   Hash_MCFM_qqVVqq(164,:) = (/ AChm_ , Str_ , AUp_ , Dn_ /)
 end subroutine
 
 subroutine Init_Hash_MCFM_qqVVqq_Generation()
diff --git a/MELA/fortran/mod_JHUGen.F90 b/MELA/fortran/mod_JHUGen.F90
index d239b04a..a31844ae 100755
--- a/MELA/fortran/mod_JHUGen.F90
+++ b/MELA/fortran/mod_JHUGen.F90
@@ -34,7 +34,7 @@ SUBROUTINE InitFirstTime(pdftable,pdfstrlength,pdfmember)
    lenLHAPDFString = pdfstrlength
 
 !---------------------------
-   call PrintLogo(io_stdout)
+   call PrintLogo(io_stdout, "JHUGen MELA")
 !---------------------------
    call ResetMubarHGabarH()
 !---------------------------
@@ -249,31 +249,4 @@ subroutine ResetPDFs(pdftable,pdfstrlength,pdfmember,pdfid)
 end subroutine
 
 
-SUBROUTINE PrintLogo(TheUnit)
-use modParameters
-use modMisc
-implicit none
-integer :: TheUnit
-integer, parameter :: linelength = 87
-
-    write(TheUnit, *) " "
-    write(TheUnit, *) " ", repeat("*", linelength)
-    write(TheUnit, *) " ", CenterWithStars("JHUGen MELA "//trim(JHUGen_Version), linelength)
-    write(TheUnit, *) " ", repeat("*", linelength)
-    write(TheUnit, *) " ", CenterWithStars("", linelength)
-    write(TheUnit, *) " ", CenterWithStars("Spin and parity determination of single-produced resonances at hadron colliders", linelength)
-    write(TheUnit, *) " ", CenterWithStars("", linelength)
-    write(TheUnit, *) " ", CenterWithStars("I. Anderson, S. Bolognesi, F. Caola, Y. Gao, A. Gritsan, Z. Guo,", linelength)
-    write(TheUnit, *) " ", CenterWithStars("C. Martin, K. Melnikov, R. Rontsch, H. Roskes, U. Sarica, M. Schulze,", linelength)
-    write(TheUnit, *) " ", CenterWithStars("N. Tran, A. Whitbeck, M. Xiao, C. You, Y. Zhou", linelength)
-    write(TheUnit, *) " ", CenterWithStars("Phys.Rev. D81 (2010) 075022;  arXiv:1001.3396  [hep-ph],", linelength)
-    write(TheUnit, *) " ", CenterWithStars("Phys.Rev. D86 (2012) 095031;  arXiv:1208.4018  [hep-ph],", linelength)
-    write(TheUnit, *) " ", CenterWithStars("Phys.Rev. D89 (2014) 035007;  arXiv:1309.4819  [hep-ph],", linelength)
-    write(TheUnit, *) " ", CenterWithStars("Phys.Rev. D94 (2016) 055023;  arXiv:1606.03107 [hep-ph].", linelength)
-    write(TheUnit, *) " ", CenterWithStars("", linelength)
-    write(TheUnit, *) " ", repeat("*", linelength)
-    write(TheUnit, *) " "
-return
-END SUBROUTINE
-
 END MODULE ModJHUGen
diff --git a/MELA/fortran/mod_Misc.F90 b/MELA/fortran/mod_Misc.F90
index d30a69bd..bcfea0da 100755
--- a/MELA/fortran/mod_Misc.F90
+++ b/MELA/fortran/mod_Misc.F90
@@ -1373,24 +1373,72 @@ END SUBROUTINE HTO_Seval3SingleHt
 
 END SUBROUTINE EvaluateSpline
 
-function CenterWithStars(string, totallength)
+function CenterWithStars(string, totallength, align, padleft, padright)
 implicit none
 character(len=*) :: string
 integer :: totallength, nspaces, nleftspaces, nrightspaces
+integer, optional :: align, padleft, padright
+integer :: align2, padleft2, padright2
 character(len=totallength) CenterWithStars
 
-    if (len(trim(string)) .gt. totallength-2) then
-        call Error("len(trim(string)) > totallength-2!")
+    align2 = 2
+    padleft2 = 0
+    padright2 = 0
+    if (present(align)) align2 = align
+    if (present(padleft)) padleft2 = padleft
+    if (present(padright)) padright2 = padright
+
+    if (len(trim(string))+padleft2+padright2 .gt. totallength-2) then
+        call Error("len(trim(string))+padleft+padright > totallength-2!")
     endif
 
     nspaces = totallength - len(trim(string)) - 2
-    nleftspaces = nspaces/2
-    nrightspaces = nspaces-nleftspaces
 
-    CenterWithStars = "*" // repeat(" ", nleftspaces) // string // repeat(" ", nrightspaces) // "*"
+    if (align2.eq.1) then !left
+      nleftspaces = padleft2
+      nrightspaces = nspaces-nleftspaces
+    elseif (align2.eq.2) then !center
+      nleftspaces = (nspaces+padleft2-padright2)/2
+      nrightspaces = nspaces-nleftspaces
+    elseif (align2.eq.3) then !right
+      nrightspaces = padright2
+      nleftspaces = nspaces-nrightspaces
+    else
+      print *, "Unknown align value", align2
+    endif
+
+    CenterWithStars = "*" // repeat(" ", nleftspaces) // trim(string) // repeat(" ", nrightspaces) // "*"
 
 end function
 
+SUBROUTINE PrintLogo(TheUnit, title)
+use modParameters
+implicit none
+integer :: TheUnit
+integer, parameter :: linelength = 87
+character(len=*) :: title
+
+    write(TheUnit, *) " "
+    write(TheUnit, *) " ", repeat("*", linelength)
+    write(TheUnit, *) " ", CenterWithStars(trim(title)//" "//trim(JHUGen_Version), linelength)
+    write(TheUnit, *) " ", repeat("*", linelength)
+    write(TheUnit, *) " ", CenterWithStars("", linelength)
+    write(TheUnit, *) " ", CenterWithStars("Spin and parity determination of single-produced resonances at hadron colliders", linelength)
+    write(TheUnit, *) " ", CenterWithStars("", linelength)
+    write(TheUnit, *) " ", CenterWithStars("I. Anderson, S. Bolognesi, F. Caola, Y. Gao, A. Gritsan,", linelength)
+    write(TheUnit, *) " ", CenterWithStars("Z. Guo, C. Martin, K. Melnikov, R. Rontsch, H. Roskes, U. Sarica,", linelength)
+    write(TheUnit, *) " ", CenterWithStars("M. Schulze, N. Tran, A. Whitbeck, M. Xiao, Y. Zhou", linelength)
+    write(TheUnit, *) " ", CenterWithStars("Phys.Rev. D81 (2010) 075022;  arXiv:1001.3396  [hep-ph],", linelength)
+    write(TheUnit, *) " ", CenterWithStars("Phys.Rev. D86 (2012) 095031;  arXiv:1208.4018  [hep-ph],", linelength)
+    write(TheUnit, *) " ", CenterWithStars("Phys.Rev. D89 (2014) 035007;  arXiv:1309.4819  [hep-ph],", linelength)
+    write(TheUnit, *) " ", CenterWithStars("Phys.Rev. D94 (2016) 055023;  arXiv:1606.03107 [hep-ph].", linelength)
+    write(TheUnit, *) " ", CenterWithStars("", linelength)
+    write(TheUnit, *) " ", repeat("*", linelength)
+    write(TheUnit, *) " "
+return
+END SUBROUTINE
+
+
 real function infinity()
   implicit none
   real :: x
diff --git a/MELA/fortran/mod_Parameters.F90 b/MELA/fortran/mod_Parameters.F90
index eb5a53d8..717d3485 100755
--- a/MELA/fortran/mod_Parameters.F90
+++ b/MELA/fortran/mod_Parameters.F90
@@ -3,7 +3,7 @@ MODULE ModParameters
 save
 !
 !
-character(len=6),parameter :: JHUGen_Version="v7.0.9"
+character(len=*),parameter :: JHUGen_Version="v7.0.11"
 !
 !
 !=====================================================
@@ -125,9 +125,6 @@ MODULE ModParameters
 logical, public :: UseUnformattedRead = .false.  !Set this to true if the regular reading fails for whatever reason
 
 logical, public :: H_DK =.false.                 ! default to false so H in V* > VH (Process = 50) does not decay to bbbar
-
-!leave this one as a parameter, no reason to ever turn it off
-logical, public, parameter :: importExternal_LHEinit = .true.
 !=====================================================
 
 
diff --git a/MELA/interface/Mela.h b/MELA/interface/Mela.h
index d72f20ad..780a735e 100755
--- a/MELA/interface/Mela.h
+++ b/MELA/interface/Mela.h
@@ -110,6 +110,31 @@ class Mela{
     float& costhetastar,
     float& Phi1
     );
+  void computeVBFAngles(
+    float& Q2V1,
+    float& Q2V2,
+    float& costheta1,
+    float& costheta2,
+    float& Phi,
+    float& costhetastar,
+    float& Phi1
+  );
+  void computeVBFAngles_ComplexBoost(
+    float& Q2V1,
+    float& Q2V2,
+    float& costheta1_real, float& costheta1_imag,
+    float& costheta2_real, float& costheta2_imag,
+    float& Phi,
+    float& costhetastar,
+    float& Phi1
+  );
+  void computeVHAngles(
+    float& costheta1,
+    float& costheta2,
+    float& Phi,
+    float& costhetastar,
+    float& Phi1
+  );
 
   void computeP_selfDspin0(
     double selfDHvvcoupl_input[nSupportedHiggses][SIZE_HVV][2],
diff --git a/MELA/interface/TNumericUtil.hh b/MELA/interface/TNumericUtil.hh
index 0dd937a5..3c367860 100755
--- a/MELA/interface/TNumericUtil.hh
+++ b/MELA/interface/TNumericUtil.hh
@@ -21,6 +21,8 @@ namespace TNumericUtil{
     triplet(){}
     T& operator[](std::size_t ipos){ return value[ipos]; } // Return by reference
     const T& operator[](std::size_t ipos)const{ return value[ipos]; } // Return by const reference
+    bool operator == (const triplet<T>& other)const{ bool res = true; for (std::size_t i=0; i<3; i++) res &= (*this)[i]==other[i]; return res; }
+    bool operator != (const triplet<T>& other)const{ return !(*this==other);  }
   };
   template<typename T> struct quadruplet{
     T value[4];
@@ -34,6 +36,8 @@ namespace TNumericUtil{
     quadruplet(){}
     T& operator[](std::size_t ipos){ return value[ipos]; } // Return by reference
     const T& operator[](std::size_t ipos)const{ return value[ipos]; } // Return by const reference
+    bool operator == (const quadruplet<T>& other)const{ bool res = true; for (std::size_t i=0; i<4; i++) res &= (*this)[i]==other[i]; return res; }
+    bool operator != (const quadruplet<T>& other)const{ return !(*this==other); }
   };
 
   typedef triplet<int> intTriplet_t;
diff --git a/MELA/interface/TUtil.hh b/MELA/interface/TUtil.hh
index 535dd8c2..5bf596bb 100755
--- a/MELA/interface/TUtil.hh
+++ b/MELA/interface/TUtil.hh
@@ -80,30 +80,30 @@ namespace TUtil{
   /// Leptons are re-ordered internally according to a standard convention:
   /// lept1 = negative-charged lepton (for OS pairs).
   void computeAngles(
-    TLorentzVector Z1_lept1, int Z1_lept1Id,
-    TLorentzVector Z1_lept2, int Z1_lept2Id,
-    TLorentzVector Z2_lept1, int Z2_lept1Id,
-    TLorentzVector Z2_lept2, int Z2_lept2Id,
     float& costhetastar,
     float& costheta1,
     float& costheta2,
     float& Phi,
-    float& Phi1
-    );
-  void computeAnglesCS(
+    float& Phi1,
     TLorentzVector Z1_lept1, int Z1_lept1Id,
     TLorentzVector Z1_lept2, int Z1_lept2Id,
     TLorentzVector Z2_lept1, int Z2_lept1Id,
-    TLorentzVector Z2_lept2, int Z2_lept2Id,
+    TLorentzVector Z2_lept2, int Z2_lept2Id
+    );
+  void computeAnglesCS(
     float pbeam,
     float& costhetastar,
     float& costheta1,
     float& costheta2,
     float& Phi,
-    float& Phi1
+    float& Phi1,
+    TLorentzVector Z1_lept1, int Z1_lept1Id,
+    TLorentzVector Z1_lept2, int Z1_lept2Id,
+    TLorentzVector Z2_lept1, int Z2_lept1Id,
+    TLorentzVector Z2_lept2, int Z2_lept2Id
     );
   // Angles of associated production
-  void computeVBFangles(
+  void computeVBFAngles(
     float& costhetastar,
     float& costheta1,
     float& costheta2,
@@ -120,7 +120,7 @@ namespace TUtil{
     TLorentzVector* injet1=0, int injet1Id=0, // Gen. partons in lab frame
     TLorentzVector* injet2=0, int injet2Id=0
     );
-  void computeVBFangles_ComplexBoost(
+  void computeVBFAngles_ComplexBoost(
     float& costhetastar,
     float& costheta1_real, float& costheta1_imag,
     float& costheta2_real, float& costheta2_imag,
@@ -137,7 +137,7 @@ namespace TUtil{
     TLorentzVector* injet1=0, int injet1Id=0, // Gen. partons in lab frame
     TLorentzVector* injet2=0, int injet2Id=0
     );
-  void computeVHangles(
+  void computeVHAngles(
     float& costhetastar,
     float& costheta1,
     float& costheta2,
diff --git a/MELA/python/mela.py b/MELA/python/mela.py
index 309b8796..3697b07c 100644
--- a/MELA/python/mela.py
+++ b/MELA/python/mela.py
@@ -63,6 +63,46 @@ class Mela(object):
       );
       return result;
     }
+    vector<float> computeVBFAngles(Mela& mela) {
+      vector<float> result(7);
+      mela.computeVBFAngles(
+        result[0],
+        result[1],
+        result[2],
+        result[3],
+        result[4],
+        result[5],
+        result[6]
+      );
+      return result;
+    }
+    vector<float> computeVBFAngles_ComplexBoost(Mela& mela) {
+      vector<float> result(9);
+      mela.computeVBFAngles_ComplexBoost(
+        result[0],
+        result[1],
+        result[2],
+        result[3],
+        result[4],
+        result[5],
+        result[6],
+        result[7],
+        result[8]
+      );
+      return result;
+    }
+    vector<float> computeVHAngles(Mela& mela, TVar::Production prod) {
+      mela.setProcess(TVar::HSMHiggs, TVar::JHUGen, prod);
+      vector<float> result(5);
+      mela.computeVHAngles(
+        result[0],
+        result[1],
+        result[2],
+        result[3],
+        result[4]
+      );
+      return result;
+    }
     //not implementing the computeP_selfD* functions here
     //would be easier to do in pure python but not worth it anyway
     float computeP(Mela& mela, bool useConstant) {
@@ -232,8 +272,17 @@ def setInputEvent_fromLHE(self, event, isgen):
     self.setInputEvent(SimpleParticleCollection_t(daughters), SimpleParticleCollection_t(associated), SimpleParticleCollection_t(mothers), isgen)
 
   def getPAux(self): return ROOT.getPAux(self.__mela)
+
   DecayAngles = namedtuple("DecayAngles", "qH m1 m2 costheta1 costheta2 Phi costhetastar Phi1")
   def computeDecayAngles(self): return self.DecayAngles(*ROOT.computeDecayAngles(self.__mela))
+  VBFAngles = namedtuple("VBFAngles", "Q2V1 Q2V2 costheta1 costheta2 Phi costhetastar Phi1")
+  def computeVBFAngles(self): return self.VBFAngles(*ROOT.computeVBFAngles(self.__mela))
+  def computeVBFAngles_ComplexBoost(self):
+    result = ROOT.computeVBFAngles_ComplexBoost(self.__mela)
+    return self.VBFAngles(result[0], result[1], result[2]+1j*result[3], result[4]+1j*result[5], *result[6:])
+  VHAngles = namedtuple("VHAngles", "costheta1 costheta2 Phi costhetastar Phi1")
+  def computeVHAngles(self, prod): return self.VHAngles(*ROOT.computeVHAngles(self.__mela, prod))
+
   def computeP(self, useConstant=True): return ROOT.computeP(self.__mela, useConstant)
   def computeD_CP(self, myME, myType): return ROOT.computeD_CP(self.__mela, myME, myType)
   def computeProdP(self, useConstant=True): return ROOT.computeProdP(self.__mela, useConstant)
@@ -537,8 +586,8 @@ def SimpleParticle_t(lineorid, pxortlv=None, py=None, pz=None, e=None):
   -11 -93.72924763700000028 39.45060783929999815 -92.98363978320000456 137.79506373300000632
   """
   associated = """
-  -11 211.33318543799998679 -14.90577872979999974 3.74371777679000006 211.89127619999999297
-  12 31.22409920730000010 -37.83127789369999761 1.23465418111000003 49.06805813689999951
+  -1 211.33318543799998679 -14.90577872979999974 3.74371777679000006 211.89127619999999297
+  2 31.22409920730000010 -37.83127789369999761 1.23465418111000003 49.06805813689999951
   """
   mothers = """
   -1 0.00000000000000000 0.00000000000000000 192.71975508899998886 192.71975508899998886
@@ -565,7 +614,7 @@ def SimpleParticle_t(lineorid, pxortlv=None, py=None, pz=None, e=None):
               )
   for _ in couplings:
     m.ghz1, m.ghz2, m.ghz4, m.ghz1_prime2 = _
-    m.setProcess(TVar.SelfDefine_spin0, TVar.JHUGen, TVar.Lep_WH)
+    m.setProcess(TVar.SelfDefine_spin0, TVar.JHUGen, TVar.Had_WH)
     prod = m.computeProdP(False)
     m.ghz1, m.ghz2, m.ghz4, m.ghz1_prime2 = _
     m.setProcess(TVar.SelfDefine_spin0, TVar.JHUGen, TVar.ZZINDEPENDENT)
@@ -573,6 +622,9 @@ def SimpleParticle_t(lineorid, pxortlv=None, py=None, pz=None, e=None):
     print prod, dec, prod*dec
 
   print m.computeDecayAngles()
+  print m.computeVBFAngles()
+  print m.computeVBFAngles_ComplexBoost()
+  print m.computeVHAngles(TVar.Had_WH)
   print "propagator:"
   print "   BW:", m.getXPropagator(TVar.FixedWidth)
   print "  CPS:", m.getXPropagator(TVar.CPS)
diff --git a/MELA/src/LinkDef.h b/MELA/src/LinkDef.h
index a677340e..6c998409 100755
--- a/MELA/src/LinkDef.h
+++ b/MELA/src/LinkDef.h
@@ -35,8 +35,9 @@
 #pragma link C++ namespace TUtil;
 #pragma link C++ function TUtil::computeAngles;
 #pragma link C++ function TUtil::computeAnglesCS;
-#pragma link C++ function TUtil::computeVBFangles;
-#pragma link C++ function TUtil::computeVHangles;
+#pragma link C++ function TUtil::computeVBFAngles;
+#pragma link C++ function TUtil::computeVBFAngles_ComplexBoost;
+#pragma link C++ function TUtil::computeVHAngles;
 #pragma link C++ function TUtil::scaleMomentumToEnergy;
 #pragma link C++ function TUtil::constrainedRemovePairMass;
 #pragma link C++ function TUtil::removeMassFromPair;
diff --git a/MELA/src/Mela.cc b/MELA/src/Mela.cc
index f9c32826..90d000c9 100755
--- a/MELA/src/Mela.cc
+++ b/MELA/src/Mela.cc
@@ -222,8 +222,34 @@ void Mela::build(double mh_){
 }
 
 void Mela::printLogo() const{
-  // To be written
-  // e.g. MELAout << '*'; MELAout.writeCentered("Hello world!", ' ', 51); MELAout << '*' << endl;
+  vector<string> logolines;
+  logolines.push_back("MELA (Matrix Element Likelihood Approach)");
+  logolines.push_back("");
+  logolines.push_back("Data analysis and Monte Carlo weights package");
+  logolines.push_back("for analyses of resonances produced at pp, ppbar, and e+e- colliders, featuring:");
+  logolines.push_back("");
+  logolines.push_back("* JHUGenMELA *");
+  logolines.push_back("Signal calculations based on analytical pdf.s, and JHU Generator (JHUGen) matrix elements");
+  logolines.push_back("(See JHUGen credits below)");
+  logolines.push_back("");
+  logolines.push_back("* MCFM *");
+  logolines.push_back("Signal, background, and interference calculations, modified based on JHUGen matrix elements");
+  logolines.push_back("(See MCFM credits below)");
+  logolines.push_back("");
+  logolines.push_back("For more details: http://spin.pha.jhu.edu");
+  logolines.push_back("");
+  size_t maxlinesize = 0;
+  for (auto const& l:logolines) maxlinesize = std::max(maxlinesize, l.length());
+  MELAout.writeCentered("", '*', maxlinesize+10); MELAout << endl;
+  unsigned int iline=0;
+  for (auto const& l:logolines){
+    MELAout << '*';
+    MELAout.writeCentered(l, ' ', maxlinesize+8);
+    MELAout << '*' << endl;
+    if (iline==0){ MELAout.writeCentered("", '*', maxlinesize+10); MELAout << endl; }
+    iline++;
+  }
+  MELAout.writeCentered("", '*', maxlinesize+10); MELAout << endl;
 }
 
 // Set-functions
@@ -406,31 +432,44 @@ void Mela::computeDecayAngles(
   float& Phi,
   float& costhetastar,
   float& Phi1
-  ){
+){
+  using TVar::simple_event_record;
+
   qH=0; m1=0; m2=0; costheta1=0; costheta2=0; Phi=0; costhetastar=0; Phi1=0;
 
   if (melaCand==0) melaCand = getCurrentCandidate();
   if (melaCand!=0){
     TLorentzVector nullVector(0, 0, 0, 0);
 
-    qH = melaCand->m();
-    m1 = melaCand->getSortedV(0)->m();
-    m2 = melaCand->getSortedV(1)->m();
+    int partIncCode=TVar::kNoAssociated; // Only use associated partons in the pT=0 frame boost
+    simple_event_record mela_event;
+    mela_event.AssociationCode=partIncCode;
+    TUtil::GetBoostedParticleVectors(melaCand, mela_event, myVerbosity_);
+    SimpleParticleCollection_t& daughters = mela_event.pDaughters;
 
-    if (melaCand->getSortedV(0)->getNDaughters()>=1 && melaCand->getSortedV(1)->getNDaughters()>=1){
-      MELAParticle* dau[2][2]={ { 0 } };
-      for (int vv=0; vv<2; vv++){
-        MELAParticle* Vi = melaCand->getSortedV(vv);
-        for (int dd=0; dd<Vi->getNDaughters(); dd++) dau[vv][dd] = Vi->getDaughter(dd);
-      }
-      TUtil::computeAngles(
-        (dau[0][0]!=0 ? dau[0][0]->p4 : nullVector), (dau[0][0]!=0 ? dau[0][0]->id : -9000),
-        (dau[0][1]!=0 ? dau[0][1]->p4 : nullVector), (dau[0][1]!=0 ? dau[0][1]->id : -9000),
-        (dau[1][0]!=0 ? dau[1][0]->p4 : nullVector), (dau[1][0]!=0 ? dau[1][0]->id : -9000),
-        (dau[1][1]!=0 ? dau[1][1]->p4 : nullVector), (dau[1][1]!=0 ? dau[1][1]->id : -9000),
-        costhetastar, costheta1, costheta2, Phi, Phi1
-        );
+    if (daughters.size()<2 || daughters.size()>4 || mela_event.intermediateVid.size()!=2){
+      if (myVerbosity_>=TVar::ERROR) MELAerr << "Mela::computeDecayAngles: Number of daughters " << daughters.size() << " or number of intermediate Vs " << mela_event.intermediateVid << " not supported!" << endl;
+      return;
     }
+
+    // Make sure there are exactly 4 daughters, null or not
+    if (daughters.size()%2==1){ for (unsigned int ipar=daughters.size(); ipar<4; ipar++) daughters.push_back(SimpleParticle_t(-9000, nullVector)); }
+    else if (daughters.size()==2){
+      daughters.push_back(SimpleParticle_t(-9000, nullVector));
+      daughters.insert(daughters.begin()+1, SimpleParticle_t(-9000, nullVector));
+    }
+    qH = (daughters.at(0).second+daughters.at(1).second+daughters.at(2).second+daughters.at(3).second).M();
+    m1 = (daughters.at(0).second+daughters.at(1).second).M();
+    m2 = (daughters.at(2).second+daughters.at(3).second).M();
+
+    TUtil::computeAngles(
+      costhetastar, costheta1, costheta2, Phi, Phi1,
+      daughters.at(0).second, daughters.at(0).first,
+      daughters.at(1).second, daughters.at(1).first,
+      daughters.at(2).second, daughters.at(2).first,
+      daughters.at(3).second, daughters.at(3).first
+    );
+
     // Protect against NaN
     if (!(costhetastar==costhetastar)) costhetastar=0;
     if (!(costheta1==costheta1)) costheta1=0;
@@ -441,6 +480,233 @@ void Mela::computeDecayAngles(
   else if (myVerbosity_>=TVar::DEBUG) MELAerr << "Mela::computeDecayAngles: No possible melaCand in TEvtProb to compute angles." << endl;
 }
 
+// VBF angles computation script of Mela to convert MELACandidates to m1, m2 etc.
+void Mela::computeVBFAngles(
+  float& Q2V1,
+  float& Q2V2,
+  float& costheta1,
+  float& costheta2,
+  float& Phi,
+  float& costhetastar,
+  float& Phi1
+){
+  using TVar::simple_event_record;
+
+  Q2V1=0; Q2V2=0; costheta1=0; costheta2=0; Phi=0; costhetastar=0; Phi1=0;
+
+  if (melaCand==0) melaCand = getCurrentCandidate();
+  if (melaCand!=0){
+    TLorentzVector nullVector(0, 0, 0, 0);
+
+    int nRequested_AssociatedJets=2;
+    int partIncCode=TVar::kUseAssociated_Jets; // Only use associated partons in the pT=0 frame boost
+    simple_event_record mela_event;
+    mela_event.AssociationCode=partIncCode;
+    mela_event.nRequested_AssociatedJets=nRequested_AssociatedJets;
+    TUtil::GetBoostedParticleVectors(melaCand, mela_event, myVerbosity_);
+    SimpleParticleCollection_t& mothers = mela_event.pMothers;
+    SimpleParticleCollection_t& aparts = mela_event.pAssociated;
+    SimpleParticleCollection_t& daughters = mela_event.pDaughters;
+
+    if ((int)aparts.size()!=nRequested_AssociatedJets){ if (myVerbosity_>=TVar::ERROR) MELAerr << "Mela::computeVBFAngles: Number of associated particles is not 2!" << endl; return; }
+
+    // Make sure there are exactly 4 daughters, null or not
+    if (daughters.size()>4){ // Unsupported size, default to undecayed Higgs
+      SimpleParticle_t& firstPart = daughters.at(0);
+      firstPart.first=25;
+      for (auto it=daughters.cbegin()+1; it!=daughters.cend(); it++){ firstPart.second = firstPart.second + it->second; }
+      daughters.erase(daughters.begin()+4, daughters.end());
+    }
+    if (daughters.size()%2==1){ for (unsigned int ipar=daughters.size(); ipar<4; ipar++) daughters.push_back(SimpleParticle_t(-9000, nullVector)); }
+    else if (daughters.size()==2){
+      daughters.push_back(SimpleParticle_t(-9000, nullVector));
+      daughters.insert(daughters.begin()+1, SimpleParticle_t(-9000, nullVector));
+    }
+
+    TUtil::computeVBFAngles(
+      costhetastar, costheta1, costheta2, Phi, Phi1, Q2V1, Q2V2,
+      daughters.at(0).second, daughters.at(0).first,
+      daughters.at(1).second, daughters.at(1).first,
+      daughters.at(2).second, daughters.at(2).first,
+      daughters.at(3).second, daughters.at(3).first,
+      aparts.at(0).second, aparts.at(0).first,
+      aparts.at(1).second, aparts.at(1).first,
+      &(mothers.at(0).second), mothers.at(0).first,
+      &(mothers.at(1).second), mothers.at(1).first
+    );
+
+    // Protect against NaN
+    if (!(costhetastar==costhetastar)) costhetastar=0;
+    if (!(costheta1==costheta1)) costheta1=0;
+    if (!(costheta2==costheta2)) costheta2=0;
+    if (!(Phi==Phi)) Phi=0;
+    if (!(Phi1==Phi1)) Phi1=0;
+  }
+  else if (myVerbosity_>=TVar::DEBUG) MELAerr << "Mela::computeVBFAngles: No possible melaCand in TEvtProb to compute angles." << endl;
+}
+void Mela::computeVBFAngles_ComplexBoost(
+  float& Q2V1,
+  float& Q2V2,
+  float& costheta1_real, float& costheta1_imag,
+  float& costheta2_real, float& costheta2_imag,
+  float& Phi,
+  float& costhetastar,
+  float& Phi1
+){
+  using TVar::simple_event_record;
+
+  Q2V1=0; Q2V2=0; costheta1_real=0; costheta2_real=0; costheta1_imag=0; costheta2_imag=0; Phi=0; costhetastar=0; Phi1=0;
+
+  if (melaCand==0) melaCand = getCurrentCandidate();
+  if (melaCand!=0){
+    TLorentzVector nullVector(0, 0, 0, 0);
+
+    int nRequested_AssociatedJets=2;
+    int partIncCode=TVar::kUseAssociated_Jets; // Only use associated partons in the pT=0 frame boost
+    simple_event_record mela_event;
+    mela_event.AssociationCode=partIncCode;
+    mela_event.nRequested_AssociatedJets=nRequested_AssociatedJets;
+    TUtil::GetBoostedParticleVectors(melaCand, mela_event, myVerbosity_);
+    SimpleParticleCollection_t& mothers = mela_event.pMothers;
+    SimpleParticleCollection_t& aparts = mela_event.pAssociated;
+    SimpleParticleCollection_t& daughters = mela_event.pDaughters;
+
+    if ((int) aparts.size()!=nRequested_AssociatedJets){ if (myVerbosity_>=TVar::ERROR) MELAerr << "Mela::computeVBFAngles_ComplexBoost: Number of associated particles is not 2!" << endl; return; }
+
+    // Make sure there are exactly 4 daughters, null or not
+    if (daughters.size()>4){ // Unsupported size, default to undecayed Higgs
+      SimpleParticle_t& firstPart = daughters.at(0);
+      firstPart.first=25;
+      for (auto it=daughters.cbegin()+1; it!=daughters.cend(); it++){ firstPart.second = firstPart.second + it->second; }
+      daughters.erase(daughters.begin()+4, daughters.end());
+    }
+    if (daughters.size()%2==1){ for (unsigned int ipar=daughters.size(); ipar<4; ipar++) daughters.push_back(SimpleParticle_t(-9000, nullVector)); }
+    else if (daughters.size()==2){
+      daughters.push_back(SimpleParticle_t(-9000, nullVector));
+      daughters.insert(daughters.begin()+1, SimpleParticle_t(-9000, nullVector));
+    }
+
+    TUtil::computeVBFAngles_ComplexBoost(
+      costhetastar, costheta1_real, costheta1_imag, costheta2_real, costheta2_imag, Phi, Phi1, Q2V1, Q2V2,
+      daughters.at(0).second, daughters.at(0).first,
+      daughters.at(1).second, daughters.at(1).first,
+      daughters.at(2).second, daughters.at(2).first,
+      daughters.at(3).second, daughters.at(3).first,
+      aparts.at(0).second, aparts.at(0).first,
+      aparts.at(1).second, aparts.at(1).first,
+      &(mothers.at(0).second), mothers.at(0).first,
+      &(mothers.at(1).second), mothers.at(1).first
+    );
+
+    // Protect against NaN
+    if (!(costhetastar==costhetastar)) costhetastar=0;
+    if (!(costheta1_real==costheta1_real)) costheta1_real=0;
+    if (!(costheta2_real==costheta2_real)) costheta2_real=0;
+    if (!(costheta1_imag==costheta1_imag)) costheta1_imag=0;
+    if (!(costheta2_imag==costheta2_imag)) costheta2_imag=0;
+    if (!(Phi==Phi)) Phi=0;
+    if (!(Phi1==Phi1)) Phi1=0;
+  }
+  else if (myVerbosity_>=TVar::DEBUG) MELAerr << "Mela::computeVBFAngles_ComplexBoost: No possible melaCand in TEvtProb to compute angles." << endl;
+}
+
+// VH angles computation script of Mela to convert MELACandidates to production angles.
+void Mela::computeVHAngles(
+  float& costheta1,
+  float& costheta2,
+  float& Phi,
+  float& costhetastar,
+  float& Phi1
+){
+  using TVar::simple_event_record;
+
+  costheta1=0; costheta2=0; Phi=0; costhetastar=0; Phi1=0;
+
+  if (melaCand==0) melaCand = getCurrentCandidate();
+  if (melaCand!=0){
+    TLorentzVector nullVector(0, 0, 0, 0);
+
+    if (!(myProduction_ == TVar::Lep_ZH || myProduction_ == TVar::Lep_WH || myProduction_ == TVar::Had_ZH || myProduction_ == TVar::Had_WH || myProduction_ == TVar::GammaH)){
+      if (myVerbosity_>=TVar::ERROR) MELAerr << "Mela::computeVHAngles: Production is not supported!" << endl;
+      return;
+    }
+
+    int nRequested_AssociatedJets=0;
+    int nRequested_AssociatedLeptons=0;
+    int nRequested_AssociatedPhotons=0;
+    int AssociationVCompatibility=0;
+    int partIncCode=TVar::kNoAssociated; // Just to avoid warnings
+    if (myProduction_ == TVar::Had_ZH || myProduction_ == TVar::Had_WH){ // Only use associated partons
+      partIncCode=TVar::kUseAssociated_Jets;
+      nRequested_AssociatedJets=2;
+    }
+    else if (myProduction_ == TVar::Lep_ZH || myProduction_ == TVar::Lep_WH){ // Only use associated leptons(+)neutrinos
+      partIncCode=TVar::kUseAssociated_Leptons;
+      nRequested_AssociatedLeptons=2;
+    }
+    else if (myProduction_ == TVar::GammaH){ // Only use associated photon
+      partIncCode=TVar::kUseAssociated_Photons;
+      nRequested_AssociatedPhotons=1;
+    }
+    if (myProduction_==TVar::Lep_WH || myProduction_==TVar::Had_WH) AssociationVCompatibility=24;
+    else if (myProduction_==TVar::Lep_ZH || myProduction_==TVar::Had_ZH) AssociationVCompatibility=23;
+    else if (myProduction_==TVar::GammaH) AssociationVCompatibility=22;
+    simple_event_record mela_event;
+    mela_event.AssociationCode=partIncCode;
+    mela_event.AssociationVCompatibility=AssociationVCompatibility;
+    mela_event.nRequested_AssociatedJets=nRequested_AssociatedJets;
+    mela_event.nRequested_AssociatedLeptons=nRequested_AssociatedLeptons;
+    mela_event.nRequested_AssociatedPhotons=nRequested_AssociatedPhotons;
+    TUtil::GetBoostedParticleVectors(melaCand, mela_event, myVerbosity_);
+    SimpleParticleCollection_t& mothers = mela_event.pMothers;
+    SimpleParticleCollection_t& aparts = mela_event.pAssociated;
+    SimpleParticleCollection_t& daughters = mela_event.pDaughters;
+
+    if ((aparts.size()<(unsigned int) (nRequested_AssociatedJets+nRequested_AssociatedLeptons) && myProduction_!=TVar::GammaH) || (aparts.size()<(unsigned int) nRequested_AssociatedPhotons && myProduction_==TVar::GammaH)){
+      if (myVerbosity_>=TVar::ERROR){
+        MELAerr << "Mela::computeVHAngles: Number of associated particles (" << aparts.size() << ") is less than ";
+        if (myProduction_!=TVar::GammaH) MELAerr << (nRequested_AssociatedJets+nRequested_AssociatedLeptons);
+        else MELAerr << nRequested_AssociatedPhotons;
+        MELAerr << endl;
+      }
+      return;
+    }
+
+    // Make sure there are exactly 4 daughters, null or not
+    if (daughters.size()>4){ // Unsupported size, default to undecayed Higgs
+      SimpleParticle_t& firstPart = daughters.at(0);
+      firstPart.first=25;
+      for (auto it=daughters.cbegin()+1; it!=daughters.cend(); it++){ firstPart.second = firstPart.second + it->second; }
+      daughters.erase(daughters.begin()+4, daughters.end());
+    }
+    if (daughters.size()%2==1){ for (unsigned int ipar=daughters.size(); ipar<4; ipar++) daughters.push_back(SimpleParticle_t(-9000, nullVector)); }
+    else if (daughters.size()==2){
+      daughters.push_back(SimpleParticle_t(-9000, nullVector));
+      daughters.insert(daughters.begin()+1, SimpleParticle_t(-9000, nullVector));
+    }
+
+    TUtil::computeVHAngles(
+      costhetastar, costheta1, costheta2, Phi, Phi1,
+      daughters.at(0).second, daughters.at(0).first,
+      daughters.at(1).second, daughters.at(1).first,
+      daughters.at(2).second, daughters.at(2).first,
+      daughters.at(3).second, daughters.at(3).first,
+      aparts.at(0).second, aparts.at(0).first,
+      aparts.at(1).second, aparts.at(1).first,
+      &(mothers.at(0).second), mothers.at(0).first,
+      &(mothers.at(1).second), mothers.at(1).first
+    );
+
+    // Protect against NaN
+    if (!(costhetastar==costhetastar)) costhetastar=0;
+    if (!(costheta1==costheta1)) costheta1=0;
+    if (!(costheta2==costheta2)) costheta2=0;
+    if (!(Phi==Phi)) Phi=0;
+    if (!(Phi1==Phi1)) Phi1=0;
+  }
+  else if (myVerbosity_>=TVar::DEBUG) MELAerr << "Mela::computeVHAngles: No possible melaCand in TEvtProb to compute angles." << endl;
+}
+
 // Regular probabilities
 void Mela::computeP_selfDspin0(
   double selfDHvvcoupl_input[nSupportedHiggses][SIZE_HVV][2],
diff --git a/MELA/src/TMCFMUtils.cc b/MELA/src/TMCFMUtils.cc
index 422d3cd8..5d52ab21 100755
--- a/MELA/src/TMCFMUtils.cc
+++ b/MELA/src/TMCFMUtils.cc
@@ -163,7 +163,9 @@ std::vector<intQuad_t> TMCFMUtils::Hash_QQVVQQ(){
         if ((idpos<2 && ipos>=2) || (idpos>=2 && ipos<2)) idAssigned = -idAssigned;
         cfg[jcfg.at(j)[ipos]] = idAssigned;
       }
-      pcfg.push_back(cfg);
+      bool exists=false;
+      for (auto const& tmpcfg:pcfg){ if (tmpcfg==cfg){ exists=true; break; } }
+      if (!exists) pcfg.push_back(cfg);
     }
   }
   // Uncommenting the lines below prints out the hash when the library is loaded.
@@ -216,7 +218,9 @@ std::vector<intQuad_t> TMCFMUtils::Hash_QQVVQQStrong(){
         if (((idpos<2 && ipos>=2) || (idpos>=2 && ipos<2)) && !PDGHelpers::isAGluon(idAssigned)) idAssigned = -idAssigned;
         cfg[jcfg.at(j)[ipos]] = idAssigned;
       }
-      pcfg.push_back(cfg);
+      bool exists=false;
+      for (auto const& tmpcfg:pcfg){ if (tmpcfg==cfg){ exists=true; break; } }
+      if (!exists) pcfg.push_back(cfg);
     }
   }
   // Uncommenting the lines below prints out the hash when the library is loaded.
diff --git a/MELA/src/TUtil.cc b/MELA/src/TUtil.cc
index 35514fa6..443df042 100755
--- a/MELA/src/TUtil.cc
+++ b/MELA/src/TUtil.cc
@@ -201,15 +201,15 @@ std::pair<TLorentzVector, TLorentzVector> TUtil::ComplexBoost(TVector3 beta, TLo
 
 /***** Decay angles *****/
 void TUtil::computeAngles(
-  TLorentzVector p4M11, int Z1_lept1Id,
-  TLorentzVector p4M12, int Z1_lept2Id,
-  TLorentzVector p4M21, int Z2_lept1Id,
-  TLorentzVector p4M22, int Z2_lept2Id,
   float& costhetastar,
   float& costheta1,
   float& costheta2,
   float& Phi,
-  float& Phi1
+  float& Phi1,
+  TLorentzVector p4M11, int Z1_lept1Id,
+  TLorentzVector p4M12, int Z1_lept2Id,
+  TLorentzVector p4M21, int Z2_lept1Id,
+  TLorentzVector p4M22, int Z2_lept2Id
   ){
   TLorentzVector nullFourVector(0, 0, 0, 0);
   if (p4M12==nullFourVector || p4M22==nullFourVector){
@@ -237,7 +237,7 @@ void TUtil::computeAngles(
 
   // Sort Z1 leptons so that:
   if (
-    Z1_lept2Id!=9000
+    Z1_lept2Id!=-9000
     &&
     (
     (Z1_lept1Id*Z1_lept2Id<0 && Z1_lept1Id<0) // for OS pairs: lep1 must be the negative one
@@ -247,7 +247,7 @@ void TUtil::computeAngles(
     ) swap(p4M11, p4M12);
   // Same for Z2 leptons
   if (
-    Z2_lept2Id!=9000
+    Z2_lept2Id!=-9000
     &&
     (
     (Z2_lept1Id*Z2_lept2Id<0 && Z2_lept1Id<0)
@@ -367,16 +367,16 @@ void TUtil::computeAngles(
   }
 }
 void TUtil::computeAnglesCS(
-  TLorentzVector p4M11, int Z1_lept1Id,
-  TLorentzVector p4M12, int Z1_lept2Id,
-  TLorentzVector p4M21, int Z2_lept1Id,
-  TLorentzVector p4M22, int Z2_lept2Id,
   float pbeam,
   float& costhetastar,
   float& costheta1,
   float& costheta2,
   float& Phi,
-  float& Phi1
+  float& Phi1,
+  TLorentzVector p4M11, int Z1_lept1Id,
+  TLorentzVector p4M12, int Z1_lept2Id,
+  TLorentzVector p4M21, int Z2_lept1Id,
+  TLorentzVector p4M22, int Z2_lept2Id
   ){
   TLorentzVector nullFourVector(0, 0, 0, 0);
   if (p4M12==nullFourVector || p4M22==nullFourVector){
@@ -588,7 +588,7 @@ void TUtil::computeAnglesCS(
 }
 
 /***** Associated production angles *****/
-void TUtil::computeVBFangles(
+void TUtil::computeVBFAngles(
   float& costhetastar,
   float& costheta1,
   float& costheta2,
@@ -740,7 +740,7 @@ void TUtil::computeVBFangles(
   costheta1 = V1.Vect().Unit().Dot(fermion1.Vect().Unit());
   costheta2 = V2.Vect().Unit().Dot(fermion2.Vect().Unit());
 }
-void TUtil::computeVBFangles_ComplexBoost(
+void TUtil::computeVBFAngles_ComplexBoost(
   float& costhetastar,
   float& costheta1_real, float& costheta1_imag,
   float& costheta2_real, float& costheta2_imag,
@@ -888,7 +888,7 @@ void TUtil::computeVBFangles_ComplexBoost(
   Q2V1 = -(V1.M2());
   Q2V2 = -(V2.M2());
 
-  // Up to here, everything has to be the same as TUtil::computeVBFangles
+  // Up to here, everything has to be the same as TUtil::computeVBFAngles
   // Computations that would have been truly in the frame of X had V1 and V2 not been virtual:
   // Use TUtil::ComplexBoost to evade imaginary gamma problems when beta**2<0
   pair<TLorentzVector, TLorentzVector> V2_BV1 = TUtil::ComplexBoost(V1.BoostVector(), V2);
@@ -901,7 +901,7 @@ void TUtil::computeVBFangles_ComplexBoost(
   costheta2_real = -(V1_BV2.first.Vect().Unit().Dot(fermion2_BV2.first.Vect().Unit()) - V1_BV2.second.Vect().Unit().Dot(fermion2_BV2.second.Vect().Unit()));
   costheta2_imag = -(V1_BV2.first.Vect().Unit().Dot(fermion2_BV2.second.Vect().Unit()) + V1_BV2.second.Vect().Unit().Dot(fermion2_BV2.first.Vect().Unit()));
 }
-void TUtil::computeVHangles(
+void TUtil::computeVHAngles(
   float& costhetastar,
   float& costheta1,
   float& costheta2,
@@ -1027,16 +1027,16 @@ void TUtil::computeVHangles(
   }
 
   TUtil::computeAngles(
-    -P1, 23, // Id is 23 to avoid an attempt to remove quark mass
-    -P2, 0, // Id is 0 to avoid swapping
-    jet1massless, 23,
-    jet2massless, 0,
     costhetastar,
     costheta1,
     costheta2,
     Phi,
-    Phi1
-    );
+    Phi1,
+    -P1, 23, // Id is 23 to avoid an attempt to remove quark mass
+    -P2, 0, // Id is 0 to avoid swapping
+    jet1massless, 23,
+    jet2massless, 0
+  );
 }
 
 
@@ -4196,6 +4196,7 @@ double TUtil::SumMatrixElementPDF(
   return msqjk;
 }
 
+// ggHdec+0J or Hdec+0J
 double TUtil::JHUGenMatEl(
   const TVar::Process& process, const TVar::Production& production, const TVar::MatrixElement& matrixElement,
   event_scales_type* event_scales, MelaIO* RcdME,
@@ -4267,7 +4268,7 @@ double TUtil::JHUGenMatEl(
     verbosity
     );
   if (mela_event.pDaughters.size()<2 || mela_event.intermediateVid.size()!=2){
-    if (verbosity>=TVar::ERROR) MELAerr << "TUtil::JHUGenMatEl: Number of daughters " << mela_event.pDaughters.size() << " or number of intermediate Vs " << mela_event.intermediateVid.size() << " not supported!" << endl;
+    if (verbosity>=TVar::ERROR) MELAerr << "TUtil::JHUGenMatEl: Number of daughters " << mela_event.pDaughters.size() << " or number of intermediate Vs " << mela_event.intermediateVid << " not supported!" << endl;
     return MatElSq;
   }
 
@@ -4568,6 +4569,7 @@ double TUtil::JHUGenMatEl(
   return MatElSq;
 }
 
+// VBF and HJJ (QCD), H undecayed, includes H BW
 double TUtil::HJJMatEl(
   const TVar::Process& process, const TVar::Production& production, const TVar::MatrixElement& matrixElement,
   event_scales_type* event_scales, MelaIO* RcdME,
@@ -6011,6 +6013,7 @@ double TUtil::HJJMatEl(
   return sum_msqjk;
 }
 
+// VH, H undecayed or Hbb
 double TUtil::VHiggsMatEl(
   const TVar::Process& process, const TVar::Production& production, const TVar::MatrixElement& matrixElement,
   event_scales_type* event_scales, MelaIO* RcdME,
@@ -6619,8 +6622,7 @@ double TUtil::VHiggsMatEl(
   return sum_msqjk;
 }
 
-
-// ttH
+// ttH, H undecayed
 double TUtil::TTHiggsMatEl(
   const TVar::Process& process, const TVar::Production& production, const TVar::MatrixElement& matrixElement,
   event_scales_type* event_scales, MelaIO* RcdME,
@@ -6995,7 +6997,7 @@ double TUtil::TTHiggsMatEl(
   return sum_msqjk;
 }
 
-// bbH
+// bbH, H undecayed
 double TUtil::BBHiggsMatEl(
   const TVar::Process& process, const TVar::Production& production, const TVar::MatrixElement& matrixElement,
   event_scales_type* event_scales, MelaIO* RcdME,
diff --git a/MELA/test/testME_v2.c b/MELA/test/testME_v2.c
index 824fc98e..da0c5be2 100755
--- a/MELA/test/testME_v2.c
+++ b/MELA/test/testME_v2.c
@@ -6,6 +6,7 @@
 #include <fstream>
 #include <cstdlib>
 #include <iomanip>
+#include <memory>
 #include "Mela.h"
 #include "TMath.h"
 #include "TLorentzVector.h"
@@ -2218,7 +2219,7 @@ void testME_ProdDec_MCFM_JHUGen_WBFZZWW_Comparison_Ping(int motherflavor=0, int
     mela.setInputEvent(&daughters, &associated, &mothers, true);
 
     if (mothers.size()>1){
-      if (vbfvhchannel==0) TUtil::computeVBFangles(
+      if (vbfvhchannel==0) TUtil::computeVBFAngles(
         costhetastar,
         costheta1,
         costheta2,
@@ -2238,7 +2239,7 @@ void testME_ProdDec_MCFM_JHUGen_WBFZZWW_Comparison_Ping(int motherflavor=0, int
         &(mothers.at(0).second), mothers.at(0).first,
         &(mothers.at(1).second), mothers.at(1).first
         );
-      else TUtil::computeVHangles(
+      else TUtil::computeVHAngles(
         costhetastar,
         costheta1,
         costheta2,
@@ -2258,7 +2259,7 @@ void testME_ProdDec_MCFM_JHUGen_WBFZZWW_Comparison_Ping(int motherflavor=0, int
         );
     }
     else{
-      if (vbfvhchannel==0) TUtil::computeVBFangles(
+      if (vbfvhchannel==0) TUtil::computeVBFAngles(
         costhetastar,
         costheta1,
         costheta2,
@@ -2275,7 +2276,7 @@ void testME_ProdDec_MCFM_JHUGen_WBFZZWW_Comparison_Ping(int motherflavor=0, int
         associated.at(0).second, associated.at(0).first,
         associated.at(1).second, associated.at(1).first
         );
-      else TUtil::computeVHangles(
+      else TUtil::computeVHAngles(
         costhetastar,
         costheta1,
         costheta2,
@@ -3929,7 +3930,7 @@ void testME_ProdDec_MCFM_JHUGen_WBFZZWW_TU_Comparison_Ping(int motherflavor=0, i
     mela.setInputEvent(&daughters_tu, &associated_tu, &mothers, true);
 
     if (mothers.size()>1){
-      if (vbfvhchannel==0) TUtil::computeVBFangles(
+      if (vbfvhchannel==0) TUtil::computeVBFAngles(
         costhetastar,
         costheta1,
         costheta2,
@@ -3949,7 +3950,7 @@ void testME_ProdDec_MCFM_JHUGen_WBFZZWW_TU_Comparison_Ping(int motherflavor=0, i
         &(mothers.at(0).second), mothers.at(0).first,
         &(mothers.at(1).second), mothers.at(1).first
         );
-      else TUtil::computeVHangles(
+      else TUtil::computeVHAngles(
         costhetastar,
         costheta1,
         costheta2,
@@ -3969,7 +3970,7 @@ void testME_ProdDec_MCFM_JHUGen_WBFZZWW_TU_Comparison_Ping(int motherflavor=0, i
         );
     }
     else{
-      if (vbfvhchannel==0) TUtil::computeVBFangles(
+      if (vbfvhchannel==0) TUtil::computeVBFAngles(
         costhetastar,
         costheta1,
         costheta2,
@@ -3986,7 +3987,7 @@ void testME_ProdDec_MCFM_JHUGen_WBFZZWW_TU_Comparison_Ping(int motherflavor=0, i
         associated.at(0).second, associated.at(0).first,
         associated.at(1).second, associated.at(1).first
         );
-      else TUtil::computeVHangles(
+      else TUtil::computeVHAngles(
         costhetastar,
         costheta1,
         costheta2,
@@ -4772,7 +4773,7 @@ void testME_ProdDec_MCFM_JHUGen_Comparison(int flavor=2, bool useBkgSample=false
       mela.setInputEvent(&daughters_ZZ, &associated, &mothers, true);
 
       if (mothers.size()>1){
-        if (vbfvhchannel==1) TUtil::computeVBFangles(
+        if (vbfvhchannel==1) TUtil::computeVBFAngles(
           costhetastar,
           costheta1,
           costheta2,
@@ -4792,7 +4793,7 @@ void testME_ProdDec_MCFM_JHUGen_Comparison(int flavor=2, bool useBkgSample=false
           &(mothers.at(0).second), mothers.at(0).first,
           &(mothers.at(1).second), mothers.at(1).first
           );
-        else TUtil::computeVHangles(
+        else TUtil::computeVHAngles(
           costhetastar,
           costheta1,
           costheta2,
@@ -4812,7 +4813,7 @@ void testME_ProdDec_MCFM_JHUGen_Comparison(int flavor=2, bool useBkgSample=false
           );
       }
       else{
-        if (vbfvhchannel==1) TUtil::computeVBFangles(
+        if (vbfvhchannel==1) TUtil::computeVBFAngles(
           costhetastar,
           costheta1,
           costheta2,
@@ -4829,7 +4830,7 @@ void testME_ProdDec_MCFM_JHUGen_Comparison(int flavor=2, bool useBkgSample=false
           associated.at(0).second, associated.at(0).first,
           associated.at(1).second, associated.at(1).first
           );
-        else TUtil::computeVHangles(
+        else TUtil::computeVHAngles(
           costhetastar,
           costheta1,
           costheta2,
diff --git a/PythonWrapper/src/MEMCalculatorsWrapper.cc b/PythonWrapper/src/MEMCalculatorsWrapper.cc
index 7dd4fb98..e5ded938 100755
--- a/PythonWrapper/src/MEMCalculatorsWrapper.cc
+++ b/PythonWrapper/src/MEMCalculatorsWrapper.cc
@@ -16,11 +16,12 @@ MEMCalculatorsWrapper::computeAngles(TLorentzVector Z1_lept1, int Z1_lept1Id,
 		   TLorentzVector Z2_lept1, int Z2_lept1Id,
 		   TLorentzVector Z2_lept2, int Z2_lept2Id) {
     Angles ret;
-    TUtil::computeAngles(Z1_lept1,Z1_lept1Id,
-                        Z1_lept2,Z1_lept2Id,
-                        Z2_lept1,Z2_lept1Id,
-                        Z2_lept2,Z2_lept2Id,
-                    ret.costhetastar,ret.costheta1,ret.costheta2,ret.phi,ret.phistar1);
+    TUtil::computeAngles(ret.costhetastar, ret.costheta1, ret.costheta2, ret.phi, ret.phistar1,
+      Z1_lept1, Z1_lept1Id,
+      Z1_lept2, Z1_lept2Id,
+      Z2_lept1, Z2_lept1Id,
+      Z2_lept2, Z2_lept2Id
+      );
     return ret;
   }