Merge pull request #39 from usarica/newVH

v2.1.8

MELA/COLLIER/setup.sh

@@ -12,7 +12,7 @@ cd $scriptdir
 pkgname="collier-1.2"
 pkgdir="COLLIER-1.2"
 tarname=$pkgname".tar.gz"
-tarweb="https://www.hepforge.org/archive/collier/"$tarname
+tarweb="https://collier.hepforge.org/downloads/"$tarname
 libname="libcollier.so"
 tmpdir="colliertmp"
 
@@ -35,14 +35,14 @@ if [[ $# > 0 ]] && [[ "$1" == *"clean"* ]];then
 else
 
   if ! [ -f "../data/$SCRAM_ARCH/$libname" ]; then
-    wget $tarweb
+    wget --no-check-certificate $tarweb
     mkdir $tmpdir
     tar -xvzf $tarname -C $tmpdir
     rm $tarname
     mv $tmpdir"/"$pkgdir"/src/"* ./
     rm -rf $tmpdir
 
-    make
+    make $@
     cp $libname "../data/$SCRAM_ARCH/$libname"
   fi
 

MELA/fortran/mod_Kinematics.F90

@@ -2,7 +2,7 @@ MODULE ModKinematics
 implicit none
 save
 
-public :: GetBWPropagator,ReweightBWPropagator,SetRunningScales,setPDFs,EvalAlphaS
+public :: GetBWPropagator,ReweightBWPropagator,SetRunningScales,setPDFs!,EvalAlphaS
 
 CONTAINS
 
@@ -1139,11 +1139,12 @@ subroutine SetRunningScales(p,id) ! p in JHU-GeV, id in JHUGen conventions
 real(dp), intent(in) :: p(1:4,4:6) ! No need to run the second index from 3 to 7: pH, pJ1, pJ2
 integer, intent(in) :: id(4:7) ! id_JJH/id_JJVV, id_J1, id_J2, id_JJ (if applicable)
 real(8) :: polemass(3:7) ! mJJH, mH, mJ1, mJ2, mJJ (if applicable)
-real(8) :: pJJHstar(4),pHstar(4),pJ(4,2),pJJ(4),pJHstar(4)
+real(8) :: pJJHstar(4),pHstar(4),pJ(4,2),pJJ(4),pJHstar(4),pTjet(5:6),maxpTjet,minpTjet
 integer idx,ip
 
    pHstar(:) = 0d0
    pJJ(:) = 0d0
+   pTjet(:) = 0d0
    polemass(3) = getMass(id(4)) ! Pole mass of the JJH system
    polemass(4) = M_Reso
    do idx=4,6
@@ -1156,8 +1157,11 @@ subroutine SetRunningScales(p,id) ! p in JHU-GeV, id in JHUGen conventions
          do ip=1,4
             pJJ(ip) = pJJ(ip) + p(ip,idx)
          enddo
+         pTjet(idx) = get_PT(p(1:4,idx))
       endif
    enddo
+   maxpTjet = maxval(pTjet)
+   minpTjet = minval(pTjet, mask=.not.all(p(1:4,5:6).eq.0d0, 1))
    polemass(7) = getMass(id(7)) ! Pole mass of the JJ system
 
    pJJHstar = pJJ + pHstar
@@ -1211,6 +1215,12 @@ subroutine SetRunningScales(p,id) ! p in JHU-GeV, id in JHUGen conventions
       Mu_Fact = Get_MInv(pJ(1:4,1))
    elseif(FacScheme .eq. -kRenFacScheme_mj) then
       Mu_Fact = polemass(5)
+   elseif(FacScheme .eq. kRenFacScheme_maxpTj) then
+      Mu_Fact = maxpTjet
+   elseif(FacScheme .eq. kRenFacScheme_minpTj) then
+      Mu_Fact = minpTjet
+   else
+      call Error("This should never be able to happen.")
    endif
 
    ! Do the same for the renormalization scale
@@ -1253,6 +1263,12 @@ subroutine SetRunningScales(p,id) ! p in JHU-GeV, id in JHUGen conventions
       Mu_Ren = Get_MInv(pJ(1:4,1))
    elseif(RenScheme .eq. -kRenFacScheme_mj) then
       Mu_Ren = polemass(5)
+   elseif(RenScheme .eq. kRenFacScheme_maxpTj) then
+      Mu_Ren = maxpTjet
+   elseif(RenScheme .eq. kRenFacScheme_minpTj) then
+      Mu_Ren = minpTjet
+   else
+      call Error("This should never be able to happen.")
    endif
 
    ! Never ever allow the scales to go negative
@@ -1493,52 +1509,52 @@ SUBROUTINE CTEQ6(X,SCALE,UPV,DNV,USEA,DSEA,STR,CHM,BOT,GLU)
 RETURN
 END SUBROUTINE
 
-! QCD scale from MCFM
-! Implementation into JHUGen by Ulascan Sarica, Dec. 2015
-subroutine EvalAlphaS()
-   use ModParameters
-   IMPLICIT NONE
-#if useLHAPDF==1
-!--- This is simply a wrapper to the LHAPDF implementation of the running coupling alphas, in the style of the native MCFM routine
-   DOUBLE PRECISION alphasPDF
-   REAL(DP) :: Q
-      Q = Mu_Ren/GeV
-      alphas=alphasPDF(Q)
-#else
-!     Evaluation of strong coupling constant alphas
-!     Original Author: R.K. Ellis
-!     q -- Scale at which alpha_s is to be evaluated
-!     alphas_mz -- ModParameters value of alpha_s at the mass of the Z-boson
-!     nloops_pdf -- ModParameters value of the number of loops (1,2, or 3) at which the beta function is evaluated to determine running.
-!     If you somehow need a more complete implementation, check everything at or before commit 28472c5bfee128dde458fd4929b4d3ece9519ab8
-   INTEGER, PARAMETER :: NF6=6
-   INTEGER, PARAMETER :: NF5=5
-   INTEGER, PARAMETER :: NF4=4
-   INTEGER, PARAMETER :: NF3=3
-   INTEGER, PARAMETER :: NF2=2
-   INTEGER, PARAMETER :: NF1=1
-
-      IF (Mu_Ren .LE. 0d0) THEN
-         WRITE(6,*) 'ModKinematics::EvalAlphaS: Mu_Ren .le. 0, Mu_Ren (GeV) = ',(Mu_Ren*GeV)
-         stop
-      ENDIF
-      IF (nQflavors_pdf .NE. NF5) THEN
-         WRITE(6,*) 'ModKinematics::EvalAlphaS: nQflavors_pdf invalid, nQflavors_pdf = ',nQflavors_pdf
-         WRITE(6,*) 'ModKinematics::EvalAlphaS: Check 28472c5bfee128dde458fd4929b4d3ece9519ab8'
-         stop
-      ENDIF
-      IF (nloops_pdf .NE. 1) THEN
-         WRITE(6,*) 'ModKinematics::EvalAlphaS: nloops_pdf invalid, nloops_pdf = ',nloops_pdf
-         WRITE(6,*) 'ModKinematics::EvalAlphaS: Check 28472c5bfee128dde458fd4929b4d3ece9519ab8'
-         stop
-      ENDIF
-
-      alphas=alphas_mz/(1.0_dp+alphas_mz*B0_PDF(NF5)*2.0_dp*dlog((Mu_Ren/zmass_pdf)))
-#endif
-      ! Calculate the derived couplings
-      call ComputeQCDVariables()
-   RETURN
-end subroutine EvalAlphaS
+!! QCD scale from MCFM
+!! Implementation into JHUGen by Ulascan Sarica, Dec. 2015
+!subroutine EvalAlphaS()
+!   use ModParameters
+!   IMPLICIT NONE
+!#if useLHAPDF==1
+!!--- This is simply a wrapper to the LHAPDF implementation of the running coupling alphas, in the style of the native MCFM routine
+!   DOUBLE PRECISION alphasPDF
+!   REAL(DP) :: Q
+!      Q = Mu_Ren/GeV
+!      alphas=alphasPDF(Q)
+!#else
+!!     Evaluation of strong coupling constant alphas
+!!     Original Author: R.K. Ellis
+!!     q -- Scale at which alpha_s is to be evaluated
+!!     alphas_mz -- ModParameters value of alpha_s at the mass of the Z-boson
+!!     nloops_pdf -- ModParameters value of the number of loops (1,2, or 3) at which the beta function is evaluated to determine running.
+!!     If you somehow need a more complete implementation, check everything at or before commit 28472c5bfee128dde458fd4929b4d3ece9519ab8
+!   INTEGER, PARAMETER :: NF6=6
+!   INTEGER, PARAMETER :: NF5=5
+!   INTEGER, PARAMETER :: NF4=4
+!   INTEGER, PARAMETER :: NF3=3
+!   INTEGER, PARAMETER :: NF2=2
+!   INTEGER, PARAMETER :: NF1=1
+!
+!      IF (Mu_Ren .LE. 0d0) THEN
+!         WRITE(6,*) 'ModKinematics::EvalAlphaS: Mu_Ren .le. 0, Mu_Ren (GeV) = ',(Mu_Ren*GeV)
+!         stop
+!      ENDIF
+!      IF (nQflavors_pdf .NE. NF5) THEN
+!         WRITE(6,*) 'ModKinematics::EvalAlphaS: nQflavors_pdf invalid, nQflavors_pdf = ',nQflavors_pdf
+!         WRITE(6,*) 'ModKinematics::EvalAlphaS: Check 28472c5bfee128dde458fd4929b4d3ece9519ab8'
+!         stop
+!      ENDIF
+!      IF (nloops_pdf .NE. 1) THEN
+!         WRITE(6,*) 'ModKinematics::EvalAlphaS: nloops_pdf invalid, nloops_pdf = ',nloops_pdf
+!         WRITE(6,*) 'ModKinematics::EvalAlphaS: Check 28472c5bfee128dde458fd4929b4d3ece9519ab8'
+!         stop
+!      ENDIF
+!
+!      alphas=alphas_mz/(1.0_dp+alphas_mz*B0_PDF(NF5)*2.0_dp*dlog((Mu_Ren/zmass_pdf)))
+!#endif
+!      ! Calculate the derived couplings
+!      call ComputeQCDVariables()
+!   RETURN
+!end subroutine EvalAlphaS
 
 
 

MELA/fortran/mod_Misc.F90

@@ -22,8 +22,14 @@ MODULE ModMisc
 END INTERFACE OPERATOR (.cross.)
 
 interface BubleSort
-   module procedure BubleSort_realinteger, BubleSort_stringlogical, BubleSort_stringinteger, BubleSort_stringreal8, BubleSort_stringcomplex8, BubleSort_stringstring
+   module procedure BubleSort_integerinteger,BubleSort_realinteger, BubleSort_stringlogical, BubleSort_stringinteger, BubleSort_stringreal8, BubleSort_stringcomplex8, BubleSort_stringstring
 end interface BubleSort
+
+interface ReOrder
+   module procedure ReOrder_integerinteger,ReOrder_integerreal
+end interface ReOrder
+
+
 contains
 
 
@@ -287,6 +293,7 @@ SUBROUTINE name(N,X, IY)                 ;\
 RETURN                                   ;\
 END SUBROUTINE
 
+MakeBubleSort(integer(8), integer, BubleSort_integerinteger)
 MakeBubleSort(real(8), integer, BubleSort_realinteger)
 MakeBubleSort(character(len=100), logical, BubleSort_stringlogical)
 MakeBubleSort(character(len=100), integer, BubleSort_stringinteger)
@@ -308,6 +315,25 @@ SUBROUTINE name(N,X, IY)                 ;\
 !     stop
 
 
+#define MakeReOrder(typey, typex, name)      \
+SUBROUTINE name(N, iy,x)                    ;\
+integer :: n,i                              ;\
+typex :: x(1:n)                             ;\
+typey :: iy(1:n)                            ;\
+typex :: temp(1:n)                          ;\
+                                            ;\
+     do i=1,N                               ;\
+         temp(i) = x( iy(i) )               ;\
+     enddo                                  ;\
+     x(1:N) = temp(1:N)                     ;\
+                                            ;\
+RETURN                                      ;\
+END SUBROUTINE                              ;\
+
+MakeReOrder(integer,integer,  ReOrder_integerinteger)
+MakeReOrder(integer, real(8), ReOrder_integerreal)
+#undef MakeReOrder
+
 
 SUBROUTINE printMom(Mom)
 implicit none
@@ -1509,6 +1535,8 @@ function CalculatesXsec(Process)
     CalculatesXsec=.true.
   elseif (Process.eq.50) then
     CalculatesXsec=.false.
+  elseif (Process.ge.51 .and. Process.le.52) then
+    CalculatesXsec=.false.
   elseif (Process.eq.60) then
     CalculatesXsec=.true.
   elseif (Process.eq.61) then