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

[Johan Alwall]

Fix the color setting in addmasses.f for multiparticle vertices, which I forgot when I switched to let the model provide color info in v. 1.3.18.

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

Concerning the dummy PDG, it should check that the number is not define in the model
(and use another dummy PDG in this case).
This should be easy to do and much safer.

Otherwise, you use a lot of substraction for checking equality.
Which is not 100% safe: on the following case for example (in fact the other case are safe).
if(icolalt(1,ida(1))+icolalt(1,ida(2))-
10 + $ icolalt(2,ida(1))-icolalt(2,ida(2)).eq.0)

We might have trouble. In the other case, it will be more clear(and faster in fact) to use direct equality.

I'm confused about this comment:
1 3/6/8 -> 3/8

Otherwise I checked all color link and they are fine.

Cheers,

Olivier

review: Needs Fixing

[Johan Alwall (johan-alwall)]

Hello Olivier,

> Concerning the dummy PDG, it should check that the number is not define in the
> model
> (and use another dummy PDG in this case).
> This should be easy to do and much safer.

Ok, I'll do that.

> Otherwise, you use a lot of substraction for checking equality.
> Which is not 100% safe: on the following case for example (in fact the other
> case are safe).
> if(icolalt(1,ida(1))+icolalt(1,ida(2))-
> 10 + $ icolalt(2,ida(1))-icolalt(2,ida(2)).eq.0)

