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

[Johan Alwall]

JA: Use model information to determine color of particles
                      for reweighting and propagator color info.
                  JA: Enable proper 4-flavor matching (such as gg>hbb~+jets)
                      Note that you need the Pythia/PGS package v. 2.1.9 or
                      later to use with 4F matching.
                      Changes include: alpha_s reweighting also for b vertices,
                      new scale treatment (mu_F for pp>hbb~ is (pT_b^max m_Th)),
                      no clustering of gluons to final-state massive particles
                      in MadEvent.
                  JA: Fixed pdf reweighting for matching, which due to a mistake
                      had never been activated
                  JA: Improved phase space integration presetting further by
                      taking into account special cuts like xpt, ht etc.
    JA: Ensure that factorization scale settings and matching works
                      also in singlet t-channel exchange processes like
                      single top and VBF.
                  JA: Introduce new convention for invariant mass cuts
                      - if max < min, exclude intermediate range
                      (allows to exclude W/Z dijet resonances in VBF processes)

[lp:~maddevelopers/mg5amcnlo/fix_4F_matching]

180. By Johan Alwall on 2011-09-23

Fixed abs(particle id) for t-channel particles in get_color

181. By Johan Alwall on 2011-09-24

Offset Pythia PARP(67)=4 for singlet t-channel exchange processes

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

This is very good as usual, this is few points open for discussion.
Even if I would strongly suggest to include some of the remarqs, I already put this as ready for merge. Since I will probably not be online for a week don't wait for me anymore for a second approval.

concerning the get_color, I have few questions:
1) Why define it in each Subprocess directory and not model basis (restricted to those used off course)
2) If this is for speed issue (avoid too long list) then I would avoid to order them. (the one comming from the leg are more probable than those of the propagators and even more than flipped sign).
3) concerning your last modifications, i.e.:
abs_ids = set([abs(i) for i in wf_ids])
1) you should add a comment in the code explaining why you are doing that
2) add a test for checking this specific problem

concerning:
JA: Introduce new convention for invariant mass cuts
                      - if max < min, exclude intermediate range
                      (allows to exclude W/Z dijet resonances in VBF processes)
Did you include a comment on this on the wiki?

concerning:
      write(*,*)'Error:Color combination',da_color(1),da_color(2),
     $ '->',mo_color,' not implemented in addmothers.f'
This might be written a HUGE number of times shouldn't we
1) Put a limit to this number of print (or why not a critical error?)
2) Put a warning at the generation level in case if we have interactions with un-supported color interactions?

concerning the removal of the following lines
c Remove propagator if "decay" to same particle (i.e. radiation of g/gamma/higgs)
...
c Always remove if daughter final-state or already removed for this reason
...
c Else remove either this resonance or daughter, which is closer to mass shell

Why did you remove those lines. I also remark that the lines removing the propagator if he is too-offshell is commented. If this normal?

As I said before, I didn't check carefully the matching modification, but just read the modifications of the code.

Also all the parralel test are not working anymore. This should be fixed. This forbids all basic checks. If you want I can work on this for PyMadEvent. (With maybe a re-definition of the parralel tests suite.

review: Approve

[Johan Alwall (johan-alwall)]

Hi Olivier,

> Even if I would strongly suggest to include some of the remarqs, I already put
> this as ready for merge. Since I will probably not be online for a week don't
> wait for me anymore for a second approval.

Thanks as usual for your excellent points. I'll answer them below, then we can check if we want to include them later.

> concerning the get_color, I have few questions:
> 1) Why define it in each Subprocess directory and not model basis (restricted
> to those used off course)

I was thinking about putting it in the MODEL directory with all particles, but there are two reasons not to: 1) some models might have huge number of particles -> slower execution, 2) I don't want to change the v4 model files (in particular we'd need to modify the makefile if nothing else). We could put it in SubProcesses, but it seemed easier (and slightly more efficient) to just put it in each subprocess directory. I checked that the time to generate the file is truly negligible.

> 2) If this is for speed issue (avoid too long list) then I would avoid to
> order them. (the one comming from the leg are more probable than those of the
> propagators and even more than flipped sign).

Good point - however most processes have a very limited number of particles, so it shouldn't make a big difference. So I preferred to stick to better readability.

> 3) concerning your last modifications, i.e.:
> abs_ids = set([abs(i) for i in wf_ids])
> 1) you should add a comment in the code explaining why you are doing that
> 2) add a test for checking this specific problem

Ok will do.

> concerning:
> JA: Introduce new convention for invariant mass cuts
> - if max < min, exclude intermediate range
> (allows to exclude W/Z dijet resonances in VBF
> processes)
> Did you include a comment on this on the wiki?

Where on the wiki? Point me and I'll be happy to add it.

> concerning:
> write(*,*)'Error:Color combination',da_color(1),da_color(2),
> $ '->',mo_color,' not implemented in addmothers.f'
> This might be written a HUGE number of times shouldn't we
> 1) Put a limit to this number of print (or why not a critical error?)
> 2) Put a warning at the generation level in case if we have interactions with
> un-supported color interactions?

Ok I'll implement that it's written into the error file.

> concerning the removal of the following lines
> c Remove propagator if "decay" to same particle (i.e. radiation of
> g/gamma/higgs)
> ...
> c Always remove if daughter final-state or already removed for this reason
> ...
> c Else remove either this resonance or daughter, which is closer to mass shell
>
> Why did you remove those lines. I also remark that the lines removing the
> propagator if he is too-offshell is commented. If this normal?

This is now all done in myamp.f (where OnBW is defined). So it shouldn't be done again in addmothers.f.

> As I said before, I didn't check carefully the matching modification, but just
> read the modifications of the code.
>
> Also all the parralel test are not working anymore. This should be fixed. This
> forbids all basic checks. If you want I can work on this for PyMadEvent. (With
...

[lp:~maddevelopers/mg5amcnlo/fix_4F_matching]

182. By Johan Alwall on 2011-09-25

Added test for abs() in get_color.f

183. By Johan Alwall on 2011-09-26

Stop execution if color combination not implemented in addmothers.f

184. By Johan Alwall on 2011-09-28

Set ptclus for non-matching particles to sqrt(s) (i.e., particles from t-channel color singlet exchange or decays), so Pythia knows not to include them in matching

185. By Johan Alwall on 2011-09-29

Fix bug for 2->1 processes

186. By Johan Alwall on 2011-09-30

Fixed some minor issues

187. By Johan Alwall on 2011-10-02

Change definition for forbidden s-channels to only exclude onshell s-channels at runtime.

188. By Johan Alwall on 2011-10-03

Updated UpdateNotes.txt

189. By Johan Alwall on 2011-10-03

It turns out that not only abs(pdg) is needed by get_color but all antipdg codes

190. By Johan Alwall on 2011-10-03

Updated UpdateNotes. Changed Pythia PS scale to max(scales) instead of average (seems to give better matching both in single top and in Higgs VBF). Also removed factor 1/2 to offset PARP(67) for t-channel singlet exchange, since it doesn't give obviously better results but slightly complicates the code.

191. By Johan Alwall on 2011-10-03

Merged with the latest version of trunk (v 1.3.17)

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

# Use from_group = None to indicate that this s-channel is forbidden
newleg = copy.copy(vert.get('legs').pop(-1))
newleg.set('from_group', None)

I don't really like to re-use from_group. This is in some sense nice because, you limit the use of RAM, but I'm worry about readability of the code (especially in the future)
from_group name is not related at all with forbiden s-channel, so this is quite confusing.
In fact, it will be nicer (but may be harder) to include this in state. Indeed, this variable is nothing else that the fact that the particle is in T or S channel. So it might be natural to modify that variable. Another Solution is to create a new object.

This is off course a point open to discussion. What's your opinion?

Cheers,

Olivier

review: Needs Information

[Johan Alwall (johan-alwall)]

Hello Olivier,

> I don't really like to re-use from_group. This is in some sense nice because, you limit the use of RAM, but I'm worry about readability of the code (especially in the future)
> from_group name is not related at all with forbiden s-channel, so this is quite confusing.
> In fact, it will be nicer (but may be harder) to include this in state. Indeed, this variable is nothing else that the fact that the particle is in T or S channel. So it might be natural to modify that variable. Another Solution is to create a new object.

The most natural would be to have an "onshell" flag for leg (just as
we have for wavefunctions). I'm just worried about memory, since legs
tend to be the most proliferous objects in the code. But if you really
prefer, I'll do that instead.

Cheers,
Johan

[lp:~maddevelopers/mg5amcnlo/fix_4F_matching]

192. By Johan Alwall on 2011-10-04

Removed the use of from_group in Leg to denote forbidden onshell/decay legs, instead use new flag onshell

193. By Johan Alwall on 2011-10-04

Included mention in UpdateNotes.txt of the special process p p > t b~ j $ w+ w- t~ + p p > t b~ j j $ w+ w- t~, where MAXJETS=1 must be set in the pythia_card.dat

194. By Johan Alwall on 2011-10-04

F**k I intruduced stupidity when trying to fix the inevitable problem due to the stupid e+e-_e+e-.pkl file before I realized what the problem was. Fixed.

[Johan Alwall (johan-alwall)]

Hi Olivier,

I have removed the reuse of from_group and instead use a new flag, "onshell" (rev. 192+194). Please let me know what you think. I'm ready to merge this tomorrow.

Cheers,
Johan

[Olivier Mattelaer (olivier-mattelaer)]

Perfect,

Just one thing, but this might be done/test later (I'm not 100 percent sure that this is possible) but couldn't we remove completely from_group

I mean something like:
for leg in all_legs:
    del leg['from_group']

With a fast grep it looks like a possible thing to do.
But I propose that for a future release.

Thanksa lot,

Olivier

review: Approve

[Johan Alwall (johan-alwall)]

Hi Olivier,

> Just one thing, but this might be done/test later (I'm not 100 percent sure that this is possible) but couldn't we remove completely from_group?

It should give an equivalent memory improvement to simply set it to
None after diagram generation, right? That would be less confusing and
less risk for problems later on I think.

Cheers,
Johan

[lp:~maddevelopers/mg5amcnlo/fix_4F_matching]

195. By Johan Alwall on 2011-10-05

Stop execution if clustering fails (since this indicates a critical bug)

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

> Hi Olivier,
>
>> Just one thing, but this might be done/test later (I'm not 100
>> percent sure that this is possible) but couldn't we remove
>> completely from_group?
>
> It should give an equivalent memory improvement to simply set it to
> None after diagram generation, right? That would be less confusing and
> less risk for problems later on I think.

This is not clear for me at all, and depend too much of the internal
optimization of python. So we might try to test it in order to be sure.

Cheers,

Olivier
On 04-oct.-11, at 22:14, Johan Alwall wrote:

> Hi Olivier,
>
>> Just one thing, but this might be done/test later (I'm not 100
>> percent sure that this is possible) but couldn't we remove
>> completely from_group?
>
> It should give an equivalent memory improvement to simply set it to
> None after diagram generation, right? That would be less confusing and
> less risk for problems later on I think.
>
> Cheers,
> Johan
>
> --
> https://code.launchpad.net/~maddevelopers/madgraph5/fix_4F_matching/
> +merge/75882
> You are reviewing the proposed merge of lp:~maddevelopers/madgraph5/
> fix_4F_matching into lp:madgraph5.

=== modified file 'Template/Source/setrun.f'
--- Template/Source/setrun.f	2011-06-01 21:40:33 +0000
+++ Template/Source/setrun.f	2011-10-05 03:29:23 +0000
@@ -58,6 +58,7 @@
       integer mothup(2,nexternal)
       integer icolup(2,nexternal,maxflow,maxsproc)
       include 'leshouche.inc'
+      data pdfwgt/.false./
 c
 c
 c
@@ -312,6 +313,7 @@
 
       if(ickkw.gt.0)then
          call get_real   (npara,param,value," alpsfact "       ,alpsfact , 1d0)
+         call get_logical   (npara,param,value," pdfwgt "       ,pdfwgt , .true.)
       endif
       if(ickkw.eq.2)then
         call get_integer(npara,param,value," highestmult "    ,nhmult, 0)

=== modified file 'Template/SubProcesses/addmothers.f'
--- Template/SubProcesses/addmothers.f	2011-09-11 13:02:35 +0000
+++ Template/SubProcesses/addmothers.f	2011-10-05 03:29:23 +0000
@@ -6,10 +6,10 @@
       include 'maxconfigs.inc'
       include 'nexternal.inc'
       include 'coupl.inc'
+      include 'maxamps.inc'
       include 'cluster.inc'
       include 'message.inc'
       include 'run.inc'
-      include 'maxamps.inc'
 
       integer jpart(7,-nexternal+3:2*nexternal-3),npart,ip,numproc
       double precision pb(0:4,-nexternal+3:2*nexternal-3)
@@ -18,6 +18,7 @@
 
       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)
       double precision qicl(-nexternal+3:2*nexternal-3), factpm
@@ -66,8 +67,9 @@
 c      common/to_colstats/ncols,ncolflow,ncolalt,icorg
 
       double precision pt
+      integer get_color
       real ran1
-      external pt,ran1
+      external pt,ran1,get_color
 
       if (first_time) then
          include 'props.inc'
@@ -121,8 +123,8 @@
          if(ic.gt.0) then
             icolalt(1,isym(i,jsym))=icolup(1,i,ic,numproc)
             icolalt(2,isym(i,jsym))=icolup(2,i,ic,numproc)
-            if (icolup(1,i,ic, numproc).gt.maxcolor) maxcolor=icolup(1,i,ic, numproc)
-            if (icolup(2,i,ic, numproc).gt.maxcolor) maxcolor=icolup(2,i,ic, numproc)
+            if (abs(icolup(1,i,ic, numproc)).gt.maxcolor) maxcolor=icolup(1,i,ic, numproc)
+            if (abs(icolup(2,i,ic, numproc)).gt.maxcolor) maxcolor=icolup(2,i,ic, numproc)
          endif
       enddo
 
@@ -182,60 +184,133 @@
 c            jpart(6,i)=3
 c            nres=nres-1
 c          endif
-c     Remove propagator if "decay" to same particle (i.e. radiation of g/gamma/higgs)
-          if(jpart(6,i).eq.2 .and. (jpart(1,i).eq.jpart(1,ida(1)).or.
-     $         jpart(1,i).eq.jpart(1,ida(2)))) then
-             if(jpart(1,i).eq.jpart(1,ida(1))) idenpart=ida(1)
-             if(jpart(1,i).eq.jpart(1,ida(2))) idenpart=ida(2)
-c     Always remove if daughter final-state or already removed for this reason
-             if(idenpart.gt.0.or.jpart(6,idenpart).eq.4) then
-                jpart(6,i)=4
-                nres=nres-1
-c     Else remove either this resonance or daughter, which is closer to mass shell
-             elseif(abs(pb(4,i)-pmass(i,iconfig)).gt.
-     $               abs(pb(4,idenpart)-pmass(i,iconfig))) then
-                jpart(6,i)=4
-                nres=nres-1
-             else if(jpart(6,idenpart).eq.2) then
-                jpart(6,idenpart)=4
-                nres=nres-1
-             endif
-          endif
 c       Set color info for all s-channels
-c       Fist set "safe" color info
-          if(icolalt(1,ida(1))+icolalt(1,ida(2))-
-     $       icolalt(2,ida(1))-icolalt(2,ida(2)).eq.0) then ! color singlet
-            icolalt(1,i) = 0
-            icolalt(2,i) = 0            
-          elseif(icolalt(1,ida(1))-icolalt(2,ida(2)).eq.0) 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))-icolalt(2,ida(1)).eq.0) 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/8 -> 3/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/8 1 -> 3/8
-            icolalt(1,i) = icolalt(1,ida(1))
-            icolalt(2,i) = icolalt(2,ida(1))
-          else if(icolalt(1,ida(2)).gt.0.and.icolalt(1,ida(1)).gt.0.and.
-     $           icolalt(2,ida(2)).le.0.and.icolalt(2,ida(1)).le.0) then         ! sextet
-            maxcolor=maxcolor+1
-            icolalt(1,i) = maxcolor
-            icolalt(2,i) = 0
-          else if(icolalt(2,ida(2)).gt.0.and.icolalt(2,ida(1)).gt.0.and.
-     $           icolalt(1,ida(2)).le.0.and.icolalt(1,ida(1)).le.0) then         ! antisextet
-            maxcolor=maxcolor+1
-            icolalt(1,i) = 0
-            icolalt(2,i) = maxcolor
-          else if(jpart(6,i).ge.3) then ! Don't need to match
-            icolalt(1,i) = icolalt(1,ida(1))+icolalt(1,ida(2))
-            icolalt(2,i) = icolalt(2,ida(1))+icolalt(2,ida(2))
+          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(2).lt.da_color(1))then
+c            Order daughters according to color
+             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          print *,'graph: ',iconfig
+c          print *,'Resonance: ',i,' daughters ',ida(1),ida(2),
+c     $         ' ids ',jpart(1,i),jpart(1,ida(1)),jpart(1,ida(2)),
+c     $         ' colors ',mo_color,da_color(1),da_color(2)
+          
+          if(mo_color.eq.1) then ! color singlet
+             icolalt(1,i) = 0
+             icolalt(2,i) = 0
+          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
+                icolalt(2,i) = icolalt(2,ida(1))
+             elseif(da_color(1).eq.-3.and.da_color(2).eq.8)then
+                icolalt(2,i) = icolalt(2,ida(2))
+             elseif(da_color(1).eq.-6.and.da_color(2).eq.3)then
+                if(-icolalt(1,ida(1)).eq.icolalt(1,ida(2)))then
+                   icolalt(2,i) = icolalt(2,ida(1))
+                else
+                   icolalt(2,i) = -icolalt(1,ida(1))
+                endif
+             elseif(da_color(1).eq.3.and.da_color(2).eq.3)then
+                maxcolor=maxcolor+1
+                icolalt(2,i) = maxcolor
+             else
+                call write_error(da_color(1), da_color(2), mo_color)
+             endif
+          elseif(mo_color.eq.3) then ! color triplet
+                icolalt(2,i) = 0
+             if(da_color(1).eq.1.and.da_color(2).eq.3)then
+                icolalt(1,i) = icolalt(1,ida(2))
+             elseif(da_color(1).eq.3.and.da_color(2).eq.8)then
+                icolalt(1,i) = icolalt(1,ida(2))
+             elseif(da_color(1).eq.-3.and.da_color(2).eq.6)then
+                if(icolalt(2,ida(1)).eq.icolalt(1,ida(2)))then
+                   icolalt(1,i) = -icolalt(2,ida(2))
+                else
+                   icolalt(1,i) = icolalt(1,ida(2))
+                endif
+             elseif(da_color(1).eq.-3.and.da_color(2).eq.-3)then
+                maxcolor=maxcolor+1
+                icolalt(1,i) = maxcolor
+             else
+                call write_error(da_color(1), da_color(2), mo_color)
+             endif
+          elseif(mo_color.eq.-6) then ! color anti-sextet
+             if(da_color(1).eq.-6.and.da_color(2).eq.1)then
+                icolalt(1,i) = icolalt(1,ida(1))
+                icolalt(2,i) = icolalt(2,ida(1))
+             elseif(da_color(1).eq.-6.and.da_color(2).eq.8)then
+                if(icolalt(2,ida(1)).eq.icolalt(1,ida(2)))then
+                   icolalt(1,i) = icolalt(1,ida(1))
+                   icolalt(2,i) = icolalt(2,ida(2))
+                else
+                   icolalt(1,i) = -icolalt(2,ida(2))
+                   icolalt(2,i) = icolalt(2,ida(1))
+                endif
+             elseif(da_color(1).eq.-3.and.da_color(2).eq.-3)then
+                icolalt(1,i) = -icolalt(2,ida(1))
+                icolalt(2,i) = icolalt(2,ida(2))
+             else
+                call write_error(da_color(1), da_color(2), mo_color)
+             endif
+          elseif(mo_color.eq.6) then ! color sextet
+             if(da_color(1).eq.1.and.da_color(2).eq.6)then
+                icolalt(1,i) = icolalt(1,ida(2))
+                icolalt(2,i) = icolalt(2,ida(2))
+             elseif(da_color(1).eq.6.and.da_color(2).eq.8)then
+                if(icolalt(1,ida(1)).eq.icolalt(2,ida(2)))then
+                   icolalt(1,i) = icolalt(1,ida(2))
+                   icolalt(2,i) = icolalt(2,ida(1))
+                else
+                   icolalt(1,i) = icolalt(1,ida(1))
+                   icolalt(2,i) = -icolalt(1,ida(2))
+                endif
+             elseif(da_color(1).eq.3.and.da_color(2).eq.3)then
+                icolalt(1,i) = icolalt(1,ida(1))
+                icolalt(2,i) = -icolalt(1,ida(2))
+             else
+                call write_error(da_color(1), da_color(2), mo_color)
+             endif
+          elseif(mo_color.eq.8) then ! color octet
+             if(da_color(1).eq.-3.and.da_color(2).eq.3)then
+                icolalt(1,i) = icolalt(1,ida(2))
+                icolalt(2,i) = icolalt(2,ida(1))             
+             elseif(da_color(1).eq.1.and.da_color(2).eq.8)then
+                icolalt(1,i) = icolalt(1,ida(2))
+                icolalt(2,i) = icolalt(2,ida(2))             
+             elseif(da_color(1).eq.8.and.da_color(2).eq.8)then
+                if(icolalt(1,ida(1)).eq.icolalt(2,ida(2)))then
+                   icolalt(1,i) = icolalt(1,ida(2))
+                   icolalt(2,i) = icolalt(2,ida(1))
+                else
+                   icolalt(1,i) = icolalt(1,ida(1))
+                   icolalt(2,i) = icolalt(2,ida(2))
+                endif
+             elseif(da_color(1).eq.-6.and.da_color(2).eq.6)then
+                if(-icolalt(1,ida(1)).eq.icolalt(1,ida(2)))then
+                   icolalt(1,i) = -icolalt(1,ida(2))
+                   icolalt(2,i) = icolalt(2,ida(1))
+                elseif(icolalt(1,ida(1)).eq.icolalt(2,ida(2)))then
+                   icolalt(1,i) = icolalt(1,ida(2))
+                   icolalt(2,i) = icolalt(2,ida(1))
+                elseif(icolalt(2,ida(1)).eq.icolalt(1,ida(2)))then
+                   icolalt(1,i) = -icolalt(2,ida(2))
+                   icolalt(2,i) = -icolalt(1,ida(1))
+                else
+                   icolalt(1,i) = icolalt(1,ida(2))
+                   icolalt(2,i) = -icolalt(1,ida(1))
+                endif
+             else
+                call write_error(da_color(1), da_color(2), mo_color)
+             endif
           else
