Merge lp:~maddevelopers/mg5amcnlo/fix_color_resonances into lp:~madteam/mg5amcnlo/trunk

[Johan Alwall]

Ensure that intermediate resonances in the presence of triplet epsilon color operators downstream from the resonances will still get the correct color indices, based on the color of the initial state particles. This is further complicated by the possibility of multiparticle vertices, and necessitated a slightly modified treatment of the color for multiparticle vertices. All modifications are in addmothers.f.

[lp:~maddevelopers/mg5amcnlo/fix_color_resonances]

199. By Johan Alwall on 2012-01-09

Fixed stupid typos in addmothers.f

200. By Johan Alwall on 2012-01-11

Fixed octet t-channel, which was not fixed correctly before

201. By Johan Alwall on 2012-01-15

Almost finished addmothers, just need to remove the check that min and max are different

202. By Johan Alwall on 2012-01-16

I think I finally nailed it (at least no problems left with the process p p > t t~ j QED=2 QCD=3, t > b~ u~ mu+, t~ > b~ w- and vice versa)

203. By Johan Alwall on 2012-01-16

Updated VERSION

[Olivier Mattelaer (olivier-mattelaer)]

Excellent work :-)

review: Approve

[Valentin Hirschi (valentin-hirschi)]

Hi Johan,

I am not too aware of the bug this revision should fix but it sounds like a good thing ;)
However, I wanted to ask the following: have you tried to merge the trunk series2.0 into it to see if the unit and parallel test are still ok with your modification?

I have not really understood what level of the code your modification affects (You mentioned addmothers.f and I do not know any fortran code with this name, but there is definitely one function called like this in MG5, so I am not sure what you are referring to) and maybe they have nothing to do with any core code of MG5 and my question is irrelevant.

See you this friday at CERN I guess, cheers.

review: Needs Information

[Johan Alwall (johan-alwall)]

Hello Valentin,

> I am not too aware of the bug this revision should fix but it sounds like a
> good thing ;)

addmothers.f is in the Template/SubProcesses directory. Its job is to assign properties to intermediate on-shell resonances. It turns out that when you have baryon violating vertices with color epsilon^{ijk}, assigning the correct color indices to the intermediate resonances is a pure hell, especially if there are also multiparticle vertices in the model. Unfortunately, it is still necessary, in order for Pythia to correctly process the events. This revision (hopefully) ensures that intermediate resonances get the right color assignment in the presence of color epsilon^{ijk} for all relevant cases.

> However, I wanted to ask the following: have you tried to merge the trunk
> series2.0 into it to see if the unit and parallel test are still ok with your
> modification?

Unit and parallel tests are not testing MadEvent, so there can't be any problem. In any case, I thought at this point the 2.0 trunk is (supposed to be) identical to the present trunk, since none of our new developments have been merged yet.

> See you this friday at CERN I guess, cheers.

I'll be there on EVO.

All the best,
Johan

=== modified file 'Template/SubProcesses/addmothers.f'
--- Template/SubProcesses/addmothers.f	2011-11-14 05:50:15 +0000
+++ Template/SubProcesses/addmothers.f	2012-01-16 09:52:10 +0000
@@ -19,8 +19,8 @@
       integer isym(nexternal,99), jsym
       integer i,j,k,ida(2),ns,nres,ires,icl,ito2,idenpart,nc,ic
       integer mo_color,da_color(2),itmp
-      integer ito(-nexternal+3:nexternal),iseed,maxcolor
-      integer icolalt(2,-nexternal+3:2*nexternal-3)
+      integer ito(-nexternal+3:nexternal),iseed,maxcolor,maxorg
+      integer icolalt(2,-nexternal+2:2*nexternal-3)
       double precision qicl(-nexternal+3:2*nexternal-3), factpm
       double precision xtarget
       data iseed/0/
@@ -35,7 +35,7 @@
       double precision pmass(-nexternal:0,lmaxconfigs)
       double precision pwidth(-nexternal:0,lmaxconfigs)
       integer pow(-nexternal:0,lmaxconfigs)
-      logical first_time
+      logical first_time,tchannel
       save pmass,pwidth,pow
       data first_time /.true./
 
@@ -72,9 +72,9 @@
 c      common/to_colstats/ncols,ncolflow,ncolalt,icorg
 
       double precision pt
-      integer get_color,elim_indices
+      integer get_color,elim_indices,set_colmp,fix_tchannel_color
       real ran1
