Add LHCB_WENU_8TEV data set

README.md

@@ -111,5 +111,11 @@ The contents of this directory are:
 * `pineappl.orders`: Output of `pineappl orders`
 * `pineappl.pdf_uncertainty`: Output of `pineappl pdf_uncertainty`
 * `results.log`: The numerical results of the run, comparing the results of the
-  grid against the results from `mg5_aMC`
+  grid against the results from `mg5_aMC`. The first column (PineAPPL) are the
+  interpolated results, which should be similar to the Monte Carlo (MC) results
+  in the second column. The third column gives the relative MC uncertainty
+  (sigma). The next column gives the differences in terms of sigma. The final
+  three columns give the per mille differences of the central, minimum, and
+  maximum scale varied results. Ideally the first two columns are the same and
+  the remaining columns are zero.
 * `time.log`: Total `time` needed for the run

nnpdf31_proc/LHCB_WENU_8TEV_WM/analysis.f

@@ -0,0 +1,61 @@
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_begin(nwgt,weights_info)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      integer nwgt
+      character*(*) weights_info(*)
+
+      call set_error_estimation(1)
+      call HwU_inithist(nwgt,weights_info)
+      call HwU_book(1,'eta', 7, 2.0d0,  3.75d0)
+      call HwU_book(2,'eta', 1, 3.75d0, 4.25d0)
+
+      return
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_end(dummy)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      double precision dummy
+      call HwU_write_file
+      return
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      include 'nexternal.inc'
+      include 'cuts.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      integer ibody  
+      integer i
+      double precision p(0:4,nexternal)
+      double precision wgts(*)
+      double precision xeta, eta_04
+      external eta_04
+
+      double precision p_reco(0:4,nexternal)
+      integer iPDG_reco(nexternal)
+
+
+
+      call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
+     $                       p, iPDG, p_reco, iPDG_reco)
+
+      do i = nincoming+1, nexternal
+        if (iPDG_reco(i).eq.11) then
+          xeta = eta_04(p_reco(0,i))
+        endif
+      enddo
+
+      call HwU_fill(1,xeta,wgts)
+      call HwU_fill(2,xeta,wgts)
+
+ 999  return
+      end

nnpdf31_proc/LHCB_WENU_8TEV_WM/launch.txt

@@ -0,0 +1,22 @@
+launch @OUTPUT@
+fixed_order = ON
+set mw @MW@
+set mz @MZ@
+set ymt @YMT@
+set ebeam1 4000
+set ebeam2 4000
+set pdlabel lhapdf
+set lhaid 324900
+set fixed_ren_scale True
+set fixed_fac_scale True
+set mur_ref_fixed @MW@
+set muf_ref_fixed @MW@
+set reweight_scale True
+set ptl = 20.0
+set etal = 4.25
+#user_defined_cut set minetal = 2.0
+set rphreco = 0.1
+set req_acc_FO 0.0002
+set pineappl True
+done
+quit

nnpdf31_proc/LHCB_WENU_8TEV_WM/metadata.txt

@@ -0,0 +1,5 @@
+arxiv=1608.01484
+description=LHCb differential Drell–Yan cross section at 8 TeV
+x1_label=eta_l
+x1_label_tex=$\eta_\ell$
+y_unit=pb

nnpdf31_proc/LHCB_WENU_8TEV_WM/output.txt

@@ -0,0 +1,7 @@
+set complex_mass_scheme True
+import model loop_qcd_qed_sm_Gmu
+define p = p b b~
+define j = p
+generate p p > e- ve~ [QCD QED]
+output @OUTPUT@
+quit

nnpdf31_proc/LHCB_WENU_8TEV_WM/postrun.sh

@@ -0,0 +1,5 @@
+#!/bin/bash
+
+# divide the grid values by two
+pineappl merge "${GRID}".tmp "${GRID}" --scale 0.5
+mv "${GRID}".tmp "${GRID}"

nnpdf31_proc/LHCB_WENU_8TEV_WP/analysis.f

@@ -0,0 +1,61 @@
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_begin(nwgt,weights_info)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      integer nwgt
+      character*(*) weights_info(*)
+
+      call set_error_estimation(1)
+      call HwU_inithist(nwgt,weights_info)
+      call HwU_book(1,'eta', 7, 2.0d0,  3.75d0)
+      call HwU_book(2,'eta', 1, 3.75d0, 4.25d0)
+
+      return
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_end(dummy)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      double precision dummy
+      call HwU_write_file
+      return
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      include 'nexternal.inc'
+      include 'cuts.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      integer ibody  
+      integer i
+      double precision p(0:4,nexternal)
+      double precision wgts(*)
+      double precision xeta, eta_04
+      external eta_04
+
+      double precision p_reco(0:4,nexternal)
+      integer iPDG_reco(nexternal)
+
+
+
+      call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
+     $                       p, iPDG, p_reco, iPDG_reco)
+
+      do i = nincoming+1, nexternal
+        if (iPDG_reco(i).eq.-11) then
+          xeta = eta_04(p_reco(0,i))
+        endif
+      enddo
+
+      call HwU_fill(1,xeta,wgts)
+      call HwU_fill(2,xeta,wgts)
+
+ 999  return
+      end

