We thank the referee for their detailed reading of the manuscript and constructive comments.

We thank the referee for their detailed reading of our manuscript and constructive comments. We have made amendments and additions within the manuscript to address the requested changes.

Weaknesses

1. We intentionally used the 1D example to introduce the MAGIC method, in order to simplify the notation and the discussion. We have also performed limited tests with multi-dimensional microscopic models and only encountered the expected issue of not being able to solve the inverse problem when there was not enough macroscopic information given. We agree with the referee that further work on the MAGIC method is warranted, and we plan to perform more dedicated tests when this will be used in more complex problems in future works. To make the reader aware of potential issues we have expanded the relevant discussion on p. 9.
2. We thank the referee for this comment. We have added the rough order of magnitude speed-up estimates for the case of our toy simulation on p.8.
3. We agree with the reviewer that a rigorous way to assess the systematic uncertainties related to the (in)correctness of the hadronization models is lacking. We did, however, explore the issues of under- and over-fitting in appendix B.5 and we have added a more detailed discussion in the main text. We believe that implementing a more rigorous set of closure tests to assess model correctness lies beyond the scope of this work, given that the model was already shown to converge properly. However, this is an important topic and deserves careful study in the future.
4. We thank the referee for pointing out the incomplete discussion. We have expanded our discussion in the main text on how such an uncertainty can be assigned by performing a model comparison between different architectures much like what was done in appendix B.5.

Requested changes

1. The referee is correct in stating that ref. [10] has similar goals of extracting microscopic information about the hadronization models using only measurable quantities. We have added two clarifying sentences in the introduction. It would also be interesting to perform a more detailed comparison between the two approaches (once they are developed enough for the comparison to be meaningful). We have mentioned this possibility on p. 9 in the main text.
2. We thank the reviewer for the suggestion and have rearranged the text accordingly.
3. We have added additional details on the multi-dimensional scenario in the second to last paragraph of section 4. We have tested and expect MAGIC training to generalize to correlated multi-dimensional distributions assuming the macroscopic distributions used during training are sensitive to the correlations exhibited by the base NF. In general, we expect multi-dimensional fine-tunings with MAGIC training to require more oversight during training to ensure proper coverage during model exploration and reweighting.
4. We have included a footnote in the second to last paragraph of the introduction and changed the text in the sixth paragraph of section 4 to clarify that the comparison and gain in computation time is meant to refer to the savings between re-simulating events versus re-weighting events, the former taking on the order of minutes per thousands of events versus seconds per thousands of events. While the NF architecture is able to generate events more quickly than the model introduced in our previous work [8], this is not the main focus of the present study and both models are still significantly slower than modern event generators. We leave attempts at producing generative-ML-based models of hadronization with sampling speeds comparable to modern event generators for future work.
5. We have added to the main text a brief discussion regarding under- and over-fitting of the model, and model correctness in general. We reference explicit techniques for model comparison in a Bayesian framework [29, 30] and the architecture scan performed in appendix B.5. In particular, we also discuss how the posterior uncertainties reflect a lack of significant under-fitting. Of course, the model is not a perfect approximation of the experimental distribution but it is not so biased towards a wrong model in such a way that the uncertainties are artificially low. A more thorough closure test can quantify the degree of compatibility between model and data. However, we believe such an analysis is not warranted in this case given the performance of the model.
6. Figure 3 shows the learned PDF as defined in eq. (20) of appendix B.1.1 obtained by feeding in an ordered list of points over the domain of the distribution. By construction, the 1D-NF is automatically normalized to unity.
7. Figure 4 is a scatter plot illustrating the relationship between the counts of bins from the left plot and the relative error in the corresponding bin. It compares the (blue) BNF estimate of the relative uncertainty, $\sigma_{\text{bin}} / \langle N_{\text{bin}} \rangle$, and the (red) relative generated uncertainty, $\sigma_{\text{gen}} / \langle N_{\text{bin}} \rangle$, for each bin. We have clarified this in the figure caption.