Exploring Theory Space with Monte Carlo Reweighting

Theories of new physics often involve a large number of unknown parameters which need to be scanned. Additionally, a putative signal in a particular channel may be due to a variety of distinct models of new physics. This makes experimental attempts to constrain the parameter space of motivated new physics models with a high degree of generality quite challenging. We describe how the reweighting of events may allow this challenge to be met, as fully simulated Monte Carlo samples generated for arbitrary benchmark models can be effectively re-used. In particular, we suggest procedures that allow more efficient collaboration between theorists and experimentalists in exploring large theory parameter spaces in a rigorous way at the LHC.


Introduction
Monte Carlo (MC) simulation is an essential and ubiquitous tool in particle physics [1][2][3][4]. For realistic studies, such simulation generally includes the modeling of the detector response as well as the underlying physical process. While theoretical studies generally require only "fast" detector simulation, such as that performed by PGS [5] or Delphes [6], the modeling of the detector response to a given physics event at the experimental level, using frameworks such as GEANT [7], is necessarily a time-consuming and computationally expensive process.
In addition to the situations above which involve many continuous parameters, there are also issues with discrete parameters; in particular there can be many choices for the spin assignments for particles in a given event topology. A well-known example of this is the fact that most if not all SUSY topologies are also produced in models of Universal Extra Dimensions (UED) [33,34]. It is therefore desirable to study, at the very least, the UED counterparts of SUSY signal benchmark points. 1 Clearly experimentalists at the CERN Large Hadron Collider (LHC) and similar experiments are faced with an overwhelming multitude of theoretical models which could be searched for. To determine the properties of a given signal in the detector generally requires detector-level MC simulation (fullsim) of that signal, e.g., with GEANT, which is a time-consuming process, as noted above. To repeat this procedure for the myriad possibilities suggested by theorists is simply not possible.
Attempts have been made to solve this problem [35][36][37][38], perhaps most notably by encouraging theorists to test the predictions of their models against suitably chosen experimental results (e.g., limits on cross-sections times branching fractions as a function of the mass spectrum). However, an actual discovery will still require the use of fullsim samples to obtain realistic distributions for the kinematic variables of interest (such as test statistics).
We propose and describe the following procedure for the experimental analysis of a large number of signal hypotheses: 1. Generate a large set of unweighted sample events at generator level, {G i } for a specific model and point in that model's parameter space, which we term (A, α) (where A describes the model and α specifies a point in that model's parameter space); a reasonable choice for the model, A, would be a simplified model [39][40][41][42][43][44][45].
2. Pass the generator-level events, {G i } through realistic detector simulation. Select desired events ("apply cuts") to obtain a set of detector level events {D i }. (Every generator-level event G i maps to some detector-level event, D i , but not every detectorlevel event passes our cuts.) 3. Note that the cross section after acceptances and efficiencies, as well as the detectorlevel predictions for all distributions in signal hypothesis (A, α) can be obtained from the resulting event sample. 4. Choose some other model, B, and parameter space point, β. Assign each detectorlevel event D i a weight w(B, β, G i )/w(A, α, G i ), which depends only on generatorlevel information and uses, in particular, the "truth" values of particles such as neutrinos or neutralinos which are not observed in the detector.
Because the reweighting involves only generator-level information, it is very fast compared with detector simulation. As we will show in the following section, where we derive and justify the method, this is not an approximation-in the limit of infinite MC samples, this method reproduces all distributions and quantities of interest exactly. Therefore, a relatively efficient method for examining the large "theory" space would be for experimentalists to provide theorists with the hard-process/ parton-level events which they have subjected to detector simulation. The theorists can then calculate the weights of these events under the original model/ parameter space point as well as the weights of these events for any desired point in model/ parameter space. The experimentalists can take these weights and obtain any desired distribution at detector level from reweighting their detector-level events. (Of course the "theorist" in this description could also be an experimentalist interested in the model being studied.) While, as we will discuss below, practical issues may arise for finite MC samples, in general this method could vastly extend the reach into theory space of LHC analyses.
We reiterate that the MC events being shared here are "truth" events. The momenta of all "invisible" final state particles (such as neutrinos or neutralinos) are fully specified. Therefore the calculation of weights in different scenarios is numerically trivial; the timeconsuming integration over undetermined momenta which is necessary in the Matrix Element Method [46][47][48][49][50][51][52][53][54], and which can be performed by tools such as MadWeight [55,56] is not necessary.
In Section 2, we will provide a brief overview of MC simulation in particle physics. It is not our intent to provide a review (see, e.g., Refs. [2,4]). We wish only to remind the reader that a set of unweighted hard-process parton-level events x i,true can be used to generate a set of unweighted events after showering, hadronization, detector simulation, object reconstruction, etc., which are given by The important points about the generalized transfer function T , which incorporates the vast amount of physics briefly summarized in Section 2, are 1. It does not depend on the new physics model or parameters.
2. It is unitary. Every hard-process event simulated corresponds to a detector-simulated event with objects reconstructed, even if, e.g., due to the failure to reconstruct the physics objects corresponding to some partons, it will not be included in our calculations of distributions for the relevant final state.
3. As T is unitary, if we take a set of unweighted events x i,true and simulate each event, x i,true , or, in equivalent language, choose a random value of x objects according to the probability distribution function (PDF), T (x i,objects , x i,true ) for each x i,true , then we will obtain an unweighted set of detector-simulated and object-reconstructed events {x i,objects }.
If the reader is willing to accept these facts about T , then Section 2 may be safely skipped.
Reweighting itself is described in Section 3. To demonstrate its use, we then provide a few examples of reweighting in action. In Section 4.1 we study a signal model with a multi-dimensional parameter space, while in Section 4.2 we use reweighting to analyze angular correlations needed to determine the spins of the new particles. In Section 4.3 we demonstrate that reweighting works even in the presence of showering, hadronization, jet clustering, and detector simulation. We present some useful notes for practitioners of our method in Section 5. Section 6 is reserved for our conclusions.

