Merge lp:~maddevelopers/mg5amcnlo/simple_cut_via_pdgcode into lp:~mg5core2/mg5amcnlo/2.6.1

[Olivier Mattelaer]

Add some basic cut (pt and invariant mass) specified by PDG code.
This is introduced both at NLO and LO.
(pseudo rapidity cut are also available at LO)

Some restriction apply to the particle that can be specified.
- No particle that can be already constraint by standard cut.
- No massless particle for NLO run

The cut is not just added in cuts.f but the phase-space integrator is also aware of such cut to improve the efficiency of the phase-space integration.
Some bug have been fix in that regards in the NLO part.

[Olivier Mattelaer (olivier-mattelaer)]

Any comment on this?

I would like to release 2.6.1 quite quickly and would like to have this include.
Any objection?

If nobody speak, I will do one run of auto-review and merge it then.

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/simple_cut_via_pdgcode]

298. By Rikkert Frederix on 2017-11-17

small typos (mainly spelling) in the (N)LO/run_card.dat

299. By Rikkert Frederix on 2017-11-17

some very mild cleaning of the NLO cuts.f and setcuts.f

300. By Rikkert Frederix on 2017-11-17

small fixes (again mostly typos)

[Rikkert Frederix (frederix)]

Hi Olivier,

I started looking more carefully through this implementation (and fixing some of the English along the way).

About this in banner.py:

        if any(pdg in pdg_to_cut for pdg in [1,2,3,4,5,21,22,11,13,15]):
            # Note that this will double check in the fortran code
            raise Exception, "Can not use PDG related cuts for massless SM particles"

Can't we directly check if they are massless? Or is this not possible at this place in the code. In particular, in the 4FS the b's are massless (and in the 3FS even the charms are massless) and it should be perfectly fine for the code to work with those.

Cheers,
Rik

[lp:~maddevelopers/mg5amcnlo/simple_cut_via_pdgcode]

301. By Olivier Mattelaer on 2017-11-18

use maxjetflavor to test if light jet are massive or not --for validity check for pdgcut--. Plus secure the mxx_only_part_antipart against not valid input

[Olivier Mattelaer (olivier-mattelaer)]

Hi Rikkert,

Thanks to look at this.

> Can't we directly check if they are massless? Or is this not possible at this place in the code.

Since that routine encapsulated only the run_card information, the only information that we can (easily) use is the “maxjetflavor”
to check if you want put cut on b/c quark.

So I have changed the line to:
if any(pdg in pdg_to_cut for pdg in [21,22,11,13,15]+ range(self['maxjetflavor']+1)):

“maxjetflavor” is not a 100% save variable (the user can change it) but
1) This is an hidden parameter (so only expert are going to change it)
2) The massive check is also done at fortran level (so not really critical if someone bypass the crash here)
3) the code also check at fortran level that the flag is_a_j, is_a_lp, is_a_lm, is_a_ph are all on false)

Note that I do not want to allow to have c/b cut at LO since dedicated cut are already existing for such situation.

Cheers,

Olivier

PS: The associate test of such changes make me discover some cases where the code was not nicely protected against invalid input. So I have fix those as well.

> On Nov 17, 2017, at 17:10, Rikkert Frederix <email address hidden> wrote:
>
> Hi Olivier,
>
> I started looking more carefully through this implementation (and fixing some of the English along the way).
>
> About this in banner.py:
>
> if any(pdg in pdg_to_cut for pdg in [1,2,3,4,5,21,22,11,13,15]):
> # Note that this will double check in the fortran code
> raise Exception, "Can not use PDG related cuts for massless SM particles"
>
> Can't we directly check if they are massless? Or is this not possible at this place in the code. In particular, in the 4FS the b's are massless (and in the 3FS even the charms are massless) and it should be perfectly fine for the code to work with those.
>
> Cheers,
> Rik
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/simple_cut_via_pdgcode/+merge/331864
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/simple_cut_via_pdgcode.

[Rikkert Frederix (frederix)]

Hi Olivier,

Thanks for implementing this. It's ready to go in my opinion.

Cheers,
Rik

review: Approve

