NLO Merging

The two NLO merging methods implemented in Pythia are called NL³ (for Nils Lavesson + Leif Lönnblad) and UNLOPS (for unitarised NLO+PS merging). Both of these schemes require Les Houches Event File input that is generated by tree-level or NLO matrix element generators. Currently, Pythia requires NLO input generated within the POWHEG framework. The generation of sensible input will be discussed below. The two NLO merging methods are illustrated in the sample main programs main87.cc (introducing NL³ ) and main88.cc (introducing UNLOPS). Before describing these programs, we would like to outline the differences between the two approaches.

NL³ is a generalisation of ";?>CKKW-L tree-level merging. The main idea of NL³ is to start from CKKW-L-reweighted multi-jet merging, and replace the α_sⁿ⁺⁰- and α_sⁿ⁺¹-terms by the NLO result of POWHEG. This "replacement" means that we subtract the α_sⁿ⁺⁰- and α_sⁿ⁺¹-terms from the CKKW-L-reweighted tree-level samples, and add another sample -- the POWHEG input. All "higher orders" are unchanged w.r.t. CKKW-L. We have implemented the "inclusive" scheme of [Lon13] in Pythia. This means that the POWHEG input will contain contributions for hard, resolved real emission jets, which are already taken care of by higher-multiplicity samples in CKKW-L. Thus, explicit phase space subtractions are also included. The sample program main87.cc, together with the input file main87.cmnd, illustrates the procedure.

UNLOPS is a generalisation of the ";?>UMEPS multi-jet merging scheme. Since UMEPS is already slightly more complicated than CKKW-L, this makes UNLOPS more complicated than NL³. The basic idea however remains the same: Start from a tree-level merging scheme (in this case UMEPS), remove all undesirable α_sⁿ⁺⁰- and α_sⁿ⁺¹-terms from this result, and add back the "correct" description via POWHEG input samples. Again, since the "inclusive" scheme of [Lon13] was implemented in Pythia, it is necessary to handle explicit phase space subtractions. Similar to UMEPS, UNLOPS further ensures that the lowest-multiplicity cross section is given by the NLO result. This means that the UMEPS philosophy of "subtract what you add" needs to be extended to multi-leg NLO inputs.

UNLOPS is a theoretically more appealing definition of NLO merging than NL³, and should thus be considered the preferred choice. However, we believe it valuable to include both methods into Pythia, so that the variation of NLO merged results due to different NLO merging schemes can be studied in situ. Furthermore, NLO merging can be outlined more pedagogically when starting from NL³. The two NLO merging methods share parts of code with CKKW-L and UMEPS, and correspondingly share many input settings with these schemes. In particular,

         The hard process (Merging:Process) needs to be defined exactly as in CKKW-L (see Defining the hard process in the ";?>CKKW-L documentation).

         The merging scale value (Merging:TMS) has to be set.

         The maximal number of additional partons in tree-level events (Merging:nJetMax) has to be set.

All settings listed under the sections "Matrix element merging and HepMC output for RIVET" and "Further variables" in the ";?>CKKW-L documentation can be accessed in NLO merging as well. Furthermore, the Merging:nRecluster switch (see the ";?>UMEPS documentation) is important. Also, all MergingHooks routines that allow for user interference in CKKW-L merging are also usable for NLO merging -- with the exception of a user-defined merging scale. The NLO merging schemes currently implemented in Pythia do not allow for a merging scale definition that differs from the parton shower evolution variable. Since this merging scale definition is not completely obvious, the NLO merging schemes also share the Merging:enforceCutOnLHE switch with CKKW-L. In this way, it is possible to use LHE files that are regularised only with weak cuts as input, while the merging machinery imposes the stronger merging scale cut automatically. This means that no merging scale implementation is required from the user side, but also means that it is the user's responsibility to ensure that the cuts used for generating input LHE files are always looser than the cut given by the merging scale value Merging:TMS. This will lead to warnings of the form "Les Houches Event fails merging scale cut. Cut by rejecting event". These warning should rather be regarded as information. An example of inclusive matrix element generation cuts would be pT_jet = 5 GeV, ΔR_{jetA jetB} = 0.01 and Q_{jetA jetB} = 5 GeV, if NLO merging with a desired merging scale value of Merging:TMS = 15 is attempted for Higgs + jets events at the LHC.