Monte Carlo Basics
The integral F of an arbitrary function f (x) over some interval (x − , x + ) is: where r i is a uniform random variate in the range (0,1), and take the average of f (x) at these points: The points x i used to evaluate this sum can be viewed as a set of "events" with weights f (x i ).
To compute the MC average of a different function, g(x), we similarly have: Importance sampling considers a change in the measure of sampling in order to reduce the variance of the MC estimate: In many applications, g(x) can be integrated analytically and has an antiderivative G with a inverse G −1 . Then x i = G −1 (r i + G(x − )). From (2.2) and (2.4) it follows that: where we emphasize that the average is with respect to the function g. It is not necessary to have a closed form for g, G or G −1 , and this will be the case for the practical application of these formulas discussed later. In fact, the selection of integration points from g(x) may be performed numerically using, for example, VEGAS [57]. As a trivial example, assume a set of N unweighted events generated according to g(x), so that each event has a weight σ/N . Then, choosing f (x) = g(x), we have: Even though the weights are uniform in this case, the values of x i associated with each event are distributed according to g(x) and can be used to calculate averages or construct histograms of other quantities based on x.

Parton-Level Event Generation
In the applications discussed here, we are interested in performing a reweighting of events at the parton-level. The parton-level is usually the first non-trivial level of approximation where physics beyond the Standard Model is needed.
In the examples considered here, we focus on hadron-hadron collisions. The basic quantities of interest are the kinematic properties of events involving some new particles or interactions. We limit ourselves to processes of the type 1 + 2 → i + j + · · · + n. Then, the partonic cross section is proportional to where the amplitude for this process is A. Imagine now, another process with the exact same initial and final state, but described by a different amplitude A . The relative probability of this process with respect to the former is given simply by |A | 2 /|A| 2 . This parton-level description is often called the "hard-process" and describes the physics occurring at energies Q Λ QCD , or equivalently length/ time scales 1/Λ QCD . The expression for the partonic cross section can be cast in the form of a probability distribution that generates parton-level "events" or configurations of kinematic variables x i . (In this heuristic description, we can think of e.g. the helicity and color of all final state partons as specified; in general, we need to generate all relevant final state helicities and color structures.) These events might be the output of a program such as MadGraph5 [58], CalcHEP [59], or CompHEP [60,61]; these programs can also provide code that can be used later in reweighting.
Since our interest is in the reweighting of new physics signal events, the hard-processes we consider will involve some new physics model, A, described by some parameters α. We will term this "theory point" (A, α). Our goal is to obtain all cross sections and distributions for a different model B, for the parameter point, β. Of course, it may be that we are scanning the parameter space of model A, in which case B = A and we are considering parameter points α 1 , α 2 , .... Also, for some processes, like Higgs boson production and decay, it may make sense for (A, α) to be the SM rather than a new physics model; we will then reweight to obtain distributions for the non-SM theory point (B, β).

