Merge pull request #31 from usarica/newVH

Prepare for v2.1.2

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
 
-}
+}

MELA/data/slc6_amd64_gcc700/download.url

@@ -0,0 +1 @@
+http://spin.pha.jhu.edu/Generator/MCFM-precompiled/v2/

MELA/data/slc7_amd64_gcc630/download.url

@@ -0,0 +1 @@
+http://spin.pha.jhu.edu/Generator/MCFM-precompiled/v2/

MELA/data/slc7_amd64_gcc700/download.url

@@ -0,0 +1 @@
+http://spin.pha.jhu.edu/Generator/MCFM-precompiled/v2/

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

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

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.
 !=====================================================
 
 

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],

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;

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,

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,14 +614,17 @@ 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)
     dec = m.computeP(False)
     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)