I want to have only one call to check for 1+1, 3+3bar, 6+6bar, 8+8 adding to a singlet. Perhaps the following is safer, and should cover all cases:
if(icolalt(1,ida(1)).eq.icolalt(2,ida(2)).and.icolalt(2,ida(1)).eq.icolalt(1,ida(2)).or.
icolalt(1,ida(1)).eq.-icolalt(1,ida(2)).and.icolalt(2,ida(1)).eq.-icolalt(2,ida(2))

> I'm confused about this comment:
> 1 3/6/8 -> 3/8

Ah - should be 1 3/6/8 -> 3/6/8

> Otherwise I checked all color link and they are fine.

Ok great!

Thanks,
Johan

[Johan Alwall (johan-alwall)]

Ok I have implemented your comments. Please approve and merge when you have a chance.

Thanks!
Johan

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

I see one modifications that was not there this afternoon.
In export_python, and which is not related to my comments.
So this require some additional tests and therefore will require some additional
time to test this properly.

So I'm not approve this for the moment.

Otherwise, I receive a bug report for one of the student.
One of the cut of the run_card is not correctly treated in the cuts.f
I'll push my fix in this version. So could you review that part.
I've send to him the fix by email, so he is suppose to check if this works and tell me tomorrow.

[lp:~maddevelopers/mg5amcnlo/fix_multiparticle_vertex_color]

186. By Olivier Mattelaer on 2011-10-29

add the cut missing in cuts.f

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

Could you explain to me why you have change export_python?
If they are a true reason, why this part is different between export_python and export_v4?

review: Needs Information

[Johan Alwall (johan-alwall)]

> Hi Johan,
>
> Could you explain to me why you have change export_python?
> If they are a true reason, why this part is different between export_python
> and export_v4?

Good point, I fixed it.

Johan

[lp:~maddevelopers/mg5amcnlo/fix_multiparticle_vertex_color]

187. By Johan Alwall on 2011-10-29

Changed my previous stupid way of checking for multiparticle vertices in export_python.py to the right way

[Olivier Mattelaer (olivier-mattelaer)]

Perfect, thanks, approved!!
Please check if you agree with my modifications with cuts.f
and make the merge.

review: Approve

[Johan Alwall (johan-alwall)]

Hello Olivier,

I honestly never understood the "lepton+neutrino" cut. What exactly is is supposed to do? It was Simon who implemented it for some particular purpose, and I never understood the point. Is it supposed to be the ET (in which case only the transverse momentum components are used) or the full invariant mass (which is not something that can be measured in the detector)?

Johan

review: Needs Information

[Olivier Mattelaer (olivier-mattelaer)]

This one was not used at all in cuts.f.
I put Simon in cc. In order that he can discuss this point with us.

Cheers,

Olivier

On 30-oct.-11, at 15:06, Johan Alwall wrote:

> Review: Needs Information
>
> Hello Olivier,
>
> I honestly never understood the "lepton+neutrino" cut. What exactly
> is is supposed to do? It was Simon who implemented it for some
> particular purpose, and I never understood the point. Is it supposed
> to be the ET (in which case only the transverse momentum components
> are used) or the full invariant mass (which is not something that
> can be measured in the detector)?
>
> Johan
> --
> https://code.launchpad.net/~maddevelopers/madgraph5/
> fix_multiparticle_vertex_color/+merge/80643
> You are reviewing the proposed merge of lp:~maddevelopers/madgraph5/
> fix_multiparticle_vertex_color into lp:madgraph5.

[lp:~maddevelopers/mg5amcnlo/fix_multiparticle_vertex_color]

188. By Olivier Mattelaer on 2011-10-30

correct install command for td script

189. By Johan Alwall on 2011-10-31

Corrected a comment about missing Et, updated UpdateNotes

=== modified file 'Template/SubProcesses/addmothers.f'
--- Template/SubProcesses/addmothers.f	2011-09-30 04:08:40 +0000
+++ Template/SubProcesses/addmothers.f	2011-10-31 13:57:25 +0000
@@ -205,6 +205,36 @@
           if(mo_color.eq.1) then ! color singlet
              icolalt(1,i) = 0
              icolalt(2,i) = 0
+          elseif(mo_color.eq.2) then ! used as dummy for multipart. vert
+             if(icolalt(1,ida(1)).eq.icolalt(2,ida(2)).and.
+     $            icolalt(2,ida(1)).eq.icolalt(1,ida(2)).or.
+     $            icolalt(1,ida(1)).eq.-icolalt(1,ida(2)).and.
+     $            icolalt(2,ida(1)).eq.-icolalt(2,ida(2))) then ! color singlet
+                icolalt(1,i) = 0
+                icolalt(2,i) = 0            
+             elseif(icolalt(1,ida(1)).eq.icolalt(2,ida(2))) then ! 3bar 3 -> 8 or 8 8 -> 8
+                icolalt(1,i) = icolalt(1,ida(2))
+                icolalt(2,i) = icolalt(2,ida(1))
+             else if(icolalt(1,ida(2)).eq.icolalt(2,ida(1))) then ! 3 3bar -> 8 or 8 8 -> 8
+                icolalt(1,i) = icolalt(1,ida(1))
+                icolalt(2,i) = icolalt(2,ida(2))
+             else if(icolalt(1,ida(1)).eq.0.and.icolalt(2,ida(1)).eq.0) then ! 1 3/6/8 -> 3/6/8
+                icolalt(1,i) = icolalt(1,ida(2))
+                icolalt(2,i) = icolalt(2,ida(2))
+             else if(icolalt(1,ida(2)).eq.0.and.icolalt(2,ida(2)).eq.0) then ! 3/6/8 1 -> 3/6/8
+                icolalt(1,i) = icolalt(1,ida(1))
+                icolalt(2,i) = icolalt(2,ida(1))
+             elseif(icolalt(1,ida(1)).eq.-icolalt(1,ida(2))) then ! antisextet
+                icolalt(1,i) = -icolalt(1,ida(2))
+                icolalt(2,i) = icolalt(2,ida(1))
+             elseif(icolalt(2,ida(1)).eq.-icolalt(2,ida(2))) then ! sextet
+                icolalt(1,i) = icolalt(1,ida(2))
+                icolalt(2,i) = -icolalt(1,ida(1))
+             else
+c               This might be due to multiparticle vertex with color combinations
+c               that can't be captured with just two sets of color indices
+                call write_error(da_color(1), da_color(2), mo_color)
+             endif
           elseif(mo_color.eq.-3) then ! color anti-triplet
                 icolalt(1,i) = 0
              if(da_color(1).eq.-3.and.da_color(2).eq.1)then

=== modified file 'Template/SubProcesses/cuts.f'
--- Template/SubProcesses/cuts.f	2011-09-16 05:33:32 +0000
+++ Template/SubProcesses/cuts.f	2011-10-31 13:57:25 +0000
@@ -53,7 +53,7 @@
       LOGICAL DEBUG
       integer i,j,njets,nheavyjets,hardj1,hardj2
       REAL*8 XVAR,ptmax1,ptmax2,htj,tmp,inclht
-      real*8 ptemp(0:3)
+      real*8 ptemp(0:3), ptemp2(0:3)
       character*20 formstr
 C
 C     PARAMETERS
@@ -262,12 +262,13 @@
          endif
       enddo
 c
-c    missing ET min & max cut 
-c    nb: missing et simply defined as the sum over the neutrino's 4 momenta
+c    missing ET min & max cut + Invariant mass of leptons and neutrino 
+c    nb: missing Et defined as the vector sum over the neutrino's pt
 c
 c-- reset ptemp(0:3)
       do j=0,3
-         ptemp(j)=0
+         ptemp(j)=0 ! for the neutrino
+         ptemp2(j)=0 ! for the leptons
       enddo
 c-  sum over the momenta
       do i=nincoming+1,nexternal
@@ -276,11 +277,19 @@
             do j=0,3
                ptemp(j)=ptemp(j)+p(j,i)
             enddo
+         elseif(is_a_l(i)) then            
+         if(debug) write (*,*) i,' -> lepton '
+            do j=0,3
+               ptemp2(j)=ptemp2(j)+p(j,i)
+            enddo
          endif
+
       enddo
 c-  check the et
       if(debug.and.ptemp(0).eq.0d0) write (*,*) 'No et miss in event'
       if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Et miss =',pt(ptemp(0)),'   ',misset,':',missetmax
+      if(debug.and.ptemp2(0).eq.0d0) write (*,*) 'No leptons in event'
+      if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Energy of leptons =',pt(ptemp2(0))
       if(ptemp(0).gt.0d0) then
          notgood=(pt(ptemp(0)) .lt. misset).or.
      &        (pt(ptemp(0)) .gt. missetmax)
@@ -290,6 +299,13 @@
             return
          endif
       endif
+      if (mmnl.gt.0d0.or.mmnlmax.lt.1d5)then
+         if(SumDot(ptemp,ptemp2,1d0).lt.mmnl.or.SumDot(ptemp, ptemp2,1d0).gt.mmnlmax) then
+            if(debug) write (*,*) 'lepton invariant mass -> fails'
+            passcuts=.false.
+            return
+         endif
+      endif
 c
 c     pt cut on heavy particles
 c     gives min(pt) for (at least) one heavy particle

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2011-10-27 03:18:04 +0000
+++ UpdateNotes.txt	2011-10-31 13:57:25 +0000
@@ -1,5 +1,12 @@
 Update notes for MadGraph 5 (in reverse time order)
 
+1.3.26 (31/10/11) JA: Fix color setting in MadEvent event file for
+                      multiparticle vertices, which was not taken into
+                      account in the upgrade in v. 1.3.18
+                  OM: Fixed mmnl cut (inv. mass of all leptons and neutrinos)
+                      which was never active.
+                  OM: Fix td install in Linux were a chmod was missing
+
 1.3.25 (27/10/11) JA: Ensure that the correct intermediate resonance
                       is always written in the event file, even when we
                       have resonances with identical properties.

=== modified file 'madgraph/VERSION'
--- madgraph/VERSION	2011-10-27 03:18:04 +0000
+++ madgraph/VERSION	2011-10-31 13:57:25 +0000
@@ -1,2 +1,2 @@
-version = 1.3.25
-date = 2011-10-27
+version = 1.3.26
+date = 2011-10-31

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2011-10-04 04:47:51 +0000
+++ madgraph/core/base_objects.py	2011-10-31 13:57:25 +0000
@@ -959,6 +959,11 @@
                     part.set('antiname', default[pdg])
                     if antipart:
                         antipart.set('antiname', default[pdg])
+
+    def get_first_non_pdg(self):
+        """Return the first positive number that is not a valid PDG code"""
+        return [c for c in range(1, len(self.get('particles')) + 1) if \
+                c not in self.get('particle_dict').keys()][0]
                 
     def write_param_card(self):
         """Write out the param_card, and return as string."""

=== modified file 'madgraph/core/helas_objects.py'
--- madgraph/core/helas_objects.py	2011-10-26 13:10:34 +0000
+++ madgraph/core/helas_objects.py	2011-10-31 13:57:25 +0000
@@ -1776,7 +1776,7 @@
             'id': self.get('interaction_id'),
             'legs': legs})
 