-      external pt,ran1,get_color,elim_indices
+      external pt,ran1,get_color,elim_indices,set_colmp,fix_tchannel_color
 
       if (first_time) then
          include 'props.inc'
@@ -142,6 +142,9 @@
 
       endif ! nc.gt.0
 
+c     Store original maxcolor to know if we have epsilon vertices
+        maxorg=maxcolor
+
 c     
 c     Get mother information from chosen graph
 c     
@@ -149,22 +152,70 @@
 c     First check number of resonant s-channel propagators
         ns=0
         nres=0
-
+        tchannel=.false.
 c     Loop over propagators to find mother-daughter information
-        do i=-1,-nexternal+3,-1
+        do i=-1,-nexternal+2,-1
 c       Daughters
-          ida(1)=iforest(1,i,iconfig)
-          ida(2)=iforest(2,i,iconfig)
-          do j=1,2
-            if(ida(j).gt.0) ida(j)=isym(ida(j),jsym)
-          enddo
+          if(i.gt.-nexternal+2)then
+             ida(1)=iforest(1,i,iconfig)
+             ida(2)=iforest(2,i,iconfig)
+             do j=1,2
+                if(ida(j).gt.0) ida(j)=isym(ida(j),jsym)
+             enddo
+          endif
 c       Decide s- or t-channel
-          if(iabs(sprop(numproc,i,iconfig)).gt.0) then ! s-channel propagator
+          if(i.gt.-nexternal+2.and.
+     $         iabs(sprop(numproc,i,iconfig)).gt.0) then ! s-channel propagator
             jpart(1,i)=sprop(numproc,i,iconfig)
             ns=ns+1
+          else if(nres.gt.0.and.maxcolor.gt.maxorg) then
+c         For t-channel propagators, just check that the colors are ok
+             if(i.eq.-nexternal+2) then
+c            This is the final t-channel, combining with leg 2
+                mo_color=0
+                if(.not.tchannel)then
+c                  There are no previous t-channels, so this is a combination of
+c                  2, 1 and the last s-channel
+                   ida(1)=1
+                   ida(2)=i+1
+                else
+c                  The daughter info is in iforest
+                   ida(1)=iforest(1,i,iconfig)
+                   ida(2)=iforest(2,i,iconfig)
+                endif
+c            Reverse colors of t-channels to get right color ordering
+                ncolmp=0
+                ncolmp=set_colmp(ncolmp,icolmp,2,jpart,
+     $               iforest(1,-max_branch,iconfig),icolalt,
+     $               icolalt(2,2),icolalt(1,2))
+             else
+                jpart(1,i)=tprid(i,iconfig)
+                mo_color=get_color(jpart(1,i))
+                ncolmp=0
+             endif
+             if(mo_color.gt.1.and.
+     $            mo_color.ne.3.and.mo_color.ne.8)then
+                da_color(1)=get_color(jpart(1,ida(1)))
+                da_color(2)=get_color(jpart(1,ida(2)))
+                call write_error(da_color(1), da_color(2), mo_color)
+             endif
+c            Set icolmp for daughters
+             ncolmp=set_colmp(ncolmp,icolmp,ida(2),jpart,
+     $            iforest(1,-max_branch,iconfig),icolalt,
+     $            icolalt(1,ida(2)),icolalt(2,ida(2)))
+c            Reverse colors of t-channels to get right color ordering
+             ncolmp=set_colmp(ncolmp,icolmp,ida(1),jpart,
+     $            iforest(1,-max_branch,iconfig),icolalt,
+     $            icolalt(2,ida(1)),icolalt(1,ida(1)))
+c            Fix t-channel color
+c             print *,'t-channel: ',i,ida(1),ida(2),mo_color
+c             print *,'colors: ',((icolmp(j,k),j=1,2),k=1,ncolmp)
+             maxcolor=fix_tchannel_color(mo_color,maxcolor,
+     $                                   ncolmp,icolmp,i,icolalt)
+             tchannel=.true.
+             cycle
           else