=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2017-07-28 13:11:54 +0000
+++ Template/LO/Cards/run_card.dat	2017-11-18 00:45:50 +0000
@@ -119,6 +119,8 @@
  %(ptamax)s  = ptamax    ! maximum pt for the photons
  %(ptlmax)s  = ptlmax    ! maximum pt for the charged leptons
  %(missetmax)s  = missetmax ! maximum missing Et (sum of neutrino's momenta)
+ %(pt_min_pdg)s = pt_min_pdg ! pt cut for other particles (use pdg code). Applied on particle and anti-particle
+ %(pt_max_pdg)s	= pt_max_pdg ! pt cut for other particles (syntax e.g. {6: 100, 25: 50}) 
 #*********************************************************************
 # Minimum and maximum E's (in the center of mass frame)              *
 #*********************************************************************
@@ -130,6 +132,8 @@
  %(ebmax)s   = ebmax ! maximum E for the b
  %(eamax)s   = eamax ! maximum E for the photons
  %(elmax)s   = elmax ! maximum E for the charged leptons
+ %(e_min_pdg)s = e_min_pdg ! E cut for other particles (use pdg code). Applied on particle and anti-particle
+ %(e_max_pdg)s = e_max_pdg ! E cut for other particles (syntax e.g. {6: 100, 25: 50}) 
 #*********************************************************************
 # Maximum and minimum absolute rapidity (for max, -1 means no cut)   *
 #*********************************************************************
@@ -141,6 +145,8 @@
  %(etabmin)s  = etabmin ! min rap for the b
  %(etaamin)s  = etaamin ! min rap for the photons
  %(etalmin)s  = etalmin ! main rap for the charged leptons
+ %(eta_min_pdg)s = eta_min_pdg ! rap cut for other particles (use pdg code). Applied on particle and anti-particle
+ %(eta_max_pdg)s = eta_max_pdg ! rap cut for other particles (syntax e.g. {6: 2.5, 23: 5})
 #*********************************************************************
 # Minimum and maximum DeltaR distance                                *
 #*********************************************************************
@@ -177,6 +183,9 @@
  %(mmbbmax)s  = mmbbmax ! max invariant mass of a b pair
  %(mmaamax)s  = mmaamax ! max invariant mass of gamma gamma pair
  %(mmllmax)s  = mmllmax ! max invariant mass of l+l- (same flavour) lepton pair
+ %(mxx_min_pdg)s = mxx_min_pdg ! min invariant mass of a pair of particles X/X~ (e.g. {6:250})
+ %(mxx_only_part_antipart)s = mxx_only_part_antipart ! if True the invariant mass is applied only 
+                       ! to pairs of particle/antiparticle and not to pairs of the same pdg codes.  
 #*********************************************************************
 # Minimum and maximum invariant mass for all letpons                 *
 #*********************************************************************

=== modified file 'Template/LO/Source/run.inc'
--- Template/LO/Source/run.inc	2016-08-20 06:27:26 +0000
+++ Template/LO/Source/run.inc	2017-11-18 00:45:50 +0000
@@ -70,4 +70,17 @@
 C
       integer pdgs_for_merging_cut(0:1000)
       common/TO_MERGE/pdgs_for_merging_cut
-
+c
+c
+c
+      integer pdg_cut(0:25)
+      double precision ptmin4pdg(0:25) 
+      double precision ptmax4pdg(0:25) 
+      double precision Emin4pdg(0:25) 
+      double precision Emax4pdg(0:25) 
+      double precision etamin4pdg(0:25) 
+      double precision etamax4pdg(0:25) 
+      double precision mxxmin4pdg(0:25) 
+      logical mxxpart_antipart(1:25)
+      common/TO_PDG_SPECIFIC_CUT/pdg_cut, ptmin4pdg,ptmax4pdg, Emin4pdg, Emax4pdg, etamin4pdg,
+     &etamax4pdg, mxxmin4pdg,mxxpart_antipart
\ No newline at end of file

=== modified file 'Template/LO/SubProcesses/setcuts.f'
--- Template/LO/SubProcesses/setcuts.f	2016-08-30 23:56:34 +0000
+++ Template/LO/SubProcesses/setcuts.f	2017-11-18 00:45:50 +0000
@@ -23,7 +23,7 @@
 c
 c     LOCAL
 c
-      integer i,j
+      integer i,j,k
       integer icollider,detail_level
       integer ncheck
       logical done,fopened
@@ -114,7 +114,7 @@
 c
 c     count the number of j/bjet/photon/lepton
 c
-      integer nb_j, nb_b, nb_a, nb_l, nb_nocut
+      integer nb_j, nb_b, nb_a, nb_l, nb_nocut,nb_pdg
       double precision smin_p, smin_m ! local variable to compute smin
 c
 c     setup masses for the final-state particles
@@ -204,7 +204,7 @@
          is_a_nu(i)=.false.
 
 
-c-do not apply cuts to these
+c-do not apply cuts to these. CAREFULL: PDG_CUT do not consider do_cuts (they simply check the cut_decays)
          if (pmass(i).gt.20d0)     do_cuts(i)=.false.  ! no cuts on top,W,Z,H
          if (abs(idup(i,1,iproc)).eq.12) do_cuts(i)=.false.  ! no cuts on ve ve~
          if (abs(idup(i,1,iproc)).eq.14) do_cuts(i)=.false.  ! no cuts on vm vm~
@@ -308,8 +308,27 @@
 c            endif
          endif
       enddo
-
-
+c
+c     check for pdg specific cut
+c
+      if (pdg_cut(1).ne.0)then
+
+         do j=1,pdg_cut(0)
+            do i=nincoming+1, nexternal
+               if(.not.cut_decays.and.from_decay(i))then
+                  cycle
+               endif
+               if (abs(idup(i, 1, iproc)).eq.pdg_cut(j))then
+                  etmin(i) = ptmin4pdg(j)
+                  etmax(i)= ptmax4pdg(j)
+                  emin(i)= Emin4pdg(j)
+                  emax(i)= Emax4pdg(j)
+                  etamin(i)= etamin4pdg(j)
+                  etamax(i)= etamax4pdg(j)
+               endif
+            enddo
+         enddo
+      endif
 
 c
 c     delta r cut
@@ -385,6 +404,35 @@
             endif
          enddo
       enddo      
+c
+c     smin cut from PDG cut 
+c
+      if (pdg_cut(1).ne.0)then
+         do k=1,pdg_cut(0)
+            do i=nincoming+1, nexternal
+               if(.not.cut_decays.and.from_decay(i))then
+                  cycle
+               endif
+               if (abs(idup(i, 1, iproc)).ne.pdg_cut(k))then
+                  cycle
+               endif
+               do j = i+1,nexternal
+                  if(.not.cut_decays.and.from_decay(j))then
+                     cycle
+                  endif
+                  if (mxxpart_antipart(k))then
+                     if  (idup(j, 1, iproc).eq.-1*idup(i, 1, iproc))then
+                        s_min(j,i) = mxxmin4pdg(k)**2
+                     endif
+                  else
+                     if  (abs(idup(j, 1, iproc)).eq.pdg_cut(k))then
+                        s_min(j,i) = mxxmin4pdg(k)**2
+                     endif
+                  endif
+               enddo
+            enddo
+         enddo
+      endif      
 
 c     
 c     ptll cut (min and max)
@@ -598,6 +646,20 @@
          endif
         smin = smin + max(smin_p**2, smin_m, mmnl**2, ptllmin**2, misset**2)
       endif
+c  check for other PDG particle
+      do i = 1,pdg_cut(0)
+         nb_pdg=0
+         do j=nincoming, nexternal
+            if (abs(idup(i,1,iproc)).eq.pdg_cut(i))then
+               k = j ! use to get the mass
+               if(do_cuts(i)) nb_pdg = nb_pdg + 1
+            endif
+         enddo
+         smin_p = nb_pdg *(pmass(k)**2 + ptmin4pdg(i)**2)
+         smin_m = nb_pdg*((2-nb_pdg)*pmass(k)**2 + (nb_pdg-1)/2.*mxxmin4pdg(i)**2)
+         smin = smin + max(smin_p, smin_m)
+      enddo
+
 c     ensure symmetry of s_min(i,j)
       do i=nincoming+1,nexternal-1
         do j=nincoming+1,nexternal-1

=== modified file 'Template/NLO/Cards/run_card.dat'
--- Template/NLO/Cards/run_card.dat	2017-07-28 07:44:57 +0000
+++ Template/NLO/Cards/run_card.dat	2017-11-18 00:45:50 +0000
@@ -71,7 +71,7 @@
 # WARNING: PYTHIA6PT works only for processes without FSR!!!!          *
 #***********************************************************************
   %(parton_shower)s   = parton_shower
-  %(shower_scale_factor)s       = shower_scale_factor ! multiply default shower starting
+  %(shower_scale_factor)s = shower_scale_factor ! multiply default shower starting
                                   ! scale by this factor
 #***********************************************************************
 # Renormalization and factorization scales                             *
@@ -156,6 +156,15 @@
   %(epsgamma)s  = epsgamma  ! epsilon_gamma parameter of eq.(3.4) in hep-ph/9801442
  %(isoem)s  = isoEM  ! isolate photons from EM energy (photons and leptons)
 #***********************************************************************
+# Cuts associated to MASSIVE particles identified by their PDG codes.  *
+# All cuts are applied to both particles and anti-particles, so use    *
+# POSITIVE PDG CODES only. Example of the syntax is {6 : 100} or       *
+# {6:100, 25:200} for multiple particles                               *
+#***********************************************************************
+  %(pt_min_pdg)s = pt_min_pdg ! Min pT for a massive particle
+  %(pt_max_pdg)s = pt_max_pdg ! Max pT for a massive particle
+  %(mxx_min_pdg)s = mxx_min_pdg ! inv. mass for any pair of (anti)particles
+#***********************************************************************
 # For aMCfast+APPLGRID use in PDF fitting (http://amcfast.hepforge.org)*
 #***********************************************************************
  %(iappl)s = iappl ! aMCfast switch (0=OFF, 1=prepare grids, 2=fill grids)

=== modified file 'Template/NLO/Source/run.inc'
--- Template/NLO/Source/run.inc	2017-06-19 17:15:00 +0000
+++ Template/NLO/Source/run.inc	2017-11-18 00:45:50 +0000
@@ -78,4 +78,13 @@
 c    For FO run (with lhe type of analysis
 c
       double precision FO_LHE_weight_ratio
-      common /FO_ANALYSIS_LHW/FO_LHE_weight_ratio
\ No newline at end of file
+      common /FO_ANALYSIS_LHW/FO_LHE_weight_ratio
+c
+c
+c
+      integer pdg_cut(0:25)
+      double precision ptmin4pdg(0:25) 
+      double precision ptmax4pdg(0:25) 
+      double precision mxxmin4pdg(0:25) 
+      logical mxxpart_antipart(1:25) 
+      common/TO_PDG_SPECIFIC_CUT/pdg_cut, ptmin4pdg,ptmax4pdg, mxxmin4pdg, mxxpart_antipart
\ No newline at end of file

=== modified file 'Template/NLO/SubProcesses/cuts.f'
--- Template/NLO/SubProcesses/cuts.f	2017-08-08 11:31:26 +0000
+++ Template/NLO/SubProcesses/cuts.f	2017-11-18 00:45:50 +0000
@@ -40,6 +40,8 @@
       external R2_04,invm2_04,pt_04,eta_04,pt,eta
 c local integers
       integer i,j
+c temporary variable for caching locally computation
+      double precision tmpvar
 c jet cluster algorithm
       integer nQCD,NJET,JET(nexternal)
       double precision pQCD(0:3,nexternal),PJET(0:3,nexternal)
@@ -62,6 +64,11 @@
       double precision p_unlops(0:3,nexternal)
       include "run.inc" ! includes the ickkw parameter
       logical passUNLOPScuts
+c PDG specific cut
+      double precision etmin(nincoming+1:nexternal-1)
+      double precision etmax(nincoming+1:nexternal-1)
+      double precision mxxmin(nincoming+1:nexternal-1,nincoming+1:nexternal-1)
+      common /to_cuts/etmin,etmax,mxxmin
 c logicals that define if particles are leptons, jets or photons. These
 c are filled from the PDG codes (iPDG array) in this function.
       logical is_a_lp(nexternal),is_a_lm(nexternal),is_a_j(nexternal)
@@ -372,6 +379,31 @@
 c End photon isolation
       endif
 
+C
+C     PDG SPECIFIC CUTS (PT/M_IJ)
+C
+      do i=nincoming+1,nexternal-1
+         if(etmin(i).gt.0d0 .or. etmax(i).gt.0d0)then
+            tmpvar = pt_04(p(0,i))
+            if (tmpvar.lt.etmin(i)) then
+               passcuts_user=.false.
+               return
+            elseif (tmpvar.gt.etmax(i) .and. etmax(i).gt.0d0) then
+               passcuts_user=.false.
+               return
+            endif
+         endif
+         do j=i+1, nexternal-1
+            if (mxxmin(i,j).gt.0d0)then
+               if (invm2_04(p(0,i),p(0,j),1d0).lt.mxxmin(i,j)**2)then
+                  passcuts_user=.false.
+                  return
+               endif
+            endif
+         enddo
+      enddo
+
+
 C***************************************************************
 C***************************************************************
 C PUT HERE YOUR USER-DEFINED CUTS
@@ -379,6 +411,7 @@
 C***************************************************************
 C
 c$$$C EXAMPLE: cut on top quark pT
+c$$$C          Note that PDG specific cut are more optimised than simple user cut
 c$$$      do i=1,nexternal   ! loop over all external particles
 c$$$         if (istatus(i).eq.1    ! final state particle
 c$$$     &        .and. abs(ipdg(i)).eq.6) then    ! top quark

=== modified file 'Template/NLO/SubProcesses/setcuts.f'
--- Template/NLO/SubProcesses/setcuts.f	2017-06-16 08:20:04 +0000
+++ Template/NLO/SubProcesses/setcuts.f	2017-11-18 00:45:50 +0000
@@ -19,7 +19,7 @@
 c
 c     LOCAL
 c
-      integer i,j
+      integer i,j,k
 C     
 C     GLOBAL
 C
@@ -36,6 +36,11 @@
       LOGICAL  IS_A_J(NEXTERNAL),IS_A_LP(NEXTERNAL),IS_A_LM(NEXTERNAL)
       LOGICAL  IS_A_PH(NEXTERNAL)
       COMMON /TO_SPECISA/IS_A_J,IS_A_LP,IS_A_LM,IS_A_PH
+c 
+      double precision etmin(nincoming+1:nexternal-1)
+      double precision etmax(nincoming+1:nexternal-1)
+      double precision mxxmin(nincoming+1:nexternal-1,nincoming+1:nexternal-1)
+      common /to_cuts/etmin,etmax, mxxmin
 c
 c     setup masses for the final-state particles (fills the /to_mass/ common block)
 c
@@ -87,6 +92,56 @@
 c-photons
          if (idup(i,1).eq.22)  is_a_ph(i)=.true. !  photon
       enddo
+c
+c     check for pdg specific cut (pt/eta)
+c
+      do i=nincoming+1, nexternal-1 ! never include last particle
+         etmin(i) = 0d0
+         etmax(i) = -1d0
+         do j = i, nexternal-1
+            mxxmin(i,j) = 0d0
+         enddo
+      enddo
+
+      if (pdg_cut(1).ne.0)then
+         do j=1,pdg_cut(0)
+            do i=nincoming+1, nexternal-1 ! never include last particle
+               if (abs(idup(i,1)).ne.pdg_cut(j)) cycle
+c     fully ensure that only massive particles are allowed at NLO
+               if(pmass(i).eq.0d0) then
+                  write(*,*) 'Illegal use of pdg specific cut.'
+                  write(*,*) 'For NLO process, '/
+     $                 /'only massive particle can be included'
+                  stop 1
+               endif
+c     fully ensure that this is not a jet/lepton/photon
+               if(is_a_lp(i) .or. is_a_lm(i) .or. is_a_j(i) .or.
+     $              is_a_ph(i)) then
+                  write(*,*) 'Illegal use of pdg specific cut.'
+                  write(*,*) 'This can not be used for '/
+     $                 /'jet/lepton/photon/gluon'
+                  stop 1
+               endif
+               etmin(i) = ptmin4pdg(j)
+               etmax(i) = ptmax4pdg(j)
+!     add the invariant mass cut
+               if(mxxmin4pdg(j).ne.0d0)then
+                  do k=i+1, nexternal-1
+                     if (mxxpart_antipart(j))then
+                        if (idup(k, 1).eq.-1*idup(i,1))then
+                           mxxmin(i,k) = mxxmin4pdg(j)
+                        endif
+                     else
+                        if (abs(idup(k, 1)).eq.pdg_cut(j))then
+                           mxxmin(i,k) = mxxmin4pdg(j)
+                        endif
+                     endif
+                  enddo
+               endif
+            enddo
+         enddo
+      endif
+
 
       RETURN
       END
@@ -157,6 +212,14 @@
       common /c_leshouche_inc/idup,mothup,icolup,niprocs
       real*8         emass(nexternal)
       common/to_mass/emass
+c     block for the (simple) cut bsed on the pdg
+      double precision etmin(nincoming+1:nexternal-1)
+      double precision etmax(nincoming+1:nexternal-1)
+      double precision mxxmin(nincoming+1:nexternal-1,nincoming+1:nexternal-1)
+      common /to_cuts/etmin,etmax,mxxmin
+      double precision smin_update , mxx
+      integer nb_iden_pdg
+c
       logical firsttime,firsttime_chans(maxchannels)
       data firsttime /.true./
       data firsttime_chans/maxchannels*.true./
@@ -201,13 +264,13 @@
 c Add the minimal jet pTs to tau
                if(IS_A_J(i) .and. i.ne.nexternal)then
                   if  (j_fks.gt.nincoming .and. j_fks.lt.nexternal) then
-                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+max(ptj,emass(i))
-                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+max(ptj,emass(i))
-                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+max(ptj,emass(i))
+                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
+                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
+                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
                   elseif (j_fks.ge.1 .and. j_fks.le.nincoming) then
                      taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
-                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+max(ptj,emass(i))
-                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+max(ptj,emass(i))
+                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
+                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
                   elseif (j_fks.eq.nexternal) then
                      write (*,*)
      &                    'ERROR, j_fks cannot be the final parton'
@@ -234,11 +297,13 @@
                   xm(i)=emass(i)+ptgmin
                elseif (is_a_lp(i)) then
 c Add the postively charged lepton pTs to tau
-                  taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
-                  if (j_fks.gt.nincoming)
-     &                 taumin(iFKS,ichan)=taumin(iFKS,ichan)+ptl
-                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+emass(i)+ptl
-                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+emass(i)+ptl
+                  if (j_fks.gt.nincoming) then
+                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+dsqrt(ptl**2+emass(i)**2)
+                  else
+                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
+                  endif
+                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+dsqrt(emass(i)**2+ptl**2)
+                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+dsqrt(emass(i)**2+ptl**2)
                   xm(i)=emass(i)+ptl
 c Add the lepton invariant mass to tau if there is at least another
 c lepton of opposite charge. (Only add half of it, i.e. 'the part
@@ -248,23 +313,23 @@
                      if (is_a_lm(j) .and. idup(i,1).eq.-idup(j,1) .and.
      $                    (mll_sf.ne.0d0 .or. mll.ne.0d0) ) then
                         if (j_fks.gt.nincoming)
-     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-ptl-emass(i) +
-     &                              max(mll/2d0,mll_sf/2d0,ptl+emass(i))
-                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-ptl-emass(i)
-     $                       + max(mll/2d0,mll_sf/2d0,ptl+emass(i))
-                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-ptl-emass(i)
-     $                       + max(mll/2d0,mll_sf/2d0,ptl+emass(i))
+     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2) +
+     &                              max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
+                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
+     $                       + max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
+                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
+     $                       + max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
                         xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,mll_sf/2d0
      $                       ,ptl+emass(i))
                         exit
                      elseif (is_a_lm(j) .and. mll.ne.0d0) then
                         if (j_fks.gt.nincoming)
-     &                       taumin(iFKS,ichan)= taumin(iFKS,ichan)-ptl-emass(i) +
-     &                                     max(mll/2d0,ptl+emass(i))
-                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-ptl-emass(i)
-     $                       + max(mll/2d0,ptl+emass(i))
-                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-ptl-emass(i)
-     $                       + max(mll/2d0,ptl+emass(i))
+     &                       taumin(iFKS,ichan)= taumin(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2) +
+     &                                     max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
+                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
+     $                       + max(mll/2d0, dsqrt(ptl**2+emass(i)**2))
+                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
+     $                       + max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
                         xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,ptl
      $                       +emass(i))
                         exit