Particle-Level Event Generation
For now we focus on perhaps the primary use of Monte Carlo event generators, which is to produce a sample of unweighted events of particles to be input to a realistic detector simulation. Generating events is important, because we are not interested only in the cross section, subject to arbitrary event selection ("cuts"), but also in the kinematic distribution, ρ(V ), for a kinematic variable of interest, V (x true ), subject to certain cuts: Here D G is the space of events (x i ) which pass cuts. These are not the final, detector-level cuts. However, cuts are often applied at this stage, either to provide a cut-off for infrared divergent quantities or simply to reduce the time spent on detector-simulation of events that will not pass the triggers or the final, detector-level cuts. If we generate a sample of N unweighted events {x i } in the allowed region, and, in the course of generating the events determine the total cross section for the specified process to be σ G,MC , then, of course, we can approximate the cross section after cuts by σ G ≈ σ G,MC , while the normalized distribution of an arbitrary kinematic variable, ρ (as in Eq. (2.9)) can be approximated by forming the histogram The generalization of Eq. (2.9) and Eq. (2.10) to multidimensional observables is obvious. The cross section in Eq. (2.8) and the distributions ρ(V ) in Eq. (2.9) and Eq. (2.10) have been calculated considering only the physics from the hard process. To obtain accurate distributions for quantities measured in the detector, we must utilize general purpose MC generators and/or showering tools like ARIADNE [62], HERWIG [63], HERWIG++ [64], PYTHIA 6 [65], PYTHIA 8 [66], SHERPA [67], or VINCIA [68][69][70]. With these tools we must simulate initial and final state radiation (ISR and FSR), decay resonances with sufficiently short life times, and hadronize colored objects. While modeling the hadronic physics we must also consider, e.g., the physics of the underlying event [71][72][73][74][75].
Formally, we can think of events at each level of the simulation as living in different spaces. Thus the initial, unweighted, hard-process parton-level event undergoes a transformation: (2.11) We note that the dimensions of each x i vector are in general quite different. Showering adds particles to the event via ISR and FSR, and decays replace resonances with two or more daughter particles. Also, colored partons, either before or after showering and decays, are obviously not in one-to-one correspondence with final state hadrons. The mappings S and H are unitary (provided we do not impose additional cuts at this stage). In the case of these QCD processes, this ability to separate physics at different length scales (and to consider the composition of probabilities rather than amplitudes) results both from factorization and from more specific results like the KLN theorem [76,77] (see e.g., Refs. [2] and [78] and references therein for more discussion of these points). Thus, for a particular event x i,true , (2.12) Likewise, for some particular showered event x i, showered we would have, (2.13) Obviously, this implies that It is clear that if we start with a set of unweighted hard-process generator-level events {x i,true }, and then, for each x i,true choose a random value for x i,showered according to the PDF S(x showered , x i,true ) corresponding to x i,true , then we will obtain a set of unweighted, showered events {x i,showered }. We can obviously continue this procedure by choosing, x i,hadron according the the PDFs H(x hadron , x i,showered ), thereby obtaining a set of unweighted, hadron-level events {x i,hadron }. That we perform this random selection according to the PDFs S(x showered , x i,true ) and H(x hadron , x i,showered ) using a simulation program such as HERWIG or PYTHIA does not affect the argument. We remind the reader that only the generation of the initial set of hard-process parton-level events {x i,true } had anything to do with the new physics theory point (A, α). The subsequent showering and hadronization depends only on the hard process event x true (by which we include both the four momenta of the particles and the particular color structure which was generated) and SM parameters, such as α S (Q).