-    def get_s_and_t_channels(self, ninitial):
+    def get_s_and_t_channels(self, ninitial, new_pdg):
         """Returns two lists of vertices corresponding to the s- and
         t-channels of this amplitude/diagram, ordered from the outermost
         s-channel and in/down towards the highest number initial state
@@ -1885,7 +1885,8 @@
                 popped_legs = \
                            base_objects.LegList([vertex.get('legs').pop(0) \
                                                     for i in [0,1]])
-                popped_legs.append(base_objects.Leg({'id': 21,
+                popped_legs.append(base_objects.Leg({\
+                    'id': new_pdg,
                     'number': min([l.get('number') for l in popped_legs]),
                     'state': True,
                     'onshell': None}))

=== modified file 'madgraph/interface/cmd_interface.py'
--- madgraph/interface/cmd_interface.py	2011-10-27 01:33:41 +0000
+++ madgraph/interface/cmd_interface.py	2011-10-31 13:57:25 +0000
@@ -2596,7 +2596,7 @@
                 logger.info('Downloading TD for Linux 32 bit')
                 target = 'http://cp3wks05.fynu.ucl.ac.be/twiki/pub/Software/TopDrawer/td'
                 subprocess.call(['wget', target], cwd=os.path.join(MG5DIR,'td'))      
-            
+                os.chmod(os.path.join(MG5DIR,'td','td'), 0775)
                 if sys.maxsize > 2**32:
                     logger.warning('''td program (needed by MadAnalysis) is not compile for 64 bit computer
                 Please follow instruction in http://cp3wks05.fynu.ucl.ac.be/twiki/bin/view/Software/TopDrawer.''')

=== modified file 'madgraph/iolibs/export_cpp.py'
--- madgraph/iolibs/export_cpp.py	2011-07-08 09:54:01 +0000
+++ madgraph/iolibs/export_cpp.py	2011-10-31 13:57:25 +0000
@@ -1184,19 +1184,22 @@
         """Return the PIDs for any resonances in 2->2 and 2->3 processes."""
 
         resonances = []
-
+        model = self.matrix_elements[0].get('processes')[0].get('model')
+        new_pdg = model.get_first_non_pdg()
         # Get a list of all resonant s-channel contributions
         diagrams = sum([me.get('diagrams') for me in self.matrix_elements], [])
         for diagram in diagrams:
             schannels, tchannels = diagram.get('amplitudes')[0].\
-                                   get_s_and_t_channels(self.ninitial)
+                                   get_s_and_t_channels(self.ninitial, new_pdg)
 
             for schannel in schannels:
                 sid = schannel.get('legs')[-1].get('id')
-                width = self.model.get_particle(sid).get('width')
-                if width.lower() != 'zero':
-                    resonances.append(sid)
-
+                try:
+                    width = self.model.get_particle(sid).get('width')
+                    if width.lower() != 'zero':
+                        resonances.append(sid)
+                except KeyError:
+                    pass
         resonance_set = set(resonances)
 
         singleres = 0
@@ -1210,7 +1213,7 @@
         # schannel is True if all diagrams are s-channel and there are
         # no QCD vertices
         schannel = not any([\
-            len(d.get('amplitudes')[0].get_s_and_t_channels(self.ninitial)[0])\
+            len(d.get('amplitudes')[0].get_s_and_t_channels(self.ninitial, new_pdg)[0])\
                  == 0 for d in diagrams]) and \
                    not any(['QCD' in d.calculate_orders() for d in diagrams])
 

=== modified file 'madgraph/iolibs/export_python.py'
--- madgraph/iolibs/export_python.py	2011-04-17 17:11:07 +0000
+++ madgraph/iolibs/export_python.py	2011-10-31 13:57:25 +0000
@@ -262,6 +262,9 @@
         """Return the amp2(i) = sum(amp for diag(i))^2 lines"""
 
         ret_lines = []
+        # Get minimum legs in a vertex
+        minvert = min([max(diag.get_vertex_leg_numbers()) for diag in \
+                       matrix_element.get('diagrams')])
         if config_map:
             # In this case, we need to sum up all amplitudes that have
             # identical topologies, as given by the config_map (which
@@ -293,16 +296,12 @@
                 ret_lines.append(line)
             ret_lines.sort()
         else:
+            wf_dict = {}
+            vx_list = []
+            optimization = 0
             for idiag, diag in enumerate(matrix_element.get('diagrams')):
-                # Ignore any diagrams with 4-particle vertices.  The
-                # easiest way to get this info is to use the
-                # get_s_and_t_channels function, which collects all
-                # vertices corresponding to this diagram.
-                schannels, tchannels = diag.get('amplitudes')[0].\
-                                             get_s_and_t_channels(2)
-                allchannels = schannels + tchannels
-                if not allchannels or \
-                       any([len(vert.get('legs')) > 3 for vert in allchannels]):
+                # Ignore any diagrams with 4-particle vertices.
+                if max(diag.get_vertex_leg_numbers()) > minvert:
                     continue
                 # Now write out the expression for AMP2, meaning the sum of
                 # squared amplitudes belonging to the same diagram

=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2011-10-12 08:55:02 +0000
+++ madgraph/iolibs/export_v4.py	2011-10-31 13:57:25 +0000
@@ -1456,6 +1456,13 @@
             get_color=%d
             return
             """ % (part_id, model.get_particle(part_id).get_color())
