setcuts.f does not work correct when combining subprocesses

[Rikkert Frederix]

Dear Experts,

I am studying WZ + 2j -> 3 lep + 2j + MET and want to investigate the ttbar + W background. I am trying to hardcode into setcuts.f and cuts.f a cut that will reject any event that does not contain
(a) a pair of oppositely charged leptons that (b) have the same flavor and (c) are in a window around mZ.
The problem is that I am finding events in my .lhe file that do not pass this criteria.

I attempt hard code my cuts into MG5 in a two-step fashion. In setcuts.f, I create boolean objects, similar to is_a_a(i) or is_a_j(i) and pass this to cuts.f. In the latter file, I ask for any same-flavor, opposite sign leptons.

In setcuts.f I have added
(Added near the beginning of the file)
      LOGICAL IS_A_EP(NEXTERNAL),IS_A_EM(NEXTERNAL)
      LOGICAL IS_A_MP(NEXTERNAL),IS_A_MM(NEXTERNAL)
      COMMON /TO_SPECISA/ (... default stuff....), IS_A_EP,IS_A_EM,IS_A_MP,IS_A_MM

(Added at the beginning of the first i=nincoming+1,nexternal loop)
         is_a_ep(i)=.false.
         is_a_em(i)=.false.
         is_a_mp(i)=.false.
         is_a_mm(i)=.false.

(Added a little further down in the same loop)
         if (idup(i,1,1).eq.(-11) ) is_a_ep(i)=.true. ! e+
         if (idup(i,1,1).eq.(-13) ) is_a_mp(i)=.true. ! mu+
         if (idup(i,1,1).eq.( 11) ) is_a_em(i)=.true. ! e-
         if (idup(i,1,1).eq.( 13) ) is_a_mm(i)=.true. ! mu-

At this point, I believe I have assigned my boolean objects the appropriate T/F values.
In cuts.f, I have added
(Added near the beginning of the file)
      LOGICAL IS_A_EP(NEXTERNAL),IS_A_EM(NEXTERNAL)
      LOGICAL IS_A_MP(NEXTERNAL),IS_A_MM(NEXTERNAL)
      COMMON /TO_SPECISA/ (... stuff that remains unchanged from default file... ) IS_A_EP,IS_A_EM,IS_A_MP,IS_A_MM
      integer zCounter
      double precision mZRng

(Added right before s-channel min & max invariant mass cuts are made)
      mZRng = 12.50D0
      zCounter = 0

      do i=nincoming+1,nexternal-1
         do j=i+1,nexternal

(Here I am asking if I have any same-flavor and oppositely charged leptons)
(Edit: The ampersands here are fortran line continuation flags)
            if( (is_a_ep(i).and.is_a_em(j)) .or.
     & (is_a_em(i).and.is_a_ep(j)) .or.
     & (is_a_mp(i).and.is_a_mm(j)) .or.
     & (is_a_mm(i).and.is_a_mp(j)) ) then

(Here I am asking if the invariant mass is less than 91.2 + 12.5 GeV and greater than 91.2 - 12.5 GeV)
               if( (dsqrt(Sumdot(p(0,i),p(0,j),+1d0)).LT.(91.2D0+mZRng))
     & .AND.
     & (dsqrt(Sumdot(p(0,i),p(0,j),+1d0)).GT.(91.2D0-mZRng))
     & ) then

( If I have a Z candidate, then I increment my Z Counter)
                  zCounter = zCounter + 1
               endif ! in mZ Range
            endif ! (e+,e-) or (e-,e+) or (mu+,mu-) or (mu-,mu+)
         enddo
      enddo

(If I have no Z candidate, then I reject the event)
      if(zCounter.eq.0) passcuts=.false.

To clarify, in my runcard.dat file, the following line has been set:
   T = cut_decays ! Cut decay products
The minimum and maximum invariant mass for pairs have been set to 0 and -1, respectively.

Now here is the kicker. For the process
define w = w+ w-
define l = l+ l-
define v = vl vl~
generate p p > t t~ w QCD=2 QED=1, t > b l v, t~ > b~ l v, w > l v