Detector-level Monte Carlo Simulation and Reconstruction
In this subsection we briefly describe the process of performing detector-level MC simulation, in order to explain why reweighting detector-level events based on the hard-process matrix elements works. Two further stages of event evolution, detector simulation and object reconstruction, also decompose into the product of PDFs: Detector simulation refers to determining the tracks produced and energy deposited in various parts of the detector, using GEANT. So in general x i,detector does not describe the four-momenta of particles, but the properties of tracks, energy in calorimeter cells or towers, etc. In order to do physics, we must take this raw detector output and map it into "physics objects", such as "electrons", "muons", or "jets". The reconstruction of jets, of course, is especially non-trivial; significant work has gone into the theoretical and experimental understanding of jet algorithms [79,80]. For our purposes, it suffices to note that the standard simulation procedures are formally analogous to the showering and hadronization described in that they can be viewed as obtaining a set of unweighted detector-level events x i,detector corresponding to hadron-level event x i,hadron , by choosing a random value of x i,detector from the PDF D(x detector , x i,hadron ) and obtaining an event consisting of reconstructed physics objects, x i,objects from this detector event, x i,detector using the PDF R(x objects , x i,detector ).
In fact, we can compose the four PDFs, S, H, D, and R into the generalized-transfer function T from Eq. 1.1, which, as a product of normalized PDFs (as each step is unitary), is also a normalized PDF. While we did not discuss MC tuning [81][82][83][84], or pileup effects explicitly, we note that these effects contribute only to T (x objects , x true ) and neither depend on the new physics parameters (A, α) nor change the interpretation of T (x objects , x true ) as a normalized probability distribution. We note also that procedures like CKKW [85,86] or MLM [87] matching, which attempt to combine aspects of the parton shower with event generation from matrix elements do not change the overall thrust of our argument. At most, we must modify our reweighting procedure to use events x i,matched and generalized transfer function T (x i,objects , x i,matched ) in Eq. (1.1); these modifications, and the associated modifications to the reweighting procedure for these events, are relatively straightforward.

Reweighting
Given a set of unweighted hard-process events, x i, true , generated with, e.g., MadGraph, and following the facts about the generalized transfer function T from Eq. (1.1) listed in Section 1 and/or the discussion in Section 2, we can produce a set of unweighted events x i,objects through MC simulation that correspond to the hard-process events x i,true .
These events can be used to obtain detector-level cross sections and distributions, in a manner exactly analogous to that specified above (c.f. Eq. (2.10)); namely if N detector of the N generator events which we simulate satisfy our cuts on our detector-simulated reconstructed physics objects, then our cross section (times efficiencies and acceptances) is given by The content of the bin in the histogram for an arbitrary kinematic variable V containing values in the range [V min , V max ] is given by Of course, these detector-level events correspond to the generator-level events {x i,true }, which were generated using the new physics theory point (A, α). This raises the question: what if we now wish to obtain detector-level cross section and observables for a different new physics model theory point, which we label (B, β)? Obviously we could simply perform the procedure described above on a new set of unweighted, parton-level events generated for (B, β). This approach, however, can become impractical very quickly, as fullsim detector simulation is very slow (rates on the order of minutes per event are typical) relative to the rest of the simulation process.
Instead, we note that the Eqs. (3.1) and (3.2) can be thought of as the MC evaluation of cross sections or histograms using the differential cross section for x objects , i.e.
If we replace the theory point (A, α) with (B, β), we have instead When we, e.g., perform the integral in Eq. (3.2) by generating unweighted events for the hard-process and then simulating, we are really replacing some region of x true space, with volume V i , with the event x i,true . The weight of this region in the calculation of quantities, . (These weights are equal for unweighted events, up to statistical fluctuations. V i is small when the differential cross section is large; equivalently the number of events selected in a region of phase space is proportional, in the limit of large statistics, to the differential cross section in that region of phase space.) When we then perform an integral like that in Eq. (3.2) using MC methods, we therefore must replace the weights: At the level of MC events, this replacement corresponds to weighting each event by the ratio of hard-process differential cross sections (3.5). With this procedure we obtain, for cross sections and for histograms