+        # Dummy particle for multiparticle vertices with pdg given by
+        # first code not in the model
+        lines += """else if(ipdg.eq.%d)then
+c           This is dummy particle used in multiparticle vertices
+            get_color=2
+            return
+            """ % model.get_first_non_pdg()
         lines += """else
         write(*,*)'Error: No color given for pdg ',ipdg
         get_color=0        
@@ -1521,16 +1528,18 @@
 
         configs = [(i+1, d) for i,d in enumerate(matrix_element.get('diagrams'))]
         mapconfigs = [c[0] for c in configs]
+        model = matrix_element.get('processes')[0].get('model')
         return mapconfigs, self.write_configs_file_from_diagrams(writer,
                                                             [[c[1]] for c in configs],
                                                             mapconfigs,
-                                                            nexternal, ninitial)
+                                                            nexternal, ninitial,
+                                                            model)
 
     #===========================================================================
     # write_configs_file_from_diagrams
     #===========================================================================
     def write_configs_file_from_diagrams(self, writer, configs, mapconfigs,
-                                         nexternal, ninitial):
+                                         nexternal, ninitial, model):
         """Write the actual configs.inc file.
         
         configs is the diagrams corresponding to configs (each
@@ -1555,6 +1564,8 @@
 
         nconfigs = 0
 
+        new_pdg = model.get_first_non_pdg()
+
         for iconfig, helas_diags in enumerate(configs):
             if any([vert > minvert for vert in
                     [d for d in helas_diags if d][0].get_vertex_leg_numbers()]):
@@ -1571,7 +1582,7 @@
                     # get_s_and_t_channels gives vertices starting from
                     # final state external particles and working inwards
                     stchannels.append(h.get('amplitudes')[0].\
-                                      get_s_and_t_channels(ninitial))
+                                      get_s_and_t_channels(ninitial, new_pdg))
                 else:
                     stchannels.append((empty_verts, None))
 