-c         Erraneous color assignment for propagator - set color to 0
-            icolalt(1,i) = 0
-            icolalt(2,i) = 0
+             call write_error(da_color(1), da_color(2), mo_color)
           endif
 c         Set tentative mothers
           jpart(2,i) = 1
@@ -295,3 +370,17 @@
 
       return
       end
+
+      subroutine write_error(ida1,ida2,imo)
+      implicit none
+      integer ida1,ida2,imo
+
+      open(unit=26,file='../../../error',status='unknown',err=999)
+      write(26,*) 'Error: Color combination ',ida1,ida2,
+     $     '->',imo,' not implemented in addmothers.f'
+      write(*,*) 'Error: Color combination ',ida1,ida2,
+     $     '->',imo,' not implemented in addmothers.f'
+      stop
+
+ 999  write(*,*) 'error'
+      end

=== modified file 'Template/SubProcesses/cluster.f'
--- Template/SubProcesses/cluster.f	2011-08-27 21:03:09 +0000
+++ Template/SubProcesses/cluster.f	2011-10-05 03:29:23 +0000
@@ -166,21 +166,22 @@
       end
 
 
-      subroutine filprp(ignum,idij)
+      subroutine filprp(iproc,ignum,idij)
 c**************************************************************************
 c     Include graph ignum in list for propagator idij
 c**************************************************************************
       implicit none
       include 'nexternal.inc'
+      include 'maxamps.inc'
       include 'cluster.inc'
-      integer ignum, idij, i
+      integer ignum, idij, iproc, i
 
       if(idij.gt.n_max_cl) return
-      do i=1,id_cl(idij,0)
-         if (id_cl(idij,i).eq.ignum) return
+      do i=1,id_cl(iproc,idij,0)
+         if (id_cl(iproc,idij,i).eq.ignum) return
       enddo
-      id_cl(idij,0)=id_cl(idij,0)+1
-      id_cl(idij,id_cl(idij,0))=ignum
+      id_cl(iproc,idij,0)=id_cl(iproc,idij,0)+1
+      id_cl(iproc,idij,id_cl(iproc,idij,0))=ignum
 c      print *,'Adding graph ',ignum,' to prop ',idij
       return
       end
@@ -200,9 +201,10 @@
       include 'genps.inc'
       include 'maxconfigs.inc'
       include 'nexternal.inc'
+      include 'maxamps.inc'
       include 'cluster.inc'
       include 'coupl.inc'
-      include 'maxamps.inc'
+      include 'message.inc'
       integer ignum, ipnum, ipids(nexternal,4,2:nexternal)
 C $B$ IFOREST $B$ !this is a tag for MadWeight
       integer i, iforest(2,-max_branch:-1,lmaxconfigs)
@@ -214,9 +216,8 @@
       INTEGER    n_max_cl_cg
       PARAMETER (n_max_cl_cg=n_max_cl*n_max_cg)
       data resmap/n_max_cl_cg*.false./
- 
 
-      Integer j, k, l, icmp(2)
+      Integer j, k, l, icmp(2), iproc
 
       double precision ZERO
       parameter (ZERO=0d0)
@@ -228,21 +229,24 @@
       data first_time /.true./
 
       integer combid
-      external combid
+      logical isjet
+      external combid,isjet
 
       if (first_time) then
          include 'props.inc'
          first_time=.false.
       endif
 
-c      write(*,*) 'graph,level,iupdown: ',ignum,ipnum,iupdown
+      if(btest(mlevel,4))
+     $     write(*,*) 'graph,level: ',ignum,ipnum
 
       filgrp=.false.
 C   Follow diagram tree down to last clustering
       do i=1,ipnum
          do j=i+1,ipnum
-c            write(*,*)'at ids   (',ipids(i,1,ipnum),',',ipids(i,2,ipnum),'), (',
-c     &           ipids(j,1,ipnum),',',ipids(j,2,ipnum),'), ',i,j
+            if(btest(mlevel,4))
+     $           write(*,*)'at ids   (',ipids(i,1,ipnum),',',ipids(i,2,ipnum),'), (',
+     $           ipids(j,1,ipnum),',',ipids(j,2,ipnum),'), ',i,j
             do k=-nexternal+1,-1
                if ((iforest(1,k,ignum).eq.ipids(i,2,ipnum).and.
      &              iforest(2,k,ignum).eq.ipids(j,2,ipnum)).or.
@@ -254,20 +258,39 @@
                   icmp(2)=ishft(1,nexternal)-1-icmp(1)
 c     Set pdg code for propagator
                   do l=1,2
-                     ipdgcl(icmp(l),ignum)=sprop(1,k,ignum)
-                     if(ipdgcl(icmp(l),ignum).eq.0)
-     $                    ipdgcl(icmp(l),ignum)=tprid(k,ignum)
-c                  write(*,*) 'add table entry for (',ipids(i,1,ipnum),
-c     &                 ',',ipids(j,1,ipnum),',',icmp,')','pdg: ',
-c     $                 ipdgcl(icmp,ignum)
-                     call filprp(ignum,icmp(l))
-
+                     do iproc=1,maxsproc
+                     if(sprop(iproc,k,ignum).ne.0)then
+                        ipdgcl(icmp(l),ignum,iproc)=sprop(iproc,k,ignum)
+c                       If this is radiation off heavy FS particle, set heavyrad to true
+                        if(isjet(ipdgcl(ipids(i,1,ipnum),ignum,iproc)).and.
+     $                       .not.isjet(ipdgcl(ipids(j,1,ipnum),ignum,iproc)).and.
+     $                       ipdgcl(ipids(j,1,ipnum),ignum,iproc).eq.sprop(iproc,k,ignum).or.
+     $                       isjet(ipdgcl(ipids(j,1,ipnum),ignum,iproc)).and.
+     $                       .not.isjet(ipdgcl(ipids(i,1,ipnum),ignum,iproc)).and.
+     $                       ipdgcl(ipids(i,1,ipnum),ignum,iproc).eq.sprop(iproc,k,ignum))then
+                           heavyrad(ignum) = .true.
+                        endif
+                     else if(tprid(k,ignum).ne.0)then
+                        ipdgcl(icmp(l),ignum,iproc)=tprid(k,ignum)
+                     else if(ipnum.eq.3)then
+                        ipdgcl(icmp(l),ignum,iproc)=ipdgcl(2,ignum,iproc)
+                     else
+                        ipdgcl(icmp(l),ignum,iproc)=0
+                        cycle
+                     endif
+                     if(btest(mlevel,4))
+     $                    write(*,*) 'add table entry for (',ipids(i,1,ipnum),
+     &                    ',',ipids(j,1,ipnum),',',icmp(l),')',
+     $                    'proc: ',iproc,
+     $                    ', pdg: ',ipdgcl(icmp(l),ignum,iproc)
+                     call filprp(iproc,ignum,icmp(l))
 c               Insert graph in list of propagators
                      if(pwidth(k,ignum).gt.ZERO) then
 c                    write(*,*)'Adding resonance ',ignum,icmp
                         resmap(icmp(l),ignum)=.true.
                      endif
                   enddo
+                  enddo
 c     proceed w/ next table, since there is no possibility,
 c     to combine the same particle in another way in this graph
                   ipids(i,1,ipnum-1)=icmp(1)
@@ -288,11 +311,6 @@
                      ipids(l,3,ipnum)=l+1
                      ipids(l,4,ipnum)=0
                   enddo
-c                  do l=1,ipnum
-c                     write(*,*) 'new: ipids(',l,') = (',ipids(l,1,ipnum),
-c     &                    ',',ipids(l,2,ipnum),',',ipids(l,3,ipnum),',',ipids(l,4,ipnum),')',
-c     $                    ' pdg: ', ipdgcl(ipids(l,1,ipnum),ignum)
-c                  enddo
                   if(ipnum.eq.2)then
 c                 Done with this diagram
                      return
@@ -317,20 +335,21 @@
       include 'genps.inc'
       include 'maxconfigs.inc'
       include 'nexternal.inc'
+      include 'maxamps.inc'
       include 'cluster.inc'
       include 'run.inc'
-      include 'maxamps.inc'
+      include 'message.inc'
 C $B$ IFOREST $B$ !this is a tag for MadWeight
       integer mapconfig(0:lmaxconfigs), this_config
       common/to_mconfigs/mapconfig, this_config
 C $E$ IFOREST $E$ !this is a tag for MadWeight
-      integer i, j, inpids, ipids(nexternal,4,2:nexternal)
+      integer i, j, inpids, iproc, ipids(nexternal,4,2:nexternal)
       integer start_config,end_config
       integer idup(nexternal,maxproc,maxsproc)
       integer mothup(2,nexternal)
       integer icolup(2,nexternal,maxflow,maxsproc)
       include 'leshouche.inc'
-
+      
       logical filgrp
       external filgrp
 
@@ -341,22 +360,41 @@
          start_config=1
          end_config=mapconfig(0)
       endif
-      do i=1,n_max_cl
-         id_cl(i,0)=0
+      do iproc=1,maxsproc
+         do i=1,n_max_cl
+            id_cl(iproc,i,0)=0
+         enddo
       enddo
       do i=start_config,end_config
+         heavyrad(i)=.false.
 c         write (*,*) ' at graph ',i
          do j=1,nexternal
             ipids(j,1,nexternal)=ishft(1,j-1)
             ipids(j,2,nexternal)=j
             ipids(j,3,nexternal)=0
             ipids(j,4,nexternal)=0
-            ipdgcl(ipids(j,1,nexternal),i)=idup(j,1,1)
+            do iproc=1,maxsproc
+               ipdgcl(ipids(j,1,nexternal),i,iproc)=idup(j,1,iproc)
+            enddo
          enddo
          inpids=nexternal
 c         print *,'Inserting graph ',i
  10      if (filgrp(i,inpids,ipids)) goto 10
-      enddo
+         if(btest(mlevel,4).and.heavyrad(i)) then
+            write(*,*)' set heavyrad of ',i,' to T'
+         endif
+      enddo
+c     Ensure that there are some allowed clusterings
+      do i=start_config,end_config
+         if(.not.heavyrad(i)) goto 20
+      enddo
+      if(btest(mlevel,4)) then
+         write(*,*)' Reset all heavyrad to .false.'
+      endif      
+      do i=start_config,end_config
+         heavyrad(i)=.false.
+      enddo
+ 20   continue
       filmap=.true.
       return
       end
@@ -431,28 +469,35 @@
 c**************************************************************************
       implicit none
       include 'nexternal.inc'
+      include 'maxamps.inc'
       include 'cluster.inc'
 
       integer idij,nbw,ibwlist(nexternal),icgs(0:n_max_cg)
       logical foundbw
       integer i, ii, j, jj, il, igsbk(0:n_max_cg)
 
+c     IPROC has the present process number
+      INTEGER IMIRROR,IPROC
+      COMMON/TO_MIRROR/IMIRROR, IPROC
+
       findmt=.false.
 c     if first clustering, set possible graphs
       if (icgs(0).eq.0) then
          ii=0
-         do i=1,id_cl(idij,0)
+         do i=1,id_cl(iproc,idij,0)
+c        check if this diagram has radiation off heavy particle
+           if(heavyrad(id_cl(iproc,idij,i))) cycle
 c        check if we have constraint from onshell resonances
            foundbw=.true.
            do j=1,nbw
-             if(resmap(ibwlist(j),id_cl(idij,i)))then
+             if(resmap(ibwlist(j),id_cl(iproc,idij,i)))then
                cycle
              endif
              foundbw=.false.
  10        enddo
-           if(nbw.eq.0.or.foundbw)then
+           if((nbw.eq.0.or.foundbw))then
               ii=ii+1
-              icgs(ii)=id_cl(idij,i)
+              icgs(ii)=id_cl(iproc,idij,i)
            endif
          enddo
          icgs(0)=ii
@@ -465,17 +510,17 @@
          j=1
          ii=0
          do i=1,icgs(0)
-            if(j.le.id_cl(idij,0).and.icgs(i).eq.id_cl(idij,j))then
+            if(j.le.id_cl(iproc,idij,0).and.icgs(i).eq.id_cl(iproc,idij,j))then
                ii=ii+1
-               icgs(ii)=id_cl(idij,j)
+               icgs(ii)=id_cl(iproc,idij,j)
                j=j+1
-            else if(j.le.id_cl(idij,0).and.icgs(i).gt.id_cl(idij,j)) then
-               do while(icgs(i).gt.id_cl(idij,j).and.j.le.id_cl(idij,0))
+            else if(j.le.id_cl(iproc,idij,0).and.icgs(i).gt.id_cl(iproc,idij,j)) then
+               do while(icgs(i).gt.id_cl(iproc,idij,j).and.j.le.id_cl(iproc,idij,0))
                   j=j+1
                enddo
-               if(j.le.id_cl(idij,0).and.icgs(i).eq.id_cl(idij,j))then
+               if(j.le.id_cl(iproc,idij,0).and.icgs(i).eq.id_cl(iproc,idij,j))then
                   ii=ii+1
-                  icgs(ii)=id_cl(idij,j)
+                  icgs(ii)=id_cl(iproc,idij,j)
                endif
             endif
          enddo
@@ -497,9 +542,11 @@
       include 'run.inc'
       include 'genps.inc'
       include 'nexternal.inc'
+      include 'maxamps.inc'
       include 'cluster.inc'
       include 'message.inc'
       include 'maxconfigs.inc'
+
       real*8 p(0:3,nexternal), pcmsp(0:3), p1(0:3)
       real*8 pi(0:3), nr(0:3), pz(0:3)
       integer i, j, k, n, idi, idj, idij, icgs(0:n_max_cg)
@@ -511,7 +558,6 @@
       common/to_mconfigs/mapconfig, this_config
 
       integer nbw,ibwlist(nexternal)
-      logical isbw(n_max_cl)
       data isbw/n_max_cl*.false./
 
       data (pz(i),i=0,3)/1d0,0d0,0d0,1d0/
@@ -569,7 +615,7 @@
                      if(isbw(idij))then
                        pt2ij(idij)=SumDot(pcl(0,idi),pcl(0,idj),1d0)
                        if (btest(mlevel,4))
-     $                    print *,'Mother ',idij,' has ptij ',
+     $                    write(*,*)'Mother ',idij,' has ptij ',
      $                    sqrt(pt2ij(idij))
                      else
                        if(ktscheme.eq.2)then
@@ -578,6 +624,7 @@
                          pt2ij(idij)=dj(pcl(0,idi),pcl(0,idj))
                        endif
                      endif
+                     zij(idij)=0d0
                   else
 c     initial state clustering, only if hadronic collision
 c     check whether 2->(n-1) process w/ cms energy > 0 remains
@@ -585,6 +632,7 @@
                      if(ickkw.eq.2.or.ktscheme.eq.2)then
                         pt2ij(idij)=pyjb(pcl(0,idi),
      $                    pcl(0,idj),pcl(0,iwinp),zij(idij))
+                        zij(idij)=0d0
                      else
                         pt2ij(idij)=djb(pcl(0,idi))
                         zij(idij)=zclus(pcl(0,idi),pcl(0,idj),pcl(0,iwinp))
@@ -613,11 +661,14 @@
 c     Take care of special 2 -> 1 case
       if (nexternal.eq.3.and.nincoming.eq.2) then
          n=1
-         imocl(n)=idij
-         idacl(n,1)=idi
-         idacl(n,2)=idj
-         pt2ijcl(n)=pcl(4,idi)
+c     Make sure that initial-state particles are daughters
+         idacl(n,1)=imap(1,2)
+         idacl(n,2)=imap(2,2)
+         imocl(n)=imap(3,2)
+         pt2ijcl(n)=pcl(4,imocl(n))
          zcl(n)=0.
+         igraphs(0)=1
+         igraphs(1)=this_config
          cluster=.true.
          clustered=.true.
          return
@@ -649,7 +700,7 @@
 c     is clustering
 c     Set mt2ij to m^2+pt^2 
             mt2ij(n)=djb(pcl(0,idacl(n,2)))
-            if (btest(mlevel,3)) then
+            if (btest(mlevel,1)) then
                write(*,*)'mtij(',n,') for ',idacl(n,2),' is ',sqrt(mt2ij(n)),
      $              ' (cf ',sqrt(pt2ijcl(n)),')'
             endif
@@ -699,7 +750,7 @@
            if(isbw(imocl(n)))then
              pcl(4,imocl(n))=pt2ijcl(n)
              if (btest(mlevel,4))
-     $          print *,'Mother ',imocl(n),' has mass**2 ',
+     $          write(*,*) 'Mother ',imocl(n),' has mass**2 ',
      $          pcl(4,imocl(n))
            endif
          endif
@@ -728,15 +779,10 @@
                   call boostx(p1(0),pcmsp(0),pcl(0,imap(3,2)))
                endif
             endif
-c         Make sure that final-state particle is always among daughters
+c         Make sure that initial-state particle is always among daughters
             idacl(n+1,1)=imap(1,2)
-            idacl(n+1,2)=imap(3,2)
-            imocl(n+1)=imap(2,2)
-c         If mother is initial state leg (2nd beam), chose other leg
-            if (imocl(n+1).eq.4)then
-              idacl(n+1,1)=imap(2,2)
-              imocl(n+1)=imap(1,2)
-            endif
+            idacl(n+1,2)=imap(2,2)
+            imocl(n+1)=imap(3,2)
 
 c            if(pcl(0,imocl(n)).gt.0d0)then
             pt2ijcl(n+1)=djb(pcl(0,imap(3,2)))
@@ -801,7 +847,7 @@
                            if(isbw(idij))then
                              pt2ij(idij)=SumDot(pcl(0,idi),pcl(0,idj),1d0)
                              if (btest(mlevel,4))
-     $                          print *,'Mother ',idij,' has ptij ',
+     $                          write(*,*) 'Mother ',idij,' has ptij ',
      $                          sqrt(pt2ij(idij))
                            else
                              if(ktscheme.eq.2)then

=== modified file 'Template/SubProcesses/cluster.inc'
--- Template/SubProcesses/cluster.inc	2011-03-09 18:29:27 +0000
+++ Template/SubProcesses/cluster.inc	2011-10-05 03:29:23 +0000
@@ -5,18 +5,20 @@
       include 'ngraphs.inc'
 c     parameters for clustering:
 c     id_cl gives diagrams for propagators     
-      integer id_cl(n_max_cl,0:n_max_cg)
-      common/cl_map/id_cl
+      integer id_cl(maxsproc,n_max_cl,0:n_max_cg)
+      logical heavyrad(n_max_cg)
+      common/cl_map/id_cl,heavyrad
 c     resmap gives potential resonances for diagrams
       logical resmap(n_max_cl,n_max_cg)
       common/res_map/resmap
       double precision pt2ijcl(nexternal),zcl(nexternal),mt2ij(nexternal),mt2last
       double precision pcl(0:4,n_max_cl) ! 4 is mass**2
+      logical isbw(n_max_cl)
       integer imocl(nexternal), idacl(nexternal,2), igraphs(0:n_max_cg)
-      integer ipdgcl(n_max_cl,n_max_cg)
+      integer ipdgcl(n_max_cl,n_max_cg,maxsproc)
       logical clustered
       common/cl_val/pcl,pt2ijcl,zcl,mt2ij,mt2last,imocl,idacl,igraphs,ipdgcl,
-     &    clustered
+     &    isbw,clustered
 c     parameters for qcd jets
       double precision ptclus(nexternal)
       common/cl_jets/ptclus

=== modified file 'Template/SubProcesses/cuts.f'
--- Template/SubProcesses/cuts.f	2011-08-07 19:55:37 +0000
+++ Template/SubProcesses/cuts.f	2011-10-05 03:29:23 +0000
@@ -388,11 +388,20 @@
             if(debug) write (*,*) dsqrt(s_min(j,i)),dsqrt(s_max(j,i))
             if(s_min(j,i).gt.0.or.s_max(j,i).lt.1d5) then
                tmp=SumDot(p(0,i),p(0,j),+1d0)
-               notgood=(tmp .lt. s_min(j,i).or.tmp .gt. s_max(j,i)) 
-               if (notgood) then
-                  if(debug) write (*,*) i,j,' -> fails'
-                  passcuts=.false.
-                  return
+               if(s_min(j,i).le.s_max(j,i))then
+                  notgood=(tmp .lt. s_min(j,i).or.tmp .gt. s_max(j,i)) 
+                  if (notgood) then
+                     if(debug) write (*,*) i,j,' -> fails'
+                     passcuts=.false.
+                     return
+                  endif
+               else
+                  notgood=(tmp .lt. s_min(j,i).and.tmp .gt. s_max(j,i)) 
+                  if (notgood) then
+                     if(debug) write (*,*) i,j,' -> fails'
+                     passcuts=.false.
+                     return
+                  endif
                endif
             endif
          enddo

=== modified file 'Template/SubProcesses/initcluster.f'
--- Template/SubProcesses/initcluster.f	2011-03-09 18:29:27 +0000
+++ Template/SubProcesses/initcluster.f	2011-10-05 03:29:23 +0000
@@ -5,6 +5,7 @@
       include 'message.inc'
       include 'run.inc'
       include 'nexternal.inc'
+      include 'maxamps.inc'
       include 'cluster.inc'
 C
 C     SPECIAL CUTS
@@ -13,7 +14,7 @@
       real*8 xetamin,xqcut,deltaeta
       common /to_specxpt/xptj,xptb,xpta,xptl,xmtc,xetamin,xqcut,deltaeta
 
-      integer i,j
+      integer i,j,iproc
       logical filmap, cluster
       external filmap, cluster
 
@@ -46,13 +47,16 @@
         stop
       endif
       if (btest(mlevel,3)) then
-        do i=1,n_max_cl
-          write(*,*) 'prop ',i,' in'
-          do j=1,id_cl(i,0)
-            write(*,*) '  graph ',id_cl(i,j)
-          enddo
+        do iproc=1,maxsproc
+           write(*,*)'for proc ',iproc
+           do i=1,n_max_cl
+              write(*,*) 'prop ',i,' in'
+              do j=1,id_cl(iproc,i,0)
+                 write(*,*) '  graph ',id_cl(iproc,i,j)
+              enddo
+           enddo
+           write(*,*)'ok'
         enddo
-        write(*,*)'ok'
       endif
       igraphs(0)=0
        

=== modified file 'Template/SubProcesses/makefile'
--- Template/SubProcesses/makefile	2011-07-10 22:43:49 +0000
+++ Template/SubProcesses/makefile	2011-10-05 03:29:23 +0000
@@ -13,11 +13,12 @@
 
 # Source files
 
-PROCESS= driver.o myamp.o genps.o unwgt.o setcuts.o \
+PROCESS= driver.o myamp.o genps.o unwgt.o setcuts.o get_color.o \
          cuts.o cluster.o reweight.o initcluster.o addmothers.o setscales.o \
          $(patsubst %.f,%.o,$(wildcard auto_dsig*.f)) \
          $(patsubst %.f,%.o,$(wildcard matrix*.f))
-SYMMETRY = myamp.o genps.o initcluster.o setscales.o reweight.o\
+SYMMETRY = symmetry.o genps.inc setcuts.o cuts.o cluster.o \
+	   myamp.o genps.o initcluster.o setscales.o reweight.o get_color.o \
            $(patsubst %.f,%.o,$(wildcard matrix*.f))
 
 # Binaries
@@ -25,8 +26,8 @@
 $(PROG): $(PROCESS) auto_dsig.o $(LIBS)
 	$(FC) $(FFLAGS) -o $(PROG) $(PROCESS) $(LINKLIBS)
 
-gensym: symmetry.o genps.inc setcuts.o cuts.o cluster.o $(SYMMETRY) configs.inc $(LIBS)
-	$(FC) $(FFLAGS) -o gensym symmetry.o setcuts.o cuts.o cluster.o $(SYMMETRY) $(LINKLIBS)
+gensym: $(SYMMETRY) configs.inc $(LIBS)
+	$(FC) $(FFLAGS) -o gensym $(SYMMETRY) $(LINKLIBS)
 
 # Dependencies
 

=== modified file 'Template/SubProcesses/myamp.f'
--- Template/SubProcesses/myamp.f	2011-09-05 06:55:55 +0000
+++ Template/SubProcesses/myamp.f	2011-10-05 03:29:23 +0000
@@ -110,7 +110,7 @@
 c
       double precision xp(0:3,-nexternal:nexternal)
       double precision mpole(-nexternal:0),shat,tsgn
-      integer i,j,iconfig
+      integer i,j,iconfig,iproc
 
       double precision pmass(-nexternal:0,lmaxconfigs)
       double precision pwidth(-nexternal:0,lmaxconfigs)
@@ -145,7 +145,7 @@
       integer icolup(2,nexternal,maxflow,maxsproc)
       include 'leshouche.inc'
 
-      logical gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
+      integer gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
       include 'decayBW.inc'
 c
 c     External
@@ -190,7 +190,15 @@
       enddo
       nbw = 0
       tsgn    = +1d0
-      do i=-1,-(nexternal-3),-1              !Loop over propagators
+c     Find non-zero process number
+      do iproc=1,maxsproc
+         if(sprop(iproc,-1,iconfig).ne.0) goto 10
+      enddo
+ 10   continue
+c     If no non-zero sprop, set iproc to 1
+      if(iproc.gt.maxsproc) iproc=1
+c     Start loop over propagators
+      do i=-1,-(nexternal-3),-1
          onbw(i) = .false.
          if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
          do j=0,3
@@ -209,16 +217,21 @@
             onshell = (abs(xmass - pmass(i,iconfig)) .lt.
      $           bwcutoff*pwidth(i,iconfig))
             if(onshell)then
-c           Only allow onshell if no "decay" to identical particle
+c     Remove on-shell forbidden s-channels (gForceBW=2) (JA 2/10/11)
+              if(gForceBW(i,iconfig).eq.2) then
+                 cut_bw = .true.
+                 return               
+              endif
+c           Only allow OnBW if no "decay" to identical particle
               OnBW(i) = .true.
               idenpart=0
               do j=1,2
                 ida(j)=iforest(j,i,iconfig)
                 if(ida(j).lt.0) then
-                   if(sprop(1,i,iconfig).eq.sprop(1,ida(j),iconfig))
+                   if(sprop(iproc,i,iconfig).eq.sprop(iproc,ida(j),iconfig))
      $                  idenpart=ida(j)
                 elseif (ida(j).gt.0) then
-                   if(sprop(1,i,iconfig).eq.IDUP(ida(j),1,1))
+                   if(sprop(iproc,i,iconfig).eq.IDUP(ida(j),1,1))
      $                  idenpart=ida(j)
                 endif
               enddo
@@ -227,10 +240,10 @@
                 OnBW(i)=.false.
 c           Else remove if daughter forced to be onshell
               elseif(idenpart.lt.0)then
-                 if(gForceBW(idenpart, iconfig)) then
+                 if(gForceBW(idenpart, iconfig).eq.1) then
                     OnBW(i)=.false.
 c           Else remove daughter if forced to be onshell
-                 elseif(gForceBW(i, iconfig)) then
+                 elseif(gForceBW(i, iconfig).eq.1) then
                     OnBW(idenpart)=.false.
 c           Else remove either this resonance or daughter, which is closer to mass shell
                  elseif(abs(xmass-pmass(i,iconfig)).gt.
@@ -242,13 +255,10 @@
                     OnBW(idenpart)=.false.
                  endif
               endif
-            endif
-c
-c     Check if we are supposed to cut forced bw (JA 4/8/11)
-c
-            if (gForceBW(i, iconfig) .and. .not. onshell)then
-               cut_bw = .true.
-               return
+            else if (gForceBW(i, iconfig).eq.1) then ! .not. onshell
+c             Check if we are supposed to cut forced bw (JA 4/8/11)
+              cut_bw = .true.
+              return
             endif
 c
 c     Here we set onshell for phase space integration (JA 4/8/11)
@@ -263,7 +273,6 @@
                return
             endif
          endif
-
       enddo
       end
 
@@ -293,7 +302,7 @@
       double precision tsgn, xo, a
       double precision x1,x2,xk(nexternal)
       double precision dr,mtot,etot,stot,xqfact
-      integer i, iconfig, l1, l2, j, nt, nbw
+      integer i, iconfig, l1, l2, j, nt, nbw, iproc
       integer iden_part(-max_branch:-1)
 
       double precision pmass(-nexternal:0,lmaxconfigs)
@@ -305,7 +314,7 @@
       integer icolup(2,nexternal,maxflow,maxsproc)
       include 'leshouche.inc'
 
-      logical gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
+      integer gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
       include 'decayBW.inc'
 c
 c     Global
@@ -341,7 +350,7 @@
       common /to_BW/ lbw
 
       include 'coupl.inc'
-
+      include 'cuts.inc'
 c
 c     External
 c
@@ -377,30 +386,39 @@
          spole(i)=0
          swidth(i)=0
       enddo
-      do i=-1,-(nexternal-3),-1              !Find all the propagotors
+c     Find non-zero process number
+      do iproc=1,maxsproc
+         if(sprop(iproc,-1,iconfig).gt.0) goto 10
+      enddo
+ 10   continue
+c     If no non-zero sprop, set iproc to 1
+      if(iproc.ge.maxsproc.and.sprop(maxsproc,-1,iconfig).eq.0)
+     $     iproc=1
+c     Start loop over propagators
+      do i=-1,-(nexternal-3),-1
 c     JA 3/31/11 Keep track of identical particles (i.e., radiation vertices)
 c     by tracing the particle identity from the external particle.
          if(iforest(1,i,iconfig).gt.0) then
-             if (sprop(1,i,iconfig).eq.idup(iforest(1,i,iconfig),1,1))
-     $        iden_part(i) = sprop(1,i,iconfig)
+             if (sprop(iproc,i,iconfig).eq.idup(iforest(1,i,iconfig),1,1))
+     $        iden_part(i) = sprop(iproc,i,iconfig)
           endif
          if(iforest(2,i,iconfig).gt.0) then
-            if(sprop(1,i,iconfig).eq.idup(iforest(2,i,iconfig),1,1))
-     $           iden_part(i) = sprop(1,i,iconfig)
+            if(sprop(iproc,i,iconfig).eq.idup(iforest(2,i,iconfig),1,1))
+     $           iden_part(i) = sprop(iproc,i,iconfig)
          endif
          if(iforest(1,i,iconfig).lt.0) then
             if((iden_part(iforest(1,i,iconfig)).ne.0.and.
-     $        sprop(1,i,iconfig).eq.iden_part(iforest(1,i,iconfig)) .or.
-     $        gforcebw(iforest(1,i,iconfig),iconfig).and.
-     $        sprop(1,i,iconfig).eq.sprop(1,iforest(1,i,iconfig),iconfig)))
-     $       iden_part(i) = sprop(1,i,iconfig)
+     $        sprop(iproc,i,iconfig).eq.iden_part(iforest(1,i,iconfig)) .or.
+     $        gforcebw(iforest(1,i,iconfig),iconfig).eq.1.and.
+     $        sprop(iproc,i,iconfig).eq.sprop(iproc,iforest(1,i,iconfig),iconfig)))
+     $       iden_part(i) = sprop(iproc,i,iconfig)
          endif
          if(iforest(2,i,iconfig).lt.0) then
             if((iden_part(iforest(2,i,iconfig)).ne.0.and.
-     $        sprop(1,i,iconfig).eq.iden_part(iforest(2,i,iconfig)).or.
-     $        gforcebw(iforest(2,i,iconfig),iconfig).and.
-     $        sprop(1,i,iconfig).eq.sprop(1,iforest(2,i,iconfig),iconfig)))
-     $           iden_part(i) = sprop(1,i,iconfig)
+     $        sprop(iproc,i,iconfig).eq.iden_part(iforest(2,i,iconfig)).or.
+     $        gforcebw(iforest(2,i,iconfig),iconfig).eq.1.and.
+     $        sprop(iproc,i,iconfig).eq.sprop(iproc,iforest(2,i,iconfig),iconfig)))
+     $           iden_part(i) = sprop(iproc,i,iconfig)
          endif
          if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
          if (tsgn .eq. 1d0) then                         !s channel