In the following, we will first describe the generation of NLO input samples, and list input settings for NLO merging in Pythia. Then, we will examine the sample main programs main87.cc and main88.cc, which implement NL³ and UNLOPS merging, respectively.

Inputs for NLO merging

         LHE files for 0... M additional partons, taken from a tree-level matrix element generator, and

         LHE files for 0... N additional partons, taken from a POWHEG NLO generator.

All input files need to be regularised, if they contain additional partons. Large files with fairly inclusive (i.e. loose) cuts are recommended. The input LHE files should further be generated with fixed renormalisation and factorisation scales. (In the POWHEG-BOX program, this means using the settings runningscales 0, btlscalereal 1, btlscalect 1, ckkwscalup 0. Some older processes in the POWHEG-BOX program need the input runningscale 0 instead of runningscales 0.)

When attempting NLO merging, the following Pythia settings are relevant.

Merging:nJetMaxNLO

(default = 0; minimum = 0)

The maximal number of additional jets for which NLO event samples are supplied by the user.

Merging:muFac

(default = -1.0)

The fixed factorisation scale used in the hard process cross section, as needed to generate the leading-order weight, in case the factorisation scale cannot be inferred from Les Houches event input. (This is the case for files that have been generated with the POWHEG-BOX program, since this program prints the transverse momentum scale of the real emission into the LH events.). If the value is not set, the SCALUP variable of the current LH event will be used instead. If wimpy showers (see ";?>Timelike Showers and ";?>Spacelike Showers) are used together with multi-jet merging, then this scale further sets the parton shower starting scale (μ_Q) for the core hard process.

Merging:muRen

(default = -1.0)

The fixed renormalisation scale used in the hard process cross section, as needed to generate the leading-order weight, in case the renormalisation scale cannot be inferred from Les Houches event input. (As mentioned above, this is the case for files generated with the POWHEG-BOX program.) If the value is not set, the SCALUP variable of the current LH event will be used instead.

Merging:muFacInME

(default = -1.0)

The fixed factorisation scale used in the matrix element calculation. This information is needed if factorisation scale variations in NLO merged results are attempted. Depending on the matrix element generator, it might not be possible to infer the factorisation scale from Les Houches event input, and thus, setting an explicit value is required. (As mentioned above, this is the case for files generated with the POWHEG-BOX program.) If the value is not set, the SCALUP variable of the current LH event will be used instead.

Merging:muRenInME

(default = -1.0)

The fixed renormalisation scale used in the matrix element calculation. This information is needed if renormalisation scale variations in NLO merged results are attempted, for the same reason as factorisation scales might be required. (As mentioned above, this is the case for files generated with the POWHEG-BOX program.) If the value is not set, the SCALUP variable of the current LH event will be used instead.

All further settings will be discussed while examining the sample main programs.

NL³ merging with main87.cc

If the user only wants to change the number of requested events (Main:numberOfEvents), the hard process (Merging:Process), the merging scale value (Merging:TMS) and the maximal number of additional tree-level or NLO-accuracte jets (Merging:nJetMax and Merging:nJetMaxNLO, respectively), and HepMC output is desired, then there is no need to change the main87.cc code. The input LHE files are also part of the (command line) input for main87.exe. The default settings in main87.cmnd are intended to work with the (very short) sample LHEF inputs (w_production_tree_0.lhe, w_production_tree_1.lhe, w_production_tree_2.lhe and w_production_powheg_0.lhe, w_production_powheg_1.lhe). For these input files, the main87.exe executable can be run with the command

./main87.exe main87.cmnd w_production myhepmc.hepmc

to produce a file myhepmc.hepmc of NLO merged HepMC event output. All mandatory Pythia input settngs have been outlined earlier. Please refrain from adding input switches than invoke any other merging scheme (e.g. e.g. Merging:doKTMerging) into the input file that you want to use in conjunction with main87.cc.

