Implement perturbative uncertainties

pineappl/src/grid.rs

@@ -286,22 +286,20 @@ impl Grid {
         order_mask: &[bool],
         bin_indices: &[usize],
         lumi_mask: &[bool],
-        xi: &(f64, f64),
+        xi: &[(f64, f64)],
     ) -> Vec<f64> {
         let bin_indices = if bin_indices.is_empty() {
             (0..self.bin_limits.bins()).collect()
         } else {
             bin_indices.to_vec()
         };
-        let mut bins: Vec<f64> = vec![0.0; bin_indices.len()];
+        let mut bins: Vec<f64> = vec![0.0; bin_indices.len() * xi.len()];
         let bin_sizes = self.bin_limits.bin_sizes();
 
         for ((i, j, k), subgrid) in self.subgrids.indexed_iter() {
             let order = &self.orders[i];
 
             if (!order_mask.is_empty() && !order_mask[i])
-                || ((order.logxir > 0) && (xi.0 == 1.0))
-                || ((order.logxif > 0) && (xi.1 == 1.0))
                 || (!lumi_mask.is_empty() && !lumi_mask[k])
             {
                 continue;
@@ -314,27 +312,33 @@ impl Grid {
 
             let lumi_entry = &self.lumi[k];
 
-            let mut value = subgrid.convolute(&|x1, x2, q2| {
-                let mut lumi = 0.0;
-                let q2 = xi.1 * q2;
-
-                for entry in lumi_entry.entry() {
-                    lumi += xfx1(entry.0, x1, q2) * xfx2(entry.1, x2, q2) * entry.2 / (x1 * x2);
+            for (l, &(xir, xif)) in xi.iter().enumerate() {
+                if ((order.logxir > 0) && (xir == 1.0)) || ((order.logxif > 0) && (xif == 1.0)) {
+                    continue;
                 }
 
-                lumi *= alphas(q2).powi(order.alphas as i32);
-                lumi
-            });
+                let mut value = subgrid.convolute(&|x1, x2, q2| {
+                    let mut lumi = 0.0;
+                    let q2 = xif * q2;
 
-            if order.logxir > 0 {
-                value *= xi.0.ln().powi(order.logxir as i32);
-            }
+                    for entry in lumi_entry.entry() {
+                        lumi += xfx1(entry.0, x1, q2) * xfx2(entry.1, x2, q2) * entry.2 / (x1 * x2);
+                    }
 
-            if order.logxif > 0 {
-                value *= xi.1.ln().powi(order.logxif as i32);
-            }
+                    lumi *= alphas(q2).powi(order.alphas as i32);
+                    lumi
+                });
+
+                if order.logxir > 0 {
+                    value *= xir.ln().powi(order.logxir as i32);
+                }
+
+                if order.logxif > 0 {
+                    value *= xif.ln().powi(order.logxif as i32);
+                }
 
-            bins[bin_index] += value / bin_sizes[j];
+                bins[l + xi.len() * bin_index] += value / bin_sizes[j];
+            }
         }
 
         bins

pineappl_cli/src/bin/pineappl.rs

@@ -67,6 +67,7 @@ fn convolute(
     pdfset: &str,
     other_pdfsets: &[&str],
     show_bins: &[usize],
+    scales: usize,
 ) -> Result<(), Box<dyn Error>> {
     let grid = Grid::read(BufReader::new(File::open(input)?))?;
     let show_bins = if show_bins.is_empty() {
@@ -77,6 +78,17 @@ fn convolute(
     let pdf = str::parse::<i32>(pdfset)
         .map(Pdf::with_lhaid)
         .unwrap_or_else(|_| Pdf::with_setname_and_member(pdfset, 0));
+    let scales_vector = vec![
+        (1.0, 1.0),
+        (2.0, 2.0),
+        (0.5, 0.5),
+        (2.0, 1.0),
+        (1.0, 2.0),
+        (0.5, 1.0),
+        (1.0, 0.5),
+        (2.0, 0.5),
+        (0.5, 2.0),
+    ];
 
     let results = grid.convolute(
         &|id, x1, q2| pdf.xfx_q2(id, x1, q2),
@@ -85,7 +97,7 @@ fn convolute(
         &[],
         &show_bins,
         &[],
-        &[(1.0, 1.0)],
+        &scales_vector[0..scales],
     );
 
     let other_results: Vec<f64> = other_pdfsets
@@ -101,24 +113,35 @@ fn convolute(
                 &[],
                 &show_bins,
                 &[],
-                &(1.0, 1.0),
+                &scales_vector[0..scales],
             )
         })
         .flatten()
         .collect();
 
     let bin_sizes = grid.bin_limits().bin_sizes();
 
-    for (bin, value) in results.iter().enumerate() {
+    for (bin, values) in results.chunks_exact(scales).enumerate() {
+        let min_value = values
+            .iter()
+            .min_by(|left, right| left.partial_cmp(right).unwrap())
+            .unwrap();
+        let max_value = values
+            .iter()
+            .max_by(|left, right| left.partial_cmp(right).unwrap())
+            .unwrap();
+
         print!(
-            "{:<3}  {:>12.7e}  {:>12.7e}",
+            "{:<3}  {:>12.7e}  {:>12.7e}  {:+5.2}% {:+5.2}%",
             show_bins[bin],
-            value,
-            value * bin_sizes[show_bins[bin]],
+            values[0],
+            values[0] * bin_sizes[show_bins[bin]],
+            (min_value / values[0] - 1.0) * 100.0,
+            (max_value / values[0] - 1.0) * 100.0,
         );
 
         for other in other_results.iter().skip(bin).step_by(show_bins.len()) {
-            print!("  {:+6.2}%", (other / value - 1.0) * 100.0);
+            print!("  {:+6.2}%", (other / values[0] - 1.0) * 100.0);
         }
 
         println!();
@@ -146,6 +169,8 @@ fn main() -> Result<(), Box<dyn Error>> {
             (@arg input: +required "Path of the input grid")
             (@arg pdfset: ... +required "LHAPDF id(s) or name(s) of the PDF set(s)")
             (@arg bins: -b --bins +takes_value "Selects a subset of bins")
+            (@arg scales: -s --scales default_value("7") possible_values(&["1", "3", "7", "9"])
+                "Set the number of scale variations")
         )
     )
     .get_matches();
@@ -173,8 +198,10 @@ fn main() -> Result<(), Box<dyn Error>> {
         let input = matches.value_of("input").unwrap();
         let pdfset: Vec<_> = matches.values_of("pdfset").unwrap().collect();
         let bins = parse_integer_list(matches.value_of("bins").unwrap_or(""))?;
+        let scales = str::parse::<usize>(matches.value_of("scales").unwrap())
+            .expect("Could not convert string to integer");
 
-        return convolute(input, pdfset.first().unwrap(), &pdfset[1..], &bins);
+        return convolute(input, pdfset.first().unwrap(), &pdfset[1..], &bins, scales);
     }
 
     Ok(())