nnpdf31_proc/LHCB_WENU_8TEV_WP/launch.txt

@@ -0,0 +1,22 @@
+launch @OUTPUT@
+fixed_order = ON
+set mw @MW@
+set mz @MZ@
+set ymt @YMT@
+set ebeam1 4000
+set ebeam2 4000
+set pdlabel lhapdf
+set lhaid 324900
+set fixed_ren_scale True
+set fixed_fac_scale True
+set mur_ref_fixed @MW@
+set muf_ref_fixed @MW@
+set reweight_scale True
+set ptl = 20.0
+set etal = 4.25
+#user_defined_cut set minetal = 2.0
+set rphreco = 0.1
+set req_acc_FO 0.0002
+set pineappl True
+done
+quit

nnpdf31_proc/LHCB_WENU_8TEV_WP/metadata.txt

@@ -0,0 +1,5 @@
+arxiv=1608.01484
+description=LHCb differential Drell–Yan cross section at 8 TeV
+x1_label=eta_l
+x1_label_tex=$\eta_{\bar{\ell}}$
+y_unit=pb

nnpdf31_proc/LHCB_WENU_8TEV_WP/output.txt

@@ -0,0 +1,7 @@
+set complex_mass_scheme True
+import model loop_qcd_qed_sm_Gmu
+define p = p b b~
+define j = p
+generate p p > e+ ve [QCD QED]
+output @OUTPUT@
+quit

nnpdf31_proc/LHCB_WENU_8TEV_WP/postrun.sh

@@ -0,0 +1,5 @@
+#!/bin/bash
+
+# divide the grid values by two
+pineappl merge "${GRID}".tmp "${GRID}" --scale 0.5
+mv "${GRID}".tmp "${GRID}"

run.sh

@@ -130,48 +130,97 @@ EOF
     grid="${dataset}".pineappl
 
     # merge the final bins