@@ -427,7 +445,7 @@
 c              JA 6/8/2011 Set xe(i) for resonances
                if (lbw(nbw).eq.1) then
                   xm(i) = max(xm(i), pmass(i,iconfig)-5d0*pwidth(i,iconfig))
-               else if (gforcebw(i,iconfig)) then
+               else if (gforcebw(i,iconfig).eq.1) then
                   xm(i) = max(xm(i), pmass(i,iconfig)-bwcutoff*pwidth(i,iconfig))
                endif
             endif
@@ -529,6 +547,22 @@
 c     Set minimum based on: 1) required energy 2) resonances 3) 1/10000 of sqrt(s)
          i = 3*(nexternal-2) - 4 + 1
          xo = max(min(etot**2/stot, 1d0-1d-8),1d0/stot)
+c        Take into account special cuts
+         xo = max(xo, xptj*dabs(xptj)/stot)
+         xo = max(xo, xptb*dabs(xptb)/stot)
+         xo = max(xo, xpta*dabs(xpta)/stot)
+         xo = max(xo, xptl*dabs(xptl)/stot)
+         xo = max(xo, xmtc*dabs(xmtc)/stot)
+         xo = max(xo, htjmin**2/stot)
+         xo = max(xo, ptj1min**2/stot)
+         xo = max(xo, ptj2min**2/stot)
+         xo = max(xo, ptj3min**2/stot)
+         xo = max(xo, ptj4min**2/stot)
+         xo = max(xo, ht2min**2/stot)
+         xo = max(xo, ht3min**2/stot)
+         xo = max(xo, ht4min**2/stot)
+         xo = max(xo, misset**2/stot)
+         xo = max(xo, ptllmin**2/stot)
          if (swidth(i).eq.0.and.xo.eq.1d0/stot) then
             write(*,*) 'Warning: No minimum found for integration'
             write(*,*) '         Setting minimum to ',1d0/stot

=== modified file 'Template/SubProcesses/reweight.f'
--- Template/SubProcesses/reweight.f	2011-08-27 21:03:09 +0000
+++ Template/SubProcesses/reweight.f	2011-10-05 03:29:23 +0000
@@ -5,6 +5,7 @@
       implicit none
       include 'nexternal.inc'
       include 'message.inc'
+      include 'maxamps.inc'
       include 'cluster.inc'
       include 'sudakov.inc'
       include 'run.inc'
@@ -76,6 +77,7 @@
       implicit none
       include 'message.inc'
       include 'nexternal.inc'
+      include 'maxamps.inc'
       include 'cluster.inc'      
       integer ipdg,imode
       double precision q0, Q11
@@ -135,14 +137,9 @@
 c**************************************************
       implicit none
       integer ipdg, irfl
-
-      isqcd=.true.
-
-c     Assume that QCD particles have pdg codes that are (multiples of 1M) +
-c     1-10 or 21
-      irfl=mod(abs(ipdg),1000000)
-      if (irfl.ge.11.and.irfl.ne.21) isqcd=.false.
-c      write(*,*)'iqcd? pdg = ',ipdg,' -> ',irfl,' -> ',isqcd
+      integer get_color
+
+      isqcd=(iabs(get_color(ipdg)).gt.1)
 
       return
       end
@@ -200,42 +197,33 @@
       integer imo,ida1,ida2,i,idmo,idda1,idda2
       integer ipdg(n_max_cl),ipart(2,n_max_cl)
 
-      do i=1,2
-        ipart(i,imo)=0
-      enddo
-
       idmo=ipdg(imo)
       idda1=ipdg(ida1)
       idda2=ipdg(ida2)
 
-      if (btest(mlevel,1)) then
-        write(*,*) ' updating ipart for: ',ida1,ida2,imo
+      if (btest(mlevel,4)) then
+        write(*,*) ' updating ipart for: ',ida1,ida2,' -> ',imo
       endif
 
-        if (btest(mlevel,1)) then
-          write(*,*) ' daughters: ',(ipart(i,ida1),i=1,2),(ipart(i,ida2),i=1,2)
-        endif
+      if (btest(mlevel,4)) then
+         write(*,*) ' daughters: ',(ipart(i,ida1),i=1,2),(ipart(i,ida2),i=1,2)
+      endif
 
 c     IS clustering - just transmit info on incoming line
       if((ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2).or.
      $   (ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2))then
-        if(ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2)
-     $     ipart(1,imo)=ipart(1,ida2)        
-        if(ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2)
-     $     ipart(1,imo)=ipart(1,ida1)
-        if (btest(mlevel,1)) then
-          write(*,*) ' -> ',(ipart(i,imo),i=1,2)
-c     Set intermediate particle identity
-          if(iabs(idmo).lt.6)then
-            if(iabs(idda1).lt.6) ipdg(imo)=-idda1
-            if(iabs(idda2).lt.6) ipdg(imo)=-idda2
-            idmo=ipdg(imo)
-            if (btest(mlevel,1)) then
-              write(*,*) ' particle identities: ',idda1,idda2,idmo
-            endif
-          endif
-        endif
-        return
+         ipart(2,imo)=0
+         if(ipart(1,ida1).le.2.and.ipart(1,ida2).le.2)then
+c           This is last clustering - keep mother ipart
+            ipart(1,imo)=ipart(1,imo)
+         elseif(ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2)then
+            ipart(1,imo)=ipart(1,ida2)        
+         elseif(ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2)then
+            ipart(1,imo)=ipart(1,ida1)
+         endif
+         if (btest(mlevel,4))
+     $        write(*,*) ' -> ',(ipart(i,imo),i=1,2)
+         return
       endif        
 
 c     FS clustering
@@ -262,33 +250,25 @@
       else if(idmo.eq.idda1.or.idmo.eq.idda1+sign(1,idda2))then
 c     quark -> quark-gluon or quark-Z or quark-h or quark-W
         ipart(1,imo)=ipart(1,ida1)
+        ipart(2,imo)=0
       else if(idmo.eq.idda2.or.idmo.eq.idda2+sign(1,idda1))then
 c     quark -> gluon-quark or Z-quark or h-quark or W-quark
         ipart(1,imo)=ipart(1,ida2)
+        ipart(2,imo)=0
       else
 c     Color singlet
          ipart(1,imo)=ipart(1,ida1)
          ipart(2,imo)=ipart(1,ida2)
       endif
       
-      if (btest(mlevel,1)) then
+      if (btest(mlevel,4)) then
         write(*,*) ' -> ',(ipart(i,imo),i=1,2)
       endif
 
-c     Set intermediate particle identity
-      if(iabs(idmo).lt.6)then
-        if(iabs(idda1).lt.6) ipdg(imo)=idda1
-        if(iabs(idda2).lt.6) ipdg(imo)=idda2
-        idmo=ipdg(imo)
-        if (btest(mlevel,1)) then
-          write(*,*) ' particle identities: ',idda1,idda2,idmo
-        endif
-      endif
-
       return
       end
       
-      logical function isjetvx(imo,ida1,ida2,ipdg,ipart)
+      logical function isjetvx(imo,ida1,ida2,ipdg,ipart,islast)
 c***************************************************
 c   Checks if a qcd vertex generates a jet
 c***************************************************
@@ -299,7 +279,7 @@
 
       integer imo,ida1,ida2,idmo,idda1,idda2,i
       integer ipdg(n_max_cl),ipart(2,n_max_cl)
-      logical isqcd,isjet
+      logical isqcd,isjet,islast
       external isqcd,isjet
 
       idmo=ipdg(imo)
@@ -307,7 +287,7 @@
       idda2=ipdg(ida2)
 
 c     Check QCD vertex