Applications
In this section, we demonstrate the use of reweighting methods through three example analyses. We also discuss the bin-by-bin errors in histograms obtained by reweighting. In Subsection 4.1 we will show how to reweight events to study different coupling structures for Higgs to four-lepton events. This will be followed, in Subsection 4.2, by an investigation of reweighting to obtain distributions for a UED model from MC samples generated for a SUSY model. Finally, in Subsection 4.3 we will study the effects of showering, hadronization, (fast) detector simulation, and jet reconstruction. Especially for the benefit of readers who may have skipped Section 2, in Fig. 1 we present a cartoon which describes the reweighting procedure. Essentially we can obtain distributions in some physics model J by reweighting MC events generated for a model I. The reweighting factor (3.5) is given by the ratio of differential cross sections evaluated for the "truth"-level parton-level MC event in each model. In the limit where model I and model J are produced from the same initial partons and have the same masses for final state particles, the parton distribution functions and phase space factors in this ratio cancel, and we are left with the ratio of squared matrix elements, as described in the caption to Fig. 1. The different colorations of the cartoons representing model I and model J in Fig. 1 represent the different distributions of the squared matrix elements over the common phase space in the two models.
In the examples provided, the factors used for reweighting were obtained using the "standalone" matrix element calculating code which can be generated automatically, in either Fortran or C++, from MadGraph5 [58]. This code requires the external parton momenta, which should be provided by a short, user-supplied, "wrapper" code. One can also generate similar standalone code in C or Mathematica using CalcHEP or CompHEP [59][60][61]. We encourage the authors of these and other similar automatic matrix element calculators to further increase the user-friendliness of their tools for the purpose of reweighting.
Before proceeding, we pause to present a caveat about the sorts of models one can study using reweighting. Obviously we can only use a reweighting procedure if the event being reweighted is possible in both models. So the final state particles must have the same masses in both models. Additionally, when one or both models contain intermediate resonances, it is important that those resonances, if sufficiently narrow, have the same masses in each model. Otherwise, the resulting extreme differences in the density of events in the phase space of the two models can lead to undersampling of important regions of phase space. However, this is not, in principle, an insurmountable difficulty, and several practical approaches for dealing with undersampling will be discussed in Section 5.

Changing the Coupling Structure: Higgs to Four Leptons
The "golden" Higgs to four-lepton channel (see Fig. 2) is quite sensitive, both for the discovery and for the subsequent measurement of the spin and parity properties of the Higgs boson , and therefore plays an important role in the experimental study of the Higgs at the LHC [145][146][147][148][149][150]. Assuming that the putative Higgs is spin-zero, we can then characterize its couplings to Z bosons, following [143], using the Lagrangian Therefore, in general, the HZZ couplings involve five parameters (up to lowest non-trivial dimension). Following the "Geolocating" approach [129], we can remove the degree of freedom associated with the partial width, leaving us with a "sphere" of four-dimensions. To explore this space fully, we may need to employ reweighting procedures, such as those described here.
As an example of such analyses, in Fig. 3 we study the distribution of the smaller dilepton invariant mass 2 M Z 2 = min (M e + e − , M µ + µ − ) in the 2e2µ final state. In each of the three panels in Fig. 3, the solid blue histogram is always the normalized M Z 2 distribution for the pure κ 5 case, i.e. when κ 5 = 0 and all other κ i (i = 1, . . . , 4) are zero. On the other hand, the solid magenta histograms correspond to the case when only one of the other κ i couplings is turned on: κ 1 = 0 in panel (a), κ 2 = 0 in panel (b), and κ 3 = 0 in panel (c). Our reweighting procedure now allows us to obtain the shape of the solid blue histogram by reweighting the corresponding magenta plot. The results are shown by the dotted red lines in Fig. 3. By comparing the solid blue and dotted red distributions, we observe a very good match, which validates our procedure.   1, 0, 0, 0), and c) κ i = (0, 0, 1, 0, 0), while the solid blue distribution is obtained by directly generating events for a model with κ i = (0, 0, 0, 0, 1). Alternatively, we can obtain the same κ 5 = 0 distribution, by recycling the events used for the solid magenta histograms, and reweighting them with the corresponding ratio of matrix elements squared. The results from such reweighting are shown by the dotted red distributions.
However, it is not sufficient to show that the central bin values for each histogram are reproduced by the reweighting procedure. To perform statistical analyses using the reweighted histograms, we need to understand the error (or uncertainty) on the histogram bin values. Following Ref. [151], if there are N events in a bin, each of which has weight w k , then the value of that bin in the reweighted histogram is The error on T is given by Clearly, in the special case of an unweighted histogram, δ yields √ N which reproduces the well-known expression from Poisson statistics. Combining Eqs. (4.2) and (4.3), and writing w 2 in terms of the variance on the weights, we find the fractional error on the bin value in a weighted histogram to be given by For the case of interest to us, the w k are the reweighting factors R(A, α, B, β) used to reweight events generated for model A, in order to obtain a histogram for model B. (Note that the value of N is the same for the two models.) Since in general the reweighting factor varies from event to event, the variance (4.4) is nonzero. Thus, when we reweight, the statistical error increases, and Eq. (4.5) quantifies this effect.
To provide a concrete demonstration of these effects, in Fig. 4(a-c) we show the fractional error, as calculated using Eq. (4.5), on each bin of the corresponding histograms    Fig. 3(a-c). In Fig. 4, solid magenta lines denote the fractional error (which scales as 1/ √ N ) in the original model with unweighted events, while the red dotted lines show the error on the corresponding histogram obtained from the same events using reweighting. In accordance with Eq. (4.5), the error is in general greater for the reweighted histogram, but not always -e.g., in Fig. 4(a), the errors before and after reweighting are essentially the same. To better illustrate the origin of the errors after reweighting, in Fig. 5 we provide two dimensional temperature plots of the reweighting factor, R, and the corresponding value of M Z2 for each event. We see a clear correlation between the spread in the values of R and the magnitude of the increase in the fractional errors in Fig. 4. At this point, one might wonder how the errors on histograms which were reweighted from an initial model A to another model B compare to the errors on histograms built from unweighted events generated directly in model B. Since the error on the bin values in the unweighted histogram is 1/ √ T , from Eqs. (4.2) and (4.5) we get where we assume sufficiently large statistics, so that the number of unweighted events is close to the expected value T . We see that the second factor in Eq. (4.6) always leads to an increase in the error when reweighting, in agreement with Fig. 4. On the other hand, the first factor, w , can be either larger or smaller than 1, depending on the relative weights between models A and B. Using Eqs. (4.5) and (4.6), one can quantify the error from reweighting in the region of interest. If that error turns out to be unacceptably large for the task at hand, one can apply the procedures discussed below in Sec. 5 to further reduce those errors.
In Fig. 3 we presented distributions for the M Z 2 variable, which we studied further in Figs. 4 and 5. However, we wish to emphasize that the variable whose distribution is obtained via reweighting could also be an optimized multivariate analysis-based variable such as the MELA KD [103,145], MEKD [124], or the output of a boosted decision tree or neural network analysis. Hence reweighting can be used in tandem with sophisticated and powerful multivariate analysis methods [152].