In the following, we will explain main87.cc in depth. Users who are willing to accept the default choices do not need to know all details, but are still encouraged to read further.

Program flow

         1. Estimate the cross section for tree-level and NLO samples after the merging scale cut.

         2. Produce reweighted tree-level events, which do not contain α_s⁰- and α_s¹-terms.

         3. Add POWHEG NLO events.

         4. Subtract phase space points with an extra (real-emission) jet above the merging scale from the POWHEG result, since such configurations have already been taken into account by processing other samples.

The first step is necessary to produce the correct weights for HepMC output events. The estimation of tree-level cross sections after the merging scale cut is generated by invoking the switch Merging:doXSectionEstimate together with Merging:doNL3Tree. In this configuration, the latter switch will only act to define the merging scale. After the tree-level cross sections have been estimated, main87.cc estimates the NLO cross sections after application of the merging scale cut, by inferring Merging:doXSectionEstimate together with Merging:doNL3Loop. Again, in this configuration, the latter switch only acts as the merging scale definition. When generating the estimates, all showering, multiparton interactions and hadronisation is turned off to not unnecessarily waste processor time. For all estimates, is further mandatory to set the value of Merging:nRequested to the jet multiplicity of the current event sample (e.g. to "2" for a sample containing W + 2 jet events). This is necessary in order to correctly apply the merging scale cut. POWHEG NLO input files for W + 1 jet e.g. contain W + 1 jet and W + 2 jet (i.e. real emission) kinematics. However, the merging scale cut aims at regularising the "underlying Born" configuration (i.e. the W + 1 states in our example). Setting Merging:nRequested = 1 for the W + 1 jet POWHEG sample ensures that even for real-emission (W + 2 jet) kinematics, the merging scale cut is applied to W + 1 jet states.

After the cross section estimation step, main87.cc proceeds to perform the actual merging. Before explaining this part, we would like to make some comments about K-factors.

main87.cc is prepared to use fixed K-factors to rescale the weight of tree-level events. This rescaling does not affect the NLO accuracy of the method, and was investigated in [Lon13]. By default, main87.cc does not use K-factors. However, if the user wants to include K-factors, this can be done by using the following input settings.

Merging:kFactor0j

(default = 1.0)

The k-Factor used to rescale the tree-level (i.e. CKKW-L or UMEPS) part of zero-jet tree-level events.

Merging:kFactor1j

(default = 1.0)

The k-Factor used to rescale the tree-level (i.e. CKKW-L or UMEPS) part of one-jet tree-level events.

Merging:kFactor2j

(default = 1.0)

The k-Factor used to rescale the tree-level (i.e. CKKW-L or UMEPS) part of two-jet tree-level events.

If the variables k0, k1, k2 in main87.cc are set to non-unity values, K-factors will be applied. The K-factor of highest jet multiplicity will then be used to also rescale tree-level samples with a number of additional jets beyond the number of the highest-multiplicity real-emission sample. If we, for example, attempt an NLO merging of W+0 jet and W+1 jet at NLO accuracy, and with W+≤4 jets at tree-level accuracy, then Merging:kFactor2j is used to rescale the W+2 jet, W+3 jets and W+4 jets tree-level samples. We recommend to not include a K-factor rescaling of the tree-level samples.

Let us turn to the production of NLO merged events. The first step in the procedure is to generate reweighted tree-level samples. This is implemented by using the following switch.

Merging:doNL3Tree On Off   (default = off)
This switch will allow the generation of the weight that should be applied to tree-level events in the NL³ merging scheme. Please note that, in order for this to work smoothly, the switch Merging:doNL3Loop and the switch Merging:doNL3Subt have to be turned off. As for the estimation of cross sections, it is mandatory to set the correct value of Merging:nRequested.