-      if(.not.isqcd(idmo).or..not.isqcd(idda1).or.
+      if(islast.or..not.isqcd(idmo).or..not.isqcd(idda1).or.
      &     .not.isqcd(idda2)) then
          isjetvx = .false.
          return
@@ -317,7 +297,7 @@
       if((ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2).or.
      $   (ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2))then
 c     Check if ida1 is outgoing parton or ida2 is outgoing parton
-        if(ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2.and.isjet(idda1).or.
+         if(ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2.and.isjet(idda1).or.
      $        ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2.and.isjet(idda2))then
            isjetvx=.true.
         else
@@ -336,6 +316,54 @@
       return
       end
 
+      logical function ispartonvx(imo,ida1,ida2,ipdg,ipart,islast)
+c***************************************************
+c   Checks if a qcd vertex generates a jet
+c***************************************************
+      implicit none
+
+      include 'ncombs.inc'
+      include 'nexternal.inc'
+
+      integer imo,ida1,ida2,idmo,idda1,idda2,i
+      integer ipdg(n_max_cl),ipart(2,n_max_cl)
+      logical isqcd,isparton,islast
+      external isqcd,isparton
+
+      idmo=ipdg(imo)
+      idda1=ipdg(ida1)
+      idda2=ipdg(ida2)
+
+c     Check QCD vertex
+      if(.not.isqcd(idmo).or..not.isqcd(idda1).or.
+     &     .not.isqcd(idda2)) then
+         ispartonvx = .false.
+         return
+      endif
+
+c     IS clustering
+      if((ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2).or.
+     $   (ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2))then
+c     Check if ida1 is outgoing parton or ida2 is outgoing parton
+         if(.not.islast.and.ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2.and.isparton(idda1).or.
+     $        ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2.and.isparton(idda2))then
+           ispartonvx=.true.
+        else
+           ispartonvx=.false.
+        endif
+        return
+      endif        
+
+c     FS clustering
+      if(isparton(idda1).or.isparton(idda2))then
+         ispartonvx=.true.
+      else
+         ispartonvx=.false.
+      endif
+      
+      return
+      end
+
       logical function setclscales(p)
 c**************************************************
 c   reweight the hard me according to ckkw
@@ -347,6 +375,7 @@
       include 'genps.inc'
       include 'maxconfigs.inc'
       include 'nexternal.inc'
+      include 'maxamps.inc'
       include 'cluster.inc'
       include 'run.inc'
       include 'coupl.inc'
@@ -355,6 +384,11 @@
 C   
       DOUBLE PRECISION P(0:3,NEXTERNAL)
 
+C   global variables
+C     Present process number
+      INTEGER IMIRROR,IPROC
+      COMMON/TO_MIRROR/IMIRROR, IPROC
+
 C   local variables
       integer i, j, idi, idj
       real*8 PI
@@ -365,23 +399,23 @@
       real*8 xptj,xptb,xpta,xptl,xmtc
       real*8 xetamin,xqcut,deltaeta
       common /to_specxpt/xptj,xptb,xpta,xptl,xmtc,xetamin,xqcut,deltaeta
+c     q2bck holds the central q2fact scales
       real*8 q2bck(2)
-      save q2bck
-      include 'maxamps.inc'
+      common /to_q2bck/q2bck
       double precision asref, pt2prev(n_max_cl),pt2min
       integer n, ibeam(2), iqcd(0:2)!, ilast(0:nexternal)
       integer idfl, idmap(-nexternal:nexternal)
       integer ipart(2,n_max_cl)
       double precision xnow(2)
-      integer jlast(2),jfirst(2),nwarning
-      logical qcdline(2),qcdrad(2)
+      integer jlast(2),jfirst(2),jcentral(2),nwarning
+      logical qcdline(2),partonline(2)
       logical failed,first
       data first/.true./
       data nwarning/0/
 
-      logical isqcd,isjet,isparton,isjetvx,cluster
+      logical isqcd,isjet,isparton,cluster
       double precision alphas
-      external isqcd, isjet, isparton, isjetvx, cluster, alphas
+      external isqcd, isjet, isparton, cluster, alphas
 
       setclscales=.true.
 
@@ -391,61 +425,50 @@
 c   Cluster the configuration
 c   
       
-      if (.not.cluster(p(0,1))) then
-c        if (xqcut.gt.0d0) then
-c          failed=.false.          
-cc          if(pt2ijcl(1).lt.xqcut**2) failed=.true.
-c          if(failed) then
-c            if (btest(mlevel,3)) then
-c              write(*,*)'q_min = ',pt2ijcl(1),' < ',xqcut**2
-c            endif
-c            setclscales=.false.
-c            return
-c          endif
-c        endif
-c      else
-        write(*,*)'setclscales: Error. Clustering failed.'
-        setclscales=.false.
-        return
-      endif
+c      if (.not.clustered) then
+         clustered = cluster(p(0,1))
+         if(.not.clustered) then
+            open(unit=26,file='../../../error',status='unknown',err=999)
+            write(26,*) 'Error: Clustering failed in cluster.f.'
+            write(*,*) 'Error: Clustering failed in cluster.f.'
+            stop
+ 999        write(*,*) 'error'
+            setclscales=.false.
+            clustered = .false.
+            return
+         endif
+c      endif
 
       if (btest(mlevel,1)) then
         write(*,*)'setclscales: identified tree {'
         do i=1,nexternal-2
-          write(*,*)'  ',i,': ',idacl(i,1),'(',ipdgcl(idacl(i,1),igraphs(1)),')',
-     $       '&',idacl(i,2),'(',ipdgcl(idacl(i,2),igraphs(1)),')',
-     $       ' -> ',imocl(i),', ptij = ',dsqrt(pt2ijcl(i)),
-     $       '(',ipdgcl(imocl(i),igraphs(1)),')'
+          write(*,*)'  ',i,': ',idacl(i,1),'(',ipdgcl(idacl(i,1),igraphs(1),iproc),')',
+     $       '&',idacl(i,2),'(',ipdgcl(idacl(i,2),igraphs(1),iproc),')',
+     $       ' -> ',imocl(i),'(',ipdgcl(imocl(i),igraphs(1),iproc),')',
+     $       ', ptij = ',dsqrt(pt2ijcl(i))
         enddo
+        write(*,*)'  process: ',iproc
         write(*,*)'  graphs (',igraphs(0),'):',(igraphs(i),i=1,igraphs(0))
         write(*,*)'}'
       endif
 
-cc
-cc   Set factorization scale as for the MLM case
-cc
-c      if(xqcut.gt.0) then
-cc     Using last clustering value
-c        if(pt2ijcl(nexternal-2).lt.max(4d0,xqcut**2))then
-c           setclscales=.false.
-c           return
-c        endif
-
 c     If last clustering is s-channel QCD (e.g. ttbar) use mt2last instead
 c     (i.e. geom. average of transverse mass of t and t~)
-        if(mt2last.gt.4d0 .and. nexternal.gt.3 .and. isqcd(ipdgcl(idacl(nexternal-3,1),igraphs(1)))
-     $      .and. isqcd(ipdgcl(idacl(nexternal-3,2),igraphs(1)))
-     $      .and. isqcd(ipdgcl(imocl(nexternal-3),igraphs(1))))then
-           mt2ij(nexternal-2)=mt2last
-           mt2ij(nexternal-3)=mt2last
-           if (btest(mlevel,3)) then
-              write(*,*)' setclscales: set last vertices to mtlast: ',sqrt(mt2last)
+        if(mt2last.gt.4d0 .and. nexternal.gt.3) then
+           if(isqcd(ipdgcl(idacl(nexternal-3,1),igraphs(1),iproc))
+     $      .and. isqcd(ipdgcl(idacl(nexternal-3,2),igraphs(1),iproc))
+     $      .and. isqcd(ipdgcl(imocl(nexternal-3),igraphs(1),iproc)))then
+              mt2ij(nexternal-2)=mt2last
+              mt2ij(nexternal-3)=mt2last
+              if (btest(mlevel,3)) then
+                 write(*,*)' setclscales: set last vertices to mtlast: ',sqrt(mt2last)
+              endif
            endif
         endif
 
 C   If we have fixed factorization scale, for ickkw>0 means central
 C   scale, i.e. last two scales (ren. scale for these vertices are
-C   anyway already set by "scale" above
+C   anyway already set by "scale" above)
       if(ickkw.gt.0) then
          if(fixed_fac_scale.and.first)then
             q2bck(1)=q2fact(1)
@@ -461,97 +484,142 @@
 
       ibeam(1)=ishft(1,0)
       ibeam(2)=ishft(1,1)
-      jlast(1)=0
-      jlast(2)=0
-      qcdline(1)=.false.
-      qcdline(2)=.false.
-      qcdrad(1)=.true.
-      qcdrad(2)=.true.
+      jlast(1)=1
+      jlast(2)=1
+      jcentral(1)=1
+      jcentral(2)=1
+      partonline(1)=isjet(ipdgcl(ibeam(1),igraphs(1),iproc))
+      partonline(2)=isjet(ipdgcl(ibeam(2),igraphs(1),iproc))
+      qcdline(1)=isqcd(ipdgcl(ibeam(1),igraphs(1),iproc))
+      qcdline(2)=isqcd(ipdgcl(ibeam(2),igraphs(1),iproc))
 
 c   Go through clusterings and set factorization scales for use in dsig
+      if (nexternal.eq.3) goto 10
       do n=1,nexternal-2
         do i=1,2
-          if (isqcd(ipdgcl(idacl(n,i),igraphs(1)))) then
             do j=1,2
-              if ((isparton(ipdgcl(idacl(n,i),igraphs(1))).and.idacl(n,i).eq.ibeam(j).or.
-     $           isparton(ipdgcl(imocl(n),igraphs(1))).and.imocl(n).eq.ibeam(j))
-     $           .and.qcdrad(j)) then
+              if (isqcd(ipdgcl(idacl(n,i),igraphs(1),iproc)).and.
+     $            idacl(n,i).eq.ibeam(j).and.qcdline(j)) then
 c             is emission - this is what we want
 c             Total pdf weight is f1(x1,pt2E)*fj(x1*z,Q)/fj(x1*z,pt2E)
 c             f1(x1,pt2E) is given by DSIG, just need to set scale.
                  ibeam(j)=imocl(n)
-                 if(ickkw.eq.0.or.jfirst(j).eq.0) jfirst(j)=n
-                 jlast(j)=n
-                 qcdline(j)=isqcd(ipdgcl(imocl(n),igraphs(1)))
-                 if(n.lt.nexternal-2)then
-                    qcdrad(j)=isqcd(ipdgcl(idacl(n,3-i),igraphs(1)))
-                 endif
-              else
-                 qcdline(j)=isqcd(ipdgcl(imocl(n),igraphs(1)))
+                 if(jfirst(j).eq.0)then
+                    if(isjet(ipdgcl(imocl(n),igraphs(1),iproc))) then
+                       jfirst(j)=n
+                    else
+                       jfirst(j)=-1
+                    endif
+                 endif
+                 if(partonline(j))then
+c                   Stop fact scale where parton line stops
+                    jlast(j)=n
+                    partonline(j)=isjet(ipdgcl(imocl(n),igraphs(1),iproc))
+                 endif
+c                Trace QCD line through event
+                 jcentral(j)=n
+                 qcdline(j)=isqcd(ipdgcl(imocl(n),igraphs(1),iproc))
               endif
             enddo
-          endif
         enddo
       enddo
+      
+ 10   if(jfirst(1).le.0) jfirst(1)=jlast(1)
+      if(jfirst(2).le.0) jfirst(2)=jlast(2)
 
       if (btest(mlevel,3))
      $     write(*,*) 'jfirst is ',jfirst(1),jfirst(2),
-     $     ' and jlast is ',jlast(1),jlast(2)
+     $     ' jlast is ',jlast(1),jlast(2),
+     $     ' and jcentral is ',jcentral(1),jcentral(2)
 
-c     Set central scale to mT2 and multiply with scalefact
-      if(jlast(1).gt.0.and.mt2ij(jlast(1)).gt.0d0)
-     $     pt2ijcl(jlast(1))=mt2ij(jlast(1))
-      if(jlast(2).gt.0.and.mt2ij(jlast(2)).gt.0d0)
-     $     pt2ijcl(jlast(2))=mt2ij(jlast(2))
+c     Set central scale to mT2
+      if(jcentral(1).gt.0.and.mt2ij(jcentral(1)).gt.0d0)
+     $     pt2ijcl(jcentral(1))=mt2ij(jcentral(1))
+      if(jcentral(2).gt.0.and.mt2ij(jcentral(2)).gt.0d0)
+     $     pt2ijcl(jcentral(2))=mt2ij(jcentral(2))
       if(btest(mlevel,4))
-     $     print *,'pt2ijcl is: ',jlast(1), sqrt(pt2ijcl(jlast(1))),
-     $     jlast(2), sqrt(pt2ijcl(jlast(2)))
-      if(qcdline(1).and.qcdline(2).and.jlast(1).ne.jlast(2)) then
-c     If not WBF or similar, set uniform scale to be maximum
-         pt2ijcl(jlast(1))=max(pt2ijcl(jlast(1)),pt2ijcl(jlast(2)))
-         pt2ijcl(jlast(2))=pt2ijcl(jlast(1))
-      endif
-      if(jlast(1).gt.0) pt2ijcl(jlast(1))=scalefact**2*pt2ijcl(jlast(1))
-      if(jlast(2).gt.0) pt2ijcl(jlast(2))=scalefact**2*pt2ijcl(jlast(2))
-      if(lpp(1).eq.0.and.lpp(2).eq.0)then
-         pt2ijcl(nexternal-2)=scalefact**2*pt2ijcl(nexternal-2)
-         pt2ijcl(nexternal-3)=pt2ijcl(nexternal-2)
-      endif
+     $     write(*,*)'pt2ijcl is: ',jlast(1), sqrt(pt2ijcl(jlast(1))),
+     $     jlast(2), sqrt(pt2ijcl(jlast(2))),
+     $     jcentral(1), sqrt(pt2ijcl(jcentral(1))),
+     $     jcentral(2), sqrt(pt2ijcl(jcentral(2)))
 
 c     Check xqcut for vertices with jet daughters only
+      ibeam(1)=ishft(1,0)
+      ibeam(2)=ishft(1,1)
       if(xqcut.gt.0) then
          do n=1,nexternal-2
-            if (n.lt.nexternal-2.and.n.ne.jlast(1).and.n.ne.jlast(2).and.
-     $           (isjet(ipdgcl(idacl(n,1),igraphs(1))).or.isjet(ipdgcl(idacl(n,2),igraphs(1)))).and.
+            do j=1,2
+               if (btest(mlevel,4))
+     $              write(*,*) 'ibeam(1,2), daughter, mother:',
+     $              ibeam(1),ibeam(2),idacl(n,j),imocl(n)
+               if (n.lt.nexternal-2) then
+                  if(n.ne.jlast(1).and.n.ne.jlast(2).and.
+     $              isjet(ipdgcl(idacl(n,j),igraphs(1),iproc)).and.
+     $              (idacl(n,3-j).eq.ibeam(1).or.
+     $              idacl(n,3-j).eq.ibeam(2)).and.
+     $              sqrt(pt2ijcl(n)).lt.xqcut)then
+c                   ISR
+                     if (btest(mlevel,3))
+     $                    write(*,*) 'Failed xqcut: ',n, ipdgcl(idacl(n,1),igraphs(1),iproc),
+     $                    ipdgcl(idacl(n,2),igraphs(1),iproc), xqcut
+                     setclscales=.false.
+                     clustered = .false.
+                     return
+                  endif
+               endif
+               if (idacl(n,j).eq.ibeam(1).and.imocl(n).ne.ibeam(2))
+     $              ibeam(1)=imocl(n)
+               if (idacl(n,j).eq.ibeam(2).and.imocl(n).ne.ibeam(1))
+     $              ibeam(2)=imocl(n)
+            enddo
+            if (n.lt.nexternal-2) then
+               if(n.ne.jlast(1).and.n.ne.jlast(2).and.
+     $           isjet(ipdgcl(idacl(n,1),igraphs(1),iproc)).and.
+     $           isjet(ipdgcl(idacl(n,2),igraphs(1),iproc)).and.
+     $           idacl(n,1).ne.ibeam(1).and.idacl(n,1).ne.ibeam(2).and.
+     $           idacl(n,2).ne.ibeam(1).and.idacl(n,2).ne.ibeam(2).and.
      $           sqrt(pt2ijcl(n)).lt.xqcut)then
-               setclscales=.false.
-               return
+c           FSR
+                  if (btest(mlevel,3))
+     $                 write(*,*) 'Failed xqcut: ',n, ipdgcl(idacl(n,1),igraphs(1),iproc),
+     $                 ipdgcl(idacl(n,2),igraphs(1),iproc), xqcut
+                  setclscales=.false.
+                  clustered = .false.
+                  return
+               endif
             endif
          enddo
       endif
 
 c     JA: Check xmtc cut for central process
-      if(pt2ijcl(jlast(1)).lt.xmtc**2.or.pt2ijcl(jlast(2)).lt.xmtc**2)then
+      if(pt2ijcl(jcentral(1)).lt.xmtc**2.or.pt2ijcl(jcentral(2)).lt.xmtc**2)then
          setclscales=.false.
+         clustered = .false.
          if(btest(mlevel,3)) write(*,*)'Failed xmtc cut ',
-     $        sqrt(pt2ijcl(jlast(1))),sqrt(pt2ijcl(jlast(1))),' < ',xmtc
+     $        sqrt(pt2ijcl(jcentral(1))),sqrt(pt2ijcl(jcentral(1))),
+     $        ' < ',xmtc
          return
       endif
       
       if(ickkw.eq.0.and.(fixed_fac_scale.or.q2fact(1).gt.0).and.
      $     (fixed_ren_scale.or.scale.gt.0)) return
 
-c     Set renormalization scale to geom. aver. of factorization scales
+c     Ensure that last scales are at least as big as first scales
+      pt2ijcl(jlast(1))=max(pt2ijcl(jlast(1)),pt2ijcl(jfirst(1)))
+      pt2ijcl(jlast(2))=max(pt2ijcl(jlast(2)),pt2ijcl(jfirst(2)))
+
+c     Set renormalization scale to geom. aver. of central scales
       if(scale.eq.0d0) then
-         if(jlast(1).gt.0.and.jlast(2).gt.0) then
-            scale=(pt2ijcl(jlast(1))*pt2ijcl(jlast(2)))**0.25d0
-         elseif(jlast(1).gt.0) then
-            scale=sqrt(pt2ijcl(jlast(1)))
-         elseif(jlast(2).gt.0) then
-            scale=sqrt(pt2ijcl(jlast(2)))
+         if(jcentral(1).gt.0.and.jcentral(2).gt.0) then
+            scale=(pt2ijcl(jcentral(1))*pt2ijcl(jcentral(2)))**0.25d0
+         elseif(jcentral(1).gt.0) then
+            scale=sqrt(pt2ijcl(jcentral(1)))
+         elseif(jcentral(2).gt.0) then
+            scale=sqrt(pt2ijcl(jcentral(2)))
          elseif(lpp(1).eq.0.and.lpp(2).eq.0) then
             scale=sqrt(pt2ijcl(nexternal-2))
          endif
+         scale = scalefact*scale
          if(scale.gt.0)
      $        G = SQRT(4d0*PI*ALPHAS(scale))
       endif
@@ -560,10 +628,35 @@
 
       if(ickkw.gt.0.and.q2fact(1).gt.0) then
 c     Use the fixed or previously set scale for central scale
-         if(jlast(1).gt.0) pt2ijcl(jlast(1))=q2fact(1)
-         if(jlast(2).gt.0) pt2ijcl(jlast(2))=q2fact(2)
-      endif
-      
+         if(jcentral(1).gt.0) pt2ijcl(jcentral(1))=q2fact(1)
+         if(jcentral(2).gt.0.and.jcentral(2).ne.jcentral(1))
+     $        pt2ijcl(jcentral(2))=q2fact(2)
+      endif
+
+      if(nexternal.eq.3.and.nincoming.eq.2.and.q2fact(1).eq.0) then
+         q2fact(1)=pt2ijcl(nexternal-2)
+         q2fact(2)=pt2ijcl(nexternal-2)
+      endif
+
+      if(q2fact(1).eq.0d0) then
+c     Use the geom. average of central scale and first non-radiation vertex
+         if(jlast(1).gt.0) q2fact(1)=sqrt(pt2ijcl(jlast(1))*pt2ijcl(jcentral(1)))
+         if(jlast(2).gt.0) q2fact(2)=sqrt(pt2ijcl(jlast(2))*pt2ijcl(jcentral(2)))
+         if(jcentral(1).eq.jcentral(2))then
+c     We have a qcd line going through the whole event, use single scale
+            q2fact(1)=max(q2fact(1),q2fact(2))
+            q2fact(2)=q2fact(1)
+         endif
+      endif
+      if(.not. fixed_fac_scale) then
+         q2fact(1)=scalefact**2*q2fact(1)
+         q2fact(2)=scalefact**2*q2fact(2)
+         q2bck(1)=q2fact(1)
+         q2bck(2)=q2fact(2)
+         if (btest(mlevel,3))
+     $      write(*,*) 'Set central fact scales to ',sqrt(q2bck(1)),sqrt(q2bck(2))
+      endif
+         
       if(lpp(1).eq.0.and.lpp(2).eq.0)then
          if(q2fact(1).gt.0)then
             pt2ijcl(nexternal-2)=q2fact(1)
@@ -574,18 +667,10 @@
          endif
       elseif(ickkw.eq.2.or.pdfwgt)then
 c     Use the minimum scale found for fact scale in ME
-         if(jlast(1).gt.0) q2fact(1)=min(pt2ijcl(jfirst(1)),pt2ijcl(jlast(1)))
-         if(jlast(2).gt.0) q2fact(2)=min(pt2ijcl(jfirst(2)),pt2ijcl(jlast(2)))
-      else if(q2fact(1).eq.0d0) then
-         if(jlast(1).gt.0) q2fact(1)=pt2ijcl(jlast(1))
-         if(jlast(2).gt.0) q2fact(2)=pt2ijcl(jlast(2))
-      endif
-
-      if(nexternal.eq.3.and.nincoming.eq.2) then
-         if(q2fact(1).eq.0)
-     $        q2fact(1)=pt2ijcl(nexternal-2)
-         if(q2fact(2).eq.0)
-     $        q2fact(2)=pt2ijcl(nexternal-2)
+         if(jlast(1).gt.0.and.jfirst(1).lt.jlast(1))
+     $        q2fact(1)=min(pt2ijcl(jfirst(1)),q2fact(1))
+         if(jlast(2).gt.0.and.jfirst(2).lt.jlast(2))
+     $        q2fact(2)=min(pt2ijcl(jfirst(2)),q2fact(2))
       endif
 
 c     Check that factorization scale is >= 2 GeV
@@ -601,6 +686,7 @@
              write(*,*) 'No more warnings written out this run.'
           endif
          setclscales=.false.
+         clustered = .false.
          return
       endif
 
@@ -620,19 +706,26 @@
       include 'genps.inc'
       include 'maxconfigs.inc'
       include 'nexternal.inc'
+      include 'maxamps.inc'
       include 'cluster.inc'
       include 'run.inc'
       include 'coupl.inc'
-      include 'maxamps.inc'
 C   
 C   ARGUMENTS 
 C   
       DOUBLE PRECISION P(0:3,NEXTERNAL)
 
 C   global variables
-      integer              IPROC 
+C   global variables
+C     Present process number
+      INTEGER IMIRROR,IPROC
+      COMMON/TO_MIRROR/IMIRROR, IPROC
+      integer              ISPROC 
       DOUBLE PRECISION PD(0:MAXPROC)
-      COMMON /SubProc/ PD, IPROC
+      COMMON /SubProc/ PD, ISPROC
+c     q2bck holds the central q2fact scales
+      real*8 q2bck(2)
+      common /to_q2bck/q2bck
 
 C   local variables
       integer i, j, idi, idj
@@ -653,27 +746,28 @@
 c     ipart gives external particle number chain
       integer ipart(2,n_max_cl)
       double precision xnow(2)
-      double precision xtarget, tmp
+      double precision xtarget,tmp,pdfj1,pdfj2,q2now,etot
       integer iseed,np
       data iseed/0/
       logical isvx
 
-      logical isqcd,isjet,isparton,isjetvx
+      logical isqcd,isjet,isparton,isjetvx,ispartonvx
       double precision alphas,getissud,pdg2pdf, sudwgt
       real xran1
-      external isqcd,isjet,isparton
+      external isqcd,isjet,isparton,ispartonvx
       external alphas, isjetvx, getissud, pdg2pdf, xran1,  sudwgt
 
       rewgt=1.0d0
 
-      if(ickkw.le.0) return
-
-      if(.not.clustered)then
+      if((ickkw.gt.0.or..not.fixed_fac_scale.or..not.fixed_ren_scale)
+     $     .and..not.clustered)then
         write(*,*)'Error: No clustering done when calling rewgt!'
         stop
       endif
       clustered=.false.
 
+      if(ickkw.le.0) return
+
 c   Set mimimum kt scale, depending on highest mult or not
       if(hmult.or.ickkw.eq.1)then
         pt2min=0
@@ -683,9 +777,12 @@
       if (btest(mlevel,3))
      $     write(*,*) 'pt2min set to ',pt2min
 
+c   Set etot, used for non-radiating partons
+      etot=2d0*sqrt(ebeam(1)*ebeam(2))
+
 c   Since we use pdf reweighting, need to know particle identities
       iprocset=1
-      np = iproc
+      np = isproc
       xtarget=xran1(iseed)*pd(np)
       iprocset = 1
       do while (pd(iprocset) .lt. xtarget .and. iprocset .lt. np)
@@ -699,15 +796,24 @@
 c   external particle clustering scales)
       do i=1,nexternal
 c        ilast(i)=ishft(1,i)
-         if(pt2min.gt.0)then
-            pt2prev(ishft(1,i-1))=max(pt2min,p(0,i)**2-p(1,i)**2-p(2,i)**2-p(3,i)**2)
-         else
-            pt2prev(ishft(1,i-1))=0d0
+         pt2prev(ishft(1,i-1))=0d0
+         if (ickkw.eq.2) then
+            if(pt2min.gt.0)then
+               pt2prev(ishft(1,i-1))=
+     $              max(pt2min,p(0,i)**2-p(1,i)**2-p(2,i)**2-p(3,i)**2)
+            endif
+            pt2pdf(ishft(1,i-1))=pt2prev(ishft(1,i-1))
+         else if(pdfwgt) then
+            pt2pdf(ishft(1,i-1))=0d0
          endif
-         pt2pdf(ishft(1,i-1))=pt2prev(ishft(1,i-1))
          ptclus(i)=sqrt(pt2prev(ishft(1,i-1)))
+         if (btest(mlevel,3))
+     $        write(*,*) 'Set ptclus for ',i,' to ', ptclus(i)
          ipart(1,ishft(1,i-1))=i
          ipart(2,ishft(1,i-1))=0
+         if (btest(mlevel,4))
+     $        write(*,*) 'Set ipart for ',ishft(1,i-1),' to ',
+     $        ipart(1,ishft(1,i-1)),ipart(2,ishft(1,i-1))
       enddo
 c      ilast(0)=nexternal
       ibeam(1)=ishft(1,0)
@@ -715,8 +821,10 @@
       if (btest(mlevel,1)) then
         write(*,*)'rewgt: identified tree {'
         do i=1,nexternal-2
-          write(*,*)'  ',i,': ',idacl(i,1),'&',idacl(i,2),
-     &       ' -> ',imocl(i),', ptij = ',dsqrt(pt2ijcl(i)) 
+          write(*,*)'  ',i,': ',idacl(i,1),'(',ipdgcl(idacl(i,1),igraphs(1),iproc),')',
+     $       '&',idacl(i,2),'(',ipdgcl(idacl(i,2),igraphs(1),iproc),')',
+     $       ' -> ',imocl(i),'(',ipdgcl(imocl(i),igraphs(1),iproc),')',
+     $       ', ptij = ',dsqrt(pt2ijcl(i))
         enddo
         write(*,*)'  graphs (',igraphs(0),'):',(igraphs(i),i=1,igraphs(0))
         write(*,*)'}'
@@ -730,7 +838,7 @@
       endif
 
 c     Prepare for resetting q2fact based on PDF reweighting
-      if(ickkw.eq.2.or.pdfwgt)then
+      if(ickkw.eq.2)then
          q2fact(1)=0d0
          q2fact(2)=0d0
       endif
@@ -741,102 +849,181 @@
 
 c   Perform alpha_s reweighting based on type of vertex
       do n=1,nexternal-2
+c       scale for alpha_s reweighting
+        q2now=pt2ijcl(n)
+        if(n.eq.nexternal-2) then
+           q2now = scale**2
+        endif
         if (btest(mlevel,3)) then
-          write(*,*)'  ',n,': ',idacl(n,1),'(',ipdgcl(idacl(n,1),igraphs(1)),
-     &       ')&',idacl(n,2),'(',ipdgcl(idacl(n,2),igraphs(1)),') -> ',
-     &       imocl(n),'(',ipdgcl(imocl(n),igraphs(1)),'), ptij = ',
-     &       dsqrt(pt2ijcl(n)) 
+          write(*,*)'  ',n,': ',idacl(n,1),'(',ipdgcl(idacl(n,1),igraphs(1),iproc),
+     &       ')&',idacl(n,2),'(',ipdgcl(idacl(n,2),igraphs(1),iproc),
+     &       ') -> ',imocl(n),'(',ipdgcl(imocl(n),igraphs(1),iproc),
+     &       '), ptij = ',dsqrt(q2now) 
         endif
-c     perform alpha_s reweighting only for vertices where a jet is produced
+c     perform alpha_s reweighting only for vertices where a parton is produced
 c     and not for the last clustering (use non-fixed ren. scale for these)
-        if (n.lt.nexternal-2.and.
-     $     isjetvx(imocl(n),idacl(n,1),idacl(n,2),ipdgcl(1,igraphs(1)),ipart)) then
+        if (n.lt.nexternal-2)then
+           if(ispartonvx(imocl(n),idacl(n,1),idacl(n,2),
+     $       ipdgcl(1,igraphs(1),iproc),ipart,.false.)) then
 c       alpha_s weight
-          rewgt=rewgt*alphas(alpsfact*sqrt(pt2ijcl(n)))/asref
-          if (btest(mlevel,3)) then
-             write(*,*)' reweight vertex: ',ipdgcl(imocl(n),igraphs(1)),ipdgcl(idacl(n,1),igraphs(1)),ipdgcl(idacl(n,2),igraphs(1))
-            write(*,*)'       as: ',alphas(alpsfact*dsqrt(pt2ijcl(n))),
-     &         '/',asref,' -> ',alphas(alpsfact*dsqrt(pt2ijcl(n)))/asref
-            write(*,*)' and G=',SQRT(4d0*PI*ALPHAS(scale))
-          endif
+              rewgt=rewgt*alphas(alpsfact*sqrt(q2now))/asref
+              if (btest(mlevel,3)) then
+                 write(*,*)' reweight vertex: ',ipdgcl(imocl(n),igraphs(1),iproc),
+     $                ipdgcl(idacl(n,1),igraphs(1),iproc),ipdgcl(idacl(n,2),igraphs(1),iproc)
+                 write(*,*)'       as: ',alphas(alpsfact*dsqrt(q2now)),
+     &                '/',asref,' -> ',alphas(alpsfact*dsqrt(q2now))/asref
+                 write(*,*)' and G=',SQRT(4d0*PI*ALPHAS(scale))
+              endif
+           endif
         endif
 c   Update starting values for FS parton showering
         do i=1,2
           do j=1,2
-            if(ipart(j,idacl(n,i)).gt.0)then
-              ptclus(ipart(j,idacl(n,i)))=dsqrt(pt2ijcl(n))
+            if(ipart(j,idacl(n,i)).gt.0.and.ipart(j,idacl(n,i)).gt.2)then
+              ptclus(ipart(j,idacl(n,i)))=
+     $              max(ptclus(ipart(j,idacl(n,i))),dsqrt(pt2ijcl(n)))
+              if(ickkw.ne.2.and.
+     $             (.not.isqcd(ipdgcl(imocl(n),igraphs(1),iproc)).or.
+     $             ipart(1,idacl(n,3-i)).le.2.and.
+     $             .not.isqcd(ipdgcl(idacl(n,3-i),igraphs(1),iproc)).or.
+     $             isbw(imocl(n))))then
+c             For particles originating in non-qcd t-channel vertices or decay vertices,
+c             set origination scale to machine energy since we don't want these
+c             to be included in matching.
+                 ptclus(ipart(j,idacl(n,i)))=etot
+              endif
+              if (btest(mlevel,3))
+     $             write(*,*) 'Set ptclus for ',ipart(j,idacl(n,i)),
+     $             ' to ', ptclus(ipart(j,idacl(n,i))),
+     $             ipdgcl(imocl(n),igraphs(1),iproc),
+     $             isqcd(ipdgcl(imocl(n),igraphs(1),iproc)),isbw(imocl(n))
             endif
           enddo
         enddo
+c     Special case for last 1,2->i vertex
+        if(n.eq.nexternal-2)then
+           ptclus(ipart(1,imocl(n)))=
+     $              max(ptclus(ipart(1,imocl(n))),dsqrt(pt2ijcl(n)))
+              if(ickkw.ne.2.and.
+     $             (.not.isqcd(ipdgcl(idacl(n,1),igraphs(1),iproc)).or.
+     $             .not.isqcd(ipdgcl(idacl(n,2),igraphs(1),iproc))))then
+c             For particles originating in non-qcd vertices or decay vertices,
+c             set origination scale to machine energy since we don't want these
+c             to be included in matching.
+                 ptclus(ipart(1,imocl(n)))=etot
+              endif
+              if (btest(mlevel,3))
+     $             write(*,*) 'Set ptclus for ',ipart(1,imocl(n)),
+     $             ' to ', ptclus(ipart(1,imocl(n))),
+     $             ipdgcl(idacl(n,1),igraphs(1),iproc),
+     $             ipdgcl(idacl(n,2),igraphs(1),iproc)
+        endif
 c   Update particle tree map
-        call ipartupdate(p,imocl(n),idacl(n,1),idacl(n,2),ipdgcl(1,igraphs(1)),ipart)
+        call ipartupdate(p,imocl(n),idacl(n,1),idacl(n,2),
+     $       ipdgcl(1,igraphs(1),iproc),ipart)
         if(ickkw.eq.2.or.pdfwgt) then
 c       Perform PDF and, if ickkw=2, Sudakov reweighting
           isvx=.false.
           do i=1,2
 c         write(*,*)'weight ',idacl(n,i),', ptij=',pt2prev(idacl(n,i))
-            if (isqcd(ipdgcl(idacl(n,i),igraphs(1)))) then
-               if(pt2min.eq.0d0) then
+            if (isqcd(ipdgcl(idacl(n,i),igraphs(1),iproc))) then
+               if(ickkw.eq.2.and.pt2min.eq.0d0) then
                   pt2min=pt2ijcl(n)
                   if (btest(mlevel,3))
      $                 write(*,*) 'pt2min set to ',pt2min
                endif
-               if(pt2prev(idacl(n,i)).eq.0d0) pt2prev(idacl(n,i))=
+               if(ickkw.eq.2.and.pt2prev(idacl(n,i)).eq.0d0)
+     $              pt2prev(idacl(n,i))=
      $              max(pt2min,p(0,i)**2-p(1,i)**2-p(2,i)**2-p(3,i)**2)
                do j=1,2
-                  if (isparton(ipdgcl(idacl(n,i),igraphs(1))).and
+                  if (isparton(ipdgcl(idacl(n,i),igraphs(1),iproc)).and
      $                 .idacl(n,i).eq.ibeam(j)) then
 c               is sudakov weight - calculate only once for each parton
 c               line where parton line ends with change of parton id or
 c               non-radiation vertex
                      isvx=.true.
                      ibeam(j)=imocl(n)
-                     if(pt2pdf(idacl(n,i)).eq.0d0) pt2pdf(idacl(n,i))=pt2prev(idacl(n,i))
-                     if(ickkw.eq.2.and.(ipdgcl(idacl(n,i),igraphs(1)).ne.
-     $                    ipdgcl(imocl(n),igraphs(1)).or.
+c                    Perform Sudakov reweighting if ickkw=2
+                     if(ickkw.eq.2.and.(ipdgcl(idacl(n,i),igraphs(1),iproc).ne.
+     $                    ipdgcl(imocl(n),igraphs(1),iproc).or.
      $                    .not.isjetvx(imocl(n),idacl(n,1),idacl(n,2),
-     $                    ipdgcl(1,igraphs(1)),ipart)).and.
-     $                    pt2prev(idacl(n,i)).lt.pt2ijcl(n).and.zcl(n).gt.1d-20)then
+     $                    ipdgcl(1,igraphs(1),iproc),ipart,n.eq.nexternal-2)).and.
+     $                    pt2prev(idacl(n,i)).lt.pt2ijcl(n))then
                         tmp=min(1d0,max(getissud(ibeam(j),ipdgcl(idacl(n,i),
-     $                       igraphs(1)),xnow(j),xnow(3-j),pt2ijcl(n)),1d-20)/
+     $                       igraphs(1),iproc),xnow(j),xnow(3-j),pt2ijcl(n)),1d-20)/
      $                       max(getissud(ibeam(j),ipdgcl(idacl(n,i),
-     $                       igraphs(1)),xnow(j),xnow(3-j),pt2prev(idacl(n,i))),1d-20))
+     $                       igraphs(1),iproc),xnow(j),xnow(3-j),pt2prev(idacl(n,i))),1d-20))
                         rewgt=rewgt*tmp
                         pt2prev(imocl(n))=pt2ijcl(n)
                         if (btest(mlevel,3)) then