my resulting .lhe file has events like the one below where there is one one e+ and two mu-. Already this is a bad sign since there are no same-flavor, oppositely charged pairs. The invariant masses are m(e+mu-1) = 292 GeV and m(e+mu-2)= 109 GeV, which is also outside my mZ window (91.2+/-12.5 GeV). My data is mostly populated with my desired data, but I cannot seem to rid myself of these events.

<event>
15 0 0.2041000E-07 0.6854600E+03 0.7546771E-02 0.1028105E+00
       -2 -1 0 0 0 501 0.00000000000E+00 0.00000000000E+00 0.24501702010E+03 0.24501702010E+03 0.00000000000E+00 0. 1.
        1 -1 0 0 502 0 0.00000000000E+00 0.00000000000E+00 -0.14051334493E+04 0.14051334493E+04 0.00000000000E+00 0. -1.
      -24 2 1 2 0 0 0.92025752585E+02 -0.11709766887E+03 -0.90199089798E+03 0.91769962143E+03 0.80027569917E+02 0. 0.
        6 2 1 2 502 0 -0.16519660079E+03 0.97615337913E+02 0.13859680161E+03 0.29270001876E+03 0.17217885020E+03 0. 0.
       24 2 4 4 0 0 -0.78710286332E+02 0.21054620753E+02 0.10659082841E+01 0.11303849640E+03 0.78345129253E+02 0. 0.
       -6 2 1 2 0 501 0.73170848209E+02 0.19482330959E+02 -0.39672233282E+03 0.43975082919E+03 0.17395013126E+03 0. 0.
      -24 2 6 6 0 0 0.45405874961E+02 0.23388338823E+02 -0.39998715766E+03 0.41113440529E+03 0.80206390323E+02 0. 0.
        5 1 4 4 502 0 -0.86486314462E+02 0.76560717159E+02 0.13753089333E+03 0.17966152237E+03 0.46999998093E+01 0. -1.
      -11 1 5 5 0 0 -0.64694693043E+02 0.84043659919E+01 0.35127366336E+02 0.74094321924E+02 0.00000000000E+00 0. 1.
       12 1 5 5 0 0 -0.14015593289E+02 0.12650254761E+02 -0.34061458052E+02 0.38944174473E+02 0.00000000000E+00 0. -1.
       -5 1 6 6 0 501 0.27764973248E+02 -0.39060078637E+01 0.32648248436E+01 0.28616423892E+02 0.46999998093E+01 0. 1.
       13 1 7 7 0 0 -0.18411721052E+02 0.18288404331E+02 -0.63274953876E+02 0.68389889553E+02 0.00000000000E+00 0. -1.
      -14 1 7 7 0 0 0.63817596012E+02 0.50999344918E+01 -0.33671220379E+03 0.34274451574E+03 0.00000000000E+00 0. 1.
       13 1 3 3 0 0 0.77625165146E+01 -0.73737860294E+02 -0.37617999987E+03 0.38341742397E+03 0.00000000000E+00 0. -1.
      -14 1 3 3 0 0 0.84263236070E+02 -0.43359808577E+02 -0.52581089811E+03 0.53428219746E+03 0.00000000000E+00 0. 1.
</event>

I appreciate any and all advice.

Best Wishes,
Richard

[Rikkert Frederix (frederix)]

All the cuts set in the setcuts.f file, are only correct for the the final subprocess when subprocesses are combined. The setcuts.f routines is called for each of the subprocesses, but the cuts information is each time overwritten. Usually, this is not very harmful, because the cuts that are applied are identical for each of the subprocesses. But in special cases (I think it might even happen with the default set of allowed cuts) this is no longer true. I think this file need some updating.

Cheers,
Rikkert

[Olivier Mattelaer (olivier-mattelaer)]

Actually, that's the reason why we keep to maintain the v4 like output where the subprocesses are not merged.
You have to have
set group_subprocesses False
in your card.

Cheers,

Olivier

[Rikkert Frederix (frederix)]

It seems rather difficult to fix this. Indeed, the way to circumvent it is to set the
set group_subprocesses False
to prevent the grouping of the subprocesses. What will become available in the next release is an automatic stop of the code if such a problem is encountered.