The weight of tree-level events can be accessed by calling the function double Info::mergingWeightNLO(). When printing (or histogramming) NLO merged events, this weight, multiplied with the estimated cross section of the current event sample, should be used as event weight (or weight of histogram bins). For Merging:doNL3Tree = on, the weight double Info::mergingWeightNLO() contains the CKKW-L weight, subtracted, if necessary, by α_s⁰- and α_s¹-terms. This weight can become negative. As an example, imagine we attempt an NLO merging of W + 0 jet and W + 1 jet at NLO accuracy, and with W + 2 jets at tree-level accuracy. This weight will then be

  Info::mergingWeightNLO() = CKKW-L-weight for zero jets - α_s⁰-terms - α_s¹-terms for events in the zero-jet sample,

  Info::mergingWeightNLO() = CKKW-L-weight for one jet - α_s⁰-terms - α_s¹-terms for events in the one-jet sample, and

  Info::mergingWeightNLO() = CKKW-L-weight for two jets for events in the two-jet sample.

After the tree-level events have been reweighted, main87.cc will move on to process the POWHEG NLO input. This is done by switching to the following flag.

Merging:doNL3Loop On Off   (default = off)
This switch will allow the processing of POWHEG NLO events in the NL³ merging scheme. Please note that, in order for this to work smoothly, the switch Merging:doNL3Tree and the switch Merging:doNL3Subt have to be turned off. As for the estimation of cross sections, it is mandatory to set the correct value of Merging:nRequested.

Also in this case, the NLO merging weight of the events can be accessed by calling the function double Info::mergingWeightNLO(). This weight should also be used when printing (or histogramming) events. For Merging:doNL3Loop = on, the weight double Info::mergingWeightNLO() is either one or zero (see Appendix E in [Lon13]). After the processing of POWHEG NLO events, main87.cc continues by generating explicit phase space subtractions. This is facilitated by the following switch.

Merging:doNL3Subt On Off   (default = off)
This switch will allow the processing of tree-level events, to produce explicit phase space subtractions in the NL³ merging scheme. Please note that, in order for this to work smoothly, the switch Merging:doNL3Tree and the switch Merging:doNL3Loop have to be turned off. As for the estimation of cross sections, it is mandatory to set the correct value of Merging:nRequested. Furthermore, it is necessary to set the value of Merging:nRecluster to one.

These contributions are necessary because we have implemented the "inclusive scheme" of [Lon13] in Pythia. The benefit of this scheme is the user does not have to intrusively change the POWHEG-BOX program to implement very particular cuts. Let us explain this comment with an example (a more detailed explanation of the idea is given in Appendix A.2 of [Lon13]). When generating W + 0 jet events with the POWHEG-BOX program, the output LHE files will contain W + 1 jet real emission events. Some of these events will contain a jet above the merging scale. However, in NLO merging methods, such configurations have already been included by a separate W + 1 jet sample. Thus, to avoid counting such events twice, we have to remove the configurations from the POWHEG-BOX output. We choose to remove such events by explicit subtraction.

As always, the NLO merging weight of the events can be accessed by calling the function double Info::mergingWeightNLO(). This weight should also be used when printing (or histogramming) events. For Merging:doNL3Subt = on, the weight double Info::mergingWeightNLO() is either one or zero (see Appendix E in [Lon13]).

After these steps, all necessary events for NL³ merging have been produced. main87.cc finishes by returning the NL³-merged total cross section.

UNLOPS merging with main88.cc

As for NL³, it is not necessary to change main88.cc if the user is only interested in changing standard settings. Thus, if the user only wants to change the number of requested events (Main:numberOfEvents), the hard process (Merging:Process), the merging scale value (Merging:TMS) and the maximal number of additional tree-level or NLO-accuracte jets (Merging:nJetMax and Merging:nJetMaxNLO, respectively), and HepMC output is desired, then there is no need to change the main88.cc code. The input LHE files are also part of the (command line) input for main88.exe. The default settings in main88.cmnd are intended to work with the (very short) sample LHEF inputs (w_production_tree_0.lhe, w_production_tree_1.lhe, w_production_tree_2.lhe and w_production_powheg_0.lhe, w_production_powheg_1.lhe). For these input files, the main88.exe executable can be run with the command

./main88.exe main88.cmnd w_production myhepmc.hepmc