@@ -1877,7 +1888,7 @@
         for iconf, configs in enumerate(s_and_t_channels):
             for vertex in configs[0] + configs[1][:-1]:
                 leg = vertex.get('legs')[-1]
-                if leg.get('id') == 21 and 21 not in particle_dict:
+                if leg.get('id') not in particle_dict:
                     # Fake propagator used in multiparticle vertices
                     mass = 'zero'
                     width = 'zero'
@@ -2414,6 +2425,7 @@
         configuration."""
 
         matrix_elements = subproc_group.get('matrix_elements')
+        model = matrix_elements[0].get('processes')[0].get('model')
 
         diagrams = []
         config_numbers = []
@@ -2437,7 +2449,8 @@
         return len(diagrams), \
                self.write_configs_file_from_diagrams(writer, diagrams,
                                                 config_numbers,
-                                                nexternal, ninitial)
+                                                nexternal, ninitial,
+                                                     model)
 
     #===========================================================================
     # write_leshouche_file

=== modified file 'tests/unit_tests/core/test_helas_objects.py'
--- tests/unit_tests/core/test_helas_objects.py	2011-10-26 13:10:34 +0000
+++ tests/unit_tests/core/test_helas_objects.py	2011-10-31 13:57:25 +0000
@@ -1444,7 +1444,7 @@
             for idiag in range(len(diagrams)):
                 if idiag in photon_none[ngluons]:
                     vertices, tchannels = \
-                          diagrams[idiag].get('amplitudes')[0].get_s_and_t_channels(2)
+                       diagrams[idiag].get('amplitudes')[0].get_s_and_t_channels(2, 20)
                     for ivert in range(len(vertices)):
                         if ivert in photon_none[ngluons][idiag]:
                             self.assertEqual(False,
@@ -1485,7 +1485,7 @@
             for idiag in range(len(diagrams)):
                 if idiag in quark_none[ngluons]:
                     vertices, tchannels = \
-                          diagrams[idiag].get('amplitudes')[0].get_s_and_t_channels(2)
+                      diagrams[idiag].get('amplitudes')[0].get_s_and_t_channels(2, 20)
                     for ivert in range(len(vertices)):
                         if ivert in quark_none[ngluons][idiag]:
                             self.assertEqual(False,
@@ -2264,7 +2264,7 @@
         for idiag in range(len(diagrams)):
             if idiag in quark_none:
                 vertices, tchannels = \
-                      diagrams[idiag].get('amplitudes')[0].get_s_and_t_channels(2)
+                      diagrams[idiag].get('amplitudes')[0].get_s_and_t_channels(2, 20)
                 for ivert in range(len(vertices)):
                     if ivert in quark_none[idiag]:
                         # This is forbidden leg
@@ -2324,7 +2324,7 @@
         for idiag in range(len(diagrams)):
             if idiag in quark_none:
                 vertices, tchannels = \
-                      diagrams[idiag].get('amplitudes')[0].get_s_and_t_channels(2)
+                      diagrams[idiag].get('amplitudes')[0].get_s_and_t_channels(2, 20)
                 for ivert in range(len(vertices)):
                     if ivert in quark_none[idiag]:
                         # This is forbidden leg

=== modified file 'tests/unit_tests/iolibs/test_export_v4.py'
--- tests/unit_tests/iolibs/test_export_v4.py	2011-10-12 16:09:26 +0000
+++ tests/unit_tests/iolibs/test_export_v4.py	2011-10-31 13:57:25 +0000
@@ -3477,6 +3477,10 @@
       ELSE IF(IPDG.EQ.21)THEN
         GET_COLOR=8
         RETURN
+      ELSE IF(IPDG.EQ.1)THEN
+C       This is dummy particle used in multiparticle vertices
+        GET_COLOR=2
+        RETURN
       ELSE
         WRITE(*,*)'Error: No color given for pdg ',IPDG
         GET_COLOR=0
@@ -5236,7 +5240,7 @@
 """C     Diagram 1
       DATA MAPCONFIG(1)/1/
       DATA (IFOREST(I,-1,1),I=1,2)/4,3/
