Merge lp:~maddevelopers/mg5amcnlo/bias_function into lp:~mg5core1/mg5amcnlo/2.5.6

[Rikkert Frederix]

Including the bias_weight_function to have biased event generation at NLO(+PS) accuracy.

[Rikkert Frederix (frederix)]

Stefano already performed some tests, showing that it works as expected with 'normal' running. Things that need to be tested are including external modules, like MadSpin, or the REWEIGHT. Also, if you have a good suggestion for a unit/acceptance test, please consider implementing it.

[Valentin Hirschi (valentin-hirschi)]

Hi Rik,

Great to see this functionality, thanks!

It would be nice to amend the bias wiki page:

https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/LOEventGenerationBias

with this new NLO functionality.

Also, how does it work exactly? Can you easily use user-defined biasing weight formulae?

A very quick look reveals that you intend to steer it via the 'event_norm' run_card entry.
Wouldn't it be handy to adopt the already existing structure at LO with:

 None = bias_module ! Bias type of bias, [None, ptj_bias, -custom_folder-]
 {} = bias_parameters ! Specifies the parameters of the module.

Where the user can create his own bias module and parameters (according to the conventions laid down on the wiki page) which will be moved automatically to the Source/BIAS directory and be compiled/link against the code?

[Rikkert Frederix (frederix)]

Hi Valentin,

Stefano and I discussed this last week with Olivier, and decided not to use the LO syntax. CMS has been asking this for some time and we want to get it out asap (essentially the coming week).
At NLO, the bias function is not as easy as at LO, e.g. the function needs to be IR safe, and therefore, at least for the moment, keep it somewhat hidden for normal users. So, we don't want that users have a function that works at LO, and then simply try to use that at NLO without thinking it over. Hence, we should not use the same syntax. Anyway, also the current way of doing it, was already available at fNLO, so, if there were some users already using that, they can effectively use exactly the same as before.

Cheers,
Rik

[Valentin Hirschi (valentin-hirschi)]

Ok, as you might suspect I have a different take on such matters, but I
understand your concerns and resulting choice.

In any case, amending the bias wiki page with the details of how to setup
bias at NLO and what the pitfalls/limirations is important I believe.

On Sat, Jul 29, 2017 at 12:44 Rikkert Frederix <email address hidden>
wrote:

> Hi Valentin,
>
> Stefano and I discussed this last week with Olivier, and decided not to
> use the LO syntax. CMS has been asking this for some time and we want to
> get it out asap (essentially the coming week).
> At NLO, the bias function is not as easy as at LO, e.g. the function needs
> to be IR safe, and therefore, at least for the moment, keep it somewhat
> hidden for normal users. So, we don't want that users have a function that
> works at LO, and then simply try to use that at NLO without thinking it
> over. Hence, we should not use the same syntax. Anyway, also the current
> way of doing it, was already available at fNLO, so, if there were some
> users already using that, they can effectively use exactly the same as
> before.
>
> Cheers,
> Rik
>
> --
>
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/bias_function/+merge/328265
> Your team mg5core1 is subscribed to branch lp:~mg5core1/mg5amcnlo/2.5.6.
>
--
Valentin

[Stefano Frixione (stefano-frixione)]

Hi Valentin,
to add to what Rikkert said.

I've taken a look at the LO bias wiki page (which, I agree, should be
amended to include information on the NLO stuff).

It is all very beautiful, but also fairly discouraging for the
average user. It is *way* too complicated for someone who simply
doesn't want to bother (which is the majority of the users and,
maybe, myself).

Notably (which is the truly important thing here) the average
ATLAS user fall into the latter category: he/she is not even able
to reproduce some very basic behaviour (as many "bug" reports can
testify).

If one wants to take a positive outlook on this, such an attitude
is not unjustified: the code should simplify life for those who want
to extract physics from it: if the technical infrastucture becomes too
heavy, its flexibility and power will simply become irrelevant: users
are put off -- trust me on that, since I talk to a very large number
of experimentalists.

Finally, in this specific case: NLO gives one many more possibilities
than LO to screw things up big time. Furthermore, the NLO and LO
structures are already different, which is why a difference in this
context is perfectly acceptable.

Cheers, Stefano.

On Sat, 29 Jul 2017, Valentin Hirschi wrote:

> Hi Rik,
>
> Great to see this functionality, thanks!
>
> It would be nice to amend the bias wiki page:
>
> https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/LOEventGenerationBias
>
> with this new NLO functionality.
>
> Also, how does it work exactly? Can you easily use user-defined biasing weight formulae?
>
> A very quick look reveals that you intend to steer it via the 'event_norm' run_card entry.
> Wouldn't it be handy to adopt the already existing structure at LO with:
>
> None = bias_module ! Bias type of bias, [None, ptj_bias, -custom_folder-]
> {} = bias_parameters ! Specifies the parameters of the module.
>
> Where the user can create his own bias module and parameters (according to the conventions laid down on the wiki page) which will be moved automatically to the Source/BIAS directory and be compiled/link against the code?
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/bias_function/+merge/328265
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/bias_function into lp:~mg5core1/mg5amcnlo/2.5.6.
>

[lp:~maddevelopers/mg5amcnlo/bias_function]

308. By Rikkert Frederix on 2017-08-03

fixed the definition of H_T in the bias-function

[Olivier Mattelaer (olivier-mattelaer)]

Hi rikkert,

I have taken a look at the code and this is my comment

1) I do not like very much to use "event_norm" to trigger the bias. The main reason is that this flag is used to control how the weight of the event should be written. How would I set the normalization of the event to "Sum" in presence of bias?

Do we really need a flag in the run_card? Can't we have the default bias function to be "1"?

Like this, this will creates problem with the re-weighting/systematics module (but trivial to fix) since they are all using that variable to know how to handle the weights.

2) the HTML/printout results, only display the un-physical number (with a very nice warning).
Should not we automatically display the correct result? I guess that this is fine like this since we want to have this only for very advanced user.

3) So For the moment, the only way to get the physical cross-section is trough the systematics printout. However for the "example" bias [pt_t**3]. This returns a wrong number (for p p > t t~):
It returns:
   --------------------------------------------------------------
      Scale variation (computed from LHE events):
          Dynamical_scale_choice -1 (envelope of 9 values):
              8.060e+02 pb +11.0% -12.0%
instead of
   --------------------------------------------------------------
      Scale variation (computed from LHE events):
          Dynamical_scale_choice -1 (envelope of 9 values):
              6.875e+02 pb +9.6% -11.7%

It seems to be a numerical precision issue, since if I change the bias to pt_top:
   --------------------------------------------------------------
      Scale variation (computed from LHE events):
          Dynamical_scale_choice -1 (envelope of 9 values):
              6.742e+02 pb +9.5% -11.6%
   --------------------------------------------------------------
However already using pt_top^2 seems bad:
   --------------------------------------------------------------
      Scale variation (computed from LHE events):
          Dynamical_scale_choice -1 (envelope of 9 values):
              7.187e+02 pb +8.8% -11.3%
   --------------------------------------------------------------
If for this bias pt_top^2, I increase the statistic to 250k, then I have:
   --------------------------------------------------------------
      Scale variation (computed from LHE events):
          Dynamical_scale_choice -1 (envelope of 9 values):
              6.886e+01 pb +8.8% -10.3%
   --------------------------------------------------------------

So maybe we should a comment about this "problem" in cuts.f (looks like we are much more sensitive than at LO)

4) After fixing systematics.py in order to assume that bias normalisation is equivalent to average.
Then I have the following results
4.1)from the online version:
   --------------------------------------------------------------
      Summary:
      Process p p > t t~ [QCD]
      Run at p-p collider (6500.0 + 6500.0 GeV)
      Number of events generated: 10000
      Total cross section, incl. bias (DO NOT USE): 8.236e+04 +- 5.0e+02 pb
   --------------------------------------------------------------
      Scale variation (computed from LHE events):
 ...

review: Needs Fixing

[Stefano Frixione (stefano-frixione)]

Hi Olivier,
some comments on your comments.

> 1) I do not like very much to use "event_norm" to trigger the bias. The
> main reason is that this flag is used to control how the weight of
> the event should be written. How would I set the normalization of the
> event to "Sum" in presence of bias?
well, tastes aside (I personally like "bias") this option has a rather
direct bearing on the nature of the weights written in the LHEF, and
this includes their "norm". Furthermore, "bias" (or equivalent)
conveys immediately the message that these, at variance with the
other weights produced out of the integration, are not unweighted.

> Do we really need a flag in the run_card? Can't we have the default bias
> function to be "1"?
I strongly oppose this. The current choice may indeed be redundant,
but it is the only safe one. The possibility of screwing with this
stuff is significant, so we do need an extra layer of protection.

What can be done (no strong feelings here) is to have a check that,
if the bias function is identically equal to one, the warning printouts
are removed from the output.

> Like this, this will creates problem with the re-weighting/systematics
> module (but trivial to fix) since they are all using that variable to
> know how to handle the weights.
I think we shouldn't do reverse engineering here: it's a new option,
so the systematics module must cope with it.

> 2) the HTML/printout results, only display the un-physical number (with
> a very nice warning). Should not we automatically display the
> correct result? I guess that this is fine like this since we want to
> have this only for very advanced user.
Well, I thought that the number printed out resulted from the integration,
not from the sum of event weights (since the former is more precise).
That being the case, the currect choice has no alternative, I'm afraid
(which is OK with me).

> 3) So For the moment, the only way to get the physical cross-section is
> trough the systematics printout. However for the "example" bias
......
> It seems to be a numerical precision issue, since if I change the bias
......
> So maybe we should a comment about this "problem" in cuts.f (looks like
> we are much more sensitive than at LO)
Well, this is expected, isn't it? By biasing, one shouldn't obtain
a particularly reliable number for inclusive rates, but smoother
high-pt tails. A comment, if any, should stress only this fact.

> 4) After fixing systematics.py in order to assume that bias
>
> So this indicates a problem in the way the weight is written in the file.
....
Note: in my tests, I've only looked at distributions. Although I've
mostly concentrated on the central weights, a quick looks at results
for non central scale values didn't show anything strange.
I'm using standard analyses as provided in ./MCatNLO subdirectories.