@@ -272,11 +337,13 @@
                   enddo
                elseif (is_a_lm(i)) then
 c Add the negatively charged lepton pTs to tau
-                  taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
-                  if (j_fks.gt.nincoming)
-     &                 taumin(iFKS,ichan)=taumin(iFKS,ichan)+ptl
-                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+emass(i)+ptl
-                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+emass(i)+ptl
+                  if (j_fks.gt.nincoming) then
+                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+dsqrt(ptl**2+emass(i)**2)
+                  else
+                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
+                  endif
+                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+dsqrt(ptl**2+emass(i)**2)
+                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+dsqrt(ptl**2+emass(i)**2)
                   xm(i)=emass(i)+ptl
 c Add the lepton invariant mass to tau if there is at least another
 c lepton of opposite charge. (Only add half of it, i.e. 'the part
@@ -286,33 +353,65 @@
                      if (is_a_lp(j) .and. idup(i,1).eq.-idup(j,1) .and.
      $                    (mll_sf.ne.0d0 .or. mll.ne.0d0) ) then
                         if (j_fks.gt.nincoming)
-     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-ptl-emass(i) +
-     &                              max(mll/2d0,mll_sf/2d0,ptl+emass(i))
-                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-ptl-emass(i)
-     $                       + max(mll/2d0,mll_sf/2d0,ptl+emass(i))
-                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-ptl-emass(i)
-     $                       + max(mll/2d0,mll_sf/2d0,ptl+emass(i))
+     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2) +
+     &                              max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
+                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
+     $                       + max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
+                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
+     $                       + max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
                         xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,mll_sf/2d0
      $                       ,ptl+emass(i))
                         exit
                      elseif (is_a_lp(j) .and. mll.ne.0d0) then
                         if (j_fks.gt.nincoming)