-      DATA (SPROP(I,-1,1),I=1,1)/21/
+      DATA (SPROP(I,-1,1),I=1,1)/1/
       DATA TPRID(-1,1)/0/
       DATA (IFOREST(I,-2,1),I=1,2)/1,5/
       DATA TPRID(-2,1)/2/
@@ -5248,7 +5252,7 @@
       DATA (SPROP(I,-1,2),I=1,1)/6/
       DATA TPRID(-1,2)/0/
       DATA (IFOREST(I,-2,2),I=1,2)/-1,3/
-      DATA (SPROP(I,-2,2),I=1,1)/21/
+      DATA (SPROP(I,-2,2),I=1,1)/1/
       DATA TPRID(-2,2)/0/
 C     Diagram 3
       DATA MAPCONFIG(3)/3/
@@ -5256,12 +5260,12 @@
       DATA (SPROP(I,-1,3),I=1,1)/6/
       DATA TPRID(-1,3)/0/
       DATA (IFOREST(I,-2,3),I=1,2)/4,-1/
-      DATA (SPROP(I,-2,3),I=1,1)/21/
+      DATA (SPROP(I,-2,3),I=1,1)/1/
       DATA TPRID(-2,3)/0/
 C     Diagram 4
       DATA MAPCONFIG(4)/4/
       DATA (IFOREST(I,-1,4),I=1,2)/4,3/
