@@ -25,15 +25,16 @@ std::vector<std::size_t> unravel_index(std::size_t flat_index, const std::vector
2525}
2626
2727std::vector<double > generate_fake_ekos (
28- std::vector<int > pids0 ,
29- std::vector<double > x0 ,
30- std::vector<int > pids1 ,
31- std::vector<double > x1
28+ std::vector<int > pids_in ,
29+ std::vector<double > x_in ,
30+ std::vector<int > pids_out ,
31+ std::vector<double > x_out
3232) {
33- std::size_t flat_len = x0 .size () * x1 .size () * pids0 .size () * pids1 .size ();
33+ std::size_t flat_len = x_out .size () * x_in .size () * pids_out .size () * pids_in .size ();
3434 std::vector<double > ops (flat_len);
3535
36- std::vector<std::size_t > shape = {pids0.size (), x0.size (), pids1.size (), x1.size ()};
36+ // NOTE: The EKO has to have as shape: (pids_in, x_in, pids_out, x_out)
37+ std::vector<std::size_t > shape = {pids_in.size (), x_in.size (), pids_out.size (), x_out.size ()};
3738 for (std::size_t i = 0 ; i != flat_len; i++) {
3839 std::vector<std::size_t > coords = unravel_index (i, shape);
3940
@@ -46,6 +47,35 @@ std::vector<double> generate_fake_ekos(
4647 return ops;
4748}
4849
50+ void print_results (std::vector<double > dxsec_grid, std::vector<double > dxsec_fktable) {
51+ const int idx_width = 6 ;
52+ const int num_width = 15 ;
53+ const int dif_width = 15 ;
54+
55+ // Print headers
56+ std::cout << std::setw (idx_width) << " Bin"
57+ << std::setw (num_width) << " Grid"
58+ << std::setw (num_width) << " FkTable"
59+ << std::setw (dif_width) << " reldiff" << std::endl;
60+
61+ // Print dashed lines
62+ std::cout << std::setw (idx_width) << std::string (idx_width - 2 , ' -' )
63+ << std::setw (num_width) << std::string (num_width - 2 , ' -' )
64+ << std::setw (num_width) << std::string (num_width - 2 , ' -' )
65+ << std::setw (dif_width) << std::string (dif_width - 2 , ' -' ) << std::endl;
66+
67+ // Print the data
68+ std::cout << std::scientific << std::setprecision (6 );
69+ for (size_t i = 0 ; i < dxsec_grid.size (); ++i) {
70+ double reldiff = (dxsec_fktable[i] - dxsec_grid[i]) / dxsec_grid[i];
71+ std::cout << std::setw (idx_width) << i
72+ << std::setw (num_width) << dxsec_grid[i]
73+ << std::setw (num_width) << dxsec_fktable[i]
74+ << std::setw (dif_width) << reldiff
75+ << std::endl;
76+ }
77+ }
78+
4979int main () {
5080 // TODO: How to get a Grid that can be evolved??
5181 std::string filename = " LHCB_WP_7TEV_opt.pineappl.lz4" ;
@@ -70,6 +100,7 @@ int main() {
70100
71101 // Get the PID basis representation
72102 pineappl_pid_basis pid_basis = pineappl_grid_pid_basis (grid);
103+ assert (pid_basis == PINEAPPL_PID_BASIS_PDG);
73104
74105 // Get the number of convolutions
75106 std::size_t n_convs = pineappl_grid_convolutions_len (grid);
@@ -86,19 +117,19 @@ int main() {
86117
87118 // Get the shape of the evolve info objects
88119 std::vector<std::size_t > evinfo_shape (5 );
89- std::vector<uint8_t > max_orders = {3 , 2 };
120+ std::vector<uint8_t > max_orders = {2 , 3 };
90121 pineappl_grid_evolve_info_shape (grid, max_orders.data (), evinfo_shape.data ());
91122
92123 // Get the values of the evolve info parameters. These contain, for example, the
93124 // information on the `x`-grid and `PID` used to interpolate the Grid.
94125 // NOTE: These are used to construct the Evolution Operator
95126 std::vector<double > fac1 (evinfo_shape[0 ]);
96127 std::vector<double > frg1 (evinfo_shape[1 ]);
97- std::vector<int > pids1 (evinfo_shape[2 ]);
98- std::vector<double > x1 (evinfo_shape[3 ]);
128+ std::vector<int > pids_in (evinfo_shape[2 ]);
129+ std::vector<double > x_in (evinfo_shape[3 ]);
99130 std::vector<double > ren1 (evinfo_shape[4 ]);
100131 pineappl_grid_evolve_info (grid, max_orders.data (), fac1.data (),
101- frg1.data (), pids1 .data (), x1 .data (), ren1.data ());
132+ frg1.data (), pids_in .data (), x_in .data (), ren1.data ());
102133
103134 // ------------------ Construct the Operator Info ------------------
104135 // The Operator Info is a vector with length `N_conv * N_Q2_slices` whose
@@ -117,12 +148,16 @@ int main() {
117148 }
118149
119150 // ------------------ Construct the Evolution Operator ------------------
151+ // Choose a different PID basis for the FK table
152+ // std::vector<int> pids_out = {-6, -5, -4, -3, -2, -1, 1, 2, 3, 4, 5, 6, 21, 22};
153+ std::vector<int > pids_out = pids_in;
154+
120155 // The Evolution Operator is a vector with length `N_conv * N_Q2_slices * Σ product(OP shape)`
121156 std::vector<double > op_slices;
122- std::size_t flat_len = x1 .size () * x1 .size () * pids1 .size () * pids1 .size ();
157+ std::size_t flat_len = x_in .size () * x_in .size () * pids_in .size () * pids_out .size ();
123158 for (std::size_t _i = 0 ; _i != conv_types.size (); _i++) {
124159 for (std::size_t j = 0 ; j != fac1.size (); j++) {
125- std::vector<double > eko = generate_fake_ekos (pids1, x1, pids1, x1 );
160+ std::vector<double > eko = generate_fake_ekos (pids_in, x_in, pids_out, x_in );
126161 for (std::size_t k = 0 ; k != flat_len; k++) {
127162 op_slices.push_back (eko[k]);
128163 }
@@ -137,11 +172,25 @@ int main() {
137172 }
138173
139174 std::vector<double > xi = {1.0 , 1.0 , 1.0 };
140- std::vector<std::size_t > tensor_shape = {pids1.size (), x1.size (), pids1.size (), x1.size ()};
141-
175+ // NOTE: The EKO has to have as shape: (pids_in, x_in, pids_out, x_out)
176+ std::vector<std::size_t > tensor_shape = {pids_in.size (), x_in.size (), pids_out.size (), x_in.size ()};
177+
178+ // NOTE: The arguments of `pineappl_grid_evolve` must follow the following orders:
179+ // - `grid`: PineAPPL Grid
180+ // - `op_info`: operator info
181+ // - `max_orders`: max orders to apply the evolution
182+ // - `operators`: evolution operator
183+ // - `x_in`: x-grid of the Grid
184+ // - `x_out`: x-grid of the FK table
185+ // - `pids_in`: PIDs basis representation of the Grid
186+ // - `pids_out`: PIDs basis representation of the FK table
187+ // - `eko_shape`: shape of the evolution operators
188+ // - `xi`: scale variation
189+ // - `ren1`: values of the renormalization scales
190+ // - `alphas_table`: values of alphas for each renormalization scales
142191 pineappl_fk_table* fktable = pineappl_grid_evolve (grid, opinfo_slices.data (),
143- max_orders.data (), op_slices.data (), x1 .data (),
144- x1 .data (), pids1 .data (), pids1 .data (),
192+ max_orders.data (), op_slices.data (), x_in .data (),
193+ x_in .data (), pids_in .data (), pids_out .data (),
145194 tensor_shape.data (), xi.data (), ren1.data (), alphas_table.data ());
146195
147196 // ------------------ Compare Grid & FK after convolution ------------------
@@ -161,32 +210,7 @@ int main() {
161210 nullptr , dxsec_fktable.data ());
162211
163212 // Print the results
164- const int idx_width = 6 ;
165- const int num_width = 15 ;
166- const int dif_width = 15 ;
167-
168- // Print headers
169- std::cout << std::setw (idx_width) << " Bin"
170- << std::setw (num_width) << " Grid"
171- << std::setw (num_width) << " FkTable"
172- << std::setw (dif_width) << " reldiff" << std::endl;
173-
174- // Print dashed lines
175- std::cout << std::setw (idx_width) << std::string (idx_width - 2 , ' -' )
176- << std::setw (num_width) << std::string (num_width - 2 , ' -' )
177- << std::setw (num_width) << std::string (num_width - 2 , ' -' )
178- << std::setw (dif_width) << std::string (dif_width - 2 , ' -' ) << std::endl;
179-
180- // Print the data
181- std::cout << std::scientific << std::setprecision (6 );
182- for (size_t i = 0 ; i < dxsec_grid.size (); ++i) {
183- double reldiff = (dxsec_fktable[i] - dxsec_grid[i]) / dxsec_grid[i];
184- std::cout << std::setw (idx_width) << i
185- << std::setw (num_width) << dxsec_grid[i]
186- << std::setw (num_width) << dxsec_fktable[i]
187- << std::setw (dif_width) << reldiff
188- << std::endl;
189- }
213+ print_results (dxsec_grid, dxsec_fktable);
190214
191215 pineappl_fktable_write (fktable, " evolved-grid.pineappl.lz4" );
192216
0 commit comments