Changing Spin Assignments: The Antler Topology
Extensive simplified model searches [39][40][41] have been performed by the ATLAS and CMS collaborations in various channels [42][43][44][45]. When performing a simplified model analysis, one generally fixes the spins of the particles involved in the simplified event topology. One common choice is to use SUSY spin assignments, though another possibility is to decay all particles by pure phase space. It is important to be able to interpret the results from these searches in the context of other models, with different spin assignments. A particularly well motivated example of an alternative spin assignment is provided by the minimal UED model (MUED) [153].
To perform reinterpretations for different spin assignments, we need to recalculate cross sections, branching ratios, and efficiencies for the given study point in the new model. We can obtain cross sections and branching ratios for the new model point via theoretical formulae, but efficiencies generally must be obtained through MC simulation. We emphasize that we can recycle existing MC samples to determine these efficiencies for different new physics models in an efficient manner. In this subsection, we demonstrate the utilization of MC event samples generated for a SUSY model to perform a collider analysis of the MUED model, in the context of the "antler" topology [154,156] depicted in Fig. 6. In particular  Figure 7. Distributions of various parton-level quantities which might be used to distinguish SUSY from MUED [154]. The SUSY parameter point used is the CMS LM6 study point [155]; in the UED scenario the Kaluza-Klein partners have the same mass as the corresponding superpartners.
• M − + : The invariant mass of the two leptons.
• √ŝ min : An estimator of the mass scale of the hard scattering, proposed in [164,165], calculated for the antler event topology assuming a particular value for theχ 0 1 mass. Since in practice we would scan over this "trial" mass, for definiteness we use the correct or "truth" value of this quantity in the analyses presented here. The definition of √ŝ min is √ŝ (4.10) • √ŝ : The actual partonic center-of-mass energy.
Distributions of the above observables are shown in Fig. 7. The MC samples used to generate these distributions consist of 100, 000 parton level events at the 14 TeV LHC. Fig. 7 suggests that MUED analyses can in fact be performed by reweighting SUSY MC events, allowing the efficient probing of non-standard spin assignments at the LHC. One could extend this technique beyond MUED to all other possible spin assignments for a given event topology with given masses. In particular, one could set robust and model independent limits on new physics in the following way. Among all possible spin assignments (or more generally, parameter values), identify the case which is most difficult to test experimentally, and use it to present limits from experimental searches. The conservative bounds set in this way will be valid for all other spin assignments as well.