to produce a file myhepmc.hepmc of UNLOPS merged HepMC event output. Please refrain from adding input switches than invoke any other merging scheme (e.g. Merging:doKTMerging) into the input file that you want to use in conjunction with main88.cc.

In the following, we will explain main88.cc in depth. To not be overly repetitive, we will at times refer to the relevant parts in the discussion of main87.cc. Users who are willing to accept the default choices do not need to know all details, but are still encouraged to read further.

Program flow

         1. Estimate the cross section for tree-level and NLO samples after the merging scale cut.

         2. Produce reweighted tree-level events, which do not contain α_s⁰- and α_s¹-terms.

         3. Add POWHEG NLO events.

         4. Subtract integrated, reweighted tree-level events, to ensure that the inclusive NLO cross section remains intact upon inclusion of multi-jet tree-level events.

         5. Subtract integrated POWHEG NLO events, to ensure that the inclusive NLO cross section remains intact upon inclusion of multi-jet tree-level events.

The estimation of cross sections after the application of the merging scale cut is nearly identical to the first step in main87.cc, and we refer to the first paragraph of the "Program flow" discussion for main87.cc for details. For main88.cc, the flags Merging:doUNLOPSTree or Merging:doUNLOPSLoop supply the merging scale definition used in the cross section estimation.

After the cross section estimation step, main88.cc proceeds to perform the actual NLO merging. The discussion of K-factors given in the NL³ section (i.e. of Merging:kFactor0j, Merging:kFactor1j and Merging:kFactor2j) also applies to main88.cc. Although UNLOPS is considerably more stable than NL³ upon changing the K-factors, we do not recommend the use of K-factors.

The production of UNLOPS-merged events with main88.cc starts by generating reweighted tree-level events. The processing of tree-level events can be invoked by setting the following flag.

Merging:doUNLOPSTree On Off   (default = off)
This switch will allow the generation of the weight that should be applied to tree-level events in the UNLOPS merging scheme. Please note that, in order for this to work smoothly, the switches Merging:doUNLOPSLoop, Merging:doUNLOPSSubt and Merging:doUNLOPSSubtNLO have to be turned off. As for the estimation of cross sections, it is mandatory to set the correct value of Merging:nRequested.

The weight of tree-level events is returned by the function double Info::mergingWeightNLO(). When printing (or histogramming) NLO merged events, this weight, multiplied with the estimated cross section of the current event sample, should be used as event weight (or weight of histogram bins). For Merging:doUNLOPSTree = on, the weight double Info::mergingWeightNLO() contains the UMEPS weight, subtracted, if necessary, by α_s⁰- and α_s¹-terms. This weight can become negative. As an example, assume that we attempt an UNLOPS merging of W + 0 jet and W + 1 jet at NLO accuracy, and with W + 2 jets at tree-level accuracy. This weight will then be

  Info::mergingWeightNLO() = UMEPS-weight for one jet - α_s⁰-terms - α_s¹-terms for events in the one-jet sample, and

  Info::mergingWeightNLO() = UMEPS-weight for two jets for events in the two-jet sample.

After reweighted tree-level events have been generated, main88.cc processes the POWHEG NLO input files. This is facilitated by the following switch.

Merging:doUNLOPSLoop On Off   (default = off)
This switch will allow the processing of POWHEG NLO events in the UNLOPS merging scheme. Please note that, in order for this to work smoothly, the switches Merging:doUNLOPSTree, Merging:doUNLOPSSubt and Merging:doUNLOPSSubtNLO have to be turned off. As for the estimation of cross sections, it is mandatory to set the correct value of Merging:nRequested.

The NLO merging weight of the events can be accessed by calling the function double Info::mergingWeightNLO(). This weight should also be used when printing (or histogramming) events. For Merging:doUNLOPSLoop = on, the weight double Info::mergingWeightNLO() is either one or zero (see Appendix E in [Lon13]).

After processing the POWHEG NLO events, main88.cc continues by generating the reweighted subtraction terms of UMEPS. This part is implemented by setting the following flag.