-     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-ptl-emass(i) +
-     &                                      max(mll/2d0,ptl+emass(i))
-                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-ptl-emass(i)
-     $                       + max(mll/2d0,ptl+emass(i))
-                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-ptl-emass(i)
-     $                       + max(mll/2d0,ptl+emass(i))
+     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2) +
+     &                                      max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
+                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
+     $                       + max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
+                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
+     $                       + max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
                         xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,ptl
      $                       +emass(i))
                         exit
                      endif
                   enddo
                else
-                  taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
-                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+emass(i)
-                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+emass(i)
-                  xm(i)=emass(i)
+                  if (i.eq.nexternal)then
+                        taumin(iFKS,ichan)=taumin(iFKS,ichan) + emass(i)
+                        taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan) +  emass(i)
+                        taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan) + emass(i)
+                        xm(i) = emass(i)
+                  else
+                     smin_update = 0
+                     nb_iden_pdg = 1
+                     mxx = 0d0
+c                    assume smin apply always on the same set of particle
+                     do j=nincoming+1,nexternal-1
+                        if (mxxmin(i,j).ne.0d0.or.mxxmin(j,i).ne.0d0) then
+                           nb_iden_pdg = nb_iden_pdg +1
+                           if (mxx.eq.0d0) mxx = max(mxxmin(i,j), mxxmin(j,i))
+                        endif
+                     enddo
+                     ! S >= (2*N-N^2)*M1^2 + (N^2-N)/2 * Mxx^2
+                     smin_update = nb_iden_pdg*((2-nb_iden_pdg)*emass(i)**2 + (nb_iden_pdg-1)/2.*mxx**2)
+                     ! compare with the update from pt cut
+                     if (smin_update.lt.nb_iden_pdg**2*(etmin(i)**2 + emass(i)**2))then
+                        ! the pt is more restrictive
+                        smin_update = dsqrt(etmin(i)**2 + emass(i)**2)
+                     else
+                        smin_update = dsqrt(smin_update)/nb_iden_pdg ! share over N particle, and change dimension
+                     endif
+                     smin_update = emass(i)
+                     ! update in sqrt(s) so take the 
+                     if  (j_fks.gt.nincoming) then
+                        taumin(iFKS,ichan)=taumin(iFKS,ichan) + smin_update
+                     else
+                        taumin(iFKS,ichan)=taumin(iFKS,ichan) + emass(i)
+                     endif
+                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan) + smin_update
+                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan) + smin_update
+                     xm(i) = smin_update
+                  endif
                endif
                xw(i)=0d0
             enddo