-    "${pineappl}" merge "${grid}" $(ls -v "${dataset}"/Events/run_01*/amcblast_obs_*.root)
-
-    lz4=$(which lz4 2> /dev/null || true)
-
-    # compress the grid with `lz4` if it's available
-    if [[ -x ${lz4} ]]; then
-        lz4 -9 "${grid}"
-        rm "${grid}"
-        grid="${grid}.lz4"
+    "${pineappl}" merge "${grid}" $(ls -v "${dataset}"/Events/run_01*/amcblast_obs_*.pineappl)
+
+    # optimize the grids
+    "${pineappl}" optimize "${grid}" "${grid}".tmp
+    mv "${grid}".tmp "${grid}"
+
+    # add metadata
+    runcard="${dataset}"/Events/run_01*/run_01_tag_1_banner.txt
+    if [[ -f ../nnpdf31_proc/"${dataset}"/metadata.txt ]]; then
+        eval $(awk -F= "BEGIN { printf \"pineappl set ${grid} ${grid}.tmp \" }
+                              { printf \"--entry %s '%s' \", \$1, \$2 }
+                        END   { printf \"--entry_from_file runcard ${runcard}\\n\" }" \
+            ../nnpdf31_proc/"${dataset}"/metadata.txt)
+    else
+        "${pineappl}" set "${grid}" "${grid}".tmp --entry_from_file runcard ${runcard}
     fi
+    mv "${grid}".tmp "${grid}"
 
     # find out which PDF set was used to generate the predictions
     pdfstring=$(grep "set lhaid" "${launch_file}" | sed 's/set lhaid \([0-9]\+\)/\1/')
 
     # (re-)produce predictions
-    "${pineappl}" convolute "${grid}" "${pdfstring}" --scales 9 > pineappl.convolute
+    "${pineappl}" convolute "${grid}" "${pdfstring}" --scales 9 --absolute > pineappl.convolute
     "${pineappl}" orders "${grid}" "${pdfstring}" --absolute > pineappl.orders
-    "${pineappl}" pdf_uncertainty "${grid}" "${pdfstring}" > pineappl.pdf_uncertainty
+    "${pineappl}" pdf_uncertainty --threads=1 "${grid}" "${pdfstring}" > pineappl.pdf_uncertainty
 
     # extract the numerical results from mg5_aMC
     sed '/^  [+-]/!d' "${dataset}"/Events/run_01*/MADatNLO.HwU > results.mg5_aMC
 
     # extract the integrated results from the PineAPPL grid
-    cat pineappl.convolute | head -n -2 | tail -n +5 | awk '{ print $5 }' > results.grid
+    cat pineappl.convolute | head -n -2 | tail -n +5 | \
+        awk '{ print $5, $6, $7, $8, $9, $10, $11, $12, $13, $14 }' > results.grid
 
     # compare the results from the grid and from mg5_aMC
     paste -d ' ' results.grid results.mg5_aMC | awk \
         'function abs(x) { return x < 0.0 ? -x : x; }
-         BEGIN { print "----------------------------------------------------------------------------------------"
-                 print "   PineAPPL       mg5_aMC   mg5_aMC unc.   sigmas   per mille   9-point scale var. (MC)"
-                 print "----------------------------------------------------------------------------------------" }
-         { printf "% e % e %e %8.3f %10.4f % e % e\n",
-                  $1, $4, $5, $5 != 0.0 ? abs($1-$4)/$5 : 0.0, $4 != 0.0 ? abs($1-$4)/$4*1000 : 0.0, $7, $8 }' | \
-        tee results.log
+         function max(x1,x2,x3,x4,x5,x6,x7,x8,x9) {
+             result=x1;
+             if (x2 > result) { result = x2; }
+             if (x3 > result) { result = x3; }
+             if (x4 > result) { result = x4; }
+             if (x5 > result) { result = x5; }
+             if (x6 > result) { result = x6; }
+             if (x7 > result) { result = x7; }
+             if (x8 > result) { result = x8; }
+             if (x9 > result) { result = x9; }
+             return result;
+         }
+         function min(x1,x2,x3,x4,x5,x6,x7,x8,x9) {
+             result=x1;
+             if (x2 < result) { result = x2; }
+             if (x3 < result) { result = x3; }
+             if (x4 < result) { result = x4; }
+             if (x5 < result) { result = x5; }
+             if (x6 < result) { result = x6; }
+             if (x7 < result) { result = x7; }
+             if (x8 < result) { result = x8; }
+             if (x9 < result) { result = x9; }
+             return result;
+         }
+         BEGIN { print "-----------------------------------------------------------------------"
+                 print "   PineAPPL         MC         sigma   diff.   central    min      max "
+                 print "                                       sigma   1/1000   1/1000   1/1000"
+                 print "-----------------------------------------------------------------------" }
+         { printf "% e % e %7.3f%% %7.3f %8.4f %8.4f %8.4f\n",
+                  $1,
+                  $13,
+                  $13 != 0.0 ? $14/$13*100 : 0.0,
+                  $14 != 0.0 ? abs($1-$13)/$14 : 0.0,
+                  $13 != 0.0 ? abs($1-$13)/$13*1000 : 0.0,
+                  $16 != 0.0 ? abs(min($2, $3, $4, $5, $6, $7, $8, $9, $10)-$16)/$16*1000 : 0.0,
+                  $17 != 0.0 ? abs(max($2, $3, $4, $5, $6, $7, $8, $9, $10)-$17)/$17*1000 : 0.0 }' | tee results.log
 
     rm results.mg5_aMC results.grid
 
+    "${pineappl}" set "${grid}" "${grid}".tmp --entry_from_file results results.log
+    mv "${grid}".tmp "${grid}"
+
     # if there is anything to do after the run, do it!
     if [[ -x ../nnpdf31_proc/"${dataset}"/postrun.sh ]]; then
         cp ../nnpdf31_proc/"${dataset}"/postrun.sh .
         GRID=$grid ./postrun.sh
     fi
+
+    lz4=$(which lz4 2> /dev/null || true)
+
+    # compress the grid with `lz4` if it's available
+    if [[ -x ${lz4} ]]; then
+        lz4 -9 "${grid}"
+        rm "${grid}"
+    fi
 }
 
 check_args_and_cd_output "$@"

run_implement_user_defined_cuts.py

@@ -111,6 +111,34 @@
         endif
       enddo
 
+''',
+    'minetal': '''c     cut on the minimum pseudorapidity of leptons
+      do i=1,nexternal-1
+        if (is_a_lm(i) .or. is_a_lp(i)) then
+          if (abs(atanh(p_reco(3,i)/sqrt(p_reco(1,i)**2+p_reco(2,i)**2+
+     &                  p_reco(3,i)**2))) .lt. {}) then
+            passcuts_user=.false.
+            return
+          endif
+        endif
+      enddo
+
+''',
+    'yzmin': '''c     cut on the rapidity of SFOS lepton pairs
+      do i=1,nexternal-1
+        if (is_a_lm(i) .or. is_a_lp(i)) then
+          do j=i+1,nexternal
+            if (ipdg_reco(i) .eq. -ipdg_reco(j)) then
+              if (abs(atanh((p_reco(3,i)+p_reco(3,j))
+     &            /(p_reco(0,i)+p_reco(0,j)))) .lt. {}) then
+                passcuts_user=.false.
+                return
+              endif
+            endif
+          enddo
+        endif
+      enddo
+
 ''',
 }