Fits to validate the decay angle distributions in polarized same sign WWjj EWK events produced with Madgraph

Dear experts,

we are using MadGraph 5 version 2.7.3 to produce polarized same sign WWjj and inclusive WZ samples. We are running Madgraph with the following process:

For the same-sign WW fully leptonic final state:

TT samples (tried 2 reference frames):

generate p p > w+{T} w+{T} j j QCD=0, w+ > e+ ve, w+ > mu+ vm @0
add process p p > w+{T} w+{T} j j QCD=0, w+ > e+ ve @1
add process p p > w+{T} w+{T} j j QCD=0, w+ > mu+ vm @1
and set reference frame:
1, 2 = me_frame for the parton parton rest frame (PPRFr) and 3, 4, 5, 6 = me frame for WW rest frame (WWRFr)

LL sample (only in the WW rest frame):

generate p p > w+{0} w+{0} j j QCD=0, w+ > e+ ve, w+ > mu+ vm @0 add process p p > w+{0} w+{0} j j QCD=0, w+ > e+ ve @1
add process p p > w+{0} w+{0} j j QCD=0, w+ > mu+ vm @1
and set reference frame:
3,4,5,6 = me_frame parton parton rest frame (WWRFr)

Here the generation happened by generating a GridPack and then producing the lhef via ./run.sh since the generation fails when just using ./bin/generate_events. (Link to other launchpad https://answers.launchpad.net/mg5amcnlo/+question/693332)

For the WZ fully leptonic final state:

TT and LL samples both in WW rest frame - tested for MG version 2.7.2 and 2.7.3
[1] define l+ = e+ mu+ ta+
[2] define l- = e- mu- ta-
[3] define vl = ve vm vt
[4] define vl~ = ve~ vm~ vt~
[5] define wpm = w+ w-
[6] generate p p > wpm{0} z{0}, w+ > l+ vl, w- > l- vl~, z > l+ l- @0
[7] add process p p > wpm{0} z{0} j, w+ > l+ vl, w- > l- vl~, z > l+ l- @1
We used this process definition, this one without line [7] with cuts as described in the slides. We also used the following definition with just the default MadGraph run_card.dat.
define l+ = e+ mu+ ta+
define l- = e- mu- ta-
define vl = ve vm vt
define vl~ = ve~ vm~ vt~
generate p p > w+{T} z{T}, w+ > l+ vl, z > l+ l- @0
add process p p > w-{T} z{T}, w- > l- vl~, z > l+ l- @1
Here the generation happened using ./bin/generate_events
In both cases for the generated events we have plotted the cosine of the W decay angle (The W is defined as the pair lepton neutrino with the mass closer to the PDG value). The problem is that when we fit that distribution using the expected analytical formula from theory (formula 2.1 here https://link.springer.com/article/10.1007%2FJHEP07%282012%29124):

For transversal y(x) = f ∗3./(8.)∗(1+x)^2 +(1−f)∗3./(8.)∗(1−x)^2
For longitudinal y(x)=3∗(1−x^2)

the agreement does not look good, and the formula does not reproduce the distribution we see from the generated events. You can see the formula and the results of the fits in the attached slides for ssWWjj and WZ inclusive.

WZ inclusive: https://cloudstore.zih.tu-dresden.de/index.php/s/t4w4s7mJYkSzxpe
WWss: https://cloudstore.zih.tu-dresden.de/index.php/s/Q2ZZ3bDCyXmsEQp

Is there something wrong with the generate lines we have written above? We are wondering if we are doing something wrong on the definition of the reference frame so that somehow we are not able to reproduce it and perform the fits anymore.

Thanks
Best