Merge pull request #12 from scarlehoff/jacobian_factors

Add factors

new_grids/compare_two.py

@@ -16,15 +16,19 @@ def run_comparison(
     dataset_vp,
     dataset_pin,
     pdf_name="NNPDF40_nnlo_as_01180",
-    ratio=1.0,
+    E906_factor=False,
     verbose=True,
-    mass_divide=False,
     remove_points=False,
 ):
     dat = API.dataset(dataset_input={"dataset": dataset_vp}, theoryid=200, use_cuts="internal")
     pdf = API.pdf(pdf=pdf_name)
     lpdf = lhapdf.mkPDF(pdf_name)
     res_vp = central_predictions(dat, pdf)
+
+    if E906_factor:
+        # Fix the result computed by vp for E906 with (Energy_605 / Energy_906)
+        res_vp *= (38.8/15.06)**2
+
     if isinstance(dataset_pin, str):
         if "*" in dataset_pin:
             dataset_pin = [
@@ -36,18 +40,14 @@ def run_comparison(
     res_pine_partial = []
     for d in dataset_pin:
         grid = Grid.read(f"{d}.pineappl.lz4")
-        if mass_divide:
-            # For some reason not all grids need that
-            ratio_mass = grid.raw.bin_left(0) ** 3 / 38.8 ** 3
-        else:
-            ratio_mass = 1.0
 
+        remove_factor = 1.0
         if remove_points:
-            ratio_mass *= np.less(grid.raw.bin_left(0), 7.52) | np.greater(
+            remove_factor *= np.less(grid.raw.bin_left(0), 7.52) | np.greater(
                 grid.raw.bin_left(0), 7.54
             )
         res_pine_partial.append(
-            grid.convolute_with_one(2212, lpdf.xfxQ2, lpdf.alphasQ2) * ratio_mass
+            grid.convolute_with_one(2212, lpdf.xfxQ2, lpdf.alphasQ2) * remove_factor
         )
 
     res_pine_raw = pd.DataFrame(res_pine_partial[0])
@@ -68,31 +68,19 @@ def run_comparison(
     res_pine = res_pine_raw.loc[res_vp.index]
     all_res = pd.concat([res_vp, res_pine, res_vp / res_pine], axis=1)
     all_res.to_csv(f"out_csv/{dataset_vp}.csv", sep="\t", header=["vp", "pine", "ratio"])
-    if ratio != 1.0:
-        all_res_corrected = pd.concat([res_vp, res_pine * ratio, res_vp / res_pine / ratio], axis=1)
-        all_res_corrected.to_csv(
-            f"out_csv/{dataset_vp}_corr.csv", sep="\t", header=["vp", f"pine*{ratio}", "ratio"]
-        )
-    if verbose and ratio != 1.0:
-        print(all_res_corrected)
-    elif verbose:
+    if verbose:
         print(all_res)
 
-    if np.allclose(res_vp, res_pine * ratio, rtol=1e-2):
-        print(f"The pienappl and NNPDF version are compatible!")
+    if np.allclose(res_vp, res_pine, rtol=1.6e-2):
+        print(f"The pineappl and NNPDF version for {dataset_vp} are compatible!")
     else:
-        print(f"XXX The pineappl and NNPDF version are NOT compatible!")
+        print(f"XXX The pineappl and NNPDF version for {dataset_vp} are NOT compatible!")
     return res_vp, res_pine
 
 
 if __name__ == "__main__":
-    #     a, b = run_comparison("DYE605", "E605nlo", ratio=54.35, verbose=False)
-    #     a, b = run_comparison("DYE886P", "rescaled_E886Pnlo", ratio=54.35, mass_divide=False, verbose=True)
-    a, b = run_comparison("DYE886P", "E886Pnlo", ratio=54.35, mass_divide=True, verbose=True)
-#     a, b = run_comparison("DYE886R", ["E886deutRnlo", "E886Rnlo"], verbose=False)
-#     a, b = run_comparison("DYE906_D", "E906deutnlo_bin_*", ratio=54.35, mass_divide=True, remove_points=False)
-#     a, b = run_comparison("DYE906R", "E906*nlo_bin_*", mass_divide=False, remove_points=False)
-#     i = 9
-#     a, b = run_comparison("DYE906_D", f"E906deutnlo_bin_0{i}", ratio=54.35, mass_divide=True, remove_points=False)
-# problems at: i=3,5
-#         a, b = run_comparison("DYE906_D", f"E906nlo_bin_0{i}", ratio=54.35, mass_divide=True, remove_points=False)
+    a, b = run_comparison("DYE605", "E605nlo", verbose=False)
+    a, b = run_comparison("DYE886P", "E866nlo", verbose=False)
+    a, b = run_comparison("DYE886R", ["E866deutRnlo", "E866Rnlo"], verbose=False)
+#     a, b = run_comparison("DYE906_D", "E906deutnlo_bin_*", E906_factor=True, verbose=False)
+    a, b = run_comparison("DYE906R", "E906*nlo_bin_*", verbose=False)

new_grids/inputE886deutnlo.dat → new_grids/inputE866deutnlo.dat

@@ -13,6 +13,7 @@
 
 Collider piso
 E_CM 38.8
+jacobian866 Yes
 
 ####################################################
 # Physical Parameters

new_grids/inputE886nlo.dat → new_grids/inputE866nlo.dat

@@ -13,6 +13,7 @@
 
 Collider pp
 E_CM 38.8
+jacobian866 Yes
 
 ####################################################
 # Physical Parameters

new_grids/inputE906deutnlo.dat

@@ -13,6 +13,7 @@
 
 Collider piso
 E_CM 15.063
+jacobian866 Yes
 
 ####################################################
 # Physical Parameters

new_grids/inputE906nlo.dat

@@ -13,6 +13,7 @@
 
 Collider pp
 E_CM 15.06
+jacobian866 Yes
 
 ####################################################
 # Physical Parameters

new_grids/run_them_all.sh

@@ -0,0 +1,28 @@
+#!/bin/bash
+
+# Script to generate all pineappl tables for NNPDF Theory 400
+
+# NOTE: pineappl optimize will fail if the target grid already exists
+# that's a feature and not a bug, either remove it yourself or use a different name
+
+exe=../build/Vrap
+kinfolder=input_kinematics
+output_placeholder=test.pineappl.lz4
+
+${exe} inputE605nlo.dat input_kinematics/E605.dat
+pineappl optimize ${output_placeholder} E605nlo.pineappl.lz4
+
+## E886P
+${exe} inputE866nlo.dat  input_kinematics/E866P.dat
+pineappl optimize ${output_placeholder} E866nlo.pineappl.lz4
+
+## E886R
+${exe} inputE866nlo.dat  input_kinematics/E866R.dat
+pineappl optimize ${output_placeholder} E866Rnlo.pineappl.lz4
+
+${exe} inputE866deutnlo.dat  input_kinematics/E866R.dat
+pineappl optimize ${output_placeholder} E866deutRnlo.pineappl.lz4
+
+## E906
+./run_E906.sh inputE906nlo.dat
+./run_E906.sh inputE906deutnlo.dat

readme.md

@@ -6,7 +6,23 @@ The code in this repository modifies version 0.9 of Vrap and it is redistributed
 
 ## Changes with respect to Vrap
 
-The code in this repository adds and option to use isoscalar targets as a hadron type (`piso`) in addition to the original `pp` and `ppbar` options. It also interfaces the code with [pineappl](https://n3pdf.github.io/pineappl/), a library that produces fast-interpolation grids for fitting parton distribution functions.
+#### Isoscalar targets
+
+The code in this repository adds and option to use isoscalar targets as a hadron type (`piso`) in addition to the original `pp` and `ppbar` options. 
+
+### Pineappl interface
+
+It also interfaces the code with [pineappl](https://n3pdf.github.io/pineappl/), a library that produces fast-interpolation grids for fitting parton distribution functions.
+
+
+#### Apfel-jacobian factors
+
+The calculations used in NNPDF need to match the format of the experimental data. For this a jacobian factor of $\sqrt{2s}$ is applied with respect to the vanilla vrap implementation.
+
+In addition, the data from DYE886 and DYE605 differ by a factor of $\left(\frac{\sqrt{s}}{M}\right)$, which is added as an option to the runcard:
+`jacobian886: True`.
+
+See [issue #10](https://github.com/scarlehoff/Hawaiian_vrap/issues/10) for more information.
 
 :warning: The goal of these modifications is to generate `pineappl` grids for fixed-target Drell-Yan. These grids are used by the NNPDF Collaboration to determine parton distribution functions. As such, this use-case is the only one that has been tested and benchmarked (see [regression tests](https://github.com/scarlehoff/Hawaiian_vrap/tree/main/regression_test)).
 

src/Vrap.C

@@ -54,7 +54,8 @@ VrapOptionsHandler::VrapOptionsHandler(){
 	add(new multipleValueOption<exchange>("VectorBoson",exchM,"gamma_only",gamma_only,"Zgamma_interf",Zgamma_interf,"Z_only",Z_only,"Zgamma",Zgamma,"Wplus",Wplus,"Wminus",Wminus,"Sets which vector boson to use."));
 	add(new ValueSettingOption<int>("NumberOfYPoints",nbrYPnts,"Sets the number of rapidity values at which to compute."));
 	add(new multipleValueOption<int>("PrintDirection",direction,"Forward",+1,"Reverse",-1,"Sets the order in which output is printed (for td). ") );  
-        add(new multipleValueOption<int>("OutputFormat",o_f,"TopDrawStyle",0,"ListValues",1,"Sets the style output is printed (for td or just list dsig/dy). ") );
+  add(new multipleValueOption<int>("OutputFormat",o_f,"TopDrawStyle",0,"ListValues",1,"Sets the style output is printed (for td or just list dsig/dy). ") );
+	add(new yesOrNoOption("jacobian866",jacobian866,"Sets whether to use the jacobian of (sqrt(s)/M)^3."));
 	//enableDebug();
 }
 

src/integration.h

@@ -634,6 +634,12 @@ DVector rap_y(){
         std::cout << std::setw(9) << std::setprecision(9) << " working on    y =  " << y << std::endl; 
     }
 
+    // Multiply the jacobian factor from https://github.com/NNPDF/Hawaiian_vrap/issues/10
+    if (jacobian866) {
+      prefactor *= pow(Q/E_CM,3);
+    }
+    prefactor *= sqrt(2.0)*E_CM;
+
     // Set the prefactor for the pineappl grid
     // this needs to be done at this stage since it is not included in the weights
     piner.set_prefactor(prefactor);

src/options.C

@@ -38,6 +38,7 @@ std::ostream& singleValueOption::printHelp(std::ostream& os) const {
 bool yesOrNoOption::process(std::istream& is,std::string& message) const {
 	std::string answer;
 	is >> answer;
+    std::cout << answer ;
 	if ( answer == "Yes" || answer == "On" || answer == "yes" || answer == "on" || answer == "YES" || answer == "ON"  ){
 		d_valueToSet=true;
 		message = "" ;

src/settings.h

@@ -40,6 +40,7 @@ std::string pdfFile;
 
 bool useMyAlphaRunning = false;
 bool useOtherPDF;
+bool jacobian866 = false;
 int pdfSet;
 int parton_flag; // 1: all, 2:qqbar 3: qg 4: gg 5:qq 6: qqbar_plus_qq