=== modified file 'madgraph/interface/common_run_interface.py'
--- madgraph/interface/common_run_interface.py	2017-11-17 12:35:40 +0000
+++ madgraph/interface/common_run_interface.py	2017-11-18 00:45:50 +0000
@@ -5169,8 +5169,11 @@
         self.mw_card[block][name] = value
     
     def setR(self, name, value):
-        logger.info('modify parameter %s of the run_card.dat to %s' % (name, value),'$MG:color:BLACK')
+
         self.run_card.set(name, value, user=True)
+        new_value = self.run_card.get(name)
+        logger.info('modify parameter %s of the run_card.dat to %s' % (name, new_value),'$MG:color:BLACK')        
+
 
     def setML(self, name, value, default=False):
         
@@ -5234,6 +5237,16 @@
         """This is run on quitting the class. Apply here all the self-consistency
         rule that you want. Do the modification via the set command."""
 
+        # For NLO run forbid any pdg specific cut on massless particle
+        if isinstance(self.run_card,banner_mod.RunCardNLO):
+            for pdg in set(self.run_card['pt_min_pdg'].keys()+self.run_card['pt_max_pdg'].keys()+
+                           self.run_card['mxx_min_pdg'].keys()): 
+            
+                if int(pdg)<0:
+                    raise Exception, "For PDG specific cuts, always use positive PDG codes: the cuts are applied to both particles and anti-particles"
+                if self.param_card.get_value('mass', int(pdg), default=0) ==0:
+                    raise Exception, "For NLO runs, you can use PDG specific cuts only for massive particles: (failed for %s)" % pdg
+        
         # if NLO reweighting is ON: ensure that we keep the rwgt information
         if 'reweight' in self.allow_arg and 'run' in self.allow_arg and \
             isinstance(self.run_card,banner_mod.RunCardNLO) and \