Merging:doUNLOPSSubt On Off   (default = off)
This switch will allow the processing of tree-level events, to produce UMEPS subtraction terms for the UNLOPS merging scheme. Please note that, in order for this to work smoothly, the switches Merging:doUNLOPSTree, Merging:doUNLOPSLoop and Merging:doUNLOPSSubtNLO have to be turned off. As for the estimation of cross sections, it is mandatory to set the correct value of Merging:nRequested. Furthermore, it is necessary to set the value of Merging:nRecluster to one.

By using this switch, main88.cc ensures that the inclusive cross section is preserved. At variance with UMEPS however, the event weight contains the UMEPS weight, subtracted, if necessary, by α_s⁺⁰- and α_s¹-terms. Otherwise, α_sⁿ⁺⁰- and α_sⁿ⁺¹-terms of the UMEPS procedure would be introduced, although our goal is to describe all α_sⁿ⁺⁰- and α_sⁿ⁺¹-terms by n-jet POWHEG input. The weight of these integrated, subtractive tree-level events is, as always, returned by the function double Info::mergingWeightNLO(). When printing (or histogramming) NLO merged events, this weight, multiplied with the estimated cross section of the current event sample, and with -1, should be used as event weight (or weight of histogram bins). As for the case of tree-level events in UNLOPS, this weight can become negative. For the example given before, i.e. attempting an UNLOPS merging of W + 0 jet and W + 1 jet at NLO accuracy, and with W + 2 jets at tree-level accuracy, this weight will be

  Info::mergingWeightNLO() = UMEPS-weight for the integrated one-jet sample - α_s⁰-terms - α_s¹-terms for events in the integrated one-jet sample, and

  Info::mergingWeightNLO() = UMEPS-weight for the integrated two-jet sample for events in the integrated two-jet sample.

This choice of weights already incorporates the fact that we have implemented the "inclusive scheme" of [Lon13], meaning that the "explicit phase space subtractions" of NL³ are (partially) included though these weights.

To ensure that the NLO inclusive cross section is unchanged, UNLOPS further requires the introduction of another sample. If POWHEG NLO events with one or more jets are included, it is necessary to subtract these samples in an integrated form. In main88.cc, this is done by setting the following flag.

Merging:doUNLOPSSubtNLO On Off   (default = off)
This switch will allow the processing of POWHEG NLO events, to produce NLO subtraction terms for the UNLOPS merging scheme. Please note that, in order for this to work smoothly, the switches Merging:doUNLOPSTree, Merging:doUNLOPSLoop and Merging:doUNLOPSSubt have to be turned off. As for the estimation of cross sections, it is mandatory to set the correct value of Merging:nRequested. Furthermore, it is necessary to set the value of Merging:nRecluster to one.

This sample also provides some "explicit phase space subtractions" of NL³, which are necessary because we implemented the "inclusive scheme" of [Lon13]. Let us again look at the example of UNLOPS merging of W + 0 jet and W + 1 jet at NLO accuracy. The integrated W + 1 jet NLO events, which are produced by the setting Merging:doUNLOPSSubtNLO = on, contain a tree-level part. This part will exactly cancel the real-emission events with one jet above the merging scale in the W + 0 jet NLO events.

The NLO merging weight of these "integrated" events can be accessed by calling the function double Info::mergingWeightNLO(). This weight should also be used when printing (or histogramming) events. For Merging:doUNLOPSSubtNLO = on, the weight double Info::mergingWeightNLO() is either one or zero (see Appendix E in [Lon13]).

After these five steps (estimation of cross sections, tree-level processing, POWHEG processing, integrated tree-level processing, integrated POWHEG processing) we have produced a UNLOPS-merged HepMC event file. main88.cc finishes by returning the UNLOPS-merged total cross section.

NLO merging and "exclusive" NLO inputs

We should be very clear that using exclusive NLO input is not recommended, since it requires hacking the NLO generator. Only for the expert user, we briefly summarise the necessary changes for using exclusive NLO input.

For the moment, assume that the NLO input has been produced in an "exclusive" way. This input can then be processed by some trivial changes in main87.cc: estimate the cross section for tree-level and NLO samples after the merging scale cut, still using inclusive NLO samples, remove the last part of main87.cc, i.e. the part that produces explicit phase space subtractions, and use the exclusive NLO files as input files for the processing of "POWHEG NLO files".

