How to generated weighted events according to a function (biased generation)

For LO generation, we will refer to https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/LOEventGenerationBias
or see the new option below (for 3.5.2 and later)

For NLO, for version before 3.5.0 see https://indico.cern.ch/event/770470/contributions/3200859/attachments/1752107/2839310/MG5_aMC.pdf

Since 3.5.0 (for NLO) and 3.5.2 (for LO):
1: You can prepare a file somewhere in your filesystem containing the fortran definition of your function (see below for LO and NLO template)
2: You can pass the path to that file within the run_card.dat within the custom_fcts entry of the run_card (this entry expects a list of path)
3: If your function needs extra parameter to be controlled directly via the run_card by adding a line in the run_card like "10.0 = newvar". Your Fortran function could then use the parameter newvar" as long as you keep the "include run.inc"

############
# LO Template (introduced in v3.5.2)
############

      subroutine bias_wgt_custom(p, original_weight, bias_weight)
      implicit none
C
C Parameters
C
      include 'maxparticles.inc'
      include 'nexternal.inc'
      include 'run.inc' ! include defition from the run_card (via common-block).
c for this example to work, the run_card NEED to have two non default entry defined:
c - ptj_bias_target
c - ptj_bias_enhancement_power
c be sure to enter them as float!
c in script use "edit run_card 1000. = ptj_bias_target" to add the line

C
C Arguments
C
      double precision p(0:3,nexternal)
      double precision original_weight, bias_weight
C
C local variables
C
      integer i
      double precision ptj(nexternal)
      double precision max_ptj
C
C Global variables
C
C
C Mandatory common block to be defined in bias modules
C
      double precision stored_bias_weight ! default is 1.0
      logical impact_xsec, requires_full_event_info ! default is True and False respectively
      common/bias/stored_bias_weight,impact_xsec,
     & requires_full_event_info

      logical first
      data first /.true./ ! used to change the common-block value
C
C Accessing the details of the event
C
      logical is_a_j(nexternal),is_a_l(nexternal),
     & is_a_b(nexternal),is_a_a(nexternal),
     & is_a_onium(nexternal),is_a_nu(nexternal),
     & is_heavy(nexternal),do_cuts(nexternal)
      common/to_specisa/is_a_j,is_a_a,is_a_l,is_a_b,is_a_nu,
     & is_heavy,is_a_onium,do_cuts

C -------------------------------
C Changing setup of the bias mode
C -------------------------------
      if (first)then
         impact_xsec = .false.
         first = .false.
      endif
C --------------------
C BEGIN IMPLEMENTATION
C --------------------

      do i=1,nexternal
         ptj(i)=-1.0d0
         if (is_a_j(i)) then
            ptj(i)=sqrt(p(1,i)**2+p(2,i)**2)
         endif
      enddo

      max_ptj=-1.0d0
      do i=1,nexternal
         max_ptj = max(max_ptj,ptj(i))
      enddo
      if (max_ptj.lt.0.0d0) then
         bias_weight = 1.0d0
         return
      endif

      bias_weight= (max_ptj/ptj_bias_target)**ptj_bias_enhancement_power
      return

      end subroutine bias_wgt_custom

############
# NLO Template (introduced in v3.5.0)
############

      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 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=3,nexternal
c H_T=H_T+sqrt(max(0d0,(p(0,i)+p(3,i))*(p(0,i)-p(3,i))))
c enddo
c bias_wgt=H_T**2
      return
      end