-c         Don't care about t-channel propagators
-            goto 100
+            goto 100             
           endif
 c       Set status codes for propagator
 c          if((igscl(0).ne.0.and.
@@ -191,70 +242,43 @@
 c          endif
 c       Set color info for all s-channels
           mo_color = get_color(jpart(1,i))
-          da_color(1) = get_color(jpart(1,ida(1)))
-          da_color(2) = get_color(jpart(1,ida(2)))
-          if(da_color(1).ne.2.and.da_color(2).lt.da_color(1).or.
-     $         da_color(2).eq.2)then
-c            Order daughters according to color, but always color 2 first
-             itmp=ida(1)
-             ida(1)=ida(2)
-             ida(2)=itmp
-             itmp=da_color(1)
-             da_color(1)=da_color(2)
-             da_color(2)=itmp
-          endif
-c     Reset list of color indices if not inside multipart. vertex
-c     (indicated by color 2)
-          if(da_color(1).ne.2)then
-             ncolmp=0
-          endif
+c     If inside multipart. vertex (indicated by color 2) cycle
+c         Set tentative mothers
+          jpart(2,i) = 1
+          jpart(3,i) = 2
+c         Set mother info for daughters
+          do j=1,2
+            jpart(2,ida(j)) = i
+            jpart(3,ida(j)) = i
+          enddo
+          if(mo_color.eq.2) cycle
+c     Reset list of color indices
+          ncolmp=0
 c     Add new color indices to list of color indices
-c     Note that color=2 means continued multiparticle index
           do j=1,2
-             if(da_color(j).eq.2.or.da_color(j).eq.1) cycle
-             ncolmp=ncolmp+1
-             icolmp(1,ncolmp)=icolalt(1,ida(j))
-             icolmp(2,ncolmp)=icolalt(2,ida(j))
-c            Avoid color sextet-type negative indices
-             if(icolmp(1,ncolmp).lt.0)then
-                ncolmp=ncolmp+1
-                icolmp(2,ncolmp)=-icolmp(1,ncolmp-1)
-                icolmp(1,ncolmp-1)=0
-                icolmp(1,ncolmp)=0
-             elseif(icolmp(2,ncolmp).lt.0)then
-                ncolmp=ncolmp+1
-                icolmp(1,ncolmp)=-icolmp(2,ncolmp-1)
-                icolmp(2,ncolmp-1)=0
-                icolmp(2,ncolmp)=0
-             endif
-             if(ncolmp.gt.maxcolmp)
-     $            call write_error(1000,ncolmp,maxcolmp)
+             ncolmp=set_colmp(ncolmp,icolmp,ida(j),jpart,
+     $            iforest(1,-max_branch,iconfig),icolalt,
+     $            icolalt(1,ida(j)),icolalt(2,ida(j)))
           enddo
-
+c          print *,'s-channel: ',i,mo_color,ida(1),ida(2)
+c          print *,'colors: ',((icolmp(j,k),j=1,2),k=1,ncolmp)
           if(mo_color.eq.1) then ! color singlet
-             icolalt(1,i) = 0
-             icolalt(2,i) = 0
+             maxcolor=elim_indices(0,0,ncolmp,icolmp,i,icolalt,maxcolor)
           elseif(mo_color.eq.-3) then ! color anti-triplet
-             maxcolor=elim_indices(0,1,ncolmp,icolmp,icolalt(1,i),maxcolor)
+             maxcolor=elim_indices(0,1,ncolmp,icolmp,i,icolalt,maxcolor)
           elseif(mo_color.eq.3) then ! color triplet
-             maxcolor=elim_indices(1,0,ncolmp,icolmp,icolalt(1,i),maxcolor)
+             maxcolor=elim_indices(1,0,ncolmp,icolmp,i,icolalt,maxcolor)
           elseif(mo_color.eq.-6) then ! color anti-sextet
-             maxcolor=elim_indices(0,2,ncolmp,icolmp,icolalt(1,i),maxcolor)
+             maxcolor=elim_indices(0,2,ncolmp,icolmp,i,icolalt,maxcolor)
           elseif(mo_color.eq.6) then ! color sextet
-             maxcolor=elim_indices(2,0,ncolmp,icolmp,icolalt(1,i),maxcolor)
+             maxcolor=elim_indices(2,0,ncolmp,icolmp,i,icolalt,maxcolor)
           elseif(mo_color.eq.8) then ! color octet
-             maxcolor=elim_indices(1,1,ncolmp,icolmp,icolalt(1,i),maxcolor)
-          elseif(mo_color.ne.2) then ! 2 indicates multipart. vertex
+             maxcolor=elim_indices(1,1,ncolmp,icolmp,i,icolalt,maxcolor)
+          else ! 2 indicates multipart. vertex
+             da_color(1) = get_color(jpart(1,ida(1)))
+             da_color(2) = get_color(jpart(1,ida(2)))
              call write_error(da_color(1), da_color(2), mo_color)
           endif
-c         Set tentative mothers
-          jpart(2,i) = 1
-          jpart(3,i) = 2
-c         Set mother info for daughters
-          do j=1,2
-            jpart(2,ida(j)) = i
-            jpart(3,ida(j)) = i
-          enddo
 c       Just zero helicity info for intermediate states
           jpart(7,i) = 0
         enddo                   ! do i
@@ -306,7 +330,9 @@
       return
       end
 
+c     *************************************
       subroutine write_error(ida1,ida2,imo)
+c     *************************************
       implicit none
       integer ida1,ida2,imo
 
@@ -332,8 +358,243 @@
  999  write(*,*) 'error'
       end
 
+c     *********************************************************************
+      function set_colmp(ncolmp,icolmp,npart,jpart,forest,icol,icol1,icol2)
+c     *********************************************************************
+      implicit none
+      integer maxcolmp
+      parameter(maxcolmp=20)
+      include 'nexternal.inc'
+      include 'genps.inc'
+c     Arguments
+      integer set_colmp
+      integer ncolmp,icolmp(2,*),npart,icol1,icol2
+      integer icol(2,-nexternal+2:2*nexternal-3)
+      integer jpart(7,-nexternal+3:2*nexternal-3)
+      integer forest(2,-max_branch:-1)
+c     Local
+      integer da_color(2),itmp,ida(2),icolor,ipart,i,j
+      integer get_color,set1colmp
+
+      set_colmp=ncolmp
+      icolor=get_color(jpart(1,npart))
+      if(icolor.eq.1) return
+      if(icolor.eq.2) then
+c     Multiparticle vertex - need to go through daughters and collect all colors
+         ipart=npart
+        do while(icolor.eq.2)
+          ida(1)=forest(1,ipart)
+          ida(2)=forest(2,ipart)
+          da_color(1)=get_color(jpart(1,ida(1)))
+          da_color(2)=get_color(jpart(1,ida(2)))
+c          print *,'iforest: ',ipart,ida(1),ida(2),da_color(1),da_color(2)
+          if(da_color(1).ne.2.and.da_color(2).lt.da_color(1).or.
+     $         da_color(2).eq.2)then
+c            Order daughters according to color, but always color 2 first
+             itmp=ida(1)
+             ida(1)=ida(2)
+             ida(2)=itmp
+             itmp=da_color(1)
+             da_color(1)=da_color(2)
+             da_color(2)=itmp
+          endif
+          do i=1,2
+             if(da_color(i).ne.1.and.da_color(i).ne.2)then
+                ncolmp=set1colmp(ncolmp,icolmp,icol(1,ida(i)),
+     $               icol(2,ida(i)))
+             endif
+          enddo
+          icolor=da_color(1)
+          ipart=ida(1)
+        enddo
+      else
+         ncolmp=set1colmp(ncolmp,icolmp,icol1,icol2)
+      endif
+      set_colmp=ncolmp
+      return
+      end
+
+c     ******************************************************
+      function set1colmp(ncolmp,icolmp,icol1,icol2)
+c     ******************************************************
+      implicit none
+c     Arguments
+      integer maxcolmp
+      parameter(maxcolmp=20)
+      integer set1colmp
+      integer ncolmp,icolmp(2,*),icol1,icol2
+
+      ncolmp=ncolmp+1
+      icolmp(1,ncolmp)=icol1
+      icolmp(2,ncolmp)=icol2
+c     Avoid color sextet-type negative indices
+      if(icolmp(1,ncolmp).lt.0)then
+         ncolmp=ncolmp+1
+         icolmp(2,ncolmp)=-icolmp(1,ncolmp-1)
+         icolmp(1,ncolmp-1)=0
+         icolmp(1,ncolmp)=0
+      elseif(icolmp(2,ncolmp).lt.0)then
+         ncolmp=ncolmp+1
+         icolmp(1,ncolmp)=-icolmp(2,ncolmp-1)
+         icolmp(2,ncolmp-1)=0
+         icolmp(2,ncolmp)=0
+      endif
+      if(ncolmp.gt.maxcolmp)
+     $     call write_error(1000,ncolmp,maxcolmp)
+      set1colmp=ncolmp
+      return
+      end
+
+c********************************************************************
+      function fix_tchannel_color(mo_color,maxcolor,ncolmp,icolmp,ires,
+     $                            icol)
+c********************************************************************
+c     Successively eliminate identical pairwise color indices from the
+c     icolmp list, until only (max) one triplet and one antitriplet remains
+c
+
+      implicit none
+      include 'nexternal.inc'
+      integer fix_tchannel_color
+      integer mo_color,maxcolor,ncolmp,icolmp(2,*)
+      integer ires,icol(2,-nexternal+2:2*nexternal-3)
+      integer i,j,i3,i3bar,max3,max3bar,min3,min3bar,maxcol,mincol
+
+c     Successively eliminate color indices in pairs until only the wanted
+c     indices remain
+      do i=1,ncolmp
+         do j=1,ncolmp
+            if(icolmp(1,i).ne.0.and.icolmp(1,i).eq.icolmp(2,j)) then
+               icolmp(1,i)=0
+               icolmp(2,j)=0
+            endif
+         enddo
+      enddo
+      
+      i3=0
+      i3bar=0
+      icol(1,ires)=0
+      icol(2,ires)=0
+      min3=1000
+      max3=0
+      min3bar=1000
+      max3bar=0
+      do i=1,ncolmp
+         if(icolmp(1,i).gt.0)then
+            i3=i3+1
+c           color for t-channels needs to be reversed
+            if(i3.eq.1) icol(2,ires)=icolmp(1,i)
+            if(i3.eq.2) icol(1,ires)=-icolmp(1,i)
+            if(icolmp(1,i).gt.max3) max3=icolmp(1,i)
+            if(icolmp(1,i).lt.min3) min3=icolmp(1,i)
+         endif
+         if(icolmp(2,i).gt.0)then
+            i3bar=i3bar+1
+c           color for t-channels needs to be reversed
+            if(i3bar.eq.1) icol(1,ires)=icolmp(2,i)
+            if(i3bar.eq.2) icol(2,ires)=-icolmp(2,i)
+            if(icolmp(2,i).gt.max3bar) max3bar=icolmp(2,i)
+            if(icolmp(2,i).lt.min3bar) min3bar=icolmp(2,i)
+         endif
+      enddo
+
+      if(mo_color.eq.0)then
+         icol(1,ires)=0
+         icol(2,ires)=0
+      endif
+
+      fix_tchannel_color=maxcolor
+      if(mo_color.le.1.and.i3.eq.0.and.i3bar.eq.0) return
+      if(mo_color.eq.3.and.(i3.eq.1.and.i3bar.eq.0
+     $     .or.i3bar.eq.1.and.i3.eq.0)) return
+      if(mo_color.eq.8.and.i3.eq.1.and.i3bar.eq.1) return
+
+c     Make sure that max and min don't come from the same octet
+      call clean_max_min(icolmp,ncolmp,max3,min3,max3bar,min3bar,i3,i3bar)
+
+      if(mo_color.le.1.and.i3-i3bar.eq.2.or.
+     $   mo_color.le.1.and.i3bar-i3.eq.2.or.
+     $   mo_color.le.1.and.i3.eq.1.and.i3bar.eq.1) then
+c     Replace the maximum index with the minimum one everywhere
+         maxcol=max(max3,max3bar)
+         if(maxcol.eq.max3) then
+            mincol=min3bar
+         else
+            mincol=min3
+         endif
+         do i=ires+1,-1
+            do j=1,2
+               if(icol(j,i).eq.maxcol)
+     $              icol(j,i)=mincol
+            enddo
+         enddo
+c         print *,'Replaced ',maxcol,' by ',mincol
+      elseif(mo_color.le.1.and.i3.eq.2.and.i3bar.eq.2) then
+c     Replace the maximum indices with the minimum ones everywhere
+         do i=ires+1,-1
+            do j=1,2
+               if(icol(j,i).eq.max3bar)
+     $              icol(j,i)=min3
+               if(icol(j,i).eq.max3)
+     $              icol(j,i)=min3bar
+            enddo
+         enddo
+c         print *,'Replaced ',max3bar,' by ',min3,' and ',max3,' by ',min3bar
+      elseif(mo_color.le.1.and.mod(i3,3).eq.0.and.mod(i3bar,3).eq.0)then
+c     This is epsilon index - do nothing
+         continue
+      else if(mo_color.eq.3.and.mod(i3-i3bar,3).eq.2) then
+c     This is an epsilon index
+         maxcolor=maxcolor+1
+         icol(1,ires)=maxcolor
+         icol(2,ires)=0
+c         print *,'Set mother color for ',ires,' to ',(icol(j,ires),j=1,2)
+      else if(mo_color.eq.3.and.mod(i3bar-i3,3).eq.2) then
+c     This is an epsilon index
+         maxcolor=maxcolor+1
+         icol(1,ires)=0
+         icol(2,ires)=maxcolor
+c         print *,'Set mother color for ',ires,' to ',(icol(j,ires),j=1,2)
+      else if(mo_color.eq.3.and.(i3-i3bar.eq.1.or.i3bar-i3.eq.1).or.
+     $        mo_color.eq.8.and.i3.eq.2.and.i3bar.eq.2) then
+c     Replace the maximum index with the minimum one everywhere
+c     (we don't know if we should replace i3 with i3bar or vice versa)
+         maxcol=max(max3,max3bar)
+         if(maxcol.eq.max3) then
+            mincol=min3bar
+         else
+            mincol=min3
+         endif
+         do i=ires+1,-1
+            do j=1,2
+               if(icol(j,i).eq.maxcol)
+     $              icol(j,i)=mincol
+            enddo
+         enddo
+c     Fix the color for ires (remember 3<->3bar for t-channels)
+         icol(1,ires)=0
+         icol(2,ires)=0
+c         print *,'Replaced ',maxcol,' by ',mincol
+         if(max3.eq.maxcol)then
+            if(i3-i3bar.ge.0) icol(2,ires)=min3
+            if(i3bar-i3.ge.0) icol(1,ires)=max3bar
+         else
+            if(i3-i3bar.ge.0) icol(2,ires)=max3
+            if(i3bar-i3.ge.0) icol(1,ires)=min3bar
+         endif
+c         print *,'Set mother color for ',ires,' to ',(icol(j,ires),j=1,2)
+      else
+c     Don't know how to deal with this
+         call write_error(i3,i3bar,mo_color)
+      endif
+
+      fix_tchannel_color=maxcolor
+      
+      return
+      end
+
 c*******************************************************************
-      function elim_indices(n3,n3bar,ncolmp,icolmp,icolres,maxcolor)
+      function elim_indices(n3,n3bar,ncolmp,icolmp,ires,icol,maxcolor)
 c*******************************************************************
 c     Successively eliminate identical pairwise color indices from the
 c     icolmp list, until only the wanted indices remain
@@ -346,8 +607,10 @@
 c
 
       implicit none
+      include 'nexternal.inc'
       integer elim_indices
-      integer n3,n3bar,ncolmp,icolmp(2,*),icolres(2),maxcolor
+      integer n3,n3bar,ncolmp,icolmp(2,*),maxcolor
+      integer ires,icol(2,-nexternal+2:2*nexternal-3)
       integer i,j,i3,i3bar
 
 c     Successively eliminate color indices in pairs until only the wanted
@@ -363,40 +626,47 @@
       
       i3=0
       i3bar=0
-      icolres(1)=0
-      icolres(2)=0
+      icol(1,ires)=0
+      icol(2,ires)=0
       do i=1,ncolmp
          if(icolmp(1,i).gt.0)then
             i3=i3+1
-            if(i3.eq.1) icolres(1)=icolmp(1,i)
-            if(i3.eq.2) icolres(2)=-icolmp(1,i)
+            if(i3.eq.1) icol(1,ires)=icolmp(1,i)
+            if(i3.eq.2) icol(2,ires)=-icolmp(1,i)
          endif
          if(icolmp(2,i).gt.0)then
             i3bar=i3bar+1
-            if(i3bar.eq.1) icolres(2)=icolmp(2,i)
-            if(i3bar.eq.2) icolres(1)=-icolmp(2,i)
+            if(i3bar.eq.1) icol(2,ires)=icolmp(2,i)
+            if(i3bar.eq.2) icol(1,ires)=-icolmp(2,i)
          endif
       enddo
 
+      if(n3.eq.0) icol(1,ires)=0
+      if(n3bar.eq.0) icol(2,ires)=0
+
       if(i3.ne.n3.or.i3bar.ne.n3bar) then
-         if(n3.gt.0.and.n3bar.eq.0)then
+         if(n3.gt.0.and.n3bar.eq.0.and.mod(i3bar+n3,3).eq.0)then
 c        This is an epsilon index interaction
             maxcolor=maxcolor+1
-            icolres(1)=maxcolor
-            icolres(2)=0
+            icol(1,ires)=maxcolor
             if(n3.eq.2)then
                maxcolor=maxcolor+1
-               icolres(2)=-maxcolor
+               icol(2,ires)=-maxcolor
             endif
-         elseif(n3bar.gt.0.and.n3.eq.0)then
+         elseif(n3bar.gt.0.and.n3.eq.0.and.mod(i3+n3bar,3).eq.0)then
 c        This is an epsilonbar index interaction
             maxcolor=maxcolor+1
-            icolres(1)=0
-            icolres(2)=maxcolor
+            icol(2,ires)=maxcolor
             if(n3.eq.2)then
                maxcolor=maxcolor+1
-               icolres(1)=-maxcolor
+               icol(1,ires)=-maxcolor
             endif
+         elseif(n3.gt.0.and.n3bar.eq.0.and.i3-i3bar.eq.n3.or.
+     $          n3bar.gt.0.and.n3.eq.0.and.i3bar-i3.eq.n3bar.or.
+     $          n3.eq.1.and.n3bar.eq.1.and.i3-i3bar.eq.0.or.
+     $          n3.eq.0.and.n3bar.eq.0.and.i3-i3bar.eq.0)then
+c        We have a previous epsilon which gives the wrong pop-up index
+            call fix_s_color_indices(n3,n3bar,i3,i3bar,ncolmp,icolmp,ires,icol)
          else
 c           Don't know how to deal with this
             call write_error(1001,n3,n3bar)
@@ -407,3 +677,153 @@
       
       return
       end
+
+c*******************************************************************
+      subroutine fix_s_color_indices(n3,n3bar,i3,i3bar,ncolmp,icolmp,
+     $                             ires,icol)
+c*******************************************************************
+c
+c     Fix color flow if some particle has got the wrong pop-up color
+c     due to epsilon-ijk vertices
+c
+
+      implicit none
+      include 'nexternal.inc'
+      integer n3,n3bar,ncolmp,icolmp(2,*),maxcolor
+      integer ires,icol(2,-nexternal+2:2*nexternal-3)
+      integer i,j,i3,i3bar
+      integer max_n3,max_n3bar,min_n3,min_n3bar,maxcol,mincol
+
+      max_n3=0
+      max_n3bar=0
+      min_n3=1000
+      min_n3bar=1000
+      do i=1,ncolmp
+         if(icolmp(1,i).gt.max_n3)
+     $        max_n3=icolmp(1,i)
+         if(icolmp(2,i).gt.max_n3bar)
+     $        max_n3bar=icolmp(2,i)
+         if(icolmp(1,i).gt.0.and.icolmp(1,i).lt.min_n3)
+     $        min_n3=icolmp(1,i)
+         if(icolmp(2,i).gt.0.and.icolmp(2,i).lt.min_n3bar)
+     $        min_n3bar=icolmp(2,i)
+      enddo
+
+      icol(1,ires)=0
+      icol(2,ires)=0
+
+c     Make sure that max and min don't come from the same octet
+      call clean_max_min(icolmp,ncolmp,max_n3,min_n3,
+     $                   max_n3bar,min_n3bar,i3,i3bar)
+
+      if(n3.eq.1.and.n3bar.eq.0.and.i3-i3bar.eq.n3.or.
+     $   n3bar.eq.1.and.n3.eq.0.and.i3bar-i3.eq.n3bar.or.
+     $   n3bar.eq.1.and.n3.eq.1.and.i3bar-i3.eq.0.or.
+     $   n3bar.eq.0.and.n3.eq.0.and.i3bar-i3.eq.0)then
+c     Replace the highest 3bar-index with the lowest 3-index,
+c     and vice versa
+         maxcol=max(max_n3,max_n3bar)
+         if(maxcol.eq.max_n3) then
+            mincol=min_n3bar
+         else
+            mincol=min_n3
+         endif
+         do i=ires,-1
+            do j=1,2
+               if(icol(j,i).eq.maxcol)
+     $              icol(j,i)=mincol
+            enddo
+         enddo
+c         print *,'Replaced ',maxcol,' with ',mincol
+         if(max_n3.eq.maxcol)then
+            if(n3.eq.1) icol(1,ires)=min_n3
+            if(n3bar.eq.1) icol(2,ires)=max_n3bar
+         else
+            if(n3.eq.1) icol(1,ires)=max_n3
+            if(n3bar.eq.1) icol(2,ires)=min_n3bar
+         endif
+c         print *,'Set mother color for ',ires,' to ',(icol(j,ires),j=1,2)
+      else
+c     Don't know how to deal with this
+         call write_error(1001,n3,n3bar)
+      endif
+      return
+      end
+
+c*******************************************************************************
+      subroutine clean_max_min(icolmp,ncolmp,max3,min3,max3bar,min3bar,i3,i3bar)
+c*******************************************************************************
+      implicit none
+      integer ncolmp,icolmp(2,*)
+      integer i,j,max3,max3bar,min3,min3bar,i3,i3bar
+      
+c      print *,'max/min3/3bar: ',max3,min3,max3bar,min3bar
+
+c     Make sure that max and min don't come from the same octet
+      do i=1,ncolmp
+         if(icolmp(1,i).eq.max3.and.icolmp(2,i).eq.min3bar)then
+            min3bar=1000
+            do j=1,ncolmp
+               if(j.eq.i) cycle
+               if(icolmp(2,j).lt.min3bar) min3bar=icolmp(2,j)
+            enddo
+         endif
+         if(icolmp(2,i).eq.max3bar.and.icolmp(1,i).eq.min3)then
+            min3=1000
+            do j=1,ncolmp
+               if(j.eq.i) cycle
+               if(icolmp(1,j).lt.min3) min3=icolmp(1,j)
+            enddo
+         endif
+      enddo
+
+c      print *,'max/min3/3bar: ',max3,min3,max3bar,min3bar
+
+c     ...and that max and min are different
+      if(i3.gt.1.and.max3.eq.min3.or.i3bar.gt.1.and.max3bar.eq.min3bar)then
+         if(max3.gt.max3bar)then
+c        Need to change min3, while still ensuring that max3bar is ok
+            max3bar=0
+            min3=1000
+            do i=1,ncolmp
+               if(icolmp(2,i).gt.max3bar.and.icolmp(2,i).ne.min3bar)
+     $              max3bar=icolmp(2,i)
+               if(icolmp(1,i).lt.min3.and.icolmp(1,i).ne.max3)
+     $              min3=icolmp(1,i)
+            enddo
+c     Make sure that max3bar and min3 don't come from the same octet
+            do i=1,ncolmp
+               if(icolmp(2,i).eq.max3bar.and.icolmp(1,i).eq.min3)then
+                  min3=1000
+                  do j=1,ncolmp
+                     if(j.eq.i.or.icolmp(1,j).eq.max3) cycle
+                     if(icolmp(1,j).lt.min3) min3=icolmp(1,j)
+                  enddo
+               endif
+            enddo
+         else
+c        Need to change min3bar, while still ensuring that max3 is ok
+            max3=0
+            min3bar=1000
+            do i=1,ncolmp
+               if(icolmp(1,i).gt.max3.and.icolmp(1,i).ne.min3)
+     $              max3=icolmp(1,i)
+               if(icolmp(2,i).lt.min3bar.and.icolmp(2,i).ne.max3bar)
+     $              min3bar=icolmp(2,i)
+            enddo
+c     Make sure that max3 and min3bar don't come from the same octet
+            do i=1,ncolmp
+               if(icolmp(1,i).eq.max3.and.icolmp(2,i).eq.min3bar)then
+                  min3bar=1000
+                  do j=1,ncolmp
+                     if(j.eq.i.or.icolmp(2,j).eq.max3bar) cycle
+                     if(icolmp(2,j).lt.min3bar) min3bar=icolmp(2,j)
+                  enddo
+               endif
+            enddo
+         endif
+      endif
+
+c      print *,'max/min3/3bar: ',max3,min3,max3bar,min3bar
+
+      end

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2011-12-25 12:39:55 +0000
+++ UpdateNotes.txt	2012-01-16 09:52:10 +0000
@@ -1,5 +1,10 @@
 Update notes for MadGraph 5 (in reverse time order)
 
+1.3.33 (01/01/12) JA: Revisited colors for propagators in addmothers.f
+                      to ensure that propagators get the correct color
+                      also from initial state particles (thanks to
+                      Michele Gabusi for forcing me to do this).
+
 1.3.32 (21/12/11) JA: Fixed a bug in the PDF reweighting routine,
                       which caused skewed eta distributions for
                       matched samples with pdfwgt=T. Thanks to Giulio

=== modified file 'madgraph/VERSION'
--- madgraph/VERSION	2011-12-21 01:15:14 +0000
+++ madgraph/VERSION	2012-01-16 09:52:10 +0000
@@ -1,2 +1,2 @@
-version = 1.3.32
-date = 2011-12-21
+version = 1.3.33_beta4
+date = 2012-01-16