Issue with decomposition of tHq process into quadratic and interference terms

Dear experts,

I am trying to simulate tHq production with CP-odd admixtures in the top-Yukawa using the NLO Higgs Characterisation model from http://feynrules.irmp.ucl.ac.be/wiki/HiggsCharacterisation in MadGraph version 2.9.18.

Ideally, I would like to have samples allowing for a continuous modelling of the process in the plane of CP-even and CP-odd top-Yukawa coupling modifiers (labelled `kHtt` and `kAtt` in the HC model) at NLO, by simulating the quadratic and interference terms individually and stitching samples together afterwards with according weights.

In the HC model, the HWW coupling is labeled with coupling order `QNP: 1`, while the Htt and Att couplings are not labeled.
For the SM tHq process, diagrams with top-Higgs coupling, W-Higgs coupling and their destructive interference contribute mostly.
As first test, I tried to retrieve the SM tHq result from the model when simulating the HWW and Htt quadratic terms and their destructive interference separately.
However, the results I received are inconsistent at LO and fail at NLO.

I use the following syntax to generate the SM tHq process (at LO here for tests):
```
import model HC_NLO_X0_UFO

generate p p > t x0 b~ j $$ w+ w-
add process p p > t~ x0 b j $$ w+ w-

```
Using the default run_card, I get a cross-section of 54fb.

I expected to receive the same result by simulating the quadratic Htt and HWW terms and their interference separately and adding up the cross-sections.

For the simulation of the quadratic Htt and HWW terms, I set other relevant couplings to zero.
So for the Htt quadratic term, I set `kSM` (scaling HWW) and `kHbb` (playing minor role here) to zero, the cross-section I get is 173fb.
To generate the HWW quadratic terms only, I set `kHtt` and `kHbb` to zero, resulting in a cross-section of 200fb.

Then, I use the following syntax for the Htt-HWW interference (no significant change observed when setting kHbb to zero or not):
```
import model HC_NLO_X0_UFO

generate p p > t x0 b~ j $$ w+ w- QNP^2==1
add process p p > t~ x0 b j $$ w+ w- QNP^2==1
```
The result is a cross-section of -216fb.
The sum of these three contributions is then approximately (200 + 173 - 216) fb = 157 fb, which is far away from the 54 fb I expected.
Do you know what might be wrong in this setup, causing the cross-sections not to add up?
As a cross-check, I tested how these three terms scale with the relevant coupling modifiers.
There are the quadratic and linear dependences on the modifiers as expected.

In order to also allow for explicit simulation of possibly relevant Htt-Att and Att-HWW interferences, we slightly changed the model such that the Htt and Att couplings are labelled with coupling orders `NPkth: 1` and `NPkta: 1`, respectively.
We made adaptions in the coupling_orders.py, couplings.py and CT_couplings.py, the adapted model can be found here: https://cernbox.cern.ch/s/ee7jEZ4gT4fez1w
With this modified model and the following syntax, I receive the same result for the quadratic Htt term, for example, as described above when setting kSM and kHbb to zero:

```
import model HC_NLO_X0_UFO_MOD

generate p p > t x0 b~ j $$ w+ w- NPkth==1
add process p p > t~ x0 b j $$ w+ w- NPkth==1
```
However, when requiring `NPkth^2==2`, which I expected to equivalently include the quadratic Htt-contributions only, the result changes (cross-section of 119 fb instead of 173 fb). This also appears when using the original model and doing the same with `QNP` for the HWW quadratic term.
Is this difference expected?
If so, what is the correct syntax for NLO simulation? I see that QNP==1 is not allowed at NLO.

Lastly, I would like to understand the following: going to NLO, I see errors in the pole checks for each of the interference terms.
Example syntax for Htt-HWW interference with the original model:
```
import model HC_NLO_X0_UFO

generate p p > t x0 b~ j $$ w+ w- QNP^2==1 WEIGHTED<=7 [QCD]
add process p p > t~ x0 b j $$ w+ w- QNP^2==1 WEIGHTED<=7 [QCD]
```
Here, I added `WEIGHTED<=7` because otherwise, I see this error:
Command "import /.automount/net_rw/net__scratch_cms3a/mausolf/NLO_samples/standaloneMG/CMSSW_13_2_4/src/MG5_aMC_v2_9_18/bin/cards/tHq_4f_NLO_MH125_kHtt_kHWW_HC.dat" interrupted in sub-command:
"generate p p > t x0 b~ j $$ w+ w- QNP^2==1 [QCD]" with error:
FKSProcessError : Cannot map born/real configurations between g u > x0 t b~ d [ all = QCD ] QNP^2==1 $$ w+ w- and g u > x0 t b~ d g WEIGHTED<=3 [ all = QCD ] QNP^2==1 $$ w+ w- (i,j=7,1): not same number of configurations: 8 0

Now, running the command from above, I get the pole cancellation error:
INFO: Result for check_poles:
Error detected in "launch auto "
write debug file /.automount/net_rw/net__scratch_cms3a/mausolf/NLO_samples/standaloneMG/CMSSW_13_2_4/src/MG5_aMC_v2_9_18/bin/tHq_4f_NLO_MH125_kHtt_kHWW_HC_NLO/run_01_tag_1_debug.log
If you need help with this issue please contact us on https://answers.launchpad.net/mg5amcnlo
str : Poles do not cancel, run cannot continue

Before, I see these warnings: WARNING: Some loop diagrams contributing to this process are discarded because they are not pure (QCD)-perturbation.
Make sure you did not want to include them.

I tried the suggestion from https://answers.launchpad.net/mg5amcnlo/+faq/2720 to set
#IRPoleCheckThreshold and #PrecisionVirtualAtRunTime to -1d0. With this, I get a positive cross-section (49 fb) for the destructive interference, so this solution does not appear to be valid here.
Do you have any suggestion how to fix this pole cancellation issue at NLO here which could still provide the correct physics?
Perhaps something else is wrong here since the results at LO already do not agree?
Or is this a conceptual issue that the interference terms cannot be computed separately?

I would greatly appreciate any guidance on the questions I raised here.

Many thanks in advance and best wishes.