=== modified file 'madgraph/various/banner.py'
--- madgraph/various/banner.py	2017-10-24 19:40:17 +0000
+++ madgraph/various/banner.py	2017-11-18 00:45:50 +0000
@@ -2328,8 +2328,17 @@
         for name, value in self.system_default.items():
                 self.set(name, value, changeifuserset=False)
 
+        for name in self.legacy_parameter:
+            if self[name] != self.legacy_parameter[name]:
+                logger.warning("The parameter %s is not supported anymore this parameter will be ignored." % name)
+               
     default_include_file = 'run_card.inc'
 
+    def update_system_parameter_for_include(self):
+        """update hidden system only parameter for the correct writtin in the 
+        include"""
+        return
+
     def write_include_file(self, output_dir):
         """Write the various include file in output_dir.
         The entry True of self.includepath will be written in run_card.inc
@@ -2338,6 +2347,9 @@
         # ensure that all parameter are coherent and fix those if needed
         self.check_validity()
         
+        #ensusre that system only parameter are correctly set
+        self.update_system_parameter_for_include()
+        
         for incname in self.includepath:
             if incname is True:
                 pathinc = self.default_include_file
@@ -2632,8 +2644,35 @@
         self.add_param('survey_nchannel_per_job', 2, hidden=True, include=False)
         self.add_param('refine_evt_by_job', -1, hidden=True, include=False)
         
-        # Specify what particle IDs to use for the CKKWL merging cut ktdurham
+        # parameter allowing to define simple cut via the pdg
+        # Special syntax are related to those. (can not be edit directly)
+        self.add_param('pt_min_pdg',{'__type__':0.}, include=False)
+        self.add_param('pt_max_pdg',{'__type__':0.}, include=False)
+        self.add_param('E_min_pdg',{'__type__':0.}, include=False)
+        self.add_param('E_max_pdg',{'__type__':0.}, include=False)
+        self.add_param('eta_min_pdg',{'__type__':0.}, include=False)
+        self.add_param('eta_max_pdg',{'__type__':0.}, include=False)
+        self.add_param('mxx_min_pdg',{'__type__':0.}, include=False)
+        self.add_param('mxx_only_part_antipart', {'default':False}, include=False, hidden=True)
         
+        self.add_param('pdg_cut',[0], hidden=True, system=True) # store which PDG are tracked
+        self.add_param('ptmin4pdg',[0.], hidden=True, system=True) # store pt min
+        self.add_param('ptmax4pdg',[-1.], hidden=True, system=True)
+        self.add_param('Emin4pdg',[0.], hidden=True, system=True) # store pt min
+        self.add_param('Emax4pdg',[-1.], hidden=True, system=True)  
+        self.add_param('etamin4pdg',[0.], hidden=True, system=True) # store pt min
+        self.add_param('etamax4pdg',[-1.], hidden=True, system=True)   
+        self.add_param('mxxmin4pdg',[-1.], hidden=True, system=True)
+        self.add_param('mxxpart_antipart', [False], hidden=True, system=True)
+        # Not implemetented right now (double particle cut)
+        #self.add_param('pdg_cut_2',[0], hidden=True, system=True)
+        # self.add_param('M_min_pdg',[0.], hidden=True, system=True) # store pt min
+        #self.add_param('M_max_pdg',[0.], hidden=True, system=True)               
+        # self.add_param('DR_min_pdg',[0.], hidden=True, system=True) # store pt min
+        #self.add_param('DR_max_pdg',[0.], hidden=True, system=True)               
+            
+            
+             
     def check_validity(self):
         """ """
         
@@ -2709,7 +2748,7 @@
                     self['mmjj'] = 0.0 
 
 
-
+    
         # check validity of the pdf set
         possible_set = ['lhapdf', 'mrs02nl','mrs02nn',
         'cteq4_m', 'cteq4_l','cteq4_d',
@@ -2724,13 +2763,68 @@
             #add warning if lhaid not define
             self.get_default('lhaid', log_level=20)
    
-        for name in self.legacy_parameter:
-            if self[name] != self.legacy_parameter[name]:
-                logger.warning("The parameter %s is not supported anymore this parameter will be ignored." % name)
-                
-
+    def update_system_parameter_for_include(self):
+        
+        # set the pdg_for_cut fortran parameter
+        pdg_to_cut = set(self['pt_min_pdg'].keys() +self['pt_max_pdg'].keys() + 
+                         self['e_min_pdg'].keys() +self['e_max_pdg'].keys() +
+                         self['eta_min_pdg'].keys() +self['eta_max_pdg'].keys()+
+                         self['mxx_min_pdg'].keys() + self['mxx_only_part_antipart'].keys())
+        pdg_to_cut.discard('__type__')
+        pdg_to_cut.discard('default')
+        if len(pdg_to_cut)>25:
+            raise Exception, "Maximum 25 different pdgs are allowed for pdg specific cut"
+        
+        if any(int(pdg)<0 for pdg in pdg_to_cut):
+            logger.warning('PDG specific cuts are always applied symmetrically on particle/anti-particle. Always use positve PDG codes')
+            raise MadGraph5Error, 'Some PDG specific cuts are defined with negative pdg code'
+        
+        
+        if any(pdg in pdg_to_cut for pdg in [1,2,3,4,5,21,22,11,13,15]):
+            raise Exception, "Can not use PDG related cut for light quark/b quark/lepton/gluon/photon"
+        
+        if pdg_to_cut:
+            self['pdg_cut'] = list(pdg_to_cut)
+            self['ptmin4pdg'] = []
+            self['Emin4pdg'] = []
+            self['etamin4pdg'] =[]
+            self['ptmax4pdg'] = []
+            self['Emax4pdg'] = []
+            self['etamax4pdg'] =[]
+            self['mxxmin4pdg'] =[]
+            self['mxxpart_antipart']  = []
+            for pdg in self['pdg_cut']:
+                for var in ['pt','e','eta', 'Mxx']:
+                    for minmax in ['min', 'max']:
+                        if var in ['Mxx'] and minmax =='max':
+                            continue
+                        new_var = '%s%s4pdg' % (var, minmax)
+                        old_var = '%s_%s_pdg' % (var, minmax)
+                        default = 0. if minmax=='min' else -1.
+                        self[new_var].append(self[old_var][str(pdg)] if str(pdg) in self[old_var] else default)
+                #special for mxx_part_antipart
+                old_var = 'mxx_only_part_antipart'
+                new_var = 'mxxpart_antipart'
+                if 'default' in self[old_var]:
+                    default = self[old_var]['default']
+                    self[new_var].append(self[old_var][str(pdg)] if str(pdg) in self[old_var] else default)
+                else:
+                    if str(pdg) not in self[old_var]:
+                        raise Exception("no default value defined for %s and no value defined for pdg %s" % (old_var, pdg)) 
+                    self[new_var].append(self[old_var][str(pdg)])
+        else:
+            self['pdg_cut'] = [0]
+            self['ptmin4pdg'] = [0.]
+            self['Emin4pdg'] = [0.]
+            self['etamin4pdg'] =[0.]
+            self['ptmax4pdg'] = [-1.]
+            self['Emax4pdg'] = [-1.]
+            self['etamax4pdg'] =[-1.]
+            self['mxxmin4pdg'] =[0.] 
+            self['mxx_part_antipart'] = [False]
             
-        
+                    
+           
     def create_default_for_process(self, proc_characteristic, history, proc_def):
         """Rules
           process 1->N all cut set on off.