Jet Clustering, Event Selection, and Detector Response
Parton level information can be distorted by various factors, including showering, hadronization, jet clustering, and detector response. All of these can "deform" the phase space of generated particles. But, as noted in Section 2, the relative weights in various models for a given observed event depend only on the parton level amplitudes for the event in each model. We present a conceptual picture of this reweighting procedure in Fig. 8.
To illustrate this point numerically, we simulate tt production at 14 TeV LHC, where the top-quark (anti-top quark) decays into a W + (W − ) boson or charged Higgs φ + (φ − ), see Fig. 9. Our goal is to evaluate the performance of the reweighting procedure described in Fig. 8, including effects that we did not account for in the parton-level examples described in Sections 4.1 and 4.2. Specifically, we generated 100, 000 parton level events with MadGraph5(aMC) v.2.1.0, passed the events to PYTHIA 6.4, and performed detector simulation with Delphes v.3.0.12. For jet clustering, we used the anti-k t algorithm with ∆R = 0.5. We denote actual event generations by solid lines and mark calculations for the reweighting procedure with dashed lines. The weight of an event in the given model is indicated by the shade of gray. The effects of showering, hadronization, detector simulation, and jet clustering are shown conceptually by the reshaping of the events which are rectangular at generator level and oddly shaped after detector simulation. Clearly the detector level events are the same whether or not we are reweighting (hence have the same shapes in our figure). Also, detector level events either pass cuts or fail cuts regardless of whether reweighting is employed. Figure 9. To demonstrate the use of reweighting events involving detector simulation and jet clustering, we study top quark pair production, followed by decays involving a b-jet and either (a) the W ± boson, as in the SM, or (b) a charged Higgs, φ ± . As we show, the charged Higgs scenario (b) can be analyzed using SM dilepton tt events (a) through reweighting. Due to the mass dependence of Higgs couplings, the lepton from charged Higgs decay is always a muon and never an electron.
For the reconstruction of visible particles, we employ default smearing factors and lepton tagging efficiencies. We apply basic analysis cuts with p T (e) > 10 GeV, p T (µ) > 5 GeV and p T (j) > 30 GeV and do not require b-tagging in our analysis.
When the top quark decays into a bottom quark and a charged Higgs boson, the  Figure 10. Invariant mass distributions, m bµ , of b-jets and the muon for events generated using the topologies shown in Fig. 9. The variable m bµ contains information about the spin of the intermediate particle, W ± or φ ± . In (a) we plot this quantity at parton level, pairing the muon with the b from the same top decay. In (b), we use events which have undergone realistic simulation, i.e. they are showered, hadronized, and passed through fast detector simulation and jet clustering. To separate the effects of detector simulation, etc. from the effects of combinatorics, we use truth information to select the correct jet to pair with the lepton. In (c), which was also generated using realistic simulation, we include the effects of combinatorics by not using truth information (to model the actual experimental situation) and instead using the mixed event subtraction method [166] to suppress contributions from incorrect parings of jets and muons.
invariant mass distribution of the bottom quark and the corresponding lepton from the charged Higgs boson decay will have a triangular shape (the solid blue line in Fig. 10(a)). 3 On the other hand, the invariant mass distribution of the b-quark and the corresponding lepton from a W ± decay will have a non-trivial shape (solid magenta lines in Fig. 10(a)). These two cases can be easily related by reweighting at the parton level, as shown by the red dotted histogram in Fig. 10(a).
At the detector level, the b-quark is reconstructed as a jet, thus we select the two highest p T jets to be our b jets. In order to study the effects of combinatorics separately, in Fig. 10(b) we show an intermediate (and unphysical) case where we use MC truth information to select the correct jet to be paired with the muon. We see that the reweighting procedure correctly reproduces the distributions of kinematic variables computed in terms of reconstructed objects.
Finally, in Fig. 10(c) we show the same three distributions, only this time ignoring the MC truth information and fully accounting for the combinatorics. In order to reduce the contamination from combinatorial pairing errors, we use the mixed event subtraction method [166], where we add both possible pairings in a given event, then subtract a pairing of the muon with a jet from a different event. We emphasize that the point of these analyses is that a distribution made through reweighting is identical to the actual distribution made with a MC events generated for the hypothesis in question.