-                           write(*,*)' reweight line: ',ipdgcl(idacl(n,i),igraphs(1)), idacl(n,i)
+                           write(*,*)' reweight line: ',ipdgcl(idacl(n,i),igraphs(1),iproc), idacl(n,i)
                            write(*,*)'     pt2prev, pt2new, x1, x2: ',pt2prev(idacl(n,i)),pt2ijcl(n),xnow(j),xnow(3-j)
                            write(*,*)'           Sud: ',tmp
                            write(*,*)'        -> rewgt: ',rewgt
                         endif
-                     else
+                     else if(ickkw.eq.2) then
                         pt2prev(imocl(n))=pt2prev(idacl(n,i))
                      endif
+c                 End Sudakov reweighting when we reach a non-radiation vertex
+                     if(ickkw.eq.2.and..not.
+     $                    ispartonvx(imocl(n),idacl(n,1),idacl(n,2),
+     $                    ipdgcl(1,igraphs(1),iproc),ipart,n.eq.nexternal-2)) then
+                        pt2prev(imocl(n))=1d30
+                        if (btest(mlevel,3)) then
+                          write(*,*)' rewgt: ending reweighting for vx ',
+     $                          idacl(n,1),idacl(n,2),imocl(n),
+     $                          ' with ids ',ipdgcl(idacl(n,1),igraphs(1),iproc),
+     $                          ipdgcl(idacl(n,2),igraphs(1),iproc),ipdgcl(imocl(n),igraphs(1),iproc)
+                        endif
+                     endif
+c               PDF reweighting
 c               Total pdf weight is f1(x1,pt2E)*fj(x1*z,Q)/fj(x1*z,pt2E)
 c               f1(x1,pt2E) is given by DSIG, already set scale for that
-                     xnow(j)=xnow(j)*zcl(n)
+                     if (zcl(n).gt.0d0.and.zcl(n).lt.1d0) then
+                        xnow(j)=xnow(j)*zcl(n)
+                     endif
+c                    PDF scale
+                     q2now=min(pt2ijcl(n), q2bck(j))
+c                    Set PDF scale to central factorization scale
+c                    if non-radiating vertex or last 2->2
+                     if(.not.isjetvx(imocl(n),idacl(n,1),idacl(n,2),
+     $                    ipdgcl(1,igraphs(1),iproc),ipart,n.eq.nexternal-2)) then
+                        q2now=q2bck(j)
+                     endif
+                     if (btest(mlevel,3))
+     $                    write(*,*)' set q2now for pdf to ',sqrt(q2now)
                      if(q2fact(j).eq.0d0.and.ickkw.eq.2)then
                         q2fact(j)=pt2min ! Starting scale for PS
-                        if (btest(mlevel,3)) then
-                           write(*,*)' reweight: set fact scale ',j,' for PS scale to: ',q2fact(j)
+                        pt2pdf(imocl(n))=q2now
+                        if (btest(mlevel,3))
+     $                       write(*,*)' set fact scale ',j,
+     $                          ' for PS scale to: ',sqrt(q2fact(j))
+                     else if(pt2pdf(idacl(n,i)).eq.0d0)then
+                        pt2pdf(imocl(n))=q2now
+                        if (btest(mlevel,3))
+     $                       write(*,*)' set pt2pdf for ',imocl(n),
+     $                          ' to: ',sqrt(pt2pdf(imocl(n)))
+                     else if(pt2pdf(idacl(n,i)).lt.q2now
+     $                       .and.isjet(ipdgcl(idacl(n,i),igraphs(1),iproc))) then
+                        pdfj1=pdg2pdf(abs(ibeam(j)),ipdgcl(idacl(n,i),
+     $                       igraphs(1),iproc)*sign(1,ibeam(j)),xnow(j),sqrt(q2now))
+                        pdfj2=pdg2pdf(abs(ibeam(j)),ipdgcl(idacl(n,i),
+     $                       igraphs(1),iproc)*sign(1,ibeam(j)),xnow(j),
+     $                       sqrt(pt2pdf(idacl(n,i))))
+                        if(pdfj2.lt.1d-10)then
+c                          Scale too low for heavy quark
+                           rewgt=0d0
+                           if (btest(mlevel,3))
+     $                        write(*,*) 'Too low scale for quark pdf: ',
+     $                        sqrt(pt2pdf(idacl(n,i))),pdfj2,pdfj1
+                           return
                         endif
