Large width Higgs production crashed with llqq final states

Hi Lailin,
I first paste here the answer you got from Stefano Frixione via email, so that it can be useful for other users:

By requiring a final state that contains strongly-interacting
particles, you generate a process whose diagrams are much more
involved than those relevant to a purely leptonic process.
In other words, a final state that contains jets is much
less liable to be correctly described by selecting diagrams
with an intermediate Higgs. Furthermore, from the technical
point of view, the fks_singular.f module that I've set up
is not really suited to describe that process, the reason
being that it misses completely the contributions due to
QCD branchings from the quarks in the final state.

Bottom line: the above is not only a tehcnical issue, but
chiefly a physics one: I'd be very reluctant to compute a jet
process by a such a strong selection of diagrams.

Technically, there is (I think) a dirty way out, namely that
of considering QCD corrections only to "production" (ie as
far as the ZZ pair), while keeping the rest at the LO.
One could think of doing this in the following way: generate
a leptonic process, eg p p > h > e+ e- mu+ mu-; then,
at the level of the LHE event file (or in the interface
to the MC), change artificially the mu+ and mu- identifiers
with those of a quark pair; note that, when one does this,
one has to associate a new colour line to the newly created
quark pair; this is something that can (probably) be easily
done within the MC, at the level of the LHE interface.

Then, a cleaner solution would be to generate p p > h > z z without specifying the decay of the z bosons, and ask madspin to decay the two z bosons

Take into account that, since (W and Z) off-shell effects are important
to you, set the parameter BW_cut in madspin as appropriate. For what
concerns this, BW_cut in madspin_card.dat is not measured in GeV, but rather
in unit of the width; in other words, the virtualities you generate
are in the range:
pole_mass - BW_cut*width < virtuality < pole_mass + BW_cut*width

Let us know if this helps.
Cheers,

Marco