Cheers, Stefano.

> Did you change that format? did you add one entry?
> This is important since this will also forbid the re-weighting actually.
> (We can skype about this if you want to).
>
>
>
>
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/bias_function/+merge/328265
> You are requested to review the propo...

[lp:~maddevelopers/mg5amcnlo/bias_function]

309. By Rikkert Frederix on 2017-08-07

Added a comment in the bias-function header saying the the statistical
uncertainty from computing the cross section by averaging the event
weights can be greatly enlarged when including a non-zero
bias-function

[Rikkert Frederix (frederix)]

Hi Olivier,

1. I agree with Stefano here. Also, the python code has to do things slightly differently when including this bias function: it might be slightly more efficient in the original way as compared to including the bias. Hence, I would prefer to have an access at the level of the python code (which means some sort of flag in the run_card).

2. See Stefano.

3. I've added a comment in the header of the bias-weight-function in cuts.f stating the the uncertainty on the cross section by summing all the event weights can be larger when including the bias function as compared to having all events of equal weight.

4. I'll look for you on Skype asap.

Cheers,
Rik

[Rikkert Frederix (frederix)]

Hi Olivier,

One more thing that I just realise only now. I see that the systematics also includes reweighting with muR=muF=sqrt(shat). This cannot be done at NLO, since shat is not an IR-safe quantity. (Think about collinear initial state radiation).

Hence, we should never include this option for NLO.

Cheers,
Rik

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

Well this is more tricky since you have case where it does make sense.
Like e+ e- > b b~ [QCD]

Cheers,

Olivier

[Stefano Frixione (stefano-frixione)]

Hi Olivier,
not quite: the scale in that case is MZ. Which is almost the same
as sqrt(s), but only if you're talking about QCD. Given that we're
about to release the EW stuff, let's try and keep things as
general as possible, in order to be able to apply them to a
broader spectrum of possibilities.
Cheers, Stefano

On Mon, 7 Aug 2017, Olivier Mattelaer wrote:

> Hi,
>
> Well this is more tricky since you have case where it does make sense.
> Like e+ e- > b b~ [QCD]
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/bias_function/+merge/328265
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/bias_function into lp:~mg5core1/mg5amcnlo/2.5.6.
>

[lp:~maddevelopers/mg5amcnlo/bias_function]

310. By olivier-mattelaer on 2017-08-07

fixing the normalisation issue related to the systematics

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

I have implemented a switch to remove it automatically as soon as we have some proton in the initial state.

I have also fix the systematics computation. I do not have test the re-weighting (but I have switch all the part depending of event_norm). For that one, I will merge this branch into reweight_mass branch to be able to do the test in that branch (since I need to do some change to that branch due to this one anyway).

Rikkert, did you check the potential side effect on the parton-shower?

I see the following line in the function banner_to_mcatnlo:
content += 'EVENT_NORM=%s\n' % self.banner.get_detail('run_card', 'event_norm').lower()
so the shower might have some side effect as well related to the new event_normalisation.

Those in turn are (I guess) related to those:
./NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIG6.Script:381: $EVENT_NORM ! Event weights are normalized to sun or average to the cross section
./NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6PT.Script:385: $EVENT_NORM ! Event weights are normalized to sun or average to the cross section
./NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6Q.Script:385: $EVENT_NORM ! Event weights are normalized to sun or average to the cross section
./NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script:355:Main:spareWord1 = $EVENT_NORM ! Event weights are normalized to sum

I check one of those:
cat <<EOF >> ./$ifile
  $EVENT_NORM ! Event weights are normalized to sun or average to the cross section
EOF
But after that I loose track on the potential impact of your change.

Cheers,

Olivier

[Rikkert Frederix (frederix)]

Hi Olivier,

The event normalisation after shower has been taken care of. In the drivers that start the shower, instead of the checking for 'event_norm==average', they check for 'event_norm!=sum'. See the first couple of files in the diff below.

Cheers,
Rik

[Olivier Mattelaer (olivier-mattelaer)]

Perfect then.

I approve the merge.

Cheers,

Olivier

review: Approve

[Olivier Mattelaer (olivier-mattelaer)]

Hi Rik,

I actually just found a bug when testing the re-weighting in the other branch.
If I use the nevt_jobs options with the 'bias' module then the results are crap.
I should have notice it above but in that case it was so close to a factor 10 that I miss it.

If I ask 250k events with nevt_jobs at 7500 then I have (for tt~):
   --------------------------------------------------------------
      Summary:
      Process p p > t t~ [QCD]
      Run at p-p collider (6500.0 + 6500.0 GeV)
      Number of events generated: 250000
      Total cross section, incl. bias (DO NOT USE): 1.427e+07 +- 3.6e+04 pb
   --------------------------------------------------------------
      Scale variation (computed from LHE events):
          Dynamical_scale_choice -1 (envelope of 9 values):
              9.571e+01 pb +8.4% -10.0%
   --------------------------------------------------------------

and you can see that the central value for the scale variation is everything but reliable.
at 250k events, I do not think that this is an effect of the statistical uncertainty.

I guess that this is a side effect of the "bias" flag.

Cheers,

Olivier

review: Needs Fixing

[Rikkert Frederix (frederix)]

Hi Olivier,

I've fixed it directly in 2.5.6, since this branch was already merged. It was only a problem with the cross section printed to the screen; the weights of the actual events were okay.

Cheers,
Rik

[Olivier Mattelaer (olivier-mattelaer)]

Thanks,

> I've fixed it directly in 2.5.6, since this branch was already merged. It was only a problem with the cross section printed to the screen; the weights of the actual events were okay.

Make sense since the (top mass) re-weighting was actually returning the correct value.

Thanks,

Olivier

> On 8 Aug 2017, at 14:42, Rikkert Frederix <email address hidden> wrote:
>
> Hi Olivier,
>
> I've fixed it directly in 2.5.6, since this branch was already merged. It was only a problem with the cross section printed to the screen; the weights of the actual events were okay.
>
> Cheers,
> Rik
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/bias_function/+merge/328265
> You are reviewing the proposed merge of lp:~maddevelopers/mg5amcnlo/bias_function into lp:~mg5core1/mg5amcnlo/2.5.6.

=== modified file 'Template/NLO/Cards/run_card.dat'
--- Template/NLO/Cards/run_card.dat	2016-09-08 23:15:34 +0000
+++ Template/NLO/Cards/run_card.dat	2017-08-07 21:43:30 +0000
@@ -36,7 +36,7 @@
 # Normalize the weights of LHE events such that they sum or average to *
 # the total cross section                                              *
 #***********************************************************************
- %(event_norm)s = event_norm      ! average or sum
+ %(event_norm)s = event_norm    ! valid settings: average, sum, bias
 #***********************************************************************
 # Number of points per itegration channel (ignored for aMC@NLO runs)   *
 #***********************************************************************

=== modified file 'Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f'
--- Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f	2017-07-25 08:06:10 +0000
+++ Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f	2017-08-07 21:43:30 +0000
@@ -323,7 +323,7 @@
 C---EVENTS ARE NORMALIZED TO SUM OR AVERAGE TO THE TOTAL CROSS SECTION
       WRITE(*,*)'How are the events normalized ("average" or "sum")?'
       READ(*,*)NORM_EVENT
-      if (NORM_EVENT.eq.'average')MQQ=1
+      if (NORM_EVENT(1:3).ne.'sum')MQQ=1
 
       MSFLAG=0
       IF (LHSOFT) THEN

=== modified file 'Template/NLO/MCatNLO/srcPythia/madfks_pydriver.f'
--- Template/NLO/MCatNLO/srcPythia/madfks_pydriver.f	2014-11-27 15:36:48 +0000
+++ Template/NLO/MCatNLO/srcPythia/madfks_pydriver.f	2017-08-07 21:43:30 +0000
@@ -318,7 +318,7 @@
 C---EVENTS ARE NORMALIZED TO SUM OR AVERAGE TO THE TOTAL CROSS SECTION
       WRITE(*,*)'How are the events normalized ("average" or "sum")?'
       READ(*,*)NORM_EVENT
-      if (NORM_EVENT.eq.'average')MQQ=1
+      if (NORM_EVENT(1:3).ne.'sum')MQQ=1
 
 C---USER'S INITIAL CALCULATIONS
       CALL PYABEG

=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia8.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia8.cc	2016-02-23 11:25:37 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia8.cc	2017-08-07 21:43:30 +0000
@@ -58,7 +58,7 @@
   int iEventshower=pythia.mode("Main:spareMode1");
   string evt_norm=pythia.word("Main:spareWord1");
   int iEventtot_norm=iEventtot;