@@ -3277,7 +3371,7 @@
 
 class RunCardNLO(RunCard):
     """A class object for the run_card for a (aMC@)NLO pocess"""
-
+    
     def default_setup(self):
         """define the default value"""
         
@@ -3361,6 +3455,20 @@
         self.add_param('FO_LHE_postprocessing',['grouping','random'], 
                        hidden=True, system=True, include=False)
     
+        # parameter allowing to define simple cut via the pdg
+        self.add_param('g',{'__type__':0.}, include=False)
+        self.add_param('pt_min_pdg',{'__type__':0.}, include=False)
+        self.add_param('pt_max_pdg',{'__type__':0.}, include=False)
+        self.add_param('mxx_min_pdg',{'__type__':0.}, include=False)
+        self.add_param('mxx_only_part_antipart', {'default':False}, include=False, hidden=True)
+        
+        #hidden parameter that are transfer to the fortran code
+        self.add_param('pdg_cut',[0], hidden=True, system=True) # store which PDG are tracked
+        self.add_param('ptmin4pdg',[0.], hidden=True, system=True) # store pt min
+        self.add_param('ptmax4pdg',[-1.], hidden=True, system=True)
+        self.add_param('mxxmin4pdg',[0.], hidden=True, system=True)
+        self.add_param('mxxpart_antipart', [False], hidden=True, system=True)
+        
     def check_validity(self):
         """check the validity of the various input"""
         