Possible extensions
As noted in Section 4 above, reweighting events generated using model I to obtain distributions in model J may become problematic for two main reasons. Firstly, in a phase space region where the cross section for the new model J is high, while the cross-section for the reference model I is low, the region is sampled by relatively few MC events. Secondly, the spread in the reweighting factors for events in the given bin may become large. Both of these effects are described by Eq. 4.5. We are obviously interested in minimizing the errors on histograms, but at the same time, we want to keep as low as possible the number of events which need to be processed by full MC simulation.
We therefore suggest three practical procedures that can be used to mitigate issues related both to undersampling and to the spread on reweighting factors, and possibly to help understand systematic effects.
1. Targeted Generation: In this procedure, whenever the error in a reweighted histogram obtained for model J becomes greater than a certain threshold, one generates additional new events in the undersampled region that are phase space "neighbors" of the events already in the bin. The simplest way to think of this is to divide phase space P into two non-overlapping regions, N and N , where N describes a phase space region where additional event generation is desired, and N is the rest of the phase space. Then we generate events for model J only in region N . Distributions for J will then be obtained by reweighting the events generated for model I in the N and for events generated for model J in the region N . 4 If one then wishes to obtain distributions for a third model, K, one will reweight events generated for model I in region N and the events generated for model J in region N , in each case using a reweighting factor involving the appropriate models.

Regions of Overlap:
A related approach is to generate large samples of unweighted events at several benchmarks in parameter space, which we label I 1 , I 2 , etc. One then obtains distributions D 1 , D 2 , etc. for model J by reweighting events generated for model I 1 , I 2 , etc. 5 If the benchmarks are chosen appropriately, every important region in phase space for model J can be described with low errors by reweighted events generated for some model I j , hence avoiding the undersampling issues described above. Also, the comparison of the values of the distribution obtained from different choices of unweighted events may allow for the calculation of systematic errors from this procedure.
3. Optimized Sample: A third, and quite distinct approach to optimizing reweighting procedures for large parameter space is to optimize the choice of model I to avoid undersampling in as much of the parameter space as possible. In general, this would be performed by maximizing an integral over the parameter space (possibly weighted 4 The two samples must be mixed accounting for the relative number of events generated in regions N and N . This is an example of multi-channel integration with channels I and J. 5 This is another example of multi-channel integration. by some prior or importance measure) of a function representing the errors that will be obtained on bins of the desired histograms. This appears to be a hard integral to perform in practice, so reasonable approaches will probably involve approximation to some degree. It is also unclear whether even an optimized model would be as efficient as generating events for at least the several benchmarks in parameter space described in the previous item.

Conclusions
We have presented and described in detail a procedure aimed at surmounting the experimental challenges presented by the multiplicity of theoretical models, which may each, in turn, have large parameter spaces. We realize that the procedure presented here is not totally unknown to experimentalists (particularly in the special case of parton distribution function reweighting [167]). However, given the importance of this challenge, we felt it important to highlight the potential of this method, especially for the study of large BSM parameter spaces.
Specifically, we have shown how reweighting events using "truth" information from generator-level Monte Carlo events, including the momenta of invisible particles, makes possible the detailed study of large signal parameter spaces and/or large numbers of signal models at the level of detail needed for experimental purposes. We demonstrated this in several motivated physics examples and also illustrated potential issues which can arise, including an explicit discussion of the errors on the weighted histograms generated from reweighting.
There are several advantages of our method: • Given a large existing sample for some theory model I, we gain speed in generating an effective fullsim sample for another theory model J by avoiding detector simulation, hadronization, showering and fragmentation.
• All existing fullsim samples can be resurrected for the study of new models.
• The reweighting is computationally very quick and simple -for example, there is no need to integrate over the momenta of any invisible particles.
• Reweighting also allows the interactions between theorists and experimentalists to be maximally efficient -the experimentalists are handling the fullsim MC generation, while the theorists are providing only the parton-level functions for reweighting the events.
We conclude with the simple recommendation that in order to be able to implement the procedure described here, when generating fullsim samples, one should make sure to record and store the generator-level information which is needed for reweighting.