The changes to main88.cc (implementing UNLOPS) are slightly more complicated. This is the case because the weights of integrated tree-level samples change when using exclusive input, as can be seen in Appendix D in [Lon13]. The correct weights can be produced by Pythia by using the following flag.

Merging:doUNLOPSTilde On Off   (default = off)
This flag allows the UNLOPS machinery to produce the event weights if exclusive NLO input is used for the merging. This flag should be set to "on" directly after the cross section estimates have been produced.

Then, it is necessary to add code for processing another sample to main88.cc, since when using exclusive inputs, it is also necessary to enforce two integrations on tree-level events (the "↑"-contributions in Appendix D of [Lon13]). This can be achieved by adding the following code at the end of main88.cc.

 
  cout << endl << endl << endl; 
  cout << "Shower subtractive events" << endl; 
 
  // Switch on processing of counter-events. 
  pythia.settings.flag("Merging:doUNLOPSTree",false); 
  pythia.settings.flag("Merging:doUNLOPSLoop",false); 
  pythia.settings.flag("Merging:doUNLOPSSubt",true); 
  pythia.settings.flag("Merging:doUNLOPSSubtNLO",false); 
  pythia.settings.mode("Merging:nRecluster",2); 
 
  nMaxCT         = nMaxNLO+1; 
  njetcounterCT  = nMaxCT; 
  iPathSubt      = iPath + "_tree"; 
 
  while(njetcounterCT >= 2){ 
 
    // From njet, choose LHE file 
    stringstream in; 
    in   << "_" << njetcounterCT << ".lhe"; 
    string LHEfile = iPathSubt + in.str(); 
 
    cout << endl << endl << endl 
         << "Start subtractive treatment for " << njetcounterCT << " jets\n" 
         << "Recluster at least 2 times" << endl; 
 
    pythia.readString("Beams:frameType = 4"); 
    pythia.settings.word("Beams:LHEF", LHEfile); 
    pythia.settings.mode("Merging:nRequested", njetcounterCT); 
    pythia.init(); 
    // Remember position in vector of cross section estimates. 
    int iNow = sizeLO-1-njetcounterCT; 
 
    // Start generation loop 
    for( int iEvent=0; iEvent < nEvent; ++iEvent ){ 
 
      // Generate next event 
      if( !pythia.next() ) { 
        if( pythia.info.atEndOfFile() ) break; 
        else continue; 
      } 
 
      // Get event weight(s). 
      double weightNLO  = pythia.info.mergingWeightNLO(); 
      double evtweight  = pythia.info.weight(); 
      weightNLO        *= evtweight; 
      // Do not print zero-weight events. 
      if ( weightNLO == 0. ) continue; 
 
      // Construct new empty HepMC event. 
      HepMC::GenEvent* hepmcevt = new HepMC::GenEvent(); 
      // Get correct cross section from previous estimate. 
      double normhepmc = -1*xsecLO[iNow] / nAcceptLO[iNow]; 
      // Set hepmc event weight. 
      hepmcevt->weights().push_back(weightNLO*normhepmc); 
      // Fill HepMC event. 
      ToHepMC.fill_next_event( pythia, hepmcevt ); 
      // Add the weight of the current event to the cross section. 
      sigmaTotal += weightNLO*normhepmc; 
      errorTotal += pow2(weightNLO*normhepmc); 
      // Report cross section to hepmc. 
      HepMC::GenCrossSection xsec; 
      xsec.set_cross_section( sigmaTotal*1e9, pythia.info.sigmaErr()*1e9 ); 
      hepmcevt->set_cross_section( xsec ); 
      // Write the HepMC event to file. Done with it. 
      ascii_io << hepmcevt; 
      delete hepmcevt; 
 
    } // end loop over events to generate 
 
    // print cross section, errors 
    pythia.stat(); 
 
    // Restart with ME of a reduced the number of jets 
    if( njetcounterCT > 2 ) 
      njetcounterCT--; 
    else 
      break; 
 
  } 

Further variables