-  if (evt_norm == "average"){
+  if (evt_norm != "sum"){
     iEventtot_norm=1;
   }
 
@@ -112,7 +112,7 @@
     double normhepmc;
     // Add the weight of the current event to the cross section.
     normhepmc = 1. / double(iEventshower);
-    if (evt_norm == "average") {
+    if (evt_norm != "sum") {
       sigmaTotal  += evtweight*normhepmc;
     } else {
       sigmaTotal  += evtweight*normhepmc*iEventtot;

=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia82.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia82.cc	2016-02-23 11:25:37 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia82.cc	2017-08-07 21:43:30 +0000
@@ -57,7 +57,7 @@
   int iEventshower=pythia.mode("Main:spareMode1");
   string evt_norm=pythia.word("Main:spareWord1");
   int iEventtot_norm=iEventtot;
-  if (evt_norm == "average"){
+  if (evt_norm != "sum"){
     iEventtot_norm=1;
   }
 
@@ -118,7 +118,7 @@
     double normhepmc;
     // Add the weight of the current event to the cross section.
     normhepmc = 1. / double(iEventshower);
-    if (evt_norm == "average") {
+    if (evt_norm != "sum") {
       sigmaTotal  += evtweight*normhepmc;
     } else {
       sigmaTotal  += evtweight*normhepmc*iEventtot;

=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia82_hep.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia82_hep.cc	2017-07-27 05:51:14 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia82_hep.cc	2017-08-07 21:43:30 +0000
@@ -83,7 +83,7 @@
     double evtweight = pythia.info.weight();
     double normhepmc;
     // ALWAYS NORMALISE HEPMC WEIGHTS TO SUM TO THE CROSS SECTION
-    if (evt_norm == "average") {
+    if (evt_norm != "sum") {
       normhepmc = 1. / double(iEventshower);
     } else {
       normhepmc = double(iEventtot) / double(iEventshower);

=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia8_hep.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia8_hep.cc	2017-07-27 05:51:14 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia8_hep.cc	2017-08-07 21:43:30 +0000
@@ -77,7 +77,7 @@
     double evtweight = pythia.info.weight();
     double normhepmc;
     // ALWAYS NORMALISE HEPMC WEIGHTS TO SUM TO THE CROSS SECTION
-    if (evt_norm == "average") {
+    if (evt_norm != "sum") {
       normhepmc = 1. / double(iEventshower);
     } else {
       normhepmc = double(iEventtot) / double(iEventshower);

=== modified file 'Template/NLO/Source/setrun.f'
--- Template/NLO/Source/setrun.f	2017-06-16 06:55:44 +0000
+++ Template/NLO/Source/setrun.f	2017-08-07 21:43:30 +0000
@@ -67,12 +67,14 @@
       buff = event_norm 
       call case_trap2(buff) ! requires a string of length 20 at least
       event_norm=buff 
-      if (event_norm(1:7).ne.'average' .and. event_norm(1:3).ne.'sum'
-     $     .and. event_norm(1:5).ne.'unity')then
+      if ( event_norm(1:7).ne.'average' .and.
+     $     event_norm(1:3).ne.'sum' .and.
+     $     event_norm(1:5).ne.'unity'.and.
+     $     event_norm(1:4).ne.'bias')then
          write (*,*) 'Do not understand the event_norm parameter'/
      &        /' in the run_card.dat. Possible options are'/
-     &        /' "average", "sum" or "unity". Current input is: ',
-     &        event_norm
+     &        /' "average", "sum", "unity" or "bias". '/
+     &        /'Current input is: ', event_norm
          stop 1
       endif
 c Check parameters for FxFx/UNLOPS/NNLL-veto

=== modified file 'Template/NLO/SubProcesses/collect_events.f'
--- Template/NLO/SubProcesses/collect_events.f	2017-07-29 08:02:08 +0000
+++ Template/NLO/SubProcesses/collect_events.f	2017-08-07 21:43:30 +0000
@@ -238,8 +238,6 @@
       data iseed/1/
       double precision rnd,fk88random
       external fk88random
-      logical debug
-      parameter (debug=.false.)
       integer nevents_file(80),proc_id(80)
       common /to_nevents_file/nevents_file,proc_id
       double precision xsecfrac_all(80)
@@ -256,11 +254,6 @@
       logical get_xsec_from_res1
       common/total_xsec/xsec,xerr,xsecABS,proc_id_tot,get_xsec_from_res1
 c
-      if(debug) then
-         write (*,*) ioutput,numoffiles,(junit(ii),ii=1,numoffiles)
-         write(ioutput,*)'test test test'
-         return
-      endif
       maxevt=0
       if (.not. get_xsec_from_res1) then
          do i=1,100

=== modified file 'Template/NLO/SubProcesses/cuts.f'
--- Template/NLO/SubProcesses/cuts.f	2017-05-09 12:37:07 +0000
+++ Template/NLO/SubProcesses/cuts.f	2017-08-07 21:43:30 +0000
@@ -898,32 +898,47 @@
       end
 
 
-      subroutine unweight_function(p_born,unwgtfun)
-c This is a user-defined function to which to unweight the events
-c A non-flat distribution will generate events with a certain
-c weight. This is particularly useful to generate more events
-c (with smaller weight) in tails of distributions.
-c It computes the unwgt factor from the momenta and multiplies
-c the weight that goes into MINT (or vegas) with this factor.
-c Before writing out the events (or making the plots), this factor
-c is again divided out.
-c This function should be called with the Born momenta to be sure
-c that it stays the same for the events, counter-events, etc.
-c A value different from 1 makes that MINT (or vegas) does not list
-c the correct cross section.
+      subroutine bias_weight_function(p,ipdg,bias_wgt)
+c This is a user-defined function to which to bias the event generation.
+c A non-flat distribution will generate events with a certain weight
+c inversely proportinal to the bias_wgt. This is particularly useful to
+c generate more events (with smaller weight) in tails of distributions.
+c It computes the bias_wgt factor from the momenta and multiplies the
+c weight that goes into MINT (or vegas) with this factor.  Before
+c writing out the events (or making the plots), this factor is again
+c divided out. A value different from 1 makes that MINT (or vegas) does
+c not list the correct cross section, but the cross section can still be
+c computed from summing all the weights of the events (and dividing by
+c the number of events). Since the weights of the events are no longer
+c identical for all events, the statistical uncertainty on this total
+c cross section can be much larger than without including the bias.
+c
+c The 'bias_wgt' should be a IR-safe function of the momenta.
+c      
+c For this to be used, the 'event_norm' option in the run_card should be
+c set to
+c      'bias' = event_norm      
+c
       implicit none
       include 'nexternal.inc'
-      double precision unwgtfun,p_born(0:3,nexternal-1),shat,sumdot
-      external sumdot
-
-      unwgtfun=1d0
-
-c How to enhance the tails is very process dependent. But, it is
-c probably easiest to enhance the tails using shat, e.g.:
-c      shat=sumdot(p_born(0,1),p_born(0,2),1d0)
-c      unwgtfun=max(100d0**2,shat)/100d0**2
-c      unwgtfun=unwgtfun**2
-
+      double precision bias_wgt,p(0:3,nexternal),H_T
+      integer ipdg(nexternal),i
+
+      bias_wgt=1d0
+
+c How to enhance the tails is very process dependent. For example for
+c top quark production one could use:
+c      do i=1,nexternal
+c         if (ipdg(i).eq.6) then
+c            bias_wgt=sqrt(p(1,i)**2+p(2,i)**2)**3
+c         endif
+c      enddo
+c Or to use H_T^2 one does     
+c      H_T=0d0
+c      do i=1,nexternal
+c         H_T=H_T+sqrt(max(0d0,(pp(0,i)+pp(3,i))*(pp(0,i)-pp(3,i))))
+c      enddo
+c      bias_wgt=H_T**2
       return
       end
 

=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2017-07-29 08:02:08 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2017-08-07 21:43:30 +0000
@@ -50,7 +50,9 @@
       double precision average_virtual(maxchannels),virtual_fraction(maxchannels)
       common/c_avg_virt/average_virtual,virtual_fraction
 
-      double precision weight
+      double precision weight,event_weight,inv_bias
+      character*7 event_norm
+      common /event_normalisation/event_norm
 c For MINT:
       real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
      $     ,ndimmax,maxchannels),ymax_virt(maxchannels),ans(nintegrals
@@ -311,7 +313,11 @@
          if (ickkw.eq.-1) putonshell=.false.
          unwgt=.true.
          open (unit=99,file='nevts',status='old',err=999)
-         read (99,*) nevts
+         if (event_norm(1:4).ne.'bias') then
+            read (99,*) nevts
+         else
+            read (99,*) nevts,event_weight
+         endif
          close(99)
          write(*,*) 'Generating ', nevts, ' events'
          if(nevts.eq.0) then
@@ -354,7 +360,11 @@
          absint=ans(1,1)+ans(5,1)
          uncer=unc(2,1)
 
-         weight=(ans(1,1)+ans(5,1))/ncall
+         if (event_norm(1:4).ne.'bias') then
+            weight=(ans(1,1)+ans(5,1))/ncall
+         else
+            weight=event_weight
+         endif
 
          if (abrv(1:3).ne.'all' .and. abrv(1:4).ne.'born' .and.
      $        abrv(1:4).ne.'virt') then
@@ -389,6 +399,10 @@
             call pick_unweight_contr(iFKS_picked)
             call update_fks_dir(iFKS_picked)
             call fill_rwgt_lines
+            if (event_norm(1:4).eq.'bias') then
+               call include_inverse_bias_wgt(inv_bias)
+               weight=event_weight*inv_bias
+            endif
             call finalize_event(x,weight,lunlhe,putonshell)
          enddo
          call deallocate_weight_lines
@@ -945,9 +959,11 @@
             call include_shape_in_shower_scale(p,iFKS)
          enddo
  12      continue
-
+         
 c Include PDFs and alpha_S and reweight to include the uncertainties
          call include_PDF_and_alphas
+c Include the weight from the bias_function
+         call include_bias_wgt
 c Sum the contributions that can be summed before taking the ABS value
          call sum_identical_contributions
 c Update the shower starting scale for the S-events after we have

=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2017-06-19 06:49:26 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2017-08-07 21:43:30 +0000
@@ -619,7 +619,7 @@
       include 'run.inc'
       include 'genps.inc'
       include 'timing_variables.inc'
-      double precision pi,unwgtfun,vegas_wgt,enhance,xnoborn_cnt,xtot
+      double precision pi,vegas_wgt,enhance,xnoborn_cnt,xtot
      $     ,bpower,cpower,tiny
       data xnoborn_cnt /0d0/
       integer inoborn_cnt,i
@@ -732,11 +732,10 @@
             enhance=0d0
          endif
       endif
-      call unweight_function(p_born,unwgtfun)
       call set_cms_stuff(0)
 c f_* multiplication factors for Born and nbody
       f_b=jac_cnt(0)*xinorm_ev/(min(xiimax_ev,xiBSVcut_used)*shat/(16
-     $     *pi**2))*enhance*unwgtfun *fkssymmetryfactorBorn*vegas_wgt
+     $     *pi**2))*enhance*fkssymmetryfactorBorn*vegas_wgt
       f_nb=f_b
       call cpu_time(tAfter)
       tf_nb=tf_nb+(tAfter-tBefore)
@@ -753,7 +752,7 @@
       include 'fks_powers.inc'
       include 'coupl.inc'
       include 'timing_variables.inc'
-      double precision unwgtfun,vegas_wgt,enhance,xnoborn_cnt,xtot
+      double precision vegas_wgt,enhance,xnoborn_cnt,xtot
      &     ,prefact,prefact_cnt_ssc,prefact_deg,prefact_c,prefact_coll
      &     ,jac_ev,pi,prefact_cnt_ssc_c,prefact_coll_c,prefact_deg_slxi
      &     ,prefact_deg_sxi,zero
@@ -845,11 +844,10 @@
             enhance=0d0
          endif
       endif
-      call unweight_function(p_born,unwgtfun)
       prefact=xinorm_ev/xi_i_fks_ev/(1-y_ij_fks_ev)
 
 c f_* multiplication factors for real-emission, soft counter, ... etc.       
-      f_r=prefact*jac_ev*enhance*unwgtfun*fkssymmetryfactor*vegas_wgt
+      f_r=prefact*jac_ev*enhance*fkssymmetryfactor*vegas_wgt
       f_MC_S=f_r
       f_MC_H=f_r
       if (.not.nocntevents) then
@@ -857,9 +855,9 @@
      &        log(xicut_used/min(xiimax_ev,xiScut_used))/(1
      &        -y_ij_fks_ev)
          f_s=(prefact+prefact_cnt_ssc)*jac_cnt(0)*enhance
-     $        *unwgtfun*fkssymmetryfactor*vegas_wgt
+     $        *fkssymmetryfactor*vegas_wgt
          f_s_MC_S=prefact*jac_cnt(0)*enhance
-     $        *unwgtfun*fkssymmetryfactor*vegas_wgt
+     $        *fkssymmetryfactor*vegas_wgt
          f_s_MC_H=f_s_MC_S
 
          if (pmass(j_fks).eq.0d0) then
@@ -867,9 +865,9 @@
             prefact_coll=xinorm_cnt(1)/xi_i_fks_cnt(1)*log(delta_used
      $           /deltaS)/deltaS
             f_c=(prefact_c+prefact_coll)*jac_cnt(1)
-     $           *enhance*unwgtfun*fkssymmetryfactor*vegas_wgt
+     $           *enhance*fkssymmetryfactor*vegas_wgt
             f_c_MC_S=prefact_c*jac_cnt(1)
-     $           *enhance*unwgtfun*fkssymmetryfactor*vegas_wgt
+     $           *enhance*fkssymmetryfactor*vegas_wgt
             f_c_MC_H=f_c_MC_S
 
             call set_cms_stuff(1)
@@ -881,12 +879,12 @@
      $           *log(xicut_used/min(xiimax_cnt(1),xiScut_used))
      $           *log(delta_used/deltaS)/deltaS
             f_dc=jac_cnt(1)*prefact_deg/(shat/(32*pi**2))*enhance
-     $           *unwgtfun*fkssymmetryfactorDeg*vegas_wgt
+     $           *fkssymmetryfactorDeg*vegas_wgt
             f_sc=(prefact_c+prefact_coll+prefact_cnt_ssc_c
-     &           +prefact_coll_c)*jac_cnt(2)*enhance*unwgtfun
+     &           +prefact_coll_c)*jac_cnt(2)*enhance
      &           *fkssymmetryfactorDeg*vegas_wgt
             f_sc_MC_S=prefact_c*jac_cnt(2)
-     $           *enhance*unwgtfun*fkssymmetryfactor*vegas_wgt
+     $           *enhance*fkssymmetryfactor*vegas_wgt
             f_sc_MC_H=f_sc_MC_S
 
             call set_cms_stuff(2)
@@ -898,13 +896,13 @@
      &           -log(min(xiimax_cnt(1),xiScut_used))**2 )*1/(2.d0
      &           *deltaS)
             f_dsc(1)=prefact_deg*jac_cnt(2)/(shat/(32*pi**2))*enhance
-     &           *unwgtfun*fkssymmetryfactorDeg*vegas_wgt
+     &           *fkssymmetryfactorDeg*vegas_wgt
             f_dsc(2)=prefact_deg_sxi*jac_cnt(2)/(shat/(32*pi**2))
-     &           *enhance*unwgtfun*fkssymmetryfactorDeg*vegas_wgt
+     &           *enhance*fkssymmetryfactorDeg*vegas_wgt
             f_dsc(3)=prefact_deg_slxi*jac_cnt(2)/(shat/(32*pi**2))
-     &           *enhance*unwgtfun*fkssymmetryfactorDeg*vegas_wgt
+     &           *enhance*fkssymmetryfactorDeg*vegas_wgt
             f_dsc(4)=( prefact_deg+prefact_deg_sxi )*jac_cnt(2)/(shat
-     &           /(32*pi**2))*enhance*unwgtfun*fkssymmetryfactorDeg
+     &           /(32*pi**2))*enhance*fkssymmetryfactorDeg
      &           *vegas_wgt
          else
             f_c=0d0
@@ -1214,6 +1212,117 @@
       return
       end
 
+      subroutine include_bias_wgt
+c Include the weight from the bias_wgt_function to all the contributions
+c in icontr. This only changes the weight of the central value (after
+c inclusion of alphaS and parton luminosity). Both for 'wgts(1,icontr)'
+c as well as the the 'parton_iproc(1:niproc(icontr),icontr)', since
+c these are the ones used in MINT as well as for unweighting. Also the
+c 'virt_wgt_mint' and 'born_wgt_mint' are updated. Furthermore, to
+c include the weight also in the 'wgt' array that contain the
+c coefficients for PDF and scale computations. 
+      use weight_lines
+      implicit none
+      integer i,j
+      double precision bias
+      character*7 event_norm
+      common /event_normalisation/event_norm
+      double precision           virt_wgt_mint,born_wgt_mint
+      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
+c Set the bias_wgt to 1 in case we do not have to do any biassing
+      if (event_norm(1:4).ne.'bias') then
+         do i=1,icontr
+            bias_wgt(i)=1d0
+         enddo
+         return
+      endif
+c loop over all contributions
+      do i=1,icontr
+         if (itype(i).eq.1) then
+            ! use (n+1)-body momenta for the real emission. Pick the
+            ! first IPROC for parton PDGs.
+            call bias_weight_function(momenta_m(0,1,2,i),parton_pdg(1,1
+     $           ,i),bias)
+         else
+            ! use n-body momenta for all the other contributions. Pick
+            ! the first IPROC for parton PDGs.
+            call bias_weight_function(momenta_m(0,1,1,i),parton_pdg(1,1
+     $           ,i),bias)
+         endif
+         bias_wgt(i)=bias
+c Update the weights:
+         wgts(1,i)=wgts(1,i)*bias_wgt(i)
+         do j=1,niproc(i)
+            parton_iproc(j,i)=parton_iproc(j,i)*bias_wgt(i)
+         enddo
+         do j=1,3
+            wgt(j,i)=wgt(j,i)*bias_wgt(i)
+         enddo
+         if (itype(i).eq.14) then
+            virt_wgt_mint=virt_wgt_mint*bias_wgt(i)
+            born_wgt_mint=born_wgt_mint*bias_wgt(i)
+         endif
+      enddo
+      return
+      end
+
+      subroutine include_inverse_bias_wgt(inv_bias)
+c Update the inverse of the bias in the event weight. All information in
+c the rwgt_lines is NOT updated.
+      use weight_lines
+      use extra_weights
+      implicit none
+      include 'genps.inc'
+      include 'nFKSconfigs.inc'
+      integer i,ict,ipr,ii
+      double precision wgt_num,wgt_denom,inv_bias
+      character*7 event_norm
+      common /event_normalisation/event_norm
+      integer iproc_save(fks_configs),eto(maxproc,fks_configs)
+     $     ,etoi(maxproc,fks_configs),maxproc_found
+      common/cproc_combination/iproc_save,eto,etoi,maxproc_found
+      logical         Hevents
+      common/SHevents/Hevents
+      if (event_norm(1:4).ne.'bias') then
+         inv_bias=1d0
+         return
+      endif
+      wgt_num=0d0
+      wgt_denom=0d0
+      do i=1,icontr_sum(0,icontr_picked)
+         ict=icontr_sum(i,icontr_picked)
+         if (bias_wgt(ict).eq.0d0) then
+            write (*,*) "ERROR in include_inverse_bias_wgt: "/
+     $           /"bias_wgt is equal to zero",ict,bias_wgt
+            stop 1
+         endif
+c for all the rwgt_lines, remove the bias-wgt contribution from the
+c weights there. Note that the wgtref (also written in the event file)
+c keeps its contribution from the bias_wgt.
+         if (.not. Hevents) then
+            ipr=eto(etoi(iproc_picked,nFKS(ict)),nFKS(ict))
+            do ii=1,iproc_save(nFKS(ict))
+               if (eto(ii,nFKS(ict)).ne.ipr) cycle
+               wgt_denom=wgt_denom+parton_iproc(ii,ict)
+               wgt_num=wgt_num+parton_iproc(ii,ict)/bias_wgt(ict)
+            enddo
+         else
+            ipr=iproc_picked
+            wgt_denom=wgt_denom+parton_iproc(ipr,ict)
+            wgt_num=wgt_num+parton_iproc(ipr,ict)/bias_wgt(ict)
+         endif
+      enddo
+      if (abs((wgtref-wgt_denom)/(wgtref+wgt_denom)).gt.1d-10) then
+         write (*,*) "ERROR in include_inverse_bias_wgt: "/
+     $        /"reference weight not equal to recomputed weight",wgtref
+     $        ,wgt_denom
+         stop 1
+      endif
+c update the event weight to be written in the file
+      inv_bias=wgt_num/wgt_denom
+      return
+      end
+      
 
       subroutine set_pdg_codes(iproc,pd,iFKS,ict)
       use weight_lines
@@ -1748,14 +1857,27 @@
       include 'nexternal.inc'
       include 'mint.inc'
       integer i
-      double precision f(nintegrals),sigint
+      double precision f(nintegrals),sigint,bias
       double precision           virt_wgt_mint,born_wgt_mint
       common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
       double precision virtual_over_born
       common /c_vob/   virtual_over_born
+      character*7 event_norm
+      common /event_normalisation/event_norm
       sigint=0d0
       do i=1,icontr
-         sigint=sigint+wgts(1,i)
+         if (event_norm(1:4).ne.'bias') then
+            sigint=sigint+wgts(1,i)
+         else
+            if (itype(i).eq.1) then
+               call bias_weight_function(momenta_m(0,1,2,i),parton_pdg(1
+     $              ,1,i),bias)
+            else
+               call bias_weight_function(momenta_m(0,1,1,i),parton_pdg(1
+     $              ,1,i),bias)
+            endif
+            sigint=sigint+wgts(1,i)*bias
+         endif
       enddo
       f(1)=abs(sigint)
       f(2)=sigint
@@ -2193,7 +2315,7 @@
       n_ctr_found=0
       n_mom_conf=0
 c Loop over all the contributions in the picked contribution (the latter
-c is chosen in the pick_unweight_cont() subroutine)
+c is chosen in the pick_unweight_contr() subroutine)
       do i=1,icontr_sum(0,icontr_picked)
          ict=icontr_sum(i,icontr_picked)
 c Check if the current set of momenta are already available in the
@@ -2270,8 +2392,7 @@
      &              trim(adjustl(n_ctr_str(n_ctr_found)))//' '
      &              //trim(adjustl(procid))
 
-               write (str_temp,
-     &                    '(i2,6(1x,d14.8),6(1x,i2),1x,i8,1x,d18.12)')
+               write (str_temp,30)
      &              QCDpower(ict),
      &              (bjx(j,ict),j=1,2),
      &              (scales2(j,ict),j=1,3),
@@ -2282,7 +2403,8 @@
      &              fks_i_d(nFKS(ict)),
      &              fks_j_d(nFKS(ict)),
      &              parton_pdg_uborn(fks_j_d(nFKS(ict)),ii,ict),
-     &              parton_iproc(ii,ict)
+     &              parton_iproc(ii,ict),
+     &              bias_wgt(ict)
                n_ctr_str(n_ctr_found) =
      &              trim(adjustl(n_ctr_str(n_ctr_found)))//' '
      &              //trim(adjustl(str_temp))
@@ -2321,7 +2443,7 @@
      &           trim(adjustl(n_ctr_str(n_ctr_found)))//' '
      &           //trim(adjustl(procid))
 
-            write (str_temp,'(i2,6(1x,d14.8),6(1x,i2),1x,i8,1x,d18.12)')
+            write (str_temp,30)
      &           QCDpower(ict),
      &           (bjx(j,ict),j=1,2),
      &           (scales2(j,ict),j=1,3),
@@ -2332,7 +2454,8 @@
      &           fks_i_d(nFKS(ict)),
      &           fks_j_d(nFKS(ict)),
      &           parton_pdg_uborn(fks_j_d(nFKS(ict)),ipr,ict),
-     &           parton_iproc(ipr,ict)
+     &           parton_iproc(ipr,ict),
+     &           bias_wgt(ict)
             n_ctr_str(n_ctr_found) =
      &           trim(adjustl(n_ctr_str(n_ctr_found)))//' '
      &           //trim(adjustl(str_temp))
@@ -2340,6 +2463,8 @@
 
          endif
       enddo
+      return
+ 30   format(i2,6(1x,d14.8),6(1x,i2),1x,i8,1x,d18.12,1x,d18.12)
       end
       
       

=== modified file 'Template/NLO/SubProcesses/reweight_xsec_events.f'
--- Template/NLO/SubProcesses/reweight_xsec_events.f	2017-06-15 16:22:23 +0000
+++ Template/NLO/SubProcesses/reweight_xsec_events.f	2017-08-07 21:43:30 +0000
@@ -394,7 +394,7 @@
          read(n_ctr_str(i),*)(wgt(j,i),j=1,3),(wgt_ME_tree(j,i),j=1,2)
      &        ,idum,(pdg(j,i),j=1,nexternal),QCDpower(i),(bjx(j,i),j=1
      &        ,2),(scales2(j,i),j=1,3),g_strong(i),(momenta_conf(j),j=1
-     &        ,2),itype(i),nFKS(i),idum,idum,idum,wgts(1,i)
+     &        ,2),itype(i),nFKS(i),idum,idum,idum,wgts(1,i),bias_wgt(i)
          do ii=1,2
             do j=1,nexternal
                do k=0,3

=== modified file 'Template/NLO/SubProcesses/symmetry_fks_v3.f'
--- Template/NLO/SubProcesses/symmetry_fks_v3.f	2017-06-16 14:41:09 +0000
+++ Template/NLO/SubProcesses/symmetry_fks_v3.f	2017-08-07 21:43:30 +0000
@@ -313,12 +313,13 @@
       RETURN
       END
 
-      subroutine unweight_function(p_born,unwgtfun)
+      subroutine bias_weight_function(p,ipdg,bias_wgt)
 c Dummy function. Should always retrun 1.
       implicit none
       include 'nexternal.inc'
-      double precision unwgtfun,p_born(0:3,nexternal-1)
-      unwgtfun=1d0
+      integer ipdg(nexternal)
+      double precision bias_wgt,p(0:3,nexternal)
+      bias_wgt=1d0
       return
       end
 

=== modified file 'Template/NLO/SubProcesses/weight_lines.f'
--- Template/NLO/SubProcesses/weight_lines.f	2017-06-13 17:00:25 +0000
+++ Template/NLO/SubProcesses/weight_lines.f	2017-08-07 21:43:30 +0000
@@ -11,7 +11,7 @@
          double precision, allocatable :: momenta(:,:,:),momenta_m(:,:,:
      $        ,:),wgt(:,:),wgt_ME_tree(:,:),bjx(:,:),scales2(:,:)
      $        ,g_strong(:),wgts(:,:),parton_iproc(:,:),y_bst(:)
-     $        ,plot_wgts(:,:),shower_scale(:),unwgt(:,:)
+     $        ,plot_wgts(:,:),shower_scale(:),unwgt(:,:),bias_wgt(:)
          save
       end module weight_lines
 
@@ -157,6 +157,10 @@
          allocate(temp1(n_contr))
          temp1(1:max_contr)=y_bst
          call move_alloc(temp1,y_bst)
+c bias_wgt
+         allocate(temp1(n_contr))
+         temp1(1:max_contr)=bias_wgt
+         call move_alloc(temp1,bias_wgt)
 c plot_wgts
          allocate(temp2(max_wgt,n_contr))
          temp2(1:max_wgt,1:max_contr)=plot_wgts
@@ -201,6 +205,7 @@
       allocate(wgts(1,1))
       allocate(parton_iproc(1,1))
       allocate(y_bst(1))
+      allocate(bias_wgt(1))
       allocate(plot_wgts(1,1))
       allocate(shower_scale(1))
       allocate(unwgt(1,1))
@@ -238,6 +243,7 @@
       if (allocated(wgts)) deallocate(wgts)
       if (allocated(parton_iproc)) deallocate(parton_iproc)
       if (allocated(y_bst)) deallocate(y_bst)
+      if (allocated(bias_wgt)) deallocate(bias_wgt)
       if (allocated(plot_wgts)) deallocate(plot_wgts)
       if (allocated(shower_scale)) deallocate(shower_scale)
       if (allocated(unwgt)) deallocate(unwgt)

=== modified file 'Template/NLO/SubProcesses/write_event.f'
--- Template/NLO/SubProcesses/write_event.f	2017-06-16 08:20:04 +0000
+++ Template/NLO/SubProcesses/write_event.f	2017-08-07 21:43:30 +0000
@@ -13,7 +13,7 @@
       integer i,j,lunlhe
       real*8 xx(ndimmax),weight,evnt_wgt
       logical putonshell
-      double precision wgt,unwgtfun
+      double precision wgt
       double precision x(99),p(0:3,nexternal)
       integer jpart(7,-nexternal+3:2*nexternal-3)
       double precision pb(0:4,-nexternal+3:2*nexternal-3)
@@ -82,14 +82,6 @@
       call add_write_info(p_born,p,ybst_til_tolab,iconfig,Hevents,
      &     putonshell,ndim,x,jpart,npart,pb,shower_scale)
 
-      call unweight_function(p_born,unwgtfun)
-      if (unwgtfun.ne.0d0) then
-         evnt_wgt=evnt_wgt/unwgtfun
-      else
-         write (*,*) 'ERROR in finalize_event, unwgtfun=0',unwgtfun
-         stop
-      endif
-
 c  Write-out the events
       call write_events_lhe(pb(0,1),evnt_wgt,jpart(1,1),npart,lunlhe
      $     ,shower_scale,ickkw)

=== modified file 'Template/NLO/Utilities/check_events.f'
--- Template/NLO/Utilities/check_events.f	2014-02-19 09:46:05 +0000
+++ Template/NLO/Utilities/check_events.f	2017-08-07 21:43:30 +0000
@@ -6,10 +6,12 @@
 c With some work on finalizeprocesses(), it should work also for 
 c LH files created by Herwig, assuming they are identified by a 
 c negative number of events
+      use extra_weights
       implicit none
       integer maxevt,ifile,efile,mfile,jfile,kfile,rfile,i,npart,
      # iuseres_1,iwhmass,ilepmass,idec,itempsc,itempPDF,isavesc,
      # isavePDF,itemp,ii
+      integer numscales,numPDFpairs,numPDFs
       double precision chtot,xint,xinterr,xinta,xintaerr,qtiny
       parameter (qtiny=1.d-4)
       double precision charges(-100:100),zmasses(1:100)
@@ -42,7 +44,8 @@
       double precision wgt3a,wgt3s
       double precision wgt4a,wgt4s
       double precision wgt5a,wgt5s
-      double precision saved_weight,tmp
+      double precision saved_weight,tmp,wmin,wmax,wlim
+      integer ipos,ineg
       character*80 event_file
       character*140 buff
       character*6 ch6
@@ -59,14 +62,14 @@
       integer j,k
       real*8 ecm,xmass(3*nexternal),xmom(0:3,3*nexternal),xnorm
 
-      include 'reweight0.inc'
-      integer kr,kf,kpdf
-      double precision sum_wgt_resc_scale(maxscales,maxscales),
-     # sum_wgt_resc_pdf(0:maxPDFs),
-     # xmax_wgt_resc_scale(maxscales,maxscales),
-     # xmax_wgt_resc_pdf(0:maxPDFs),
-     # xmin_wgt_resc_scale(maxscales,maxscales),
-     # xmin_wgt_resc_pdf(0:maxPDFs)
+      integer kk,kr,kf,kpdf
+      double precision 
+     # sum_wgt_resc_scale(maxscales,maxscales,maxdynscales),
+     # sum_wgt_resc_pdf(0:maxPDFs,maxPDFsets),
+     # xmax_wgt_resc_scale(maxscales,maxscales,maxdynscales),
+     # xmax_wgt_resc_pdf(0:maxPDFs,maxPDFsets),
+     # xmin_wgt_resc_scale(maxscales,maxscales,maxdynscales),
+     # xmin_wgt_resc_pdf(0:maxPDFs,maxPDFsets)
       integer istep
       double precision percentage
       include 'dbook.inc'
@@ -144,6 +147,8 @@
       unweighted=.true.
       keepevent=.true.
       shower=.false.
+      wmax=-1.d100
+      wmin=1.d100
 
       call read_lhef_header_full(ifile,maxevt,itempsc,itempPDF,
      #                           MonteCarlo)
@@ -220,7 +225,7 @@
      &     XSECUP,XERRUP,XMAXUP,LPRUP)
 
       
-      do i=1,min(10,maxevt)
+      do i=1,min(100,maxevt)
         call read_lhef_event_catch(ifile,
      &       NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
      &       IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
@@ -247,10 +252,12 @@
               write(*,*)'Inconsistency #2 in event file',i,' ',buff
               stop
             endif
-          unweighted=unweighted.and.
+            unweighted=unweighted.and.
      #               abs(1.d0-abs(XWGTUP)/saved_weight).lt.1.d-5
           endif
         endif
+        wmax=max(wmax,XWGTUP)
+        wmin=min(wmin,XWGTUP)
         if(idec.eq.0)call findmothers(NUP,IDUP,ISTUP)
 
       enddo
@@ -263,6 +270,11 @@
         write(*,*)'The events appear to be unweighted'
       else
         write(*,*)'The events appear to be weighted'
+        if(event_norm.eq.'ave'.or.event_norm.eq.'bia')then
+          wmax=wmax/maxevt
+          wmin=wmin/maxevt
+        endif
+        wlim=max(abs(wmax),abs(wmin))
       endif
 
       if(rwgtinfo)then
@@ -304,25 +316,40 @@
       itoterr=0
       mtoterr=0
       if(jwgtinfo.eq.8.or.jwgtinfo.eq.9)then
-        do kr=1,maxscales
-          do kf=1,maxscales
-            sum_wgt_resc_scale(kr,kf)=0.d0
-            xmax_wgt_resc_scale(kr,kf)=-1.d100
-            xmin_wgt_resc_scale(kr,kf)=1.d100
+        do kk=1,maxdynscales
+          do kr=1,maxscales
+            do kf=1,maxscales
+              sum_wgt_resc_scale(kr,kf,kk)=0.d0
+              xmax_wgt_resc_scale(kr,kf,kk)=-1.d100
+              xmin_wgt_resc_scale(kr,kf,kk)=1.d100
+            enddo
           enddo
         enddo
-        do kpdf=0,maxPDFs
-          sum_wgt_resc_pdf(kpdf)=0.d0
-          xmax_wgt_resc_pdf(kpdf)=-1.d100
-          xmin_wgt_resc_pdf(kpdf)=1.d100
+        do kk=1,maxPDFsets
+          do kpdf=0,maxPDFs
+            sum_wgt_resc_pdf(kpdf,kk)=0.d0
+            xmax_wgt_resc_pdf(kpdf,kk)=-1.d100
+            xmin_wgt_resc_pdf(kpdf,kk)=1.d100
+          enddo
         enddo
       endif
 
       i=0
+      ipos=0
+      ineg=0
       call inihist
       call bookup(1,'scalup',2d0,0d0,200d0)
       call bookup(2,'scalup',8d0,0d0,800d0)
       call bookup(3,'log scalup',0.1d0,0d0,4d0)
+      if(.not.unweighted)then
+        tmp=2*wlim/100.d0
+        call bookup(4,'weight',tmp,-wlim,wlim)
+        tmp=4*wlim/100.d0
+        call bookup(5,'weight',tmp,-2*wlim,2*wlim)
+        tmp=8*wlim/100.d0
+        call bookup(6,'weight',tmp,-4*wlim,4*wlim)
+      endif
+
       dowhile(i.lt.maxevt.and.keepevent)
          call read_lhef_event_catch(ifile,
      &        NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
@@ -330,6 +357,17 @@
          call mfill(1,SCALUP,XWGTUP)
          call mfill(2,SCALUP,XWGTUP)
          if(SCALUP.gt.0d0)call mfill(3,log10(SCALUP),XWGTUP)
+         if(.not.unweighted)then
+           if(event_norm.eq.'ave'.or.event_norm.eq.'bia')then
+             call mfill(4,XWGTUP/maxevt,abs(XWGTUP))
+             call mfill(5,XWGTUP/maxevt,abs(XWGTUP))
+             call mfill(6,XWGTUP/maxevt,abs(XWGTUP))
+           else
+             call mfill(4,XWGTUP,abs(XWGTUP))
+             call mfill(5,XWGTUP,abs(XWGTUP))
+             call mfill(6,XWGTUP,abs(XWGTUP))
+           endif
+         endif
          if(index(buff,'endoffile').ne.0)then
            keepevent=.false.
            goto 111
@@ -338,32 +376,41 @@
          i=i+1
          sum_wgt=sum_wgt+XWGTUP
          sum_abs_wgt=sum_abs_wgt+abs(XWGTUP)
+         if(sign(1.d0,XWGTUP).gt.0.d0)then
+           ipos=ipos+1
+         else
+           ineg=ineg+1
+         endif
 
 c Note: with pre-beta2 convention, the reweighting cross sections were
 c normalized such that one needed to compute e.g. 
 c XWGTUP*wgtxsecmu(kr,kf)/wgtref
          if(jwgtinfo.eq.8.or.jwgtinfo.eq.9)then
-           do kr=1,numscales
-             do kf=1,numscales
-               sum_wgt_resc_scale(kr,kf)=sum_wgt_resc_scale(kr,kf)+
-     #                                   wgtxsecmu(kr,kf)
-               xmax_wgt_resc_scale(kr,kf)=
-     #           max(xmax_wgt_resc_scale(kr,kf),
-     #               wgtxsecmu(kr,kf))
-               xmin_wgt_resc_scale(kr,kf)=
-     #           min(xmin_wgt_resc_scale(kr,kf),
-     #               wgtxsecmu(kr,kf))
+           do kk=1,dyn_scale(0)
+             do kr=1,numscales
+               do kf=1,numscales
+                 sum_wgt_resc_scale(kr,kf,kk)=
+     #              sum_wgt_resc_scale(kr,kf,kk)+wgtxsecmu(kr,kf,kk)
+                 xmax_wgt_resc_scale(kr,kf,kk)=
+     #             max(xmax_wgt_resc_scale(kr,kf,kk),
+     #               wgtxsecmu(kr,kf,kk))
+                 xmin_wgt_resc_scale(kr,kf,kk)=
+     #             min(xmin_wgt_resc_scale(kr,kf,kk),
+     #               wgtxsecmu(kr,kf,kk))
+               enddo
              enddo
            enddo
-           do kpdf=1,2*numPDFpairs
-             sum_wgt_resc_pdf(kpdf)=sum_wgt_resc_pdf(kpdf)+
-     #                              wgtxsecPDF(kpdf)
-             xmax_wgt_resc_pdf(kpdf)=
-     #         max(xmax_wgt_resc_pdf(kpdf),
-     #             wgtxsecPDF(kpdf))
-             xmin_wgt_resc_pdf(kpdf)=
-     #         min(xmin_wgt_resc_pdf(kpdf),
-     #             wgtxsecPDF(kpdf))
+           do kk=1,lhaPDFid(0)
+             do kpdf=1,2*numPDFpairs
+               sum_wgt_resc_pdf(kpdf,kk)=sum_wgt_resc_pdf(kpdf,kk)+
+     #                                wgtxsecPDF(kpdf,kk)
+               xmax_wgt_resc_pdf(kpdf,kk)=
+     #           max(xmax_wgt_resc_pdf(kpdf,kk),
+     #             wgtxsecPDF(kpdf,kk))
+               xmin_wgt_resc_pdf(kpdf,kk)=
+     #           min(xmin_wgt_resc_pdf(kpdf,kk),
+     #             wgtxsecPDF(kpdf,kk))
+             enddo
            enddo
          endif
 
@@ -476,18 +523,20 @@
              tmp=wgtcentral/XWGTUP
              wgt1a=wgt1a+tmp
              wgt1s=wgt1s+tmp**2
-             tmp=wgtmumin/wgtcentral
-             wgt2a=wgt2a+tmp
-             wgt2s=wgt2s+tmp**2
-             tmp=wgtmumax/wgtcentral
-             wgt3a=wgt3a+tmp
-             wgt3s=wgt3s+tmp**2
-             tmp=wgtpdfmin/wgtcentral
-             wgt4a=wgt4a+tmp
-             wgt4s=wgt4s+tmp**2
-             tmp=wgtpdfmax/wgtcentral
-             wgt5a=wgt5a+tmp
-             wgt5s=wgt5s+tmp**2
+             if(wgtcentral.ne.0.d0)then
+               tmp=wgtmumin/wgtcentral
+               wgt2a=wgt2a+tmp
+               wgt2s=wgt2s+tmp**2
+               tmp=wgtmumax/wgtcentral
+               wgt3a=wgt3a+tmp
+               wgt3s=wgt3s+tmp**2
+               tmp=wgtpdfmin/wgtcentral
+               wgt4a=wgt4a+tmp
+               wgt4s=wgt4s+tmp**2
+               tmp=wgtpdfmax/wgtcentral
+               wgt5a=wgt5a+tmp
+               wgt5s=wgt5s+tmp**2
+             endif
            endif
          endif
 
@@ -505,6 +554,11 @@
  111     continue
       enddo
 
+      if(maxevt.ne.(ipos+ineg))then
+        write(*,*)'Something wrong with counting events:',
+     #             maxevt,ipos,ineg
+      endif
+
       open(unit=99,file='SCALUP.top',status='unknown')
       call mclear
       xnorm=1d0
@@ -515,13 +569,19 @@
       call multitop(1,3,2,'scalup    ',' ','LOG')
       call multitop(2,3,2,'scalup    ',' ','LOG')
       call multitop(3,3,2,'log scalup',' ','LOG')
+      call multitop(4,3,2,'weight    ',' ','LOG')
+      call multitop(5,3,2,'weight    ',' ','LOG')
+      call multitop(6,3,2,'weight    ',' ','LOG')
       close(99)
 
-      if(event_norm.eq.'ave')sum_wgt=sum_wgt/maxevt
-      if(event_norm.eq.'ave')sum_abs_wgt=sum_abs_wgt/maxevt
+      if(event_norm.eq.'ave'.or.event_norm.eq.'bia')
+     #   sum_wgt=sum_wgt/maxevt
+      if(event_norm.eq.'ave'.or.event_norm.eq.'bia')
+     #   sum_abs_wgt=sum_abs_wgt/maxevt
       err_wgt=sum_abs_wgt/sqrt(dfloat(maxevt))
       write(*,*)'  '
       write (*,*) 'The total number of events is:',i
+      write (*,*) ' of which:',ipos,' w>0',ineg,' w<0'
       write (*,*) 'Sum of weights is    :',sum_wgt,' +-',err_wgt
       write (*,*) 'Sum of abs weights is:',sum_abs_wgt,' +-',err_wgt
 
@@ -586,56 +646,68 @@
         if(jwgtinfo.eq.8.or.jwgtinfo.eq.9)then
           write(64,*)'  '
           write(64,*)'Sums of rescaled weights'
-          do kr=1,numscales
-            do kf=1,numscales
-              if(event_norm.eq.'ave')then
-                write(64,300)'scales',kr,kf,' ->',
-     #                     sum_wgt_resc_scale(kr,kf)/maxevt
-              else
-                write(64,300)'scales',kr,kf,' ->',
-     #                     sum_wgt_resc_scale(kr,kf)
-              endif
+          do kk=1,dyn_scale(0)
+            do kr=1,numscales
+              do kf=1,numscales
+                if(event_norm.eq.'ave'.or.event_norm.eq.'bia')then
+                  write(64,300)'scales',kk,kr,kf,' ->',
+     #                       sum_wgt_resc_scale(kr,kf,kk)/maxevt
+                else
+                  write(64,300)'scales',kk,kr,kf,' ->',
+     #                       sum_wgt_resc_scale(kr,kf,kk)
+                endif
+              enddo
             enddo
           enddo
-          if(event_norm.eq.'ave')then
-            do kpdf=1,2*numPDFpairs
-              write(64,301)'PDF',kpdf,' ->',
-     #                   sum_wgt_resc_pdf(kpdf)/maxevt
+          if(event_norm.eq.'ave'.or.event_norm.eq.'bia')then
+            do kk=1,lhaPDFid(0)
+              do kpdf=1,2*numPDFpairs
+                write(64,301)'PDF',kk,kpdf,' ->',
+     #                     sum_wgt_resc_pdf(kpdf,kk)/maxevt
+              enddo
             enddo
           else
-            do kpdf=1,2*numPDFpairs
-              write(64,301)'PDF',kpdf,' ->',
-     #                   sum_wgt_resc_pdf(kpdf)
+            do kk=1,lhaPDFid(0)
+              do kpdf=1,2*numPDFpairs
+                write(64,301)'PDF',kk,kpdf,' ->',
+     #                     sum_wgt_resc_pdf(kpdf,kk)
+              enddo
             enddo
           endif
         endif
         if(jwgtinfo.eq.8.or.jwgtinfo.eq.9)then
           write(64,*)'  '
           write(64,*)'Max and min of rescaled weights'
-          do kr=1,numscales
-            do kf=1,numscales
-              if(event_norm.eq.'ave')then
-                write(64,400)'scales',kr,kf,' ->',
-     #                     xmax_wgt_resc_scale(kr,kf),
-     #                     xmin_wgt_resc_scale(kr,kf)
-              else
-                write(64,400)'scales',kr,kf,' ->',
-     #                     xmax_wgt_resc_scale(kr,kf),
-     #                     xmin_wgt_resc_scale(kr,kf)
-              endif
+          do kk=1,dyn_scale(0)
+            do kr=1,numscales
+              do kf=1,numscales
+                if(event_norm.eq.'ave'.or.event_norm.eq.'bia')then
+                  write(64,400)'scales',kk,kr,kf,' ->',
+     #                       xmax_wgt_resc_scale(kr,kf,kk)/maxevt,
+     #                       xmin_wgt_resc_scale(kr,kf,kk)/maxevt
+                else
+                  write(64,400)'scales',kk,kr,kf,' ->',
+     #                       xmax_wgt_resc_scale(kr,kf,kk),
+     #                       xmin_wgt_resc_scale(kr,kf,kk)
+                endif
+              enddo
             enddo
           enddo
-          if(event_norm.eq.'ave')then
-            do kpdf=1,2*numPDFpairs
-              write(64,401)'PDF',kpdf,' ->',
-     #                   xmax_wgt_resc_pdf(kpdf),
-     #                   xmin_wgt_resc_pdf(kpdf)
+          if(event_norm.eq.'ave'.or.event_norm.eq.'bia')then
+            do kk=1,lhaPDFid(0)
+              do kpdf=1,2*numPDFpairs
+                write(64,401)'PDF',kk,kpdf,' ->',
+     #                   xmax_wgt_resc_pdf(kpdf,kk)/maxevt,
+     #                   xmin_wgt_resc_pdf(kpdf,kk)/maxevt
+              enddo
             enddo
           else
+            do kk=1,lhaPDFid(0)
             do kpdf=1,2*numPDFpairs
-              write(64,401)'PDF',kpdf,' ->',
-     #                   xmax_wgt_resc_pdf(kpdf),
-     #                   xmin_wgt_resc_pdf(kpdf)
+              write(64,401)'PDF',kk,kpdf,' ->',
+     #                   xmax_wgt_resc_pdf(kpdf,kk),
+     #                   xmin_wgt_resc_pdf(kpdf,kk)
+              enddo
             enddo
           endif
         endif
@@ -655,10 +727,10 @@
 
       if(idec.eq.0)call finalizedecays()
 
- 300  format(a,2(1x,i2),a,(1x,f16.6))
- 301  format(a,6x,i3,a,(1x,f16.6))
- 400  format(a,2(1x,i2),a,2(1x,f16.6))
- 401  format(a,6x,i3,a,2(1x,f16.6))
+ 300  format(a,3(1x,i2),a,(1x,f16.6))
+ 301  format(a,6x,i2,i3,a,(1x,f16.6))
+ 400  format(a,3(1x,i2),a,2(1x,f16.6))
+ 401  format(a,6x,i2,i3,a,2(1x,f16.6))
 
       end
 

=== added symlink 'Template/NLO/Utilities/extra_weights.f'
=== target is u'../Source/extra_weights.f'
=== modified file 'Template/NLO/Utilities/makefile'
--- Template/NLO/Utilities/makefile	2014-10-31 10:50:31 +0000
+++ Template/NLO/Utilities/makefile	2017-08-07 21:43:30 +0000
@@ -2,8 +2,9 @@
 P01=$(shell echo $(P0) | cut -d'/' -f3)
 
 INCLUDE=-I../SubProcesses/$(P01)
-FF=gfortran -g -fno-automatic -ffixed-line-length-132 $(INCLUDE)
-FILES=handling_lhe_events.o fill_MC_mshell.o
+EXTRAW=-I../Source
+FF=gfortran -g -fno-automatic -ffixed-line-length-132 $(INCLUDE) $(EXTRAW) 
+FILES=handling_lhe_events.o fill_MC_mshell.o extra_weights.o
 FILES2=madfks_plot.o dbook.o setcuts.o analysis_td_template.o \
        open_output_files.o open_output_files_dummy.o appl_interface_dummy.o
 

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2017-08-04 23:05:05 +0000
+++ UpdateNotes.txt	2017-08-07 21:43:30 +0000
@@ -16,7 +16,11 @@
       RF: Fixed bug #1706072 related to wrong path with NLO gridpack mode
           diagrams with 1->3 decays.
       RF: Refactor of the (NLO) code related to extra weights (for PDF/scale uncertainties).
+<<<<<<< TREE
       RF: Fixed a bug (found by SF) that gave a seriously biased resonance mass when using MadSpin for decaying a 2->1 process.
+=======
+      RF+OM: Added the possibility to also have a bias-function for event generation at (f)NLO(+PS)
+>>>>>>> MERGE-SOURCE
 
 2.5.5(26/05/17)
       OM: Fixing bug in the creation of the LO gridpack introduced in 2.4.3. Since 2.4.3 the generated 

=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
--- madgraph/interface/amcatnlo_run_interface.py	2017-08-05 18:44:21 +0000
+++ madgraph/interface/amcatnlo_run_interface.py	2017-08-07 21:43:30 +0000
@@ -1755,7 +1755,10 @@
         """write the nevts files in the SubProcesses/P*/G*/ directories"""
         for job in jobs:
             with open(pjoin(job['dirname'],'nevts'),'w') as f:
-                f.write('%i\n' % job['nevents'])
+                if self.run_card['event_norm'].lower()=='bias':
+                    f.write('%i %f\n' % (job['nevents'],self.cross_sect_dict['xseca']))
+                else:
+                    f.write('%i\n' % job['nevents'])
 
     def combine_split_order_run(self,jobs_to_run):
         """Combines jobs and grids from split jobs that have been run"""
@@ -2579,6 +2582,8 @@
             self.cross_sect_dict['unit']='pb'
             self.cross_sect_dict['xsec_string']='Total cross section'
             self.cross_sect_dict['axsec_string']='Total abs(cross section)'
+        if self.run_card['event_norm'].lower()=='bias':
+            self.cross_sect_dict['xsec_string']+=', incl. bias (DO NOT USE)'
 
         if mode in ['aMC@NLO', 'aMC@LO', 'noshower', 'noshowerLO']:
             status = ['Determining the number of unweighted events per channel',
@@ -3165,6 +3170,8 @@
             p.communicate(input = '1\n')
         elif event_norm.lower() == 'unity':
             p.communicate(input = '3\n')
+        elif event_norm.lower() == 'bias':
+            p.communicate(input = '0\n')
         else:
             p.communicate(input = '2\n')
 
@@ -3576,7 +3583,7 @@
 
                     if self.banner.get('run_card', 'event_norm').lower() == 'sum':
                         norm = 1.
-                    elif self.banner.get('run_card', 'event_norm').lower() == 'average':
+                    else:
                         norm = 1./float(self.shower_card['nsplit_jobs'])
 
                     plotfiles2 = []
@@ -3981,6 +3988,8 @@
         #  number of events is not 0
         evt_files = [line.split()[0] for line in lines[:-1] if line.split()[1] != '0']
         evt_wghts = [float(line.split()[3]) for line in lines[:-1] if line.split()[1] != '0']
+        if self.run_card['event_norm'].lower()=='bias' and self.run_card['nevents'] != 0:
+            evt_wghts[:]=[wgt/float(self.run_card['nevents']) for wgt in evt_wghts]
         #prepare the job_dict
         job_dict = {}
         exe = 'reweight_xsec_events.local'

=== modified file 'madgraph/interface/common_run_interface.py'
--- madgraph/interface/common_run_interface.py	2017-08-04 23:05:05 +0000
+++ madgraph/interface/common_run_interface.py	2017-08-07 21:43:30 +0000
@@ -1733,7 +1733,7 @@
                     key = tuple(float(x) for x in split[:-1])
                     cross= float(split[-1])
                     if 'event_norm' in self.run_card and \
-                            self.run_card['event_norm'] in ['average', 'unity']:
+                            self.run_card['event_norm'] in ['average', 'unity', 'bias']:
                         cross *= (event_per_job+1 if i <nb_job_with_plus_one else event_per_job)
                     if len(all_cross) > pos:
                         all_cross[pos] += cross
@@ -1955,7 +1955,7 @@
                             raise Exception, "Some of the run failed: Please read %s_%s_debug.log" % (f, self.run_tag) 
                 
                 
-                if 'event_norm' in self.run_card and self.run_card['event_norm'] == 'average':
+                if 'event_norm' in self.run_card and self.run_card['event_norm'] in ['average','bias']:
                     for key, value in cross_sections.items():
                         cross_sections[key] = value / (nb_event+1)
                 lhe.remove()

=== modified file 'madgraph/interface/reweight_interface.py'
--- madgraph/interface/reweight_interface.py	2017-04-12 15:03:26 +0000
+++ madgraph/interface/reweight_interface.py	2017-08-07 21:43:30 +0000
@@ -717,7 +717,7 @@
                 
         # check normalisation of the events:
         if 'event_norm' in self.run_card:
-            if self.run_card['event_norm'] == 'average':
+            if self.run_card['event_norm'] in ['average','bias']:
                 for key, value in cross.items():
                     cross[key] = value / (event_nb+1)
                 

=== modified file 'madgraph/madevent/gen_crossxhtml.py'
--- madgraph/madevent/gen_crossxhtml.py	2016-09-05 10:10:45 +0000
+++ madgraph/madevent/gen_crossxhtml.py	2017-08-07 21:43:30 +0000
@@ -750,6 +750,7 @@
         # define at run_result
         self['run_name'] = run_name
         self['tag'] = run_card['run_tag']
+        self['event_norm'] = run_card['event_norm']
         self.event_path = pjoin(path,'Events')
         self.me_dir = path
         self.debug = None
@@ -1336,7 +1337,7 @@
         # Compute the text for eachsubpart
         
         sub_part_template_parton = """
-        <td rowspan=%(cross_span)s><center><a href="./HTML/%(run)s/results.html"> %(cross).4g <font face=symbol>&#177;</font> %(err).2g</a> %(syst)s </center></td>
+        <td rowspan=%(cross_span)s><center><a href="./HTML/%(run)s/results.html"> %(cross).4g <font face=symbol>&#177;</font> %(err).2g %(bias)s</a> %(syst)s </center></td>
         <td rowspan=%(cross_span)s><center> %(nb_event)s<center></td><td> %(type)s </td>
         <td> %(links)s</td>
         <td> %(action)s</td>
@@ -1411,6 +1412,10 @@
                 continue
             local_dico = {'type': ttype, 'run': self['run_name'], 'syst': '',
                           'tag': self['tag']}
+            if self['event_norm'].lower()=='bias':
+                local_dico['bias']='(biased, do not use)'
+            else:
+                local_dico['bias']=''
 
             if 'run_mode' in self.keys():
                 local_dico['run_mode'] = self['run_mode']
@@ -1545,7 +1550,8 @@
                            'links': 'banner only',
                            'action': action,
                            'run_mode': '',
-                           'syst':''
+                           'syst':'',
+                           'bias':''
                            }                                
                                   
         if self.debug is KeyboardInterrupt:

=== modified file 'madgraph/madevent/sum_html.py'
--- madgraph/madevent/sum_html.py	2017-07-07 10:47:20 +0000
+++ madgraph/madevent/sum_html.py	2017-08-07 21:43:30 +0000
@@ -337,7 +337,7 @@
         # this is for amcatnlo: the number of events has to be read from another file
         if self.nevents == 0 and self.nunwgt == 0 and isinstance(filepath, str) and \
                 os.path.exists(pjoin(os.path.split(filepath)[0], 'nevts')): 
-            nevts = int(open(pjoin(os.path.split(filepath)[0], 'nevts')).read())
+            nevts = int((open(pjoin(os.path.split(filepath)[0], 'nevts')).read()).split()[0])
             self.nevents = nevts
             self.nunwgt = nevts
         

=== modified file 'madgraph/various/lhe_parser.py'
--- madgraph/various/lhe_parser.py	2017-08-03 07:23:07 +0000
+++ madgraph/various/lhe_parser.py	2017-08-07 21:43:30 +0000
@@ -2236,8 +2236,9 @@
         return out
         
         
-    def parse(self, text):
-        """parse the line and create the related object"""
+    def parse(self, text, keep_bias=False):
+        """parse the line and create the related object.
+           keep bias allow to not systematically correct for the bias in the written information"""
         #0.546601845792D+00 0.000000000000D+00 0.000000000000D+00 0.119210435309D+02 0.000000000000D+00  5 -1 2 -11 12 21 0 0.24546101D-01 0.15706890D-02 0.12586055D+04 0.12586055D+04 0.12586055D+04  1  2  2  2  5  2  2 0.539995789976D+04
         #0.274922677249D+01 0.000000000000D+00 0.000000000000D+00 0.770516514633D+01 0.113763730192D+00  5 21 2 -11 12 1 2 0.52500539D-02 0.30205908D+00 0.45444066D+04 0.45444066D+04 0.45444066D+04 0.12520062D+01  1  2  1  3  5  1       -1 0.110944218997D+05
         # below comment are from Rik description email
@@ -2341,6 +2342,17 @@
         #     (and is passed to the integrator). It contains everything.
         #     from example: 0.110944218997D+05  
         self.ref_wgt = float(data[flag+14])
+        # 15. The bias weight. This weight is included in the self.ref_wgt, as well as in
+        #     the self.pwgt. However, it is already removed from the XWGTUP (and
+        #     scale/pdf weights). That means that in practice this weight is not used.
+        try:
+            self.bias_wgt = float(data[flag+15])
+        except KeyError:
+            self.bias_wgt = 1.0
+            
+        if not keep_bias:
+            self.ref_wgt /= self.bias_wgt
+            self.pwgt = [p/self.bias_wgt for p in self.pwgt]
 
         #check the momenta configuration linked to the event
         if self.type in self.real_type:

=== modified file 'madgraph/various/systematics.py'
--- madgraph/various/systematics.py	2017-08-03 07:23:07 +0000
+++ madgraph/various/systematics.py	2017-08-07 21:43:30 +0000
@@ -81,7 +81,8 @@
         self.force_write_banner = bool(write_banner)
         self.orig_dyn = self.banner.get('run_card', 'dynamical_scale_choice')
         self.orig_pdf = self.banner.run_card.get_lhapdf_id()
-    
+        matching_mode = self.banner.get('run_card', 'ickkw')
+
         #check for beam
         beam1, beam2 = self.banner.get_pdg_beam()
         if abs(beam1) != 2212 and abs(beam2) != 2212:
@@ -124,7 +125,13 @@
         if isinstance(dyn, str):
             dyn = dyn.split(',')
         self.dyn=[int(i) for i in dyn]
- 
+        # For FxFx only mode -1 makes sense
+        if matching_mode == 3:
+            self.dyn = [-1]
+        # avoid sqrts at NLO if ISR is possible
+        if 4 in self.dyn and self.b1 and self.b2 and not self.is_lo:
+            self.dyn.remove(4)
+
         if isinstance(together, str):
             self.together = together.split(',')
         else:
@@ -375,7 +382,7 @@
 
         if norm == 'sum':
             norm = 1
-        elif norm == 'average':
+        elif norm in ['average', 'bias']:
             norm = 1./nb_event
         elif norm == 'unity':
             norm = 1
@@ -838,6 +845,7 @@
                 tmp *= wgtpdf                
                 wgt += tmp
                 
+                
                 if __debug__ and dyn== -1 and Dmur==1 and Dmuf==1 and pdf==self.orig_pdf:
                     if not misc.equal(tmp, onewgt.ref_wgt, sig_fig=2):
                         misc.sprint(tmp, onewgt.ref_wgt, (tmp-onewgt.ref_wgt)/tmp)
@@ -846,7 +854,6 @@
                         misc.sprint(mur2,muf2)
                         raise Exception, 'not enough agreement between stored value and computed one'
                 
-                
         return wgt