-                     else if(q2fact(j).eq.0d0)then
-                        q2fact(j)=pt2ijcl(n)
-                     else if(pt2pdf(idacl(n,i)).lt.pt2ijcl(n).and.zcl(n).gt.1d-20) then
-                        if(ickkw.eq.1) q2fact(j)=pt2ijcl(n)
-                        rewgt=rewgt*max(pdg2pdf(abs(ibeam(j)),ipdgcl(idacl(n,i),
-     $                       igraphs(1))*sign(1,ibeam(j)),xnow(j),sqrt(pt2ijcl(n))),1d-20)/
-     $                       max(pdg2pdf(abs(ibeam(j)),ipdgcl(idacl(n,i),
-     $                       igraphs(1))*sign(1,ibeam(j)),xnow(j),
-     $                       sqrt(pt2pdf(idacl(n,i)))),1d-20)
+                        rewgt=rewgt*pdfj1/pdfj2
                         if (btest(mlevel,3)) then
-                           write(*,*)' reweight ',ipdgcl(idacl(n,i),igraphs(1)),' by pdfs: '
-                           write(*,*)'     x, pt2prev, ptnew: ',xnow(j),pt2pdf(idacl(n,i)),pt2ijcl(n)
-                           write(*,*)'           PDF: ',
-     $                          pdg2pdf(abs(ibeam(j)),ipdgcl(idacl(n,i),
-     $                          igraphs(1))*sign(1,ibeam(j)),xnow(j),sqrt(pt2ijcl(n))),' / ',
-     $                          pdg2pdf(abs(ibeam(j)),ipdgcl(idacl(n,i),
-     $                          igraphs(1))*sign(1,ibeam(j)),xnow(j),sqrt(pt2pdf(idacl(n,i))))
+                           write(*,*)' reweight ',n,i,ipdgcl(idacl(n,i),igraphs(1),iproc),' by pdfs: '
+                           write(*,*)'     x, ptprev, ptnew: ',xnow(j),
+     $                          sqrt(pt2pdf(idacl(n,i))),sqrt(q2now)
+                           write(*,*)'           PDF: ',pdfj1,' / ',pdfj2
                            write(*,*)'        -> rewgt: ',rewgt
 c                           write(*,*)'  (compare for glue: ',
 c     $                          pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2pdf(idacl(n,i)))),' / ',
@@ -846,35 +1033,26 @@
 c                           write(*,*)'       -> ',pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2pdf(idacl(n,i))))/
 c     $                          pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2ijcl(n)))*rewgt,' )'
                         endif
