Splitting a four fermion operator and restore the right coupling

Dear Madgraph staff,
at the end of this message you find all the model files I use (for better understanding my problems).

I am trying to work with a SM extension which includes four-fermions operators (identical particles). One such operator is of the form

$O1 = c \bar{q}_l \gamma^\mu T^A q_l \bar{q}_l \gamma_\mu T_A q_l$ where $q_l$ is the first family SU(2) doublet, $T^A$ is the SU(3) generator and $c$ is the coupling.

I know that it is not possible for madgraph to identify the correct fermionic flow when four identical particles are considered in the process. Then, I tried to split by hand this vertex.

In order to do so, I worked with FeynRules and I created an auxiliary colored field $G^c_\mu$ where $c$ is the color index. Then, I wrote the vertex $k \bar{q}_l \gamma^\mu T^A q_l G^A_\mu$ where $k$ is to be choose in such a way to match the coefficient $c$.

Then, I exported my model in a UFO output form and I set the value of auxiliary field mass as $m_{aux} = 1000 GeV$.

I am interested in contributions of O1 under processes p p > j j when c = -0.5. I have two questions:

1) Integrating out the heavy field, the matching of the couplings should give $k^2/(2*m^2_{aux}) =c/{m^2_{aux}$. Since I want c negative, I would define $k = i $ where $i$ is the imaginary part. In the file couplings.py of the exported UFO model I wrote

GC_1 = Coupling(name = 'GC_1',
                value = '-(CoeffO8q*complex(1,0))',
                order = {'NP':1})

where CoeffO8q = 1 and "NP" is a new interactionOrder I defined.

          Is this way of matching the desired coupling of the effective operator right?

2) I generated processes with my new model using

generate p p > j j
add process p p > auxiliarycolored >>jj

where "auxiliarycolored" is the name of the colored vector auxiliary field.

Is this syntax enough to simulate SM-interaction plus 4-fermions contact interaction?
Usually, I would like to substitute the auxiliary field propagator with a constant $(1/(2*m^2_{aux})$, but I don't know how Madgraph works with these kind of fields. I didn't define it as an auxiliary field in FeynRules, then I think Madgraph treats it as a physical field at all.

Doing as I told, I get that for small value of $k$, also for $k=1$, I don't see any significant deviations from the SM. Since I am trying to reproduce the results of this paper https://arxiv.org/pdf/1201.6510v2.pdf (look at fig. 1 pag. 7), I would expect to see big deviation for small value of the variable $\chi = exp(|y1-y2|)$ where $y_i$ are the pseudorapidities of the jets. And I don't see these elevations. So I suspect I am making some mistakes. Any help?

You can find the model files here https://www.dropbox.com/sh/hm77oyedcb3p055/AABzBlEPZvcPjAstjoclCR3ya?dl=0