-      DATA (SPROP(I,-1,4),I=1,1)/21/
+      DATA (SPROP(I,-1,4),I=1,1)/1/
       DATA TPRID(-1,4)/0/
       DATA (IFOREST(I,-2,4),I=1,2)/1,-1/
       DATA TPRID(-2,4)/2/
@@ -5279,7 +5283,7 @@
         self.assertFileContains('test',
 """      PMASS(-1,1)  = ZERO
       PWIDTH(-1,1) = ZERO
-      POW(-1,1) = 2
+      POW(-1,1) = 0
       PMASS(-2,1)  = ZERO
       PWIDTH(-2,1) = ZERO
       POW(-2,1) = 1
@@ -5288,16 +5292,16 @@
       POW(-1,2) = 1
       PMASS(-2,2)  = ZERO
       PWIDTH(-2,2) = ZERO
-      POW(-2,2) = 2
+      POW(-2,2) = 0
       PMASS(-1,3)  = ABS(MT)
       PWIDTH(-1,3) = ABS(WT)
       POW(-1,3) = 1
       PMASS(-2,3)  = ZERO
       PWIDTH(-2,3) = ZERO
-      POW(-2,3) = 2
+      POW(-2,3) = 0
       PMASS(-1,4)  = ZERO
       PWIDTH(-1,4) = ZERO
-      POW(-1,4) = 2
+      POW(-1,4) = 0
       PMASS(-2,4)  = ZERO
       PWIDTH(-2,4) = ZERO
       POW(-2,4) = 1
@@ -8534,7 +8538,7 @@
         for (idiag, diagram) in enumerate(diagrams):
 
             schannels, tchannels = diagram.get('amplitudes')[0].\
-                                         get_s_and_t_channels(2)
+                                         get_s_and_t_channels(2, 20)
 
             self.assertEqual([[l.get('number') for l in v.get('legs')] for v \
                               in schannels],
@@ -8736,6 +8740,10 @@
       ELSE IF(IPDG.EQ.24)THEN
         GET_COLOR=1
         RETURN
+      ELSE IF(IPDG.EQ.1)THEN
+C       This is dummy particle used in multiparticle vertices
+        GET_COLOR=2
+        RETURN
       ELSE
         WRITE(*,*)'Error: No color given for pdg ',IPDG
         GET_COLOR=0