-                     endif
-c               End both Sudakov and pdf reweighting when we reach a
-c               non-radiation vertex
-                     if(isjetvx(imocl(n),idacl(n,1),idacl(n,2),ipdgcl(1,igraphs(1)),ipart)) then
-                        pt2pdf(imocl(n))=pt2ijcl(n)
-                     else
-                        pt2pdf(imocl(n))=1d30
-                        pt2prev(imocl(n))=1d30
-                        if (btest(mlevel,3)) then
-                           write(*,*)' rewgt: for vertex ',idacl(n,1),idacl(n,2),imocl(n),
-     $                          ' with ids ',ipdgcl(idacl(n,1),igraphs(1)),
-     $                          ipdgcl(idacl(n,2),igraphs(1)),ipdgcl(imocl(n),igraphs(1))
-                           write(*,*)'    set pt2prev, pt2pdf: ',pt2prev(imocl(n)),pt2pdf(imocl(n))
-                        endif
+c                       Set scale for mother as this scale
+                        pt2pdf(imocl(n))=q2now                           
+                     else if(pt2pdf(idacl(n,i)).ge.q2now) then
+c                    If no reweighting, just copy daughter scale for mother
+                        pt2pdf(imocl(n))=pt2pdf(idacl(n,i))
                      endif
                      goto 10
                   endif
                enddo
 c           fs sudakov weight
                if(ickkw.eq.2.and.pt2prev(idacl(n,i)).lt.pt2ijcl(n).and.
-     $              (isvx.or.ipdgcl(idacl(n,i),igraphs(1)).ne.ipdgcl(imocl(n),igraphs(1)).or.
-     $              (ipdgcl(idacl(n,i),igraphs(1)).ne.
-     $              ipdgcl(idacl(n,3-i),igraphs(1)).and.
+     $              (isvx.or.ipdgcl(idacl(n,i),igraphs(1),iproc).ne.ipdgcl(imocl(n),igraphs(1),iproc).or.
+     $              (ipdgcl(idacl(n,i),igraphs(1),iproc).ne.
+     $              ipdgcl(idacl(n,3-i),igraphs(1),iproc).and.
      $              pt2prev(idacl(n,i)).gt.pt2prev(idacl(n,3-i))))) then
                   tmp=sudwgt(sqrt(pt2min),sqrt(pt2prev(idacl(n,i))),
-     &                 dsqrt(pt2ijcl(n)),ipdgcl(idacl(n,i),igraphs(1)),1)
+     &                 dsqrt(pt2ijcl(n)),ipdgcl(idacl(n,i),igraphs(1),iproc),1)
                   rewgt=rewgt*tmp
                   if (btest(mlevel,3)) then
-                     write(*,*)' reweight fs line: ',ipdgcl(idacl(n,i),igraphs(1)), idacl(n,i)
+                     write(*,*)' reweight fs line: ',ipdgcl(idacl(n,i),igraphs(1),iproc), idacl(n,i)
                      write(*,*)'     pt2prev, pt2new: ',pt2prev(idacl(n,i)),pt2ijcl(n)
                      write(*,*)'           Sud: ',tmp
                      write(*,*)'        -> rewgt: ',rewgt
@@ -886,13 +1064,13 @@
             endif
  10         continue
           enddo
-          if (ickkw.eq.2.and.n.eq.nexternal-2.and.isqcd(ipdgcl(imocl(n),igraphs(1))).and.
+          if (ickkw.eq.2.and.n.eq.nexternal-2.and.isqcd(ipdgcl(imocl(n),igraphs(1),iproc)).and.
      $         pt2prev(imocl(n)).lt.pt2ijcl(n)) then
              tmp=sudwgt(sqrt(pt2min),sqrt(pt2prev(imocl(n))),
-     &            dsqrt(pt2ijcl(n)),ipdgcl(imocl(n),igraphs(1)),1)
+     &            dsqrt(pt2ijcl(n)),ipdgcl(imocl(n),igraphs(1),iproc),1)
              rewgt=rewgt*tmp
              if (btest(mlevel,3)) then
-                write(*,*)' reweight last fs line: ',ipdgcl(imocl(n),igraphs(1)), imocl(n)
+                write(*,*)' reweight last fs line: ',ipdgcl(imocl(n),igraphs(1),iproc), imocl(n)
                 write(*,*)'     pt2prev, pt2new: ',pt2prev(imocl(n)),pt2ijcl(n)
                 write(*,*)'           Sud: ',tmp
                 write(*,*)'        -> rewgt: ',rewgt
@@ -901,9 +1079,15 @@
         endif
       enddo
 
-      if((ickkw.eq.2.or.pdfwgt).and.lpp(1).eq.0.and.lpp(2).eq.0)then
+      if(ickkw.eq.2.and.lpp(1).eq.0.and.lpp(2).eq.0)then
          q2fact(1)=pt2min
          q2fact(2)=q2fact(1)
+      else if (ickkw.eq.1.and.pdfwgt) then
+         q2fact(1)=q2bck(1)
+         q2fact(2)=q2bck(2)         
+         if (btest(mlevel,3))
+     $        write(*,*)' set fact scales for PS to ',
+     $        sqrt(q2fact(1)),sqrt(q2fact(2))
       endif
 
       if (btest(mlevel,3)) then

=== modified file 'Template/SubProcesses/setcuts.f'
--- Template/SubProcesses/setcuts.f	2011-09-04 14:11:41 +0000
+++ Template/SubProcesses/setcuts.f	2011-10-05 03:29:23 +0000
@@ -566,7 +566,7 @@
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
       common/to_sprop/sprop,tprid
-      logical gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
+      integer gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
       include 'decayBW.inc'
 
 c
@@ -580,7 +580,7 @@
 
 c     Set who comes from decay based on forced BW
       do i=-(nexternal-3),-1
-         if(tprid(i,1).eq.0.and.gForceBW(i,1).or.
+         if(tprid(i,1).eq.0.and.gForceBW(i,1).eq.1.or.
      $        from_decay(i)) then
             from_decay(i)=.true.
             from_decay(iforest(1,i,1))=.true.

=== modified file 'Template/SubProcesses/setscales.f'
--- Template/SubProcesses/setscales.f	2010-11-02 20:52:12 +0000
+++ Template/SubProcesses/setscales.f	2011-10-05 03:29:23 +0000
@@ -132,6 +132,7 @@
       
       
       q2fact(1)=0d0             !factorization scale**2 for pdf1
+      q2fact(2)=0d0             !factorization scale**2 for pdf2
 
 c      call set_ren_scale(P,q2fact(1))
 c

=== removed file 'Template/SubProcesses/status'
--- Template/SubProcesses/status	2010-10-30 03:26:37 +0000
+++ Template/SubProcesses/status	1970-01-01 00:00:00 +0000
@@ -1,1 +0,0 @@
- 0 / 1 

=== modified file 'Template/SubProcesses/symmetry.f'
--- Template/SubProcesses/symmetry.f	2011-08-27 21:03:09 +0000
+++ Template/SubProcesses/symmetry.f	2011-10-05 03:29:23 +0000
@@ -233,7 +233,7 @@
       integer iarray(imax)
       logical lconflict(-max_branch:nexternal)
       logical done
-      logical gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
+      integer gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
       include 'decayBW.inc'
 
 c-----
@@ -373,7 +373,7 @@
       integer itree(2,-max_branch:-1),iconfig
       logical lconflict(-max_branch:nexternal)
       integer sprop(maxsproc,-max_branch:-1)  ! Propagator id
-      logical forcebw(-max_branch:-1)  ! Forced BW, for identical particle conflicts
+      integer forcebw(-max_branch:-1) ! Forced BW, for identical particle conflicts
 c
 c     local
 c
@@ -428,14 +428,14 @@
             if(itree(1,-i).lt.0) then
                if(iden_part(itree(1,-i)).ne.0.and.
      $           sprop(1,-i).eq.iden_part(itree(1,-i)) .or.
-     $         forcebw(itree(1,-i)).and.
+     $         forcebw(itree(1,-i)).eq.1.and.
      $           sprop(1,-i).eq.sprop(1,itree(1,-i)))
      $              iden_part(-i) = sprop(1,-i)
             endif
             if(itree(2,-i).lt.0) then
                if(iden_part(itree(2,-i)).ne.0.and.
      $           sprop(1,-i).eq.iden_part(itree(2,-i)).or.
-     $         forcebw(itree(2,-i)).and.
+     $         forcebw(itree(2,-i)).eq.1.and.
      $           sprop(1,-i).eq.sprop(1,itree(2,-i)))
      $              iden_part(-i) = sprop(1,-i)
             endif

=== modified file 'Template/SubProcesses/unwgt.f'
--- Template/SubProcesses/unwgt.f	2011-09-11 13:02:35 +0000
+++ Template/SubProcesses/unwgt.f	2011-10-05 03:29:23 +0000
@@ -476,6 +476,9 @@
          jpart(1,isym(i,jsym)) = idup(i,ip,numproc)
          jpart(2,isym(i,jsym)) = mothup(1,i)
          jpart(3,isym(i,jsym)) = mothup(2,i)
+c        Color info is filled in mothup
+         jpart(4,isym(i,jsym)) = 0
+         jpart(5,isym(i,jsym)) = 0
          jpart(6,isym(i,jsym)) = 1
       enddo
       do i=1,nincoming
@@ -550,7 +553,7 @@
 c
 c      write(*,*) 'Writing event'
       if(q2fact(1).gt.0.and.q2fact(2).gt.0)then
-         sscale = (q2fact(1)*q2fact(2))**0.25
+         sscale = sqrt(max(q2fact(1),q2fact(2)))
       else if(q2fact(1).gt.0)then
          sscale = sqrt(q2fact(1))
       else if(q2fact(2).gt.0)then

=== modified file 'Template/bin/survey'
--- Template/bin/survey	2011-08-24 23:38:57 +0000
+++ Template/bin/survey	2011-10-05 03:29:23 +0000
@@ -122,6 +122,10 @@
 		time nice ./$job >/dev/null
 		../../bin/sum_html >/dev/null
 	    fi
+	    if [[ -e ../../error ]]; then
+		echo "File error present, stop execution"
+		exit
+	    fi
 	done
 	cd ..
 	if [[ $mode == 0  ]]; then 

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2011-09-30 13:55:43 +0000
+++ UpdateNotes.txt	2011-10-05 03:29:23 +0000
@@ -1,6 +1,50 @@
 Update notes for MadGraph 5 (in reverse time order)
 
-1.3.17 (30/10/11) OM: Fix a crash created by ALOHA when it tries to create the full
+1.3.18 (04/10/11) JA: Use model information to determine color of particles
+                      for reweighting and propagator color info.
+                  JA: Changed the definition of "forbidden s-channels"
+                      denoted by "$" to exclude on-shell s-channels while
+                      keeping all diagrams (i.e., complemetary to the decay
+                      chain formalism). This reduces the problems with 
+                      gauge invariance compared to previously.
+                      "Onshell" is as usual defined by the "bwcutoff" flag 
+                      in the run_card.dat.
+                  JA: Enable proper 4-flavor matching (such as gg>hbb~+jets)
+                      Note that you need the Pythia/PGS package v. 2.1.9 or 
+                      later to use with 4F matching.
+                      Changes include: alpha_s reweighting also for b vertices,
+                      new scale treatment (mu_F for pp>hbb~ is (pT_b^max*m_Th)),
+                      no clustering of gluons to final-state massive particles
+                      in MadEvent.
+		  JA: Ensure that factorization scale settings and matching works
+                      also in singlet t-channel exchange processes like
+                      single top and VBF. The dynamic factorization
+                      scale is given by the pT of the scattered quark
+                      (on each side of the event).
+                Note: You need the Pythia/PGS package v. 2.1.10 or later
+                      to use with VBF matching, to ensure that both radiated
+                      and scattered partons are treated correctly
+                      - scattered partons need to be excluded from the matching,
+                      since their pT can be below QCUT. An even better
+                      treatment would require to individually shower and match
+                      the two sides in Pythia, which is not presently possible.
+                Note: In the matched 4-flavor process p p > t b~ j $ w+ w- t~ +
+                      p p > t b~ j j $ w+ w- t~, there is an admixture
+                      of t-channel single top (with up to 1 radiated jet) 
+                      and s-channel single top (with up to 2 radiated jets). 
+                      In this case, the automatic determination of maximum 
+                      multiplicity sample doesn't work (since max in the file 
+                      will be 2 jets, but for t-channel max is 1 jet).
+                      So MAXJETS=1 must be specified in the pythia_card.dat.
+                  JA: Fixed pdf reweighting for matching, which due to a mistake
+                      had never been activated.
+                  JA: Improved phase space integration presetting further by 
+                      taking into account special cuts like xpt, ht etc.
+                  JA: Introduce new convention for invariant mass cuts
+                      - if max < min, exclude intermediate range
+                      (allows to exclude W/Z dijet resonances in VBF processes)
+
+1.3.17 (30/09/11) OM: Fix a crash created by ALOHA when it tries to create the full
                       set of ALOHA routines (pythia8 output only).
 
 1.3.16 (11/09/11) JA: Fixed the problem from 1.3.12.

=== modified file 'madgraph/VERSION'
--- madgraph/VERSION	2011-09-30 13:55:43 +0000
+++ madgraph/VERSION	2011-10-05 03:29:23 +0000
@@ -1,3 +1,3 @@
-version = 1.3.17
-date = 2011-09-30
+version = 1.3.18
+date = 2011-10-04
 

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2011-08-27 02:01:27 +0000
+++ madgraph/core/base_objects.py	2011-10-05 03:29:23 +0000
@@ -1051,9 +1051,12 @@
 
         self['id'] = 0
         self['number'] = 0
-        # True = final, False = initial (boolean to save memory)
+        # state: True = final, False = initial (boolean to save memory)
         self['state'] = True
+        # from_group: Used in diagram generation
         self['from_group'] = True
+        # onshell: decaying leg (True), forbidden s-channel (False), none (None)
+        self['onshell'] = None
 
     def filter(self, name, value):
         """Filter for valid leg property values."""
@@ -1075,12 +1078,18 @@
                         "%s is not a valid boolean for leg flag from_group" % \
                                                                     str(value)
 
+        if name == 'onshell':
+            if not isinstance(value, bool) and value != None:
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid boolean for leg flag onshell" % \
+                                                                    str(value)
+
         return True
 
     def get_sorted_keys(self):
         """Return particle property names as a nicely sorted list."""
 
-        return ['id', 'number', 'state', 'from_group']
+        return ['id', 'number', 'state', 'from_group', 'onshell']
 
     def is_fermion(self, model):
         """Returns True if the particle corresponding to the leg is a

=== modified file 'madgraph/core/diagram_generation.py'
--- madgraph/core/diagram_generation.py	2011-07-22 19:48:59 +0000
+++ madgraph/core/diagram_generation.py	2011-10-05 03:29:23 +0000
@@ -546,14 +546,18 @@
         if process.get('forbidden_s_channels'):
             ninitial = len(filter(lambda leg: leg.get('state') == False,
                               process.get('legs')))
-            res = base_objects.DiagramList(\
-                filter(lambda diagram: \
-                       not any([vertex.get_s_channel_id(\
-                                process.get('model'), ninitial) \
-                                in process.get('forbidden_s_channels')
-                                for vertex in diagram.get('vertices')[:-1]]),
-                       res))
-
+            verts = base_objects.VertexList(sum([[vertex for vertex \
+                                                  in diagram.get('vertices')[:-1]
+                                           if vertex.get_s_channel_id(\
+                                               process.get('model'), ninitial) \
+                                           in process.get('forbidden_s_channels')] \
+                                               for diagram in res], []))
+            for vert in verts:
+                # Use onshell = False to indicate that this s-channel is forbidden
+                newleg = copy.copy(vert.get('legs').pop(-1))
+                newleg.set('onshell', False)
+                vert.get('legs').append(newleg)
+                
         # Set diagrams to res
         self['diagrams'] = res
 
@@ -577,6 +581,9 @@
                     ntlnumber = legs[-1].get('number')
                     lastleg = filter(lambda leg: leg.get('number') != ntlnumber,
                                      lastvx.get('legs'))[0]
+                    # Reset onshell in case we have forbidden s-channels
+                    if lastleg.get('onshell') == False:
+                        lastleg.set('onshell', None)
                     # Replace the last leg of nexttolastvertex
                     legs[-1] = lastleg
                     nexttolastvertex.set('legs', legs)
@@ -908,15 +915,17 @@
         vertices = []
 
         for diagram in self.get('diagrams'):
+            # Keep track of external legs (leg numbers already used)
+            leg_external = set()
             for ivx, vertex in enumerate(diagram.get('vertices')):
                 for ileg, leg in enumerate(vertex.get('legs')):
-                    if leg.get('state') and leg.get('id') in decay_ids:
-                        # Use from_group to indicate decaying legs,
+                    # Ensure that only external legs get decay flag
+                    if leg.get('state') and leg.get('id') in decay_ids and \
+                           leg.get('number') not in leg_external:
+                        # Use onshell to indicate decaying legs,
                         # i.e. legs that have decay chains
                         leg = copy.copy(leg)
-                        leg.set('from_group', True)
-                    else:
-                        leg.set('from_group', False)
+                        leg.set('onshell', True)
                     try:
                         index = legs.index(leg)
                     except ValueError:
@@ -924,6 +933,7 @@
                         legs.append(leg)
                     else: # Found a leg
                         vertex.get('legs')[ileg] = legs[index]
+                    leg_external.add(leg.get('number'))
                 try:
                     index = vertices.index(vertex)
                     diagram.get('vertices')[ivx] = vertices[index]
@@ -986,7 +996,7 @@
                 self['decay_chains'].append(\
                     DecayChainAmplitude(process, collect_mirror_procs,
                                         ignore_six_quark_processes))
-            # Flag decaying legs in the core process by from_group = True
+            # Flag decaying legs in the core process by onshell = True
             decay_ids = sum([[a.get('process').get('legs')[0].get('id') \
                               for a in dec.get('amplitudes')] for dec in \
                              self['decay_chains']], [])

=== modified file 'madgraph/core/helas_objects.py'
--- madgraph/core/helas_objects.py	2011-09-11 13:02:35 +0000
+++ madgraph/core/helas_objects.py	2011-10-05 03:29:23 +0000
@@ -85,12 +85,13 @@
         part = model.get_particle(leg.get('id'))
 
         # For legs with decay chains defined, include leg id (don't combine)
-        if leg.get('from_group'): id = leg.get('id')
+        if leg.get('onshell'): id = leg.get('id')
         else: id = 0
         # For FS legs, don't care about number (but do for IS legs)
         if leg.get('state'): number = 0
         else: number = leg.get('number')
-        return [((number, id, part.get('spin'),
+        # Include also onshell, since this specifies forbidden s-channel
+        return [((number, id, part.get('spin'), leg.get('onshell'),
                   part.get('is_part'), part.get('self_antipart'),
                   part.get('mass'), part.get('width'), part.get('color')),
                  leg.get('number'))]
@@ -191,7 +192,7 @@
         # fermionflow = 1    fermions have +-1 for flow (bosons always +1),
         #                    -1 is used only if there is a fermion flow clash
         #                    due to a Majorana particle 
-        self['state'] = 'incoming'
+        self['state'] = 'initial'
         self['leg_state'] = True
         self['mothers'] = HelasWavefunctionList()
         self['number_external'] = 0
@@ -202,8 +203,9 @@
         self['decay'] = False
         # The onshell flag is used in processes with defined decay
         # chains, to indicate that this wavefunction is decayed and
-        # should be onshell
-        self['onshell'] = False
+        # should be onshell (True), as well as for forbidden s-channels (False).
+        # Default is None
+        self['onshell'] = None
 
     # Customized constructor
     def __init__(self, *arguments):
@@ -227,6 +229,9 @@
                 self.set('number_external', leg.get('number'))
                 self.set('number', leg.get('number'))
                 self.set('state', {False: 'initial', True: 'final'}[leg.get('state')])
+                if leg.get('onshell') == False:
+                    # Denotes forbidden s-channel
+                    self.set('onshell', leg.get('onshell'))
                 self.set('leg_state', leg.get('state'))
                 # Need to set 'decay' to True for particles which will be
                 # decayed later, in order to not combine such processes
@@ -365,11 +370,17 @@
                       "%s is not a valid list of mothers for wavefunction" % \
                       str(value)
 
-        if name in ['decay', 'onshell']:
+        if name in ['decay']:
             if not isinstance(value, bool):
                 raise self.PhysicsObjectError, \
                         "%s is not a valid bool" % str(value) + \
-                        " for decay or onshell"
+                        " for decay"
+
+        if name in ['onshell']:
+            if not isinstance(value, bool) and value != None:
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid bool" % str(value) + \
+                        " for onshell"
 
         return True
 
@@ -933,32 +944,33 @@
         # Generate last vertex
         legs = base_objects.LegList()
 
-        # We use the from_group flag to indicate whether this outgoing
-        # leg corresponds to a decaying (onshell) particle or not
+        # We use the onshell flag to indicate whether this outgoing
+        # leg corresponds to a decaying (onshell) particle, forbidden
+        # s-channel, or regular
         try:
-            lastleg = wf_dict[self.get('number')]
+            lastleg = wf_dict[(self.get('number'),self.get('onshell'))]
         except KeyError:            
             lastleg = base_objects.Leg({
                 'id': self.get_pdg_code(),
                 'number': self.get('number_external'),
                 'state': self.get('leg_state'),
-                'from_group': self.get('onshell')
+                'onshell': self.get('onshell')
                 })
             if optimization != 0:
-                wf_dict[self.get('number')] = lastleg
+                wf_dict[(self.get('number'),self.get('onshell'))] = lastleg
 
         for mother in self.get('mothers'):
             try:
-                leg = wf_dict[mother.get('number')]
+                leg = wf_dict[(mother.get('number'),False)]
             except KeyError:
                 leg = base_objects.Leg({
                     'id': mother.get_pdg_code(),
                     'number': mother.get('number_external'),
                     'state': mother.get('leg_state'),
-                    'from_group': mother.get('onshell')
+                    'onshell': None
                     })
                 if optimization != 0:
-                    wf_dict[mother.get('number')] = leg
+                    wf_dict[(mother.get('number'),False)] = leg
             legs.append(leg)
 
         legs.append(lastleg)
@@ -1025,7 +1037,7 @@
                     'id': mother.get_pdg_code(),
                     'number': mother.get('number_external'),
                     'state': mother.get('leg_state'),
-                    'from_group': False
+                    'onshell': mother.get('onshell')
                     }))
 
             if init_mothers[0].get('number_external') == 1 and \
@@ -1091,7 +1103,7 @@
                     'id': mother.get_pdg_code(),
                     'number': mother.get('number_external'),
                     'state': mother.get('leg_state'),
-                    'from_group': False
+                    'onshell': mother.get('onshell')
                     }))
             legs.insert(0, mother_leg)
 
@@ -1729,16 +1741,16 @@
         legs = base_objects.LegList()
         for mother in self.get('mothers'):
             try:
-                leg = wf_dict[mother.get('number')]
+                leg = wf_dict[(mother.get('number'),False)]
             except KeyError:
                 leg = base_objects.Leg({
                     'id': mother.get_pdg_code(),
                     'number': mother.get('number_external'),
                     'state': mother.get('leg_state'),
-                    'from_group': mother.get('onshell')
+                    'onshell': None
                     })
                 if optimization != 0:
-                    wf_dict[mother.get('number')] = leg
+                    wf_dict[(mother.get('number'),False)] = leg
             legs.append(leg)
 
         return base_objects.Vertex({
@@ -1781,7 +1793,7 @@
                     'id': mother.get_pdg_code(),
                     'number': mother.get('number_external'),
                     'state': mother.get('leg_state'),
-                    'from_group': False
+                    'onshell': mother.get('onshell')
                     }))
 
             # Renumber resulting leg according to minimum leg number
@@ -1815,7 +1827,7 @@
                     'id': mother.get_pdg_code(),
                     'number': mother.get('number_external'),
                     'state': mother.get('leg_state'),
-                    'from_group': False
+                    'onshell': mother.get('onshell')
                     }))
             # Renumber resulting leg according to minimum leg number
             legs[-1].set('number', min([l.get('number') for l in legs[:-1]]))
@@ -1857,7 +1869,7 @@
                 popped_legs.append(base_objects.Leg({'id': 21,
                     'number': min([l.get('number') for l in popped_legs]),
                     'state': True,
-                    'from_group': False}))
+                    'onshell': None}))
 
                 new_vertex = base_objects.Vertex({
                     'id': vertex.get('id'),

=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2011-09-11 13:02:35 +0000
+++ madgraph/iolibs/export_v4.py	2011-10-05 03:29:23 +0000
@@ -621,7 +621,7 @@
                         if subproc_group:
                             pdf_lines = pdf_lines + \
                                         ("%s%d=PDG2PDF(ABS(LPP(IB(%d))),%d*LP," + \
-                                         "XBK(IB(%d)),DSQRT(Q2FACT(IB(%d))))\n") % \
+                                         "XBK(IB(%d)),DSQRT(Q2FACT(%d)))\n") % \
                                          (pdf_codes[initial_state],
                                           i + 1, i + 1, pdgtopdf[initial_state],
                                           i + 1, i + 1)
@@ -1006,6 +1006,10 @@
                              mapconfigs,
                              matrix_element)
 
+        filename = 'get_color.f'
+        self.write_colors_file(writers.FortranWriter(filename),
+                               matrix_element)
+
         filename = 'decayBW.inc'
         self.write_decayBW_file(writers.FortranWriter(filename),
                            s_and_t_channels)
@@ -1410,6 +1414,58 @@
         return True
 
     #===========================================================================
+    # write_coloramps_file
+    #===========================================================================
+    def write_colors_file(self, writer, matrix_elements):
+        """Write the get_color.f file for MadEvent, which returns color
+        for all particles used in the matrix element."""
+
+        if isinstance(matrix_elements, helas_objects.HelasMatrixElement):
+            matrix_elements = [matrix_elements]
+
+        model = matrix_elements[0].get('processes')[0].get('model')
+
+        # We need the both particle and antiparticle wf_ids, since the identity
+        # depends on the direction of the wf.
+        wf_ids = set(sum([sum([sum([[wf.get_pdg_code(),wf.get_anti_pdg_code()] \
+                                    for wf in d.get('wavefunctions')],[]) \
+                               for d in me.get('diagrams')], []) \
+                          for me in matrix_elements], []))
+
+        leg_ids = set(sum([sum([[l.get('id') for l in \
+                                 p.get_legs_with_decays()] for p in \
+                                me.get('processes')], []) for me in
+                           matrix_elements], []))
+        particle_ids = sorted(list(wf_ids.union(leg_ids)))
+
+        lines = """function get_color(ipdg)
+        implicit none
+        integer get_color, ipdg
+
+        if(ipdg.eq.%d)then
+        get_color=%d
+        return
+        """ % (particle_ids[0], model.get_particle(particle_ids[0]).get_color())
+
+        for part_id in particle_ids[1:]:
+            lines += """else if(ipdg.eq.%d)then
+            get_color=%d
+            return
+            """ % (part_id, model.get_particle(part_id).get_color())
+        lines += """else
+        write(*,*)'Error: No color given for pdg ',ipdg
+        get_color=0        
+        return
+        endif
+        end
+        """
+        
+        # Write the file
+        writer.writelines(lines)
+
+        return True
+
+    #===========================================================================
     # write_maxconfigs_file
     #===========================================================================
     def write_maxconfigs_file(self, writer, matrix_elements):
@@ -1591,17 +1647,17 @@
 
         lines = []
 
-        booldict = {False: ".false.", True: ".true."}
+        booldict = {None: "0", True: "1", False: "2"}
 
         for iconf, config in enumerate(s_and_t_channels):
             schannels = config[0]
             for vertex in schannels:
                 # For the resulting leg, pick out whether it comes from
-                # decay or not, as given by the from_group flag
+                # decay or not, as given by the onshell flag
                 leg = vertex.get('legs')[-1]
                 lines.append("data gForceBW(%d,%d)/%s/" % \
                              (leg.get('number'), iconf + 1,
-                              booldict[leg.get('from_group')]))
+                              booldict[leg.get('onshell')]))
 
         # Write the file
         writer.writelines(lines)
@@ -2067,6 +2123,10 @@
                                    maxflows,
                                    matrix_elements)
 
+        filename = 'get_color.f'
+        self.write_colors_file(writers.FortranWriter(filename),
+                               matrix_elements)
+
         filename = 'config_subproc_map.inc'
         self.write_config_subproc_map_file(writers.FortranWriter(filename),
                                            subproc_diagrams_for_config)

=== modified file 'madgraph/iolibs/template_files/super_auto_dsig_group_v4.inc'
--- madgraph/iolibs/template_files/super_auto_dsig_group_v4.inc	2011-05-18 02:29:29 +0000
+++ madgraph/iolibs/template_files/super_auto_dsig_group_v4.inc	2011-10-05 03:29:23 +0000
@@ -149,14 +149,14 @@
 C     Select among the subprocesses based on PDF weight
       SUMPROB=0D0
       DO J=1,SYMCONF(0)
-         DO I=1,MAXSPROC
-            IF(CONFSUB(I,SYMCONF(J)).NE.0) THEN
-               DO K=1,2
-                  IF(K.EQ.1.OR.MIRRORPROCS(I))THEN
+         DO IPROC=1,MAXSPROC
+            IF(CONFSUB(IPROC,SYMCONF(J)).NE.0) THEN
+               DO IMIRROR=1,2
+                  IF(IMIRROR.EQ.1.OR.MIRRORPROCS(IPROC))THEN
 C                 Calculate PDF weight for all subprocesses
-                     SELPROC(K,I,J)=DSIGPROC(PP,J,I,K,SYMCONF,CONFSUB,1d0,4)
-                     SUMPROB=SUMPROB+SELPROC(K,I,J)
-                     IF(K.EQ.2)THEN
+                     SELPROC(IMIRROR,IPROC,J)=DSIGPROC(PP,J,IPROC,IMIRROR,SYMCONF,CONFSUB,1d0,4)
+                     SUMPROB=SUMPROB+SELPROC(IMIRROR,IPROC,J)
+                     IF(IMIRROR.EQ.2)THEN
 C                 Need to flip back x values
                         XDUM=XBK(1)
                         XBK(1)=XBK(2)

=== modified file 'tests/acceptance_tests/test_cmd.py'
--- tests/acceptance_tests/test_cmd.py	2011-08-07 19:55:37 +0000
+++ tests/acceptance_tests/test_cmd.py	2011-10-05 03:29:23 +0000
@@ -182,6 +182,10 @@
                                                     'SubProcesses',
                                                     'P0_epem_epem',
                                                     'maxconfigs.inc')))
+        self.assertTrue(os.path.exists(os.path.join(self.out_dir,
+                                                    'SubProcesses',
+                                                    'P0_epem_epem',
+                                                    'get_color.f')))
         self.assertFalse(os.path.exists(os.path.join(self.out_dir,
                                                     'SubProcesses',
                                                     'P0_epem_epem',
@@ -613,6 +617,10 @@
                                                     'SubProcesses',
                                                     'P2_gg_qq',
                                                     'maxconfigs.inc')))
+        self.assertTrue(os.path.exists(os.path.join(self.out_dir,
+                                                    'SubProcesses',
+                                                    'P2_gg_qq',
+                                                    'get_color.f')))
         # Check that the Source directory compiles
         status = subprocess.call(['make'],
                                  stdout=devnull, 

=== modified file 'tests/input_files/e+e-_e+e-.pkl'
--- tests/input_files/e+e-_e+e-.pkl	2011-07-21 04:45:53 +0000
+++ tests/input_files/e+e-_e+e-.pkl	2011-10-05 03:29:23 +0000
@@ -41,11982 +41,12003 @@
 p14
 g5
 (dp15
-S'state'
+S'onshell'
 p16
+NsS'state'
+p17
 I00
 sS'id'
-p17
+p18
 I-11
 sS'from_group'
-p18
+p19
 I01
 sS'number'
-p19
+p20
 I1
-stp20
-Rp21
+stp21
+Rp22
 ag0
 (g14
 g5
-(dp22
+(dp23
 g16
+Nsg17
 I00
-sg17
+sg18
 I11
-sg18
+sg19
 I01
-sg19
+sg20
 I2
-stp23
-Rp24
+stp24
+Rp25
 ag0
 (g14
 g5
-(dp25
+(dp26
 g16
+Nsg17
 I01
-sg17
+sg18
 I-11
-sg18
+sg19
 I01
-sg19
+sg20
 I3
-stp26
-Rp27
+stp27
+Rp28
 ag0
 (g14
 g5
-(dp28
+(dp29
 g16
+Nsg17
 I01
-sg17
+sg18
 I11
-sg18
-I01
 sg19
+I01
+sg20
 I4
-stp29
-Rp30
-atp31
-Rp32
+stp30
+Rp31
+atp32
+Rp33
 sS'uid'
-p33
+p34
 I0
 sS'forbidden_particles'
-p34
-(lp35
+p35
+(lp36
 sS'orders'
-p36
-(dp37
+p37
+(dp38
 S'WEIGHTED'
-p38
+p39
 I4
 ssS'forbidden_s_channels'
-p39
-(lp40
+p40
+(lp41
 sS'required_s_channels'
-p41
-(lp42
+p42
+(lp43
 sS'model'
-p43
+p44
 g0
 (cmadgraph.core.base_objects
 Model
-p44
+p45
 g5
-(dp45
+(dp46
 S'ref_dict_to0'
-p46
-(dp47
+p47
+(dp48
 (I-4
 I4
 I21
-tp48
-(lp49
+tp49
+(lp50
 I4
 as(I-24
 I24
-tp50
-(lp51
+tp51
+(lp52
 I0
 as(I-1
 I1
 I23
-tp52
-(lp53
+tp53
+(lp54
 I18
 as(I-5
 I5
 I23
-tp54
-(lp55
+tp55
+(lp56
 I22
 as(I-15
 I15
 I25
-tp56
-(lp57
+tp57
+(lp58
 I42
 as(I-24
 I-11
 I12
-tp58
-(lp59
+tp59
+(lp60
 I39
 as(I-24
 I-15
 I16
-tp60
-(lp61
+tp61
+(lp62
 I41
 as(I21
 I21
 I21
 I21
-tp62
-(lp63
+tp63
+(lp64
 I8
 as(I11
 I-11
-tp64
-(lp65
+tp65
+(lp66
 I0
 as(I-13
 I13
 I22
-tp66
-(lp67
+tp67
+(lp68
 I16
 as(I-4
 I4
-tp68
-(lp69
+tp69
+(lp70
 I0
 as(I-5
 I5
-tp70
-(lp71
+tp71
+(lp72
 I0
 as(I-1
 I1
 I21
-tp72
-(lp73
+tp73
+(lp74
 I1
 as(I-3
 I3
 I21
-tp74
-(lp75
+tp75
+(lp76
 I3
 as(I-2
 I2
-tp76
-(lp77
+tp77
+(lp78
 I0
 as(I5
 I-5
-tp78
-(lp79
+tp79
+(lp80
 I0
 as(I-6
 I5
 I24
-tp80
-(lp81
+tp81
+(lp82
 I35
 as(I-13
 I13
 I23
-tp82
-(lp83
+tp83
+(lp84
 I25
 as(I1
 I-1
-tp84
-(lp85
+tp85
+(lp86
 I0
 as(I-24
 I-3
 I4
-tp86
-(lp87
+tp87
+(lp88
 I31
 as(I21
 I21
-tp88
-(lp89
+tp89
+(lp90
 I0
 as(I-1
 I1
-tp90
-(lp91
+tp91
+(lp92
 I0
 as(I16
 I-16
-tp92
-(lp93
+tp93
+(lp94
 I0
 as(I-3
 I3
 I23
-tp94
-(lp95
+tp95
+(lp96
 I20
 as(I-2
 I1
 I24
-tp96
-(lp97
+tp97
+(lp98
 I33
 as(I12
 I-12
-tp98
-(lp99
+tp99
+(lp100
 I0
 as(I-2
 I2
 I23
-tp100
-(lp101
+tp101
+(lp102
 I19
 as(I4
 I-4
-tp102
-(lp103
+tp103
+(lp104
 I0
 as(I-16
 I16
 I23
-tp104
-(lp105
+tp105
+(lp106
 I29
 as(I-13
 I13
-tp106
-(lp107
+tp107
+(lp108
 I0
 as(I23
 I23
 I25
 I25
-tp108
-(lp109
+tp109
+(lp110
 I55
 as(I-2
 I2
 I22
-tp110
-(lp111
+tp111
+(lp112
 I10
 as(I-24
 I23
 I24
-tp112
-(lp113
+tp113
+(lp114
 I46
 as(I-24
 I23
 I23
 I24
-tp114
-(lp115
+tp115
+(lp116
 I52
 as(I-24
 I24
 I25
-tp116
-(lp117
+tp117
+(lp118
 I47
 as(I-15
 I15
 I23
-tp118
-(lp119
+tp119
+(lp120
 I26
 as(I-11
 I11
-tp120
-(lp121
+tp121
+(lp122
 I0
 as(I-2
 I2
 I21
-tp122
-(lp123
+tp123
+(lp124
 I2
 as(I23
 I23
-tp124
-(lp125
+tp125
+(lp126
 I0
 as(I-5
 I5
 I22
-tp126
-(lp127
+tp127
+(lp128
 I13
 as(I-14
 I14
-tp128
-(lp129
+tp129
+(lp130
 I0
 as(I6
 I-6
-tp130
-(lp131
+tp131
+(lp132
 I0
 as(I-5
 I5
 I25
-tp132
-(lp133
+tp133
+(lp134
 I43
 as(I23
 I23
 I25
-tp134
-(lp135
+tp135
+(lp136
 I48
 as(I14
 I-14
-tp136
-(lp137
+tp137
+(lp138
 I0
 as(I-16
 I15
 I24
-tp138
-(lp139
+tp139
+(lp140
 I38
 as(I-24
 I-13
 I14
-tp140
-(lp141
+tp141
+(lp142
 I40
 as(I-3
 I3
-tp142
-(lp143
+tp143
+(lp144
 I0
 as(I22
 I22
-tp144
-(lp145
+tp145
+(lp146
 I0
 as(I24
 I-24
-tp146
-(lp147
+tp147
+(lp148
 I0
 as(I2
 I-2
-tp148
-(lp149
+tp149
+(lp150
 I0
 as(I-6
 I6
 I25
-tp150
-(lp151
+tp151
+(lp152
 I44
 as(I25
 I25
-tp152
-(lp153
+tp153
+(lp154
 I0
 as(I-3
 I3
 I22
-tp154
-(lp155
+tp155
+(lp156
 I11
 as(I-12
 I12
 I23
-tp156
-(lp157
+tp157
+(lp158
 I27
 as(I-16
 I16
-tp158
-(lp159
+tp159
+(lp160
 I0
 as(I-5
 I5
 I21
-tp160
-(lp161
+tp161
+(lp162
 I5
 as(I-12
 I12
-tp162
-(lp163
+tp163
+(lp164
 I0
 as(I-6
 I6
 I23
-tp164
-(lp165
+tp165
+(lp166
 I23
 as(I-24
 I22
 I24
-tp166
-(lp167
+tp167
+(lp168
 I45
 as(I-12
 I11
 I24
-tp168
-(lp169
+tp169
+(lp170
 I36
 as(I21
 I21
 I21
-tp170
-(lp171
+tp171
+(lp172
 I7
 as(I-15
 I15
-tp172
-(lp173
+tp173
+(lp174
 I0
 as(I3
 I-3
-tp174
-(lp175
+tp175
+(lp176
 I0
 as(I-6
 I6
 I22
-tp176
-(lp177
+tp177
+(lp178
 I14
 as(I25
 I25
 I25
-tp178
-(lp179
+tp179
+(lp180
 I49
 as(I-24
 I-24
 I24
 I24
-tp180
-(lp181
+tp181
+(lp182
 I53
 as(I-15
 I15
 I22
-tp182
-(lp183
+tp183
+(lp184
 I17
 as(I-11
 I11
 I23
-tp184
-(lp185
+tp185
+(lp186
 I24
 as(I-24
 I24
 I25
 I25
-tp186
-(lp187
+tp187
+(lp188
 I54
 as(I-6
 I6
 I21
-tp188
-(lp189
+tp189
+(lp190
 I6
 as(I-24
 I-1
 I2
-tp190
-(lp191
+tp191
+(lp192
 I30
 as(I-11
 I11
 I22
-tp192
-(lp193
+tp193
+(lp194
 I15
 as(I-14
 I13
 I24
-tp194
-(lp195
+tp195
+(lp196
 I37
 as(I-4
 I4
 I23
-tp196
-(lp197
+tp197
+(lp198
 I21
 as(I-24
 I22
 I22
 I24
-tp198
-(lp199
+tp199
+(lp200
 I50
 as(I-24
 I22
 I23
 I24
-tp200
-(lp201
+tp201
+(lp202
 I51
 as(I-24
 I-5
 I6
-tp202
-(lp203
+tp203
+(lp204
 I32
 as(I13
 I-13
-tp204
-(lp205
+tp205
+(lp206
 I0
 as(I-4
 I4
 I22
-tp206
-(lp207
+tp207
+(lp208
 I12
 as(I-6
 I6
-tp208
-(lp209
+tp209
+(lp210
 I0
 as(I-4
 I3
 I24
-tp210
-(lp211
+tp211
+(lp212
 I34
 as(I-14
 I14
 I23
-tp212
-(lp213
+tp213
+(lp214
 I28
 as(I-1
 I1
 I22
-tp214
-(lp215
+tp215
+(lp216
 I9
 as(I15
 I-15
-tp216
-(lp217
+tp217
+(lp218
 I0
 assS'interactions'
-p218
+p219
 g0
 (cmadgraph.core.base_objects
 InteractionList
-p219
+p220
 g2
-(lp220
+(lp221
 g0
 (cmadgraph.core.base_objects
 Interaction
-p221
+p222
 g5
-(dp222
+(dp223
 S'particles'
-p223
+p224
 g0
 (cmadgraph.core.base_objects
 ParticleList
-p224
+p225
 g2
-(lp225
+(lp226
 g0
 (cmadgraph.core.base_objects
 Particle
-p226
+p227
 g5
-(dp227
+(dp228
 S'texname'
-p228
-S'd'
 p229
+S'd'
+p230
 sS'is_part'
-p230
+p231
 I00
 sS'name'
-p231
-S'd'
 p232
+S'd'
+p233
 sS'self_antipart'
-p233
+p234
 I00
 sS'color'
-p234
+p235
 I3
 sS'width'
-p235
-S'ZERO'
 p236
+S'ZERO'
+p237
 sS'charge'
-p237
+p238
 F0.0
 sS'mass'
-p238
-S'ZERO'
 p239
+S'ZERO'
+p240
 sS'antiname'
-p240
+p241
 S'd~'
-p241
+p242
 sS'line'
-p242
+p243
 S'straight'
-p243
+p244
 sS'propagating'
-p244
+p245
 I01
 sS'spin'
-p245
+p246
 I2
 sS'antitexname'
-p246
-g229
+p247
+g230
 sS'pdg_code'
-p247
+p248
 I1
-stp248
-Rp249
+stp249
+Rp250
 ag0
-(g226
+(g227
 g5
-(dp250
-g228
+(dp251
 g229
-sg230
-I01
+g230
 sg231
-g232
-sg233
+I01
+sg232
+g233
+sg234
 I00
-sg234
+sg235
 I3
-sg235
-g236
-sg237
+sg236
+g237
+sg238
 F0.0
-sg238
-g239
-sg240
-g241
-sg242
-g243
-sg244
-I01
+sg239
+g240
+sg241
+g242
+sg243
+g244
 sg245
+I01
+sg246
 I2
-sg246
-g229
 sg247
+g230
+sg248
 I1
-stp251
-Rp252
+stp252
+Rp253
 ag0
-(g226
+(g227
 g5
-(dp253
-g228
+(dp254
+g229
 S'_'
-p254
-sg230
-I01
-sg231
-S'g'
 p255
-sg233
+sg231
 I01
-sg234
-I8
-sg235
-S'ZERO'
+sg232
+S'g'
 p256
-sg237
-F0.0
-sg238
+sg234
+I01
+sg235
+I8
+sg236
 S'ZERO'
 p257
-sg240
-S'g'
+sg238
+F0.0
+sg239
+S'ZERO'
 p258
-sg242
+sg241
+S'g'
+p259
+sg243
 S'curly'
-p259
-sg244
+p260
+sg245
 I01
-sg245
+sg246
 I3
-sg246
-g254
 sg247
+g255
+sg248
 I21
-stp260
-Rp261
-atp262
-Rp263
+stp261
+Rp262
+atp263
+Rp264
 sS'lorentz'
-p264
-(lp265
+p265
+(lp266
 S''
-p266
-asg234
-(lp267
+p267
+asg235
+(lp268
 g0
 (cmadgraph.core.color_algebra
 ColorString
-p268
+p269
 g2
-(lp269
+(lp270
 cmadgraph.core.color_algebra
 T
-p270
+p271
 (I2
 I1
 I0
-tp271
-Rp272
-atp273
-Rp274
-(dp275
+tp272
+Rp273
+atp274
+Rp275
+(dp276
 S'coeff'
-p276
+p277
 cfractions
 Fraction
-p277
+p278
 (S'1'
-p278
-tp279
-Rp280
+p279
+tp280
+Rp281
 sS'canonical'
-p281
+p282
 (((S'T'
-p282
+p283
 (I1
 I2
 I3
-tp283
 tp284
 tp285
-(dp286
+tp286
+(dp287
 I0
 I3
 sI1
 I2
 sI2
 I1
-stp287
+stp288
 sS'immutable'
-p288
-((g282
+p289
+((g283
 (I2
 I1
 I0
-tp289
 tp290
 tp291
+tp292
 sS'is_imaginary'
-p292
+p293
 I00
 sS'Nc_power'
-p293
+p294
 I0
-sbasg17
+sbasg18
 I1
-sg36
-(dp294
+sg37
+(dp295
 S'QCD'
-p295
+p296
 I1
 ssS'couplings'
-p296
-(dp297
+p297
+(dp298
 (I0
 I0
-tp298
+tp299
 S'GG'
-p299
-sstp300
-Rp301
+p300
+sstp301
+Rp302
 ag0
-(g221
+(g222
 g5
-(dp302
-g223
+(dp303
+g224
 g0
-(g224
+(g225
 g2
-(lp303
+(lp304
 g0
-(g226
+(g227
 g5
-(dp304
-g228
-S'u'
-p305
-sg230
-I00
-sg231
+(dp305
+g229
 S'u'
 p306
-sg233
+sg231
 I00
-sg234
-I3
-sg235
-S'ZERO'
+sg232
+S'u'
 p307
-sg237
-F0.0
-sg238
+sg234
+I00
+sg235
+I3
+sg236
 S'ZERO'
 p308
-sg240
+sg238
+F0.0
+sg239
+S'ZERO'
+p309
+sg241
 S'u~'
-p309
-sg242
-g243
-sg244
-I01
+p310
+sg243
+g244
 sg245
-I2
+I01
 sg246
-g305
+I2
 sg247
+g306
+sg248
 I2
-stp310
-Rp311
+stp311
+Rp312
 ag0
-(g226
+(g227
 g5
-(dp312
-g228
-g305
-sg230
-I01
+(dp313
+g229
+g306
 sg231
-g306
-sg233
+I01
+sg232
+g307
+sg234
 I00
-sg234
+sg235
 I3
-sg235
-g307
-sg237
-F0.0
-sg238
+sg236
 g308
-sg240
+sg238
+F0.0
+sg239
 g309
-sg242
-g243
-sg244
-I01
+sg241
+g310
+sg243
+g244
 sg245
-I2
+I01
 sg246
-g305
+I2
 sg247
+g306
+sg248
 I2
-stp313
-Rp314
+stp314
+Rp315
 ag0
-(g226
+(g227
 g5
-(dp315
-g228
-g254
-sg230
-I01
+(dp316
+g229
+g255
 sg231
-g255
-sg233
 I01
+sg232
+g256
 sg234
+I01
+sg235
 I8
-sg235
-g256
-sg237
-F0.0
-sg238
+sg236
 g257
-sg240
+sg238
+F0.0
+sg239
 g258
-sg242
+sg241
 g259
-sg244
+sg243
+g260
+sg245
 I01
-sg245
+sg246
 I3
-sg246
-g254
 sg247
+g255
+sg248
 I21
-stp316
-Rp317
-atp318
-Rp319
-sg264
-(lp320
-g266
-asg234
+stp317
+Rp318
+atp319
+Rp320
+sg265
 (lp321
-g0
-(g268
-g2
+g267
+asg235
 (lp322
-g270
+g0
+(g269
+g2
+(lp323
+g271
 (I2
 I1
 I0
-tp323
-Rp324
-atp325
-Rp326
-(dp327
-g276
-g280
-sg281
-(((g282
+tp324
+Rp325
+atp326
+Rp327
+(dp328
+g277
+g281
+sg282
+(((g283
 (I1
 I2
 I3
-tp328
 tp329
 tp330
-(dp331
+tp331
+(dp332
 I0
 I3
 sI1
 I2
 sI2
 I1
-stp332
-sg288
-((g282
+stp333
+sg289
+((g283
 (I2
 I1
 I0
-tp333
 tp334
 tp335
-sg292
-I00
+tp336
 sg293
+I00
+sg294
 I0
-sbasg17
+sbasg18
 I2
-sg36
-(dp336
+sg37
+(dp337
 S'QCD'
-p337
+p338
 I1
-ssg296
-(dp338
-g298
+ssg297
+(dp339
+g299
 S'GG'
-p339
-sstp340
-Rp341
+p340
+sstp341
+Rp342
 ag0
-(g221
+(g222
 g5
-(dp342
-g223
+(dp343
+g224
 g0
-(g224
+(g225
 g2
-(lp343
+(lp344
 g0
-(g226
+(g227
 g5
-(dp344
-g228
-S's'
-p345
-sg230
-I00
-sg231
+(dp345
+g229
 S's'
 p346
-sg233
+sg231
 I00
-sg234
-I3
-sg235
-S'ZERO'
+sg232
+S's'
 p347
-sg237
-F0.0
-sg238
+sg234
+I00
+sg235
+I3
+sg236
 S'ZERO'
 p348
-sg240
+sg238
+F0.0
+sg239
+S'ZERO'
+p349
+sg241
 S's~'
-p349
-sg242
-g243
-sg244
-I01
+p350
+sg243
+g244
 sg245
+I01
+sg246
 I2
-sg246
-g345
 sg247
+g346
+sg248
 I3
-stp350
-Rp351
+stp351
+Rp352
 ag0
-(g226
+(g227
 g5
-(dp352
-g228
-g345
-sg230
-I01
+(dp353
+g229
+g346
 sg231
-g346
-sg233
+I01
+sg232
+g347
+sg234
 I00
-sg234
+sg235
 I3
-sg235
-g347
-sg237
-F0.0
-sg238
+sg236
 g348
-sg240
+sg238
+F0.0
+sg239
 g349
-sg242
-g243
-sg244
-I01
+sg241
+g350
+sg243
+g244
 sg245
+I01
+sg246
 I2
-sg246
-g345
 sg247
+g346
+sg248
 I3
-stp353
-Rp354
+stp354
+Rp355
 ag0
-(g226
+(g227
 g5
-(dp355
-g228
-g254
-sg230
-I01
+(dp356
+g229
+g255
 sg231
-g255
-sg233
 I01
+sg232
+g256
 sg234
+I01
+sg235
 I8
-sg235
-g256
-sg237
-F0.0
-sg238
+sg236
 g257
-sg240
+sg238
+F0.0
+sg239
 g258
-sg242
+sg241
 g259
-sg244
+sg243
+g260
+sg245
 I01
-sg245
+sg246
 I3
-sg246
-g254
 sg247
+g255
+sg248
 I21
-stp356
-Rp357
-atp358
-Rp359
-sg264
-(lp360
-g266
-asg234
+stp357
+Rp358
+atp359
+Rp360
+sg265
 (lp361
-g0
-(g268
-g2
+g267
+asg235
 (lp362
-g270
+g0
+(g269
+g2
+(lp363
+g271
 (I2
 I1
 I0
-tp363
-Rp364
-atp365
-Rp366
-(dp367
-g276
-g280
-sg281
-(((g282
+tp364
+Rp365
+atp366
+Rp367
+(dp368
+g277
+g281
+sg282
+(((g283
 (I1
 I2
 I3
-tp368
 tp369
 tp370
-(dp371
+tp371
+(dp372
 I0
 I3
 sI1
 I2
 sI2
 I1
-stp372
-sg288
-((g282
+stp373
+sg289
+((g283
 (I2
 I1
 I0
-tp373
 tp374
 tp375
-sg292
-I00
+tp376
 sg293
+I00
+sg294
 I0
-sbasg17
+sbasg18
 I3
-sg36
-(dp376
+sg37
+(dp377
 S'QCD'
-p377
+p378
 I1
-ssg296
-(dp378
-g298
+ssg297
+(dp379
+g299
 S'GG'
-p379
-sstp380
-Rp381
+p380
+sstp381
+Rp382
 ag0
-(g221
+(g222
 g5
-(dp382
-g223
+(dp383
+g224
 g0
-(g224
+(g225
 g2
-(lp383
+(lp384
 g0
-(g226
+(g227
 g5
-(dp384
-g228
-S'c'
-p385
-sg230
-I00
-sg231
+(dp385
+g229
 S'c'
 p386
-sg233
+sg231
 I00
-sg234
-I3
-sg235
-S'ZERO'
+sg232
+S'c'
 p387
-sg237
-F0.0
-sg238
+sg234
+I00
+sg235
+I3
+sg236
 S'ZERO'
 p388
-sg240
+sg238
+F0.0
+sg239
+S'ZERO'
+p389
+sg241
 S'c~'
-p389
-sg242
-g243
-sg244
-I01
+p390
+sg243
+g244
 sg245
+I01
+sg246
 I2
-sg246
-g385
 sg247
+g386
+sg248
 I4
-stp390
-Rp391
+stp391
+Rp392
 ag0
-(g226
+(g227
 g5
-(dp392
-g228
-g385
-sg230
-I01
+(dp393
+g229
+g386
 sg231
-g386
-sg233
+I01
+sg232
+g387
+sg234
 I00
-sg234
+sg235
 I3
-sg235
-g387
-sg237
-F0.0
-sg238
+sg236
 g388
-sg240
+sg238
+F0.0
+sg239
 g389
-sg242
-g243
-sg244
-I01
+sg241
+g390
+sg243
+g244
 sg245
+I01
+sg246
 I2
-sg246
-g385
 sg247
+g386
+sg248
 I4
-stp393
-Rp394
+stp394
+Rp395
 ag0
-(g226
+(g227
 g5
-(dp395
-g228
-g254
-sg230
-I01
+(dp396
+g229
+g255
 sg231
-g255
-sg233
 I01
+sg232
+g256
 sg234
+I01
+sg235
 I8
-sg235
-g256
-sg237
-F0.0
-sg238
+sg236
 g257
-sg240
+sg238
+F0.0
+sg239
 g258
-sg242
+sg241
 g259
-sg244
+sg243
+g260
+sg245
 I01
-sg245
+sg246
 I3
-sg246
-g254
 sg247
+g255
+sg248
 I21
-stp396
-Rp397
-atp398
-Rp399
-sg264
-(lp400
-g266
-asg234
+stp397
+Rp398
+atp399
+Rp400
+sg265
 (lp401
-g0
-(g268
-g2
+g267
+asg235
 (lp402
-g270
+g0
+(g269
+g2
+(lp403
+g271
 (I2
 I1
 I0
-tp403
-Rp404
-atp405
-Rp406
-(dp407
-g276
-g280
-sg281
-(((g282
+tp404
+Rp405
+atp406
+Rp407
+(dp408
+g277
+g281
+sg282
+(((g283
 (I1
 I2
 I3
-tp408
 tp409
 tp410
-(dp411
+tp411
+(dp412
 I0
 I3
 sI1
 I2
 sI2
 I1
-stp412
-sg288
-((g282
+stp413
+sg289
+((g283
 (I2
 I1
 I0
-tp413
 tp414
 tp415
-sg292
-I00
+tp416
 sg293
+I00
+sg294
 I0
-sbasg17
+sbasg18
 I4
-sg36
-(dp416
+sg37
+(dp417
 S'QCD'
-p417
+p418
 I1
-ssg296
-(dp418
-g298
+ssg297
+(dp419
+g299
 S'GG'
-p419
-sstp420
-Rp421
+p420
+sstp421
+Rp422
 ag0
-(g221
+(g222
 g5
-(dp422
-g223
+(dp423
+g224
 g0
-(g224
+(g225
 g2
-(lp423
+(lp424
 g0
-(g226
+(g227
 g5
-(dp424
-g228
-S'b'
-p425
-sg230
-I00
-sg231
+(dp425
+g229
 S'b'
 p426
-sg233
+sg231
 I00
-sg234
-I3
-sg235
-S'ZERO'
+sg232
+S'b'
 p427
-sg237
+sg234
+I00
+sg235
+I3
+sg236
+S'ZERO'
+p428
+sg238
 F0.0
-sg238
+sg239
 S'BMASS'
-p428
-sg240
+p429
+sg241
 S'b~'
-p429
-sg242
-g243
-sg244
-I01
+p430
+sg243
+g244
 sg245
+I01
+sg246
 I2
-sg246
-g425
 sg247
+g426
+sg248
 I5
-stp430
-Rp431
+stp431
+Rp432
 ag0
-(g226
+(g227
 g5
-(dp432
-g228
-g425
-sg230
-I01
+(dp433
+g229
+g426
 sg231
-g426
-sg233
+I01
+sg232
+g427
+sg234
 I00
-sg234
+sg235
 I3
-sg235
-g427
-sg237
-F0.0
-sg238
+sg236
 g428
-sg240
+sg238
+F0.0
+sg239
 g429
-sg242
-g243
-sg244
-I01
+sg241
+g430
+sg243
+g244
 sg245
+I01
+sg246
 I2
-sg246
-g425
 sg247
+g426
+sg248
 I5
-stp433
-Rp434
+stp434
+Rp435
 ag0
-(g226
+(g227
 g5
-(dp435
-g228
-g254
-sg230
-I01
+(dp436
+g229
+g255
 sg231
-g255
-sg233
 I01
+sg232
+g256
 sg234
+I01
+sg235
 I8
-sg235
-g256
-sg237
-F0.0
-sg238
+sg236
 g257
-sg240
+sg238
+F0.0
+sg239
 g258
-sg242
+sg241
 g259
-sg244
+sg243
+g260
+sg245
 I01
-sg245
+sg246
 I3
-sg246
-g254
 sg247
+g255
+sg248
 I21
-stp436
-Rp437
-atp438
-Rp439
-sg264
-(lp440
-g266
-asg234
+stp437
+Rp438
+atp439
+Rp440
+sg265
 (lp441
-g0
-(g268
-g2
+g267
+asg235
 (lp442
-g270
+g0
+(g269
+g2
+(lp443
+g271
 (I2
 I1
 I0
-tp443
-Rp444
-atp445
-Rp446
-(dp447
-g276
-g280
-sg281
-(((g282
+tp444
+Rp445
+atp446
+Rp447
+(dp448
+g277
+g281
+sg282
+(((g283
 (I1
 I2
 I3
-tp448
 tp449
 tp450
-(dp451
+tp451
+(dp452
 I0
 I3
 sI1
 I2
 sI2