@@ -3505,6 +3613,56 @@
                 raise InvalidRunCard, "'rw_fscale' has two or more identical entries. They have to be all different for the code to work correctly."
 
 
+    def update_system_parameter_for_include(self):
+        
+        # set the pdg_for_cut fortran parameter
+        pdg_to_cut = set(self['pt_min_pdg'].keys() +self['pt_max_pdg'].keys()+
+                         self['mxx_min_pdg'].keys()+ self['mxx_only_part_antipart'].keys())
+        pdg_to_cut.discard('__type__')
+        pdg_to_cut.discard('default')
+        if len(pdg_to_cut)>25:
+            raise Exception, "Maximum 25 different PDGs are allowed for PDG specific cut"
+        
+        if any(int(pdg)<0 for pdg in pdg_to_cut):
+            logger.warning('PDG specific cuts are always applied symmetrically on particle/anti-particle. Always use positve PDG codes')
+            raise MadGraph5Error, 'Some PDG specific cuts are defined with negative PDG codes'
+        
+        
+        if any(pdg in pdg_to_cut for pdg in [21,22,11,13,15]+ range(self['maxjetflavor']+1)):
+            # Note that this will double check in the fortran code
+            raise Exception, "Can not use PDG related cuts for massless SM particles/leptons"
+        if pdg_to_cut:
+            self['pdg_cut'] = list(pdg_to_cut)
+            self['ptmin4pdg'] = []
+            self['ptmax4pdg'] = []
+            self['mxxmin4pdg'] = []
+            self['mxxpart_antipart']  = []
+            for pdg in self['pdg_cut']:
+                for var in ['pt','mxx']:
+                    for minmax in ['min', 'max']:
+                        if var == 'mxx' and minmax == 'max':
+                            continue
+                        new_var = '%s%s4pdg' % (var, minmax)
+                        old_var = '%s_%s_pdg' % (var, minmax)
+                        default = 0. if minmax=='min' else -1.
+                        self[new_var].append(self[old_var][str(pdg)] if str(pdg) in self[old_var] else default)
+                #special for mxx_part_antipart
+                old_var = 'mxx_only_part_antipart'
+                new_var = 'mxxpart_antipart'
+                if 'default' in self[old_var]:
+                    default = self[old_var]['default']
+                    self[new_var].append(self[old_var][str(pdg)] if str(pdg) in self[old_var] else default)
+                else:
+                    if str(pdg) not in self[old_var]:
+                        raise Exception("no default value defined for %s and no value defined for pdg %s" % (old_var, pdg)) 
+                    self[new_var].append(self[old_var][str(pdg)])
+        else:
+            self['pdg_cut'] = [0]
+            self['ptmin4pdg'] = [0.]
+            self['ptmax4pdg'] = [-1.]
+            self['mxxmin4pdg'] = [0.]
+            self['mxx_part_antipart'] = [False]
+
     def write(self, output_file, template=None, python_template=False):
         """Write the run_card in output_file according to template 
            (a path to a valid run_card)"""