Merge lp:~maddevelopers/mg5amcnlo/3.0.2-alpha0 into lp:~maddevelopers/mg5amcnlo/3.2.1

[marco zaro]

This merge includes the possibility to use tagged particles in the final state, which is particularly relevant for EW corrections with tagged photons (see 2106.02059).
Note that cuts.f has been deeply refactored.

[lp:~maddevelopers/mg5amcnlo/3.0.2-alpha0]

1000. By marco zaro on 2021-09-16

removed process folders

[Rikkert Frederix (frederix)]

Hi Marco,

I'll have a closer look next week. But for now:

- There is a text conflict in the UpdateNotes.txt
- Why do you add the new models? Shouldn't they part of FeynRules webpage so that they are simply downloaded when needed?

Best,
Rikkert

[lp:~maddevelopers/mg5amcnlo/3.0.2-alpha0]

1001. By marco zaro on 2021-09-16

merged with 3.2.1

[marco zaro (marco-zaro)]

Hi Rikkert,
thanks.

> On 16 Sep 2021, at 16:22, Rikkert Frederix <email address hidden> wrote:
>
> Hi Marco,
>
> I'll have a closer look next week. But for now:
>
> - There is a text conflict in the UpdateNotes.txt
I solved it, putting the changes for 3.1.2 together with 3.2.1, since the former was never released
> - Why do you add the new models? Shouldn't they part of FeynRules webpage so that they are simply downloaded when needed?
I will discuss with Olivier and put the models on the web repository

Cheers,

Marco
>
> Best,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2-alpha0/+merge/408723
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.2.1.
>

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

I have given a quick look and this sounds very good.

- I agree with Rikkert that it is better to put the model online (maybe not on the Feynrules website but on our madgraph website) the branch lp:~maddevelopers/mg5amcnlo/UFOmodel_db is designed for that actually.

- I'm worried about the re-introduction of the fortran function SORTTI and SORTTC.
If I remember correctly, such function had some issue with GCC10, did you check.

I do not see in the diff any kind of test (maybe it is just truncated). Can you add at least on parralel test such that I can run it on various GCC version and check that they are no issue on that.

Did you add unittest for the new syntax? if not I will likely add those when running some tests.

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/3.0.2-alpha0]

1002. By marco zaro on 2021-09-20

fix in cuts.f to remove sortti/sorttc; unit+acceptance test added; models removed and added to the UFOmodel branch

[marco zaro (marco-zaro)]

Hello both
thanks for the reviews.
I have removed the models rom the branch and added to the UFOmodel one. I should have taken care of all the other comments in revision 1002.

Cheers,

Marco

> On 17 Sep 2021, at 09:54, Olivier Mattelaer <email address hidden> wrote:
>
> Hi,
>
> I have given a quick look and this sounds very good.
>
> - I agree with Rikkert that it is better to put the model online (maybe not on the Feynrules website but on our madgraph website) the branch lp:~maddevelopers/mg5amcnlo/UFOmodel_db is designed for that actually.
>
> - I'm worried about the re-introduction of the fortran function SORTTI and SORTTC.
> If I remember correctly, such function had some issue with GCC10, did you check.
>
> I do not see in the diff any kind of test (maybe it is just truncated). Can you add at least on parralel test such that I can run it on various GCC version and check that they are no issue on that.
>
> Did you add unittest for the new syntax? if not I will likely add those when running some tests.
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2-alpha0/+merge/408723
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.2.1.
>

[Rikkert Frederix (frederix)]

Hi Marco,

I've done the following test with this branch:

import model loop_qcd_qed_sm_Gmu-a0
generate p p > !a! !a! [real=QCD] ; output ; launch -p -f

and it crashes stating that there are "no integration channels". Looking in the relevant gensym.log, there is "ERROR, too many photons 2". I thought that the above process should have worked, no? Have I misunderstood something here?

Best,
Rikkert

[lp:~maddevelopers/mg5amcnlo/3.0.2-alpha0]

1003. By marco zaro on 2021-10-04

cuts.f updated so that recombine_momenta won't crash for too many photons (thanks Ioannis Tsinikos); the gmu<>alpha scheme conversion is written dinamically, to ensure compatibility with processes without tagged photons

[marco zaro (marco-zaro)]

Hi Rikkert,
Indeed, Ioannis had produces a modified version of cuts.f that makes disappear your error. I have updated cuts.f in the branch accordingly, and fixed some other things in the code.
Cheers,

Marco

> On 21 Sep 2021, at 09:56, Rikkert Frederix <email address hidden> wrote:
>
> Hi Marco,
>
> I've done the following test with this branch:
>
> import model loop_qcd_qed_sm_Gmu-a0
> generate p p > !a! !a! [real=QCD] ; output ; launch -p -f
>
> and it crashes stating that there are "no integration channels". Looking in the relevant gensym.log, there is "ERROR, too many photons 2". I thought that the above process should have worked, no? Have I misunderstood something here?
>
> Best,
> Rikkert
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2-alpha0/+merge/408723
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.2.1z

[Olivier Mattelaer (olivier-mattelaer)]

This is very good and basically ready to be merged.

Some minor point.
1) !z! works ... should not be this forbidden?
2) can we use tagged photon with QCD correction only?
-- In the commit I did, I have added unittest that crash related to the above command, you can remove those test if those are fine--
3) 2!a! was making the code to crash but !2a! was working fine. I have updated the parser to support both.
4) I have spotted a typo in banner.py (see the comment within the diff for that)
5) I have loose track of all the change in cuts.f but did not spot anything obviously wrong. If you are not confident enough about that part maybe I should go deeper and track that we did not "loose" any previous cut but sounds like it.

Thanks a lot and sorry for being slow on this review

review: Approve

[lp:~maddevelopers/mg5amcnlo/3.0.2-alpha0]

1004. By olivier-mattelaer on 2021-10-14

add unittest --not all passing-- and update the parser to allow duplicate tagged particcle

1005. By marco zaro on 2021-10-15

fix in banner.py according to Olivier's comment

[marco zaro (marco-zaro)]

Hi Olivier,
thanks for your comments

> On 14 Oct 2021, at 16:26, Olivier Mattelaer <email address hidden> wrote:
>
> Review: Approve
>
> This is very good and basically ready to be merged.
>
> Some minor point.
> 1) !z! works ... should not be this forbidden?
yes, but it should do nothing
> 2) can we use tagged photon with QCD correction only?
> -- In the commit I did, I have added unittest that crash related to the above command, you can remove those test if those are fine—
In principle the syntax should work, giving identical results with the case without tagging
> 3) 2!a! was making the code to crash but !2a! was working fine. I have updated the parser to support both.
I was not aware of the fact that one can use the numbers to specify multiple times the same particle
> 4) I have spotted a typo in banner.py (see the comment within the diff for that)
I have fixed it, please have a look
> 5) I have loose track of all the change in cuts.f but did not spot anything obviously wrong. If you are not confident enough about that part maybe I should go deeper and track that we did not "loose" any previous cut but sounds like it.
>
Let us see what Rik says here…

Cheers,

Marco

> Thanks a lot and sorry for being slow on this review
>
> Diff comments:
>
>>
>> === modified file 'madgraph/various/banner.py'
>> --- madgraph/various/banner.py 2021-07-29 07:14:14 +0000
>> +++ madgraph/various/banner.py 2021-10-04 14:52:23 +0000
>> @@ -4552,11 +4553,20 @@
>> if proc_characteristic['ninitial'] == 1:
>> #remove all cut
>> self.remove_all_cut()
>> +
>> + # check for tagged photons
>> + tagged_particles = set()
>>
>> # Check if need matching
>> min_particle = 99
>> max_particle = 0
>> for proc in proc_def:
>> + for leg in proc['legs']:
>> + if leg['is_tagged']:
>> + tagged_particles.add(leg['id'])
>> +
>> + if 22 in tagged_particles:
>> + self['gamma_is_j'] = False
>> min_particle = min(len(proc['legs']), min_particle)
>> max_particle = max(len(proc['legs']), max_particle)
>
> The min_particle/max_particle should not be in the if ""22" block but still in the
> for proc in proc_def: loop
>
>> matching = False
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2-alpha0/+merge/408723
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.2.1.
>

[Rikkert Frederix (frederix)]

Hi Marco,

I think this is ready for merging.

best,
Rikkert

review: Approve

=== modified file 'Template/NLO/Cards/run_card.dat'
--- Template/NLO/Cards/run_card.dat	2021-04-05 16:07:15 +0000
+++ Template/NLO/Cards/run_card.dat	2021-10-15 11:42:16 +0000
@@ -129,7 +129,10 @@
  %(bwcutoff)s  = bwcutoff
 #***********************************************************************
 # Cuts on the jets. Jet clustering is performed by FastJet.            *
-#  - If gamma_is_j, photons are also clustered                            *
+#  - If gamma_is_j, photons are also clustered with jets.              *
+#    Otherwise, they will be treated as tagged particles and photon    * 
+#    isolation will be applied. Note that photons in the real emission *
+#    will always be clustered with QCD partons.                        *
 #  - When matching to a parton shower, these generation cuts should be *
 #    considerably softer than the analysis cuts.                       *
 #  - More specific cuts can be specified in SubProcesses/cuts.f        *
 
=== modified file 'Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f'
--- Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f	2018-10-26 09:55:25 +0000
+++ Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f	2021-10-15 11:42:16 +0000
@@ -5,8 +5,12 @@
       integer nwgt
       character*(*) weights_info(*)
       integer i,l
-      character*6 cc(2)
-      data cc/'|T@NLO','|T@LO '/
+c      character*6 cc(2)
+c      data cc/'|T@NLO','|T@LO '/
+      character*13 cc(9)
+      data cc/ ' |T@NLO      ',' |T@LO       ',' |T@LO1      '
+     $        ,' |T@LO2      ',' |T@LO3      ',' |T@NLO1     '
+     $        ,' |T@NLO2     ',' |T@NLO3     ',' |T@NLO4     '/
 
 c     Also specific perturbative orders can be directly plotted, adding for examples
 c     the following further entries in the variable data
@@ -18,8 +22,8 @@
 c     
 c     See also line 376 in this file 
       call HwU_inithist(nwgt,weights_info)
-      do i=1,2
-      l=(i-1)*55
+      do i=1,9
+      l=(i-1)*59
 c transverse momenta
       call HwU_book(l+ 1,'total rate                 '//cc(i),1,0.5d0,1.5d0)
       call HwU_book(l+ 2,'1st charged lepton log pT  '//cc(i),25,-0.2d0,3.8d0)
@@ -76,6 +80,14 @@
 c HT
       call HwU_book(l+54,'log HT (partons)           '//cc(i),25,-0.2d0,3.8d0)
       call HwU_book(l+55,'log HT (reconstructed)     '//cc(i),25,-0.2d0,3.8d0)
+
+
+      call HwU_book(l+56,'1st isolated ph log pT     '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+57,'2nd isolated ph log pT     '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+58,'3rd isolated ph log pT     '//cc(i),25,-0.2d0,3.8d0)
+
+      call HwU_book(l+59,'3 isolated ph log M        '//cc(i),25,-0.2d0,3.8d0)
+
       enddo
       return
       end
@@ -110,7 +122,7 @@
      $     ,ptmpmm,ptepve,ptmmvm,ptepemmpmm,ptepvemmvm,ptt,ptat,pttt
      $     ,etmiss,pth(2),pthh,ptv(3),ptjet(3),Mepem,Mmpmm,Mepve,Mmmvm
      $     ,Mepemmpmm,Mepvemmvm,Mtt,Mhh,Mj1j2,Mj1j3,Mj2j3,Mj1j2j3,Mvvv
-     $     ,HTparton,HTreco,p_reco(0:4,nexternal)
+     $     ,HTparton,HTreco,p_reco(0:4,nexternal),ptphiso(nexternal),Mphphph
       integer nQCD,jet(nexternal),njet,itop,iatop,iem,iep,imp,imm,ive
      $     ,ivm,iv1,iv2,iv3,ih1,ih2,il,ipdg_reco(nexternal)
       double precision getptv4,getptv4_2,getptv4_4,getinvm4_2
@@ -119,25 +131,47 @@
      $     ,getinvm4_4,l10
       integer orders_tag_plot
       common /corderstagplot/ orders_tag_plot
+
+
+c Photon isolation
+      integer nph,nem,nin,nphiso
+      double precision ptg,chi_gamma_iso,iso_getdrv40
+      double precision Etsum(0:nexternal)
+      real drlist(nexternal)
+      double precision pgamma(0:3,nexternal),pgammaiso(0:3,nexternal),pem(0:3,nexternal)
+      logical alliso,isolated
+c Sort array of results: ismode>0 for real, isway=0 for ascending order
+      integer ismode,isway,izero,isorted(nexternal)
+      parameter (ismode=1)
+      parameter (isway=0)
+      parameter (izero=0)
+      integer get_n_tagged_photons      
+
+
+      logical is_a_lp(nexternal),is_a_lm(nexternal),is_a_j(nexternal)
+     $     ,is_a_ph(nexternal)
+      REAL*8 pt,eta
+      external pt,eta,chi_gamma_iso,sortzv
+c      integer iph1,iph2,iph3
 c First, try to recombine photons with leptons      
-      if (.not.quarkphreco) then
-         write (*,*) 'quark-photon recombination is turned off. '/
-     $        /'Do need it'
-         stop
-      endif
-      if (.not. lepphreco) then
-         write (*,*) 'lepton-photon recombination is turned off. '
-     $        //'Do need it.'
-         stop
-      endif
+c      if (.not.quarkphreco) then
+c         write (*,*) 'quark-photon recombination is turned off. '/
+c     $        /'Do need it'
+c         stop
+c      endif
+c      if (.not. lepphreco) then
+c         write (*,*) 'lepton-photon recombination is turned off. '
+c     $        //'Do need it.'
+c         stop
+c      endif
       call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
      $                       p, iPDG, p_reco, iPDG_reco)
 
 c Put all (light) QCD partons(+photon) in momentum array for jet clustering.
       nQCD=0
       do j=nincoming+1,nexternal
-         if (abs(ipdg_reco(j)).le.5 .or. ipdg_reco(j).eq.21 .or.
-     $        ipdg_reco(j).eq.22)then
+         if (abs(ipdg_reco(j)).le.5 .or. ipdg_reco(j).eq.21 
+     $      .or.  (ipdg_reco(j).eq.22.and.gamma_is_j)) then
             nQCD=nQCD+1
             do i=0,3
                pQCD(i,nQCD)=p_reco(i,j)
@@ -173,6 +207,161 @@
 c
 c******************************************************************************
 
+
+c PHOTON (ISOLATION) CUTS
+c
+c find the photons
+      do i=1,nexternal
+         if (istatus(i).eq.1 .and. ipdg(i).eq.22 .and. .not.gamma_is_j) then
+            is_a_ph(i)=.true.
+         else
+            is_a_ph(i)=.false.
+         endif
+      enddo
+      if (ptgmin.ne.0d0) then
+         nph=0
+         do j=nincoming+1,nexternal
+            if (is_a_ph(j)) then
+               nph=nph+1
+               do i=0,3
+                  pgamma(i,nph)=p(i,j)
+               enddo
+            endif
+         enddo
+         if(nph.eq.0)goto 444
+c         write(*,*) 'ERROR in cuts.f: photon isolation is not working'
+c     $           // ' for mixed QED-QCD corrections'
+c         stop 1
+
+         if(isoEM)then
+            nem=nph
+            do k=1,nem
+               do i=0,3
+                  pem(i,k)=pgamma(i,k)
+               enddo
+            enddo
+            do j=nincoming+1,nexternal
+               if (is_a_lp(j).or.is_a_lm(j)) then
+                  nem=nem+1
+                  do i=0,3
+                     pem(i,nem)=p(i,j)
+                  enddo
+               endif
+            enddo
+         endif
+
+c         alliso=.true.
+         nphiso=0
+
+         j=0
+c         do while(j.lt.nph.and.alliso)
+         do while(j.lt.nph)
+
+c Loop over all photons
+            j=j+1
+
+            ptg=pt(pgamma(0,j))
+            if(ptg.lt.ptgmin)then
+               cycle
+c               return
+            endif
+            if (etagamma.gt.0d0) then
+               if (abs(eta(pgamma(0,j))).gt.etagamma) then
+                  cycle
+c                  return
+               endif
+            endif
+
+c Isolate from hadronic energy
+            do i=1,nQCD
+               drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pQCD(0,i)))
+            enddo
+            call sortzv(drlist,isorted,nQCD,ismode,isway,izero)
+            Etsum(0)=0.d0
+            nin=0
+            do i=1,nQCD
+               if(dble(drlist(isorted(i))).le.R0gamma)then
+                  nin=nin+1
+                  Etsum(nin)=Etsum(nin-1)+pt(pQCD(0,isorted(i)))
+               endif
+            enddo
+            isolated=.True.
+            do i=1,nin
+c               alliso=alliso .and.
+c     $              Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
+c     $              R0gamma,xn,epsgamma,ptg)
+                if(Etsum(i).gt.chi_gamma_iso(dble(drlist(isorted(i))),
+     $              R0gamma,xn,epsgamma,ptg)) then
+                    isolated=isolated.and..False.
+                endif
+            enddo
+            if(.not.isolated)cycle
+
+c Isolate from EM energy
+            if(isoEM.and.nem.gt.1)then
+               do i=1,nem
+                  drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i)))
+               enddo
+               call sortzv(drlist,isorted,nem,ismode,isway,izero)
+c First of list must be the photon: check this, and drop it
+               if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then
+                  write(*,*)'Error #1 in photon isolation'
+                  write(*,*)j,isorted(1),drlist(isorted(1))
+                  stop
+               endif
+               Etsum(0)=0.d0
+               nin=0
+               do i=2,nem
+                  if(dble(drlist(isorted(i))).le.R0gamma)then
+                     nin=nin+1
+                     Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i)))
+                  endif
+               enddo
+               isolated=.True.
+               do i=1,nin
+c                  alliso=alliso .and.
+c     $               Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
+c     $               R0gamma,xn,epsgamma,ptg)
+                if(Etsum(i).gt.chi_gamma_iso(dble(drlist(isorted(i))),
+     $               R0gamma,xn,epsgamma,ptg)) then
+                    isolated=isolated.and..False.
+                endif
+            enddo
+            if(.not.isolated)cycle
+            endif
+c End of loop over photons
+
+         nphiso=nphiso+1
+
+         do i=0,3
+           pgammaiso(i,nphiso)=pgamma(i,j)
+         enddo
+
+
+         enddo
+         if(nphiso.lt.get_n_tagged_photons())then
+            print*,"mismatch with cuts.f"
+            stop
+         endif
+
+
+444     continue
+c End photon isolation
+       endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 c Look for the other physics objects
       itop=0
       iatop=0
@@ -189,7 +378,14 @@
       iv3=0
       ih1=0
       ih2=0
+c      iph1=0
+c      iph2=0
+c      iph3=0
+c          print*,"nell'analisi"
       do i=1,nexternal
+c          print*,"idpg di ",i,"=",ipdg(i)
+c          print*,"idpg_reco di ",i,"=",ipdg_reco(i)
+
          if (ipdg_reco(i).eq.6) then
             itop=i
          elseif(ipdg_reco(i).eq.-6) then
@@ -232,8 +428,27 @@
                   stop
                endif
             endif
+c         elseif(abs(ipdg_reco(i)).eq.22.and.i.gt.nincoming) then
+c            if (iph1.eq.0) then
+c               iph1=i
+c            else
+c               if (iph2.eq.0) then
+c                  iph2=i
+c               else
+c                  if (iph3.eq.0) then
+c                      iph3=i
+c                  else  
+c                     write (*,*) 'too many photonss'
+c                     stop
+c                  endif
+c               endif
+c            endif
          endif
+
       enddo
+c      print*,itop,iatop
+c      stop
+
       if (itop.ne.0) ptt=getptv4(p_reco(0,itop))
       if (iatop.ne.0) ptat=getptv4(p_reco(0,iatop))
       if (itop.ne.0 .and. iatop.ne.0) then
@@ -318,11 +533,21 @@
             pth(2)=tmp
          endif
       endif
+
+
       if (iv1.ne.0) ptv(1)=getptv4(p_reco(0,iv1))
       if (iv2.ne.0) ptv(2)=getptv4(p_reco(0,iv2))
       if (iv3.ne.0) ptv(3)=getptv4(p_reco(0,iv3))
       if (iv1.ne.0 .and. iv2.ne.0 .and. iv3.ne.0) then
          Mvvv=getinvm4_3(p_reco(0,iv1),p_reco(0,iv2),p_reco(0,iv3))
+      endif
+
+      do i=1,nphiso
+         ptphiso(i)=getptv4(pgammaiso(0,i))
+      enddo
+
+      if (iv1.ne.0 .and. iv2.ne.0 .and. iv3.ne.0) then
+         Mvvv=getinvm4_3(p_reco(0,iv1),p_reco(0,iv2),p_reco(0,iv3))
 c     order the vector bosons (if there are 3)
          do i=1,2
             do j=1,3-i
@@ -341,6 +566,31 @@
             ptv(2)=tmp
          endif
       endif
+
+
+
+
+      if (nphiso.eq.3) then
+         Mphphph=getinvm4_3(pgammaiso(0,1),pgammaiso(0,2),pgammaiso(0,3))
+c     order the isolated photons (if there are 3)
+         do i=1,2
+            do j=1,3-i
+               if (ptphiso(j).lt.ptphiso(j+1)) then
+                  tmp=ptphiso(j)
+                  ptphiso(j)=ptphiso(j+1)
+                  ptphiso(j+1)=tmp
+               endif
+            enddo
+         enddo
+      elseif (nphiso.eq.2) then
+c     order the isolated photons (if there are 2)
+         if (ptphiso(1).lt.ptphiso(2)) then
+            tmp=ptphiso(1)
+            ptphiso(1)=ptphiso(2)
+            ptphiso(2)=tmp
+         endif
+      endif
+      
       do i=1,njet
          ptjet(i)=getptv4(pjet(0,i))
       enddo
@@ -369,9 +619,18 @@
       enddo
       if (ive.ne.0 .or. ivm.ne.0) HTreco=HTreco+etmiss
          
-      do i=1,2
-         l=(i-1)*55
+      do i=1,9
+         l=(i-1)*59
          if (ibody.ne.3 .and.i.eq.2) cycle
+        if (i.eq. 3.and.orders_tag_plot.ne.204) cycle
+         if (i.eq. 4.and.orders_tag_plot.ne.402) cycle
+         if (i.eq. 5.and.orders_tag_plot.ne.600) cycle
+         if (i.eq. 6.and.orders_tag_plot.ne.206) cycle
+         if (i.eq. 7.and.orders_tag_plot.ne.404) cycle
+         if (i.eq. 8.and.orders_tag_plot.ne.602) cycle
+         if (i.eq. 9.and.orders_tag_plot.ne.800) cycle
+
+
 
 c         How to tag orders (QCD+QED*100)
 c
@@ -493,6 +752,18 @@
 c HT
          call HwU_fill(l+54,l10(HTparton),wgts)
          call HwU_fill(l+55,l10(HTreco),wgts)
+
+         if (nphiso.ge.1) call HwU_fill(l+56,l10(ptphiso(1)),wgts)
+         if (nphiso.ge.2) call HwU_fill(l+57,l10(ptphiso(2)),wgts)
+         if (nphiso.ge.3) call HwU_fill(l+58,l10(ptphiso(3)),wgts)
+
+         if (nphiso.ge.3) call HwU_fill(l+59,l10(Mphphph),wgts)
+
+
+
+
+
+
       enddo
 
  999  return      
 
=== modified file 'Template/NLO/SubProcesses/BinothLHA.f'
--- Template/NLO/SubProcesses/BinothLHA.f	2020-06-26 09:54:58 +0000
+++ Template/NLO/SubProcesses/BinothLHA.f	2021-10-15 11:42:16 +0000
@@ -8,6 +8,7 @@
       implicit none
       include "nexternal.inc"
       include "coupl.inc"
+      include "../../Source/MODEL/input.inc"
       include 'born_nhel.inc'
       double precision pi, zero,mone
       parameter (pi=3.1415926535897932385d0)
@@ -57,6 +58,7 @@
       integer getordpowfromindex_ml5
       integer orders_to_amp_split_pos
       logical, allocatable, save :: keep_order(:)
+      include 'a0Gmuconv.inc'
       include 'orders.inc'
       integer amp_orders(nsplitorders)
       integer split_amp_orders(nsplitorders), iamp
@@ -74,10 +76,17 @@
 
       logical updateloop
       common /to_updateloop/updateloop
+
+      integer get_n_tagged_photons
+      external get_n_tagged_photons
+      double precision get_virtual_a0Gmu_conv
+      external get_virtual_a0Gmu_conv 
+      integer ntagph, qed_pow_b
 c masses
       include 'pmass.inc'
       data nbad / 0 /
 
+
       IOErrCounter = 0
 c update the ren_scale for MadLoop and the couplings (should be the
 c Ellis-Sexton scale)
@@ -93,8 +102,7 @@
       single  = 0d0
       double  = 0d0
       born_hel_from_virt = 0d0
-C     reset the amp_split array
-      amp_split(1:amp_split_size) = 0d0
+C     reset the various arrays
       amp_split_finite_ML(1:amp_split_size) = 0d0
       amp_split_poles_ML(1:amp_split_size,1) = 0d0
       amp_split_poles_ML(1:amp_split_size,2) = 0d0
@@ -256,6 +264,32 @@
 c      endif
 c      virt_wgt=virt_wgt+conversion*born_wgt*ao2pi
 c======================================================================
+
+c======================================================================
+c If there are tagged photon and other particles in the process, 
+C one must use a mixed Gmu-alpha0 renormalisation.
+      ntagph = get_n_tagged_photons() 
+      if (ntagph.ne.0) then
+        do i = 1, AMP_SPLIT_SIZE_BORN 
+          call amp_split_pos_to_orders(i, amp_orders)
+          born_wgt = amp_split(i)
+          ! this is the number of powers of 'e' in the born
+          qed_pow_b = amp_orders(qed_pos)
+          amp_orders(qed_pos) = amp_orders(qed_pos) + 2
+          if (amp_orders(qed_pos).gt.NLO_ORDERS(qed_pos)) cycle
+
+          amp_split_poles_ML(orders_to_amp_split_pos(amp_orders),1) = 
+     $        amp_split_poles_ML(orders_to_amp_split_pos(amp_orders),1) +
+     $         get_virtual_a0Gmu_conv(qed_pow_b,ntagph,1,born_wgt)
+
+          amp_split_finite_ML(orders_to_amp_split_pos(amp_orders)) = 
+     $        amp_split_finite_ML(orders_to_amp_split_pos(amp_orders)) +
+     $         get_virtual_a0Gmu_conv(qed_pow_b,ntagph,0,born_wgt)
+
+        enddo
+      endif
+c======================================================================
+
 c
 c Check poles for the first PS points when doing MC over helicities, and
 c for all phase-space points when not doing MC over helicities. Skip
@@ -268,7 +302,7 @@
          polecheck_passed = .true.
          ! loop over the full result and each of the amp_split
          ! contribution
-         do iamp=0,amp_split_size
+         do iamp=1,amp_split_size
           ! skip 0 contributions in the amp_split array
             if (iamp.ne.0) then
                if (amp_split_poles_FKS(iamp,1).eq.0d0.and.
 
=== modified file 'Template/NLO/SubProcesses/chooser_functions.f'
--- Template/NLO/SubProcesses/chooser_functions.f	2018-04-05 17:59:21 +0000
+++ Template/NLO/SubProcesses/chooser_functions.f	2021-10-15 11:42:16 +0000
@@ -100,6 +100,10 @@
       double precision ch_i,ch_j,ch_m
       integer particle_type_born(nexternal-1)
       common /c_particle_type_born/particle_type_born
+      logical particle_tag(nexternal)
+      common /c_particle_tag/particle_tag
+      logical particle_tag_born(nexternal-1)
+      common /c_particle_tag/particle_tag_born
       logical need_color_links, need_charge_links
       common /c_need_links/need_color_links, need_charge_links
       integer extra_cnt, isplitorder_born, isplitorder_cnt
@@ -130,6 +134,7 @@
             stop
          endif
          particle_type(i)=particle_type_D(nFKSprocess,i)
+         particle_tag(i)=particle_tag_D(nFKSprocess,i)
          particle_charge(i)=particle_charge_D(nFKSprocess,i)
          pdg_type(i)=pdg_type_D(nFKSprocess,i)
          ! is_aorg is true if the particle can induce soft singularities
@@ -144,18 +149,22 @@
          if (i.lt.min(i_fks,j_fks)) then
             particle_type_born(i)=particle_type(i)
             particle_charge_born(i)=particle_charge(i)
+            particle_tag_born(i)=particle_tag(i)
          elseif (i.gt.max(i_fks,j_fks)) then
             particle_type_born(i-1)=particle_type(i)
             particle_charge_born(i-1)=particle_charge(i)
+            particle_tag_born(i-1)=particle_tag(i)
          elseif (i.eq.min(i_fks,j_fks)) then
             i_type=particle_type(i_fks)
             j_type=particle_type(j_fks)
             ch_i=particle_charge(i_fks)
             ch_j=particle_charge(j_fks)
+            particle_tag_born(i) = particle_tag(j_fks)
             call get_mother_col_charge(i_type,ch_i,j_type,ch_j,m_type,ch_m) 
             particle_type_born(i)=m_type
             particle_charge_born(i)=ch_m
          elseif (i.ne.max(i_fks,j_fks)) then
+            particle_tag_born(i) = particle_tag(i)
             particle_type_born(i)=particle_type(i)
             particle_charge_born(i)=particle_charge(i)
          endif
 
=== modified file 'Template/NLO/SubProcesses/cuts.f'
--- Template/NLO/SubProcesses/cuts.f	2021-01-08 10:23:30 +0000
+++ Template/NLO/SubProcesses/cuts.f	2021-10-15 11:42:16 +0000
@@ -40,42 +40,28 @@
       external R2_04,invm2_04,pt_04,eta_04,pt,eta
 C     recombination of photons
       double precision p_reco(0:4,nexternal), R_reco
-      integer iPDG_reco(nexternal)
+      integer iPDG_reco(nexternal),nphiso
 c local integers
       integer i,j
-c temporary variable for caching locally computation
-      double precision tmpvar
+c bare parton algorithm
+      integer nPART
+      double precision pPART(0:3,nexternal)
 c jet cluster algorithm
       integer nQCD,NJET,JET(nexternal)
       double precision pQCD(0:3,nexternal),PJET(0:3,nexternal)
-      double precision rfj,sycut,palg,amcatnlo_fastjetdmerge
       integer njet_eta
-      integer mm
-c Photon isolation
-      integer nph,nem,k,nin
-      double precision ptg,chi_gamma_iso,iso_getdrv40
-      double precision Etsum(0:nexternal)
-      real drlist(nexternal)
-      double precision pgamma(0:3,nexternal),pem(0:3,nexternal)
-      logical alliso
-c Sort array of results: ismode>0 for real, isway=0 for ascending order
-      integer ismode,isway,izero,isorted(nexternal)
-      parameter (ismode=1)
-      parameter (isway=0)
-      parameter (izero=0)
-c The UNLOPS cut
-      double precision p_unlops(0:3,nexternal)
+      double precision pgamma(0:3,nexternal),pgamma_iso(0:3,nexternal)
+      integer nph
       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)
-     $     ,is_a_ph(nexternal)
+     $,is_a_ph(nexternal),is_nph_iso(nexternal),is_nextph_iso(nexternal)
+     $,is_nextph_iso_reco(nexternal)
+      logical is_a_lp_reco(nexternal),is_a_lm_reco(nexternal)
+      logical passcuts_leptons, passcuts_unlops_jv, passcuts_photons, 
+     $        passcuts_jets, passcuts_pdgs 
+
 
       passcuts_user=.true. ! event is okay; otherwise it is changed
 
@@ -84,98 +70,329 @@
 C Cuts from the run_card.dat
 C***************************************************************
 C***************************************************************
-      !first recombine the photons and fermions
+
+      ! Find the bare QCD partons
+      ! This is used only as input to the photon isolation
+      call identify_PART_partons(p,istatus,ipdg,pPART,nPART,is_a_lp,is_a_lm)
+
+      ! Apply the Photon cuts on isolated photons based on the bare particles
+      passcuts_user = passcuts_user .and.
+     $ passcuts_photons(p,istatus,ipdg,is_a_lp,is_a_lm,pPART,nPART,
+     $ pgamma,nph,is_nph_iso,is_nextph_iso)
+      if (.not.passcuts_user) return
+
+      ! Recombine the photons and fermions
       call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
-     $                       p, iPDG, p_reco, iPDG_reco)
-
-c
-c CHARGED LEPTON CUTS
-c
-c find the charged leptons (also used in the photon isolation cuts below)
-      do i=1,nexternal
-         if(istatus(i).eq.1 .and.
-     &    (ipdg_reco(i).eq.11 .or. ipdg_reco(i).eq.13 .or. ipdg_reco(i).eq.15)) then
-            is_a_lm(i)=.true.
-         else
-            is_a_lm(i)=.false.
-         endif
-         if(istatus(i).eq.1 .and.
-     &    (ipdg_reco(i).eq.-11 .or. ipdg_reco(i).eq.-13 .or. ipdg_reco(i).eq.-15)) then
-            is_a_lp(i)=.true.
-         else
-            is_a_lp(i)=.false.
-         endif
-      enddo
-c apply the charged lepton cuts
-      do i=nincoming+1,nexternal
-         if (is_a_lp(i).or.is_a_lm(i)) then
-c transverse momentum
-            if (ptl.gt.0d0) then
-               if (pt_04(p_reco(0,i)).lt.ptl) then
-                  passcuts_user=.false.
-                  return
-               endif
-            endif
-c pseudo-rapidity
-            if (etal.gt.0d0) then
-               if (abs(eta_04(p_reco(0,i))).gt.etal) then
-                  passcuts_user=.false.
-                  return
-               endif
-            endif
-c DeltaR and invariant mass cuts
-            if (is_a_lp(i)) then
-               do j=nincoming+1,nexternal
-                  if (is_a_lm(j)) then
-                     if (drll.gt.0d0) then
-                        if (R2_04(p_reco(0,i),p_reco(0,j)).lt.drll**2) then
-                           passcuts_user=.false.
-                           return
-                        endif
-                     endif
-                     if (mll.gt.0d0) then
-                        if (invm2_04(p_reco(0,i),p_reco(0,j),1d0).lt.mll**2) then
-                           passcuts_user=.false.
-                           return
-                        endif
-                     endif
-                     if (ipdg_reco(i).eq.-ipdg_reco(j)) then
-                        if (drll_sf.gt.0d0) then
-                           if (R2_04(p_reco(0,i),p_reco(0,j)).lt.drll_sf**2) then
-                              passcuts_user=.false.
-                              return
-                           endif
-                        endif
-                        if (mll_sf.gt.0d0) then
-                           if (invm2_04(p_reco(0,i),p_reco(0,j),1d0).lt.mll_sf**2)
-     $                          then
-                              passcuts_user=.false.
-                              return
-                           endif
-                        endif
-                     endif
+     $     p, iPDG, is_nextph_iso,  p_reco, iPDG_reco, is_nextph_iso_reco)
+
+      ! Apply the reco lepton cuts
+      passcuts_user = passcuts_user .and. 
+     $                  passcuts_leptons(p_reco,istatus,ipdg_reco,is_a_lp_reco,is_a_lm_reco)
+      if (.not.passcuts_user) return
+
+      ! Find the reco QCD partons including 
+      ! A. All photons if gamma_is_j is on
+      ! B. Non-iso, non-reco photons if reco is on
+      call identify_QCD_partons(is_nextph_iso_reco,p_reco,istatus,ipdg_reco,is_a_j,pQCD,nQCD)
+
+      ! Apply the UNLOPS/JetVeto cuts
+      passcuts_user = passcuts_user .and. 
+     $                  passcuts_unlops_jv(p_reco,istatus,ipdg_reco,pQCD,nQCD,ickkw)
+      if (.not.passcuts_user) return
+
+      ! Apply the Jet cuts
+      passcuts_user = passcuts_user .and. 
+     $                  passcuts_jets(p_reco,pQCD,nQCD,pgamma,nph,is_nph_iso,ickkw)
+      if (.not.passcuts_user) return
+
+      ! Apply PDG specific cuts
+      passcuts_user = passcuts_user .and. 
+     $                  passcuts_pdgs(p_reco,istatus,ipdg_reco)
+      if (.not.passcuts_user) return
+
+C***************************************************************
+C***************************************************************
+C PUT HERE YOUR USER-DEFINED CUTS
+C***************************************************************
+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
+c$$$C apply the pT cut (pT should be large than 200 GeV for the event to
+c$$$C pass cuts)
+c$$$            if ( p(1,i)**2+p(2,i)**2 .lt. 200d0**2 ) then
+c$$$C momenta do not pass cuts. Set passcuts_user to false and return
+c$$$               passcuts_user=.false.
+c$$$               return
+c$$$            endif
+c$$$         endif
+c$$$      enddo
+c
+      return
+      end
+
+      subroutine identify_PART_partons(p,istatus,ipdg,pPART,nPART,is_a_lp,is_a_lm)
+      implicit none
+      include 'nexternal.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      double precision p(0:4,nexternal)
+      integer nPART
+      double precision pPART(0:3,nexternal)
+      logical is_a_lp(nexternal),is_a_lm(nexternal)
+      include "run.inc"
+      include "cuts.inc"
+
+      integer i, j
+c
+c Bare partons and leptons
+c
+      nPART=0
+      do j=1,nexternal
+         is_a_lp(j)=.false.
+         is_a_lm(j)=.false.
+c Partons
+        if (istatus(j).eq.1 .and.
+     &     (abs(ipdg(j)).le.maxjetflavor .or. ipdg(j).eq.21)
+     &) then
+            nPART=nPART+1
+            do i=0,3
+               pPART(i,nPART)=p(i,j)
+            enddo
+        endif
+c Leptons
+         if (ipdg(j).eq.11 .or. ipdg(j).eq.13
+     $      .or.  ipdg(j).eq.15) then
+            is_a_lm(j)=.true.
+         endif
+         if (ipdg(j).eq.-11 .or. ipdg(j).eq.-13
+     $      .or.  ipdg(j).eq.-15) then
+            is_a_lp(j)=.true.
+         endif
+      enddo
+
+      return
+      end
+
+      logical function passcuts_photons(p,istatus,ipdg,is_a_lp,is_a_lm,
+     $pPART,nPART,pgamma,nph,is_nph_iso,is_nextph_iso)
+      implicit none
+      include 'nexternal.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      double precision p(0:4,nexternal)
+      logical is_a_lp(nexternal),is_a_lm(nexternal)
+      integer nPART, nph
+      double precision pPART(0:3,nexternal), pgamma(0:3,nexternal)
+      double precision pgamma_iso(0:3,nexternal)
+      logical is_nph_iso(nexternal),is_nextph_iso(nexternal)
+      include "cuts.inc"
+      include "run.inc"
+      integer i,j,k,mu
+c Sort array of results: ismode>0 for real, isway=0 for ascending order
+      integer ismode,isway,izero,isorted(nexternal)
+      parameter (ismode=1)
+      parameter (isway=0)
+      parameter (izero=0)
+
+c Photon isolation
+      integer nem,nin,nphiso
+      double precision ptg,chi_gamma_iso,iso_getdrv40
+      double precision Etsum(0:nexternal)
+      real drlist(nexternal)
+      double precision pem(0:3,nexternal)
+
+      logical alliso,isolated
+      integer get_n_tagged_photons
+      logical is_a_ph(nexternal)
+
+      REAL*8 pt,eta
+      external pt,eta
+ 
+      passcuts_photons = .true.
+
+c
+c PHOTON (ISOLATION) CUTS
+c
+c Initialise common logical iso
+      do i=nincoming+1,nexternal
+        is_nextph_iso(i)=.False.
+      enddo
+c find the photons
+      do i=nincoming+1,nexternal
+         if (ipdg(i).eq.22 .and. .not.gamma_is_j) then
+            is_a_ph(i)=.true.
+         else
+            is_a_ph(i)=.false.
+         endif
+      enddo
+
+      if (ptgmin.ne.0d0) then
+         nph=0
+         do j=nincoming+1,nexternal
+            if (is_a_ph(j)) then
+               nph=nph+1
+               do i=0,3
+                  pgamma(i,nph)=p(i,j)
+               enddo
+            endif
+         enddo
+         if(nph.eq.0) return
+         
+         if(isoEM)then
+            nem=nph
+            do k=1,nem
+               do i=0,3
+                  pem(i,k)=pgamma(i,k)
+               enddo
+            enddo
+            do j=nincoming+1,nexternal
+               if (is_a_lp(j).or.is_a_lm(j)) then
+                  nem=nem+1
+                  do i=0,3
+                     pem(i,nem)=p(i,j)
+                  enddo
+               endif
+            enddo
+         endif
+         
+         nphiso=0
+
+         j=0
+c Loop over all photons
+         do while(j.lt.nph)
+
+            j=j+1
+            is_nph_iso(j)=.False. 
+            ptg=pt(pgamma(0,j))
+            if(ptg.lt.ptgmin)then
+               cycle
+            endif
+            if (etagamma.gt.0d0) then
+               if (abs(eta(pgamma(0,j))).gt.etagamma) then
+                  cycle
+               endif
+            endif
+         
+c Isolate from hadronic energy
+            do i=1,nPART
+               drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pPART(0,i)))
+            enddo
+            call sortzv(drlist,isorted,nPART,ismode,isway,izero)
+            Etsum(0)=0.d0
+            nin=0
+            do i=1,nPART
+               if(dble(drlist(isorted(i))).le.R0gamma)then
+                  nin=nin+1
+                  Etsum(nin)=Etsum(nin-1)+pt(pPART(0,isorted(i)))
+               endif
+            enddo
+            isolated=.True.
+            do i=1,nin
+               if(Etsum(i).gt.chi_gamma_iso(dble(drlist(isorted(i))),
+     $             R0gamma,xn,epsgamma,ptg)) then
+                   isolated=.False.
+                   exit
+               endif
+            enddo
+            if(.not.isolated)cycle
+            
+c Isolate from EM energy
+            if(isoEM.and.nem.gt.1)then
+               do i=1,nem
+                  drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i)))
+               enddo
+               call sortzv(drlist,isorted,nem,ismode,isway,izero)
+c First of list must be the photon: check this, and drop it
+                 if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then
+                    write(*,*)'Error #1 in photon isolation'
+                    write(*,*)j,isorted(1),drlist(isorted(1))
+                    stop
+                 endif
+               Etsum(0)=0.d0
+               nin=0
+               do i=2,nem
+                  if(dble(drlist(isorted(i))).le.R0gamma)then
+                     nin=nin+1
+                     Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i)))
                   endif
                enddo
+               isolated=.True.
+               do i=1,nin
+                 if(Etsum(i).gt.chi_gamma_iso(dble(drlist(isorted(i))),
+     $               R0gamma,xn,epsgamma,ptg)) then
+                    isolated=.False.
+                    exit
+                 endif
+               enddo
+               if(.not.isolated)cycle
             endif
+            is_nph_iso(j)=.True.
+            nphiso=nphiso+1
+
+           if (nphiso.gt.0) then
+             do mu=0,3
+               pgamma_iso(mu,nphiso)=pgamma(mu,j)
+             enddo
+
+             do i=nincoming+1,nexternal
+               if ( ipdg(i).eq.22 .and. 
+     $              pt(p(0,i)).eq.pt(pgamma_iso(0,nphiso)) ) then
+                 is_nextph_iso(i)=.True.
+               endif
+             enddo
+           endif
+         enddo
+c End of loop over photons
+
+         if(nphiso.lt.get_n_tagged_photons())then
+            passcuts_photons=.false.
+            return
          endif
-      enddo
+      endif
+
+      return
+      end
+
+
+      subroutine identify_QCD_partons(is_iso,p,istatus,ipdg,is_a_j,pQCD,nQCD)
+      implicit none
+      include 'nexternal.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      double precision p(0:4,nexternal)
+      logical is_a_j(nexternal)
+      integer nQCD
+      double precision pQCD(0:3,nexternal)
+      logical is_iso(nexternal)
+      REAL*8 pt,eta
+      external pt,eta
+      include "run.inc" 
+      include "cuts.inc"
+
+      integer i, j 
+
 c
 c JET CUTS
 c
 c find the jets
       do i=1,nexternal
          if (istatus(i).eq.1 .and.
-     &        (abs(ipdg_reco(i)).le.maxjetflavor .or. ipdg_reco(i).eq.21
-     &         .or.(ipdg_reco(i).eq.22.and.gamma_is_j))) then
+     &        (  abs(ipdg(i)).le.maxjetflavor .or. ipdg(i).eq.21
+     &         .or. (ipdg(i).eq.22.and.gamma_is_j) .or.
+     &         (ipdg(i).eq.22.and. .not.is_iso(i))  )
+     &) then
             is_a_j(i)=.true.
          else
             is_a_j(i)=.false.
          endif
       enddo
-
 c If we do not require a mimimum jet energy, there's no need to apply
 c jet clustering and all that.
-      if (ptj.gt.0d0.or.ptgmin.gt.0d0) then
+      if (ptj.ne.0d0.or.ptgmin.ne.0d0) then
 c Put all (light) QCD partons in momentum array for jet clustering.
 c From the run_card.dat, maxjetflavor defines if b quark should be
 c considered here (via the logical variable 'is_a_jet').  nQCD becomes
@@ -186,41 +403,43 @@
             if (is_a_j(j)) then
                nQCD=nQCD+1
                do i=0,3
-                  pQCD(i,nQCD)=p_reco(i,j)
+                  pQCD(i,nQCD)=p(i,j)
                enddo
             endif
          enddo
       endif
 
-c THE UNLOPS CUT:
-      if (ickkw.eq.4 .and. ptj.gt.0d0) then
-c Use special pythia pt cut for minimal pT
-         do i=1,nexternal
-            do j=0,3
-               p_unlops(j,i)=p_reco(j,i)
-            enddo
-         enddo
-         call pythia_UNLOPS(p_unlops,passUNLOPScuts)
-         if (.not. passUNLOPScuts) then
-            passcuts_user=.false.
-            return
-         endif
-c Bypass normal jet cuts
-         goto 122
-c THE VETO XSEC CUT:
-      elseif (ickkw.eq.-1 .and. ptj.gt.0d0) then
-c Use veto'ed Xsec for analytic NNLL resummation
-         if (nQCD.ne.1) then
-            write (*,*) 'ERROR: more than one QCD parton in '/
-     $           /'this event in cuts.f. There should only be one'
-            stop
-         endif
-         if (pt(pQCD(0,1)) .gt. ptj) then
-            passcuts_user=.false.
-            return
-         endif
-      endif
-
+      return
+      end
+
+
+      logical function passcuts_jets(p,pQCD,nQCD,pgamma,nph,is_nph_iso,ickkw)
+      implicit none
+      include 'nexternal.inc'
+      double precision p(0:4,nexternal)
+      integer nQCD, nph
+      double precision pQCD(0:3,nexternal), pgamma(0:3,nexternal)
+      logical is_nph_iso(nexternal)
+      integer ickkw
+      include "cuts.inc"
+
+      integer NJET,JET(nexternal)
+      double precision rfj,sycut,palg,amcatnlo_fastjetdmerge
+      double precision PJET(0:3,nexternal)
+      integer mm
+      integer i,j
+
+      integer get_n_tagged_photons
+
+      REAL*8 pt,eta
+      external pt,eta
+
+      passcuts_jets=.true.
+
+c JET CUTS
+
+C       do nothing if ickkw=4 (UNLOPS)
+      if (ickkw.eq.4)return
 
       if (ptj.gt.0d0.and.nQCD.gt.1) then
 
@@ -233,11 +452,10 @@
             if(abs(pQCD(0,j)/p(0,1)).lt.1.d-8) mm=mm+1
          enddo
          if(mm.gt.1)then
-            passcuts_user=.false.
+            passcuts_jets=.false.
             return
          endif
 
-
 c Define jet clustering parameters (from cuts.inc via the run_card.dat)
          palg=JETALGO           ! jet algorithm: 1.0=kt, 0.0=C/A, -1.0 = anti-kt
          rfj=JETRADIUS          ! the radius parameter
@@ -268,127 +486,178 @@
 
 c Apply the jet cuts
          if (njet .ne. nQCD .and. njet .ne. nQCD-1) then
-            passcuts_user=.false.
-            return
-         endif
-      endif
- 122  continue
-c
-c PHOTON (ISOLATION) CUTS
-c
-c find the photons
+            passcuts_jets=.false.
+            return
+         endif
+      endif
+
+      return
+      end
+
+
+
+      logical function passcuts_unlops_jv(p,istatus,ipdg,pQCD,nQCD,ickkw)
+      implicit none
+      include 'nexternal.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      double precision p(0:4,nexternal)
+      logical is_a_j(nexternal)
+      integer nQCD
+      double precision pQCD(0:3,nexternal)
+      integer ickkw
+      include "cuts.inc"
+      double precision p_unlops(0:3,nexternal)
+      logical passUNLOPScuts
+      integer i, j 
+
+      REAL*8 pt
+      external pt
+
+      passcuts_unlops_jv=.true.
+
+
+c THE UNLOPS CUT:
+      if (ickkw.eq.4 .and. ptj.gt.0d0) then
+c Use special pythia pt cut for minimal pT
+         do i=1,nexternal
+            do j=0,3
+               p_unlops(j,i)=p(j,i)
+            enddo
+         enddo
+         call pythia_UNLOPS(p_unlops,passUNLOPScuts)
+         if (.not. passUNLOPScuts) then
+            passcuts_unlops_jv=.false.
+            return
+         endif
+c THE VETO XSEC CUT:
+      elseif (ickkw.eq.-1 .and. ptj.gt.0d0) then
+c Use veto'ed Xsec for analytic NNLL resummation
+         if (nQCD.ne.1) then
+            write (*,*) 'ERROR: more than one QCD parton in '/
+     $           /'this event in cuts.f. There should only be one'
+            stop
+         endif
+         if (pt(pQCD(0,1)) .gt. ptj) then
+            passcuts_unlops_jv=.false.
+            return
+         endif
+      endif
+      return
+      end
+
+
+
+      logical function passcuts_leptons(p,istatus,ipdg,is_a_lp_reco,is_a_lm_reco)
+      implicit none
+      include 'nexternal.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      double precision p(0:4,nexternal)
+      logical is_a_lp_reco(nexternal),is_a_lm_reco(nexternal)
+
+      REAL*8 R2_04,invm2_04,pt_04,eta_04,pt,eta
+      external R2_04,invm2_04,pt_04,eta_04,pt,eta
+      integer i,j
+
+      include 'cuts.inc'
+
+      passcuts_leptons=.true.
+
+c
+c CHARGED LEPTON CUTS
+c
+c find the charged leptons (also used in the photon isolation cuts below)
       do i=1,nexternal
-         if (istatus(i).eq.1 .and. ipdg(i).eq.22 .and. .not.gamma_is_j) then
-            is_a_ph(i)=.true.
-         else
-            is_a_ph(i)=.false.
+         if(istatus(i).eq.1 .and.
+     &    (ipdg(i).eq.11 .or. ipdg(i).eq.13 .or. ipdg(i).eq.15)) then
+            is_a_lm_reco(i)=.true.
+         else
+            is_a_lm_reco(i)=.false.
+         endif
+         if(istatus(i).eq.1 .and.
+     &    (ipdg(i).eq.-11 .or. ipdg(i).eq.-13 .or. ipdg(i).eq.-15)) then
+            is_a_lp_reco(i)=.true.
+         else
+            is_a_lp_reco(i)=.false.
          endif
       enddo
-      if (ptgmin.gt.0d0) then
-         nph=0
-         do j=nincoming+1,nexternal
-            if (is_a_ph(j)) then
-               nph=nph+1
-               do i=0,3
-                  pgamma(i,nph)=p(i,j)
-               enddo
-            endif
-         enddo
-         if(nph.eq.0)goto 444
-         write(*,*) 'ERROR in cuts.f: photon isolation is not working'
-     $           // ' for mixed QED-QCD corrections'
-         stop 1
-         
-         if(isoEM)then
-            nem=nph
-            do k=1,nem
-               do i=0,3
-                  pem(i,k)=pgamma(i,k)
-               enddo
-            enddo
-            do j=nincoming+1,nexternal
-               if (is_a_lp(j).or.is_a_lm(j)) then
-                  nem=nem+1
-                  do i=0,3
-                     pem(i,nem)=p(i,j)
-                  enddo
-               endif
-            enddo
-         endif
-         
-         alliso=.true.
-
-         j=0
-         do while(j.lt.nph.and.alliso)
-c Loop over all photons
-            j=j+1
-            
-            ptg=pt(pgamma(0,j))
-            if(ptg.lt.ptgmin)then
-               passcuts_user=.false.
-               return
-            endif
-            if (etagamma.gt.0d0) then
-               if (abs(eta(pgamma(0,j))).gt.etagamma) then
-                  passcuts_user=.false.
-                  return
-               endif
-            endif
-         
-c Isolate from hadronic energy
-            do i=1,nQCD
-               drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pQCD(0,i)))
-            enddo
-            call sortzv(drlist,isorted,nQCD,ismode,isway,izero)
-            Etsum(0)=0.d0
-            nin=0
-            do i=1,nQCD
-               if(dble(drlist(isorted(i))).le.R0gamma)then
-                  nin=nin+1
-                  Etsum(nin)=Etsum(nin-1)+pt(pQCD(0,isorted(i)))
-               endif
-            enddo
-            do i=1,nin
-               alliso=alliso .and.
-     $              Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
-     $              R0gamma,xn,epsgamma,ptg)
-            enddo
-            
-c Isolate from EM energy
-            if(isoEM.and.nem.gt.1)then
-               do i=1,nem
-                  drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i)))
-               enddo
-               call sortzv(drlist,isorted,nem,ismode,isway,izero)
-c First of list must be the photon: check this, and drop it
-               if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then
-                  write(*,*)'Error #1 in photon isolation'
-                  write(*,*)j,isorted(1),drlist(isorted(1))
-                  stop
-               endif
-               Etsum(0)=0.d0
-               nin=0
-               do i=2,nem
-                  if(dble(drlist(isorted(i))).le.R0gamma)then
-                     nin=nin+1
-                     Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i)))
+c apply the charged lepton cuts
+      do i=nincoming+1,nexternal
+         if (is_a_lp_reco(i).or.is_a_lm_reco(i)) then
+c transverse momentum
+            if (ptl.gt.0d0) then
+               if (pt_04(p(0,i)).lt.ptl) then
+                  passcuts_leptons=.false.
+                  return
+               endif
+            endif
+c pseudo-rapidity
+            if (etal.gt.0d0) then
+               if (abs(eta_04(p(0,i))).gt.etal) then
+                  passcuts_leptons=.false.
+                  return
+               endif
+            endif
+c DeltaR and invariant mass cuts
+            if (is_a_lp_reco(i)) then
+               do j=nincoming+1,nexternal
+                  if (is_a_lm_reco(j)) then
+                     if (drll.gt.0d0) then
+                        if (R2_04(p(0,i),p(0,j)).lt.drll**2) then
+                           passcuts_leptons=.false.
+                           return
+                        endif
+                     endif
+                     if (mll.gt.0d0) then
+                        if (invm2_04(p(0,i),p(0,j),1d0).lt.mll**2) then
+                           passcuts_leptons=.false.
+                           return
+                        endif
+                     endif
+                     if (ipdg(i).eq.-ipdg(j)) then
+                        if (drll_sf.gt.0d0) then
+                           if (R2_04(p(0,i),p(0,j)).lt.drll_sf**2) then
+                              passcuts_leptons=.false.
+                              return
+                           endif
+                        endif
+                        if (mll_sf.gt.0d0) then
+                           if (invm2_04(p(0,i),p(0,j),1d0).lt.mll_sf**2)
+     $                          then
+                              passcuts_leptons=.false.
+                              return
+                           endif
+                        endif
+                     endif
                   endif
                enddo
-               do i=1,nin
-                  alliso=alliso .and.
-     $               Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
-     $               R0gamma,xn,epsgamma,ptg)
-               enddo
             endif
-c End of loop over photons
-         enddo
-         if(.not.alliso)then
-            passcuts_user=.false.
-            return
          endif
- 444     continue
-c End photon isolation
-      endif
+      enddo
+
+      return
+      end
+
+      logical function passcuts_pdgs(p,istatus,ipdg)
+      implicit none
+      include 'nexternal.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      double precision p(0:4,nexternal)
+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
+      REAL*8 invm2_04,pt_04
+      external invm2_04,pt_04
+c temporary variable for caching locally computation
+      double precision tmpvar
+      integer i,j
+
+      passcuts_pdgs = .true.
+
 
 C
 C     PDG SPECIFIC CUTS (PT/M_IJ)
@@ -397,49 +666,24 @@
          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.
+               passcuts_pdgs=.false.
                return
             elseif (tmpvar.gt.etmax(i) .and. etmax(i).gt.0d0) then
-               passcuts_user=.false.
+               passcuts_pdgs=.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.
+                  passcuts_pdgs=.false.
                   return
                endif
             endif
-         enddo
+         enddo  
       enddo
-
-
-C***************************************************************
-C***************************************************************
-C PUT HERE YOUR USER-DEFINED CUTS
-C***************************************************************
-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
-c$$$C apply the pT cut (pT should be large than 200 GeV for the event to
-c$$$C pass cuts)
-c$$$            if ( p(1,i)**2+p(2,i)**2 .lt. 200d0**2 ) then
-c$$$C momenta do not pass cuts. Set passcuts_user to false and return
-c$$$               passcuts_user=.false.
-c$$$               return
-c$$$            endif
-c$$$         endif
-c$$$      enddo
-c
-      return
-      end
-
-
+      return  
+      end 
 
 
 C***************************************************************
@@ -983,5 +1227,21 @@
       return
       end
 
+      integer function get_n_tagged_photons()
+      implicit none
+      integer i
+      include "nexternal.inc"
+      logical particle_tag(nexternal)
+      common /c_particle_tag/particle_tag
+      get_n_tagged_photons = 0
+
+      do i = nincoming+1, nexternal
+        if (particle_tag(i))
+     $     get_n_tagged_photons = get_n_tagged_photons+1
+      enddo
+
+      return
+      end
+
 
 
 
=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2021-09-30 16:30:15 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2021-10-15 11:42:16 +0000
@@ -1519,6 +1519,12 @@
       double precision       wgt_ME_born,wgt_ME_real
       common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
 
+      integer ntagph
+      double precision resc
+      integer get_n_tagged_photons
+      double precision get_rescale_alpha_factor
+      external get_n_tagged_photons get_rescale_alpha_factor
+
       if (wgt1.eq.0d0 .and. wgt2.eq.0d0 .and. wgt3.eq.0d0) return
 c Check for NaN's and INF's. Simply skip the contribution
       if (wgt1.ne.wgt1) return
@@ -1582,9 +1588,21 @@
       icontr=icontr+1
       call weight_lines_allocated(nexternal,icontr,max_wgt,max_iproc)
       itype(icontr)=type
-      wgt(1,icontr)=wgt1
-      wgt(2,icontr)=wgt2
-      wgt(3,icontr)=wgt3
+
+C here we rescale the contributions by the ratio of alpha's in different
+C schemes; it is needed when there are tagged photons around
+      ntagph = get_n_tagged_photons()
+      if (ntagph.eq.0) then
+        wgt(1,icontr)=wgt1
+        wgt(2,icontr)=wgt2
+        wgt(3,icontr)=wgt3
+      else if (ntagph.gt.0) then
+          resc = get_rescale_alpha_factor(ntagph, orders(qed_pos)) 
+          wgt(1,icontr) = wgt1 * resc
+          wgt(2,icontr) = wgt2 * resc
+          wgt(3,icontr) = wgt3 * resc
+      endif
+
       bjx(1,icontr)=xbk(1)
       bjx(2,icontr)=xbk(2)
       scales2(1,icontr)=QES2
@@ -5555,6 +5573,8 @@
       integer fks_j_from_i(nexternal,0:nexternal)
      &     ,particle_type(nexternal),pdg_type(nexternal)
       common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
+      logical particle_tag(nexternal)
+      common /c_particle_tag/particle_tag
       double precision particle_charge(nexternal)
       common /c_charges/particle_charge
       include "run.inc"
@@ -5787,7 +5807,7 @@
                   if (particle_charge(i).eq.0d0.and.pdg_type(i).ne.22)
      $                 cycle
 C     set charge factors
-                  if (pdg_type(i).eq.22) then
+                  if (pdg_type(i).eq.22.and..not.particle_tag(i)) then
                      c_used = 0d0
                      gamma_used = gamma_ph
                      gammap_used = gammap_ph
@@ -6067,7 +6087,7 @@
                      if (particle_charge(i).eq.0d0.and.pdg_type(i).ne.22)
      $                    cycle
 C     set charge factors
-                     if (pdg_type(i).eq.22) then
+                     if (pdg_type(i).eq.22.and..not.particle_tag(i)) then
                         c_used = 0d0
                         gamma_used = gamma_ph
                         gammap_used = gammap_ph
@@ -6396,6 +6416,8 @@
       double precision particle_charge(nexternal), particle_charge_born(nexternal-1)
       common /c_charges/particle_charge
       common /c_charges_born/particle_charge_born
+      logical particle_tag(nexternal)
+      common /c_particle_tag/particle_tag
       include 'coupl.inc'
       include 'q_es.inc'
       double precision p(0:3,nexternal),xmu2,double,single
@@ -6504,7 +6526,7 @@
             if (pdg_type(i).ne.22) then
               contr2=contr2-particle_charge(i)**2
               contr1=contr1-3d0/2d0*particle_charge(i)**2
-            else
+            elseif (.not.particle_tag(i)) then
               contr1=contr1-gamma_ph
             endif
           else
@@ -6695,6 +6717,8 @@
       double precision particle_charge(nexternal), particle_charge_born(nexternal-1)
       common /c_charges/particle_charge
       common /c_charges_born/particle_charge_born
+      logical particle_tag(nexternal)
+      common /c_particle_tag/particle_tag
       double precision zero
       parameter (zero=0d0)
 
@@ -6839,7 +6863,7 @@
          do i=nincoming+1,nexternal
             if (pdg_type(i).eq.21) ngluons_FKS(nFKSprocess)
      $           =ngluons_FKS(nFKSprocess)+1
-            if (pdg_type(i).eq.22) nphotons_FKS(nFKSprocess)
+            if (pdg_type(i).eq.22.and..not.particle_tag(i)) nphotons_FKS(nFKSprocess)
      $           =nphotons_FKS(nFKSprocess)+1
          enddo
 
 
=== modified file 'Template/NLO/SubProcesses/makefile_fks_dir'
--- Template/NLO/SubProcesses/makefile_fks_dir	2020-11-27 13:41:46 +0000
+++ Template/NLO/SubProcesses/makefile_fks_dir	2021-10-15 11:42:16 +0000
@@ -32,6 +32,7 @@
 	 $(patsubst %.f,%.o,$(wildcard b_sf_???.f)) 			\
 	 born.o sborn_sf.o extra_cnt_wrapper.o				\
 	 $(patsubst %.f,%.o,$(wildcard born_cnt_*.f))			\
+	 rescale_alpha_tagged.o						\
 	 fks_Sij.o $(fastjetfortran_madfks) fks_singular.o		\
 	 montecarlocounter.o reweight_xsec.o boostwdir2.o       	\
 	 initcluster.o cluster.o splitorders_stuff.o 			\
 
=== modified file 'Template/NLO/SubProcesses/recmom.f'
--- Template/NLO/SubProcesses/recmom.f	2020-09-16 12:55:48 +0000
+++ Template/NLO/SubProcesses/recmom.f	2021-10-15 11:42:16 +0000
@@ -1,12 +1,17 @@
-      subroutine recombine_momenta(R, etaph, reco_l, reco_q, p_in, pdg_in, p_out, pdg_out)
+      subroutine recombine_momenta(R, etaph, reco_l, reco_q, p_in, 
+     $pdg_in, is_nextph_iso, p_out, pdg_out, is_nextph_iso_reco)
       implicit none
       ! recombine photons with the closest fermion if the distance is
       ! less than R and if the rapidity of photons is < etaph (etaph < 0
       ! means no cut). Output a new set of momenta and pdgs corresponding
       ! to the recombined particles. If recombination occurs the photon
       ! disappears from the output particles
+      ! If isolated photons exist (is_nextph_iso(i)=True), they are not 
+      ! used for the recombination
       ! arguments
       include 'nexternal.inc'
+      include 'run.inc'
+      include 'cuts.inc'
       double precision R, etaph, p_in(0:4,nexternal), p_out(0:4,nexternal)
       logical reco_l, reco_q
       integer pdg_in(nexternal), pdg_out(nexternal)
@@ -17,13 +22,24 @@
       integer n_ph, i_ph
       integer i,j
       integer ifreco
-      double precision dreco, dthis
+      double precision dreco, dthis, iso_getdrv40
       integer skip
       logical is_light_charged_fermion
-      double precision R2, eta
+      logical iso_in(nexternal), iso_out(nexternal)
+      double precision R2_04, eta_04
+      ! use from photon isolation
+      logical is_nextph_iso(nexternal)
+c      integer nphiso
+c      double precision pgamma_iso(0:3,nexternal)
+c      common / photon / is_nextph_iso,nphiso,pgamma_iso
+      ! Save for reco photon isolation
+      logical is_nextph_iso_reco(nexternal)
+c      common / photon_reco / is_nextph_iso_reco
       ! 
       integer times_reco
       common/to_times_reco/ times_reco
+      REAL*8 pt,eta
+      external pt,eta
       ! reset everything
       do j=1,nexternal
         pdg_out(j)=0
@@ -46,28 +62,36 @@
           nq = 0
       endif
 
-      ! count the photons
+      ! count the photons that are not isolated
       n_ph=0
       do i=nincoming+1, nexternal
-        if (pdg_in(i).eq.id_ph.and.
-     $   (abs(eta(p_in(0,i))).lt.etaph.or.etaph.lt.0d0)) then
+        if (pdg_in(i).eq.id_ph .and. pt(p_in(0,i)).ne.0d0 .and.
+     $   (abs(eta_04(p_in(0,i))).lt.etaph.or.etaph.lt.0d0)
+     $    .and. .not.is_nextph_iso(i) ) then
             n_ph=n_ph+1
             i_ph=i
         endif
       enddo
+
+c Define iso_in as the iso input from photon isolation
+      do i=nincoming+1,nexternal
+        iso_in(i)=is_nextph_iso(i)
+      enddo
+
       if (n_ph.eq.0 .or. (nl.eq.0 .and. nq.eq.0)) then
         ! do nothing
         do j=1,nexternal
           pdg_out(j)=pdg_in(j)
+          iso_out(j)=iso_in(j)
           do i=0,4
             p_out(i,j)=p_in(i,j)
           enddo
         enddo
-        return
       elseif (n_ph.eq.1) then
         ! do nothing for initial states
         do j=1,nincoming
           pdg_out(j)=pdg_in(j)
+          iso_out(j)=iso_in(j)
           do i=0,4
             p_out(i,j)=p_in(i,j)
           enddo
@@ -78,7 +102,7 @@
         if (i_ph.gt.0) then
           do i = nincoming+1, nexternal
             if (is_light_charged_fermion(pdg_in(i),nq,nl)) then
-              dthis=dsqrt(R2(p_in(0,i_ph),p_in(0,i)))
+              dthis=dsqrt(R2_04(p_in(0,i_ph),p_in(0,i)))
               if (dthis.le.dreco) then
                 dreco=dthis
                 ifreco=i
@@ -90,6 +114,7 @@
         ! do nothing also for final states
           do j=nincoming+1,nexternal
             pdg_out(j)=pdg_in(j)
+            iso_out(j)=iso_in(j)
             do i=0,4
               p_out(i,j)=p_in(i,j)
             enddo
@@ -100,11 +125,13 @@
           do j=nincoming+1,nexternal
             if (j.ne.i_ph.and.j.ne.ifreco) then
               pdg_out(j-skip)=pdg_in(j)
+              iso_out(j-skip)=iso_in(j)
               do i=0,4
                 p_out(i,j-skip)=p_in(i,j)
               enddo
             elseif (j.eq.ifreco) then
               pdg_out(j-skip)=pdg_in(j)
+              iso_out(j-skip)=iso_in(j)
               do i=0,3
                 p_out(i,j-skip)=p_in(i,j)+p_in(i,i_ph)
               enddo
@@ -119,6 +146,35 @@
         stop 1
       endif
 
+      do i=nincoming+1,nexternal
+        is_nextph_iso_reco(i)=iso_out(i)
+      enddo
+
+      return 
+      end
+
+
+
+
+
+      subroutine recombine_momenta_notagph(R, etaph, reco_l, reco_q, p_in, pdg_in, p_out, pdg_out)
+      implicit none
+      ! wraps recombine_momenta, it provides a function interface
+      ! for the case where no tagged photon exist 
+      ! (mostly for backward complatibility)
+      ! arguments
+      include 'nexternal.inc'
+      double precision R, etaph, p_in(0:4,nexternal), p_out(0:4,nexternal)
+      logical reco_l, reco_q
+      integer pdg_in(nexternal), pdg_out(nexternal)
+
+      logical is_iso_photon_in(nexternal), is_iso_photon_out(nexternal)
+
+      is_iso_photon_in(:) = .false.
+
+      call recombine_momenta(R, etaph, reco_l, reco_q, p_in, 
+     $pdg_in, is_iso_photon_in, p_out, pdg_out, is_iso_photon_out)
+
       return 
       end
 
 
=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2021-09-08 08:59:52 +0000
+++ UpdateNotes.txt	2021-10-15 11:42:16 +0000
@@ -1,6 +1,7 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
-3.1.2 (??/??/??)
+3.2.1 (XX/XX/21)
+      DP+HS+IT+MZ: EW corrections can be computed for tagged photons
       MZ+CS: Fixes relevant to PineAPPL: 
             - minor fix on the gluon pdg code in pineappl_interface.cc
             - the integration channels are threated the same way with or without PineAPPL 
 
=== modified file 'madgraph/core/diagram_generation.py'
--- madgraph/core/diagram_generation.py	2021-07-21 19:09:13 +0000
+++ madgraph/core/diagram_generation.py	2021-10-15 11:42:16 +0000
@@ -30,6 +30,7 @@
 
 import madgraph.core.base_objects as base_objects
 import madgraph.various.misc as misc
+import madgraph.fks.fks_tag as fks_tag
 from madgraph import InvalidCmd, MadGraph5Error
 from six.moves import range
 from six.moves import zip
@@ -1713,6 +1714,15 @@
                  if leg['state'] == True]
         polids = [tuple(leg['polarization'])  for leg in process_definition['legs'] \
                  if leg['state'] == True]
+
+        # keep track of the 'is_tagged' property of the legs if needed
+        try:
+            fstags = [leg['is_tagged'] for leg in process_definition['legs'] \
+                 if 'is_tagged' in leg.keys() and leg['state'] == True]
+
+        except KeyError:
+            fstags = []
+
         # Generate all combinations for the initial state
         for prod in itertools.product(*isids):
             islegs = [\
@@ -1735,10 +1745,16 @@
                 red_fsidlist.add(tuple(tag))
                 # Generate leg list for process
                 leg_list = [copy.copy(leg) for leg in islegs]
-                leg_list.extend([\
-                        base_objects.Leg({'id':id, 'state': True, 'polarization': fslegs[i]['polarization']}) \
-                        for i,id  in enumerate(prod)])
                 
+                if not fstags:
+                    leg_list.extend([\
+                            base_objects.Leg({'id':id, 'state': True, 'polarization': fsleg['polarization']}) \
+                            for id, fsleg in zip(prod, fslegs)])
+                else:
+                    leg_list.extend([\
+                            fks_tag.TagLeg({'id':id, 'state': True, 'polarization': fsleg['polarization'], 'is_tagged': tag}) \
+                            for id, fsleg, tag in zip(prod, fslegs, fstags)])
+
                 legs = base_objects.LegList(leg_list)
 
                 # Check for crossed processes
 
=== modified file 'madgraph/fks/fks_base.py'
--- madgraph/fks/fks_base.py	2021-05-23 14:17:24 +0000
+++ madgraph/fks/fks_base.py	2021-10-15 11:42:16 +0000
@@ -416,6 +416,7 @@
         legs = [(leg.get('id'), leg) for leg in leglist]
         self.pdgs = array.array('i',[s[0] for s in legs])
         self.colors = [leg['color'] for leg in leglist]
+        self.particle_tags = [leg['is_tagged'] for leg in leglist]
         if not self.process['perturbation_couplings'] == ['QCD']:
             self.charges = [leg['charge'] for leg in leglist]
         else:
@@ -529,6 +530,13 @@
                     leg in self.born_amp['process']['legs']]                    
 
 
+    def get_is_tagged(self):
+        """return the list of the 'is_tagged' keys
+        of each leg in born_amp"""
+        return [leg.get('is_tagged') for \
+                    leg in self.born_amp['process']['legs']]                    
+
+
     def get_leglist(self):
         """return the leg list
         for the born amp"""
 
=== modified file 'madgraph/fks/fks_common.py'
--- madgraph/fks/fks_common.py	2021-03-16 15:13:03 +0000
+++ madgraph/fks/fks_common.py	2021-10-15 11:42:16 +0000
@@ -263,7 +263,6 @@
     if dict == {}:
         dict = find_pert_particles_interactions(model, pert)
     splittings = []
-#check that the leg is a qcd leg
 
     if leg.get('id') in dict['pert_particles']:
         part = model.get('particle_dict')[leg.get('id')]
@@ -284,6 +283,13 @@
                         nsoft += 1
                 if nsoft >= 1:
                     for split in split_leg(leg, parts, model):
+                        # check if the leg is tagged, that 
+                        # the same particles appear also in the two daughters
+                        if 'is_tagged' in leg.keys() and leg['is_tagged'] and \
+                           leg['id'] != split[0]['id'] and \
+                           leg['id'] != split[1]['id']:
+                            continue
+
                         # add the splitting, but check if there is 
                         # an initial-state lepton if the flag
                         # include_init_leptons is False
@@ -833,6 +839,7 @@
     -'charge', which gives the charge of the leg
     -'massless', boolean, true if leg is massless
     -'spin' which gives the spin of leg
+    -'is_tagged', boolean, true if leg is tagged in the final state
     -'is_part', boolean, true if leg is an particle
     -'self_antipart', boolean, true if leg is an self-conjugated particle
     """
@@ -846,13 +853,14 @@
         self['charge'] = 0.
         self['massless'] = True
         self['spin'] = 0
+        self['is_tagged'] = False
         self['is_part'] = True
         self['self_antipart'] = False
     
     def get_sorted_keys(self):
         """Return particle property names as a nicely sorted list."""
         keys = super(FKSLeg, self).get_sorted_keys()
-        keys += ['fks', 'color','charge', 'massless', 'spin','is_part','self_antipart']
+        keys += ['fks', 'color','charge', 'massless', 'spin','is_tagged','is_part','self_antipart',]
         return keys
 
     
@@ -868,7 +876,7 @@
                 six.reraise(self.PhysicsObjectError, "%s is not a valid leg %s flag" % \
                                                  str(value), name)
                                                  
-        if name in ['massless','self_antipart','is_part']:
+        if name in ['massless','is_tagged','self_antipart','is_part']:
             if not isinstance(value, bool):
                 six.reraise(self.PhysicsObjectError, "%s is not a valid boolean for leg flag %s" % \
                                                                     str(value), name)
@@ -879,3 +887,4 @@
         return super(FKSLeg,self).filter(name, value)
     
      
+
 
=== modified file 'madgraph/fks/fks_helas_objects.py'
--- madgraph/fks/fks_helas_objects.py	2021-02-08 22:58:11 +0000
+++ madgraph/fks/fks_helas_objects.py	2021-10-15 11:42:16 +0000
@@ -901,6 +901,7 @@
     contains:
     -- colors
     -- charges
+    -- particle_tags
     -- i/j/ij fks, ij refers to the born leglist
     -- ijglu
     -- need_color_links
@@ -919,6 +920,7 @@
         if fksrealproc != None:
             self.isfinite = False
             self.colors = fksrealproc.colors
+            self.particle_tags = fksrealproc.particle_tags
             self.charges = fksrealproc.charges
             self.fks_infos = fksrealproc.fks_infos
             self.is_to_integrate = fksrealproc.is_to_integrate
 
=== added file 'madgraph/fks/fks_tag.py'
--- madgraph/fks/fks_tag.py	1970-01-01 00:00:00 +0000
+++ madgraph/fks/fks_tag.py	2021-10-15 11:42:16 +0000
@@ -0,0 +1,75 @@
+################################################################################
+#
+# Copyright (c) 2009 The MadGraph5_aMC@NLO Development team and Contributors
+#
+# This file is a part of the MadGraph5_aMC@NLO project, an application which 
+# automatically generates Feynman diagrams and matrix elements for arbitrary
+# high-energy processes in the Standard Model and beyond.
+#
+# It is subject to the MadGraph5_aMC@NLO license which should accompany this 
+# distribution.
+#
+# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
+#
+################################################################################
+
+"""Definitions of the objects needed both for MadFKS with tagged particles"""
+
+import madgraph.core.base_objects as MG
+
+
+class MultiTagLeg(MG.MultiLeg):
+    """a daughter class of MultiLeg, with the extra possibility of specifying
+    whether a given leg is tagged or not, via the "is_tagged" key
+    """
+
+    def default_setup(self):
+        """Default values for all properties"""
+        super(MultiTagLeg, self).default_setup()
+        self['is_tagged'] = False
+    
+    def get_sorted_keys(self):
+        """Return particle property names as a nicely sorted list."""
+        keys = super(MultiTagLeg, self).get_sorted_keys()
+        keys += ['is_tagged']
+        return keys
+
+    
+    def filter(self, name, value):
+        """Filter for valid leg property values."""
+
+        if name == 'is_tagged':
+            if not isinstance(value, bool):
+                raise self.PhysicsObjectError( \
+                        "%s is not a valid string for leg 'is_tagged' flag" \
+                                                        % str(value))
+        return super(MultiTagLeg,self).filter(name, value)
+    
+     
+
+class TagLeg(MG.Leg):
+    """a daughter class of Leg, with the extra possibility of specifying
+    whether a given leg is tagged or not, via the "is_tagged" key
+    """
+
+    def default_setup(self):
+        """Default values for all properties"""
+        super(TagLeg, self).default_setup()
+        self['is_tagged'] = False
+    
+    def get_sorted_keys(self):
+        """Return particle property names as a nicely sorted list."""
+        keys = super(TagLeg, self).get_sorted_keys()
+        keys += ['is_tagged']
+        return keys
+
+    
+    def filter(self, name, value):
+        """Filter for valid leg property values."""
+
+        if name == 'is_tagged':
+            if not isinstance(value, bool):
+                raise self.PhysicsObjectError( \
+                        "%s is not a valid string for leg 'is_tagged' flag" \
+                                                        % str(value))
+        return super(TagLeg,self).filter(name, value)
 
=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py	2021-08-20 21:40:18 +0000
+++ madgraph/interface/madgraph_interface.py	2021-10-15 11:42:16 +0000
@@ -93,6 +93,8 @@
 import madgraph.various.misc as misc
 import madgraph.various.cluster as cluster
 
+import madgraph.fks.fks_tag as fks_tag
+
 import models as ufomodels
 import models.import_ufo as import_ufo
 import models.write_param_card as param_writer
@@ -4927,6 +4929,21 @@
                 state = True
                 continue
 
+            # check if the particle is tagged (!PART!)
+            if part_name.startswith('!') and part_name.endswith('!'):
+                part_name = part_name[1:-1]
+                is_tagged = True
+            elif part_name.endswith('!') and part_name.count('!') == 2 and part_name[:part_name.find('!')].isdigit():
+                part_name = part_name.replace('!','')
+                is_tagged = True
+#                misc.sprint(part_name)
+            else:
+                is_tagged = False
+
+            # check that only final-state particles are tagged
+            if is_tagged and not state:
+                raise self.InvalidCmd("initial particles cannot be tagged")
+
             mylegids = []
             polarization = []
             if '{' in part_name:
@@ -5006,6 +5023,9 @@
 
             duplicate =1
             if part_name in self._multiparticles:
+                # multiparticles cannot be tagged
+                if is_tagged:
+                    raise self.InvalidCmd("Multiparticles cannot be tagged")
                 if isinstance(self._multiparticles[part_name][0], list):
                     raise self.InvalidCmd("Multiparticle %s is or-multiparticle" % part_name + \
                           " which can be used only for required s-channels")
@@ -5036,12 +5056,26 @@
 
             if mylegids:
                 for _ in range(duplicate):
-                    myleglist.append(base_objects.MultiLeg({'ids':mylegids,
-                                                        'state':state,
-                                                        'polarization': polarization}))
+                    if LoopOption in ['virt','sqrvirt','tree','noborn']:
+                        # check that no tagged particles exist in this mode
+                        if is_tagged:
+                            raise self.InvalidCmd(
+                                "%s mode does not handle tagged particles" % LoopOption)
+
+                        myleglist.append(base_objects.MultiLeg({'ids':mylegids,
+                                                            'state':state,
+                                                            'polarization': polarization}))
+                    else:
+                        myleglist.append(fks_tag.MultiTagLeg({'ids':mylegids,
+                                                          'state':state,
+                                                          'polarization': polarization,
+                                                          'is_tagged':is_tagged}))
             else:
                 raise self.InvalidCmd("No particle %s in model" % part_name)
 
+        if any(['is_tagged' in l.keys()  and l['is_tagged'] for l in myleglist]):
+            logger.warning('The process involves tagged particles. Please consider citing arXiv:2106.02059 if relevant.')
+
         # Apply the keyword 'all' for perturbed coupling orders.
         if perturbation_couplings.lower() in ['all', 'loonly']:
             if perturbation_couplings.lower() in ['loonly']:
 
=== modified file 'madgraph/iolibs/export_fks.py'
--- madgraph/iolibs/export_fks.py	2021-07-21 19:12:36 +0000
+++ madgraph/iolibs/export_fks.py	2021-10-15 11:42:16 +0000
@@ -627,6 +627,16 @@
                             writers.FortranWriter(filename),
                             matrix_element)
 
+        filename = 'a0Gmuconv.inc'
+        startfroma0 = self.write_a0gmuconv_file(
+                            writers.FortranWriter(filename),
+                            matrix_element)
+
+        filename = 'rescale_alpha_tagged.f'
+        self.write_rescale_a0gmu_file(
+                            writers.FortranWriter(filename),
+                            startfroma0, matrix_element)
+
         filename = 'orders.h'
         self.write_orders_c_header_file(
                             writers.CPPWriter(filename),
@@ -1096,6 +1106,23 @@
         writer.writelines(text)
 
 
+    def write_a0gmuconv_file(self, writer, matrix_element):
+        """writes an include file with the informations about the 
+        alpha0 < > gmu conversion, to be used when the process has
+        tagged photons
+        """
+
+        bool_dict = {True: '.true.', False: '.false.'}
+        bornproc = matrix_element.born_me['processes'][0]
+        startfromalpha0 = False
+        if any([l['is_tagged'] and l['id'] == 22 for l in bornproc['legs']]):
+            if 'loop_qcd_qed_sm_a0' in bornproc['model'].get('modelpath'):
+                startfromalpha0 = True
+
+        text = 'logical  startfroma0\nparameter (startfroma0=%s)\n' % bool_dict[startfromalpha0]
+        writer.writelines(text)
+        return startfromalpha0
+
     def write_orders_c_header_file(self, writer, amp_split_size, amp_split_size_born):
         """writes the header file including the amp_split_size declaration for amcblast
 	"""
@@ -1105,6 +1132,57 @@
         writer.writelines(text)
 
 
+    def write_rescale_a0gmu_file(self, writer, startfroma0, matrix_element):
+        """writes the function that computes the rescaling factor needed in
+        the case of external photons
+        """
+
+        # get the model parameters
+        params = sum([v for v in self.model.get('parameters').values()], [])
+        parnames = [p.name.lower() for p in params]
+
+        bornproc = matrix_element.born_me['processes'][0]
+        # this is to ensure compatibility with standard processes
+        if not any([l['is_tagged'] and l['id'] == 22 for l in bornproc['legs']]):
+            to_check = []
+            expr = '1d0'
+            conv_pol = '0d0'
+            conv_fin = '0d0'
+        
+        elif startfroma0:
+            to_check = ['mdl_aewgmu', 'mdl_aew']
+            base = 'mdl_aewgmu/mdl_aew'
+            exp = 'qed_pow/2d0-ntag'
+            expr = '(%s)**(%s)' % (base, exp)
+            conv_fin = '(qed_pow - ntagph * 2d0) * MDL_ECOUP_DGMUA0_UV_EW_FIN_ * born_wgt'
+            conv_pol = '(qed_pow - ntagph * 2d0) * MDL_ECOUP_DGMUA0_UV_EW_1EPS_ * born_wgt'
+        else:
+            to_check = ['mdl_aew', 'mdl_aew0']
+            base = 'mdl_aew0/mdl_aew'
+            exp = 'ntag'
+            expr = '(%s)**(%s)' % (base, exp)
+            conv_fin = '- ntagph * 2d0 * MDL_ECOUP_DGMUA0_UV_EW_FIN_ * born_wgt'
+            conv_pol = '- ntagph * 2d0 * MDL_ECOUP_DGMUA0_UV_EW_1EPS_ * born_wgt'
+
+        replace_dict = {'rescale_fact': expr,
+                        'virtual_a0Gmu_conv_finite': conv_fin,
+                        'virtual_a0Gmu_conv_pole': conv_pol}
+
+        if not all(p in parnames for p in to_check):
+            raise fks_common.FKSProcessError(
+                    'Some parameters needed when there are tagged '+\
+                    'photons cannot be found in the model.\n' +\
+                    'Please load the correct model and restriction ' +\
+                    '(e.g loop_qcd_qed_sm_Gmu-a0 or loop_qcd_qed_sm_a0-Gmu)')
+
+        file = open(os.path.join(_file_path, \
+                          'iolibs/template_files/rescale_alpha_tagged.inc')).read()
+        file = file % replace_dict
+        
+        # Write the file
+        writer.writelines(file)
+
+
 
     def write_orders_file(self, writer, matrix_element):
         """writes the include file with the informations about coupling orders.
@@ -2857,6 +2935,7 @@
             col_lines = []
             pdg_lines = []
             charge_lines = []
+            tag_lines = []
             fks_j_from_i_lines = []
             split_type_lines = []
             for i, info in enumerate(fks_info_list):
@@ -2870,6 +2949,9 @@
                     'DATA (PARTICLE_CHARGE_D(%d, IPOS), IPOS=1, NEXTERNAL) / %s /'\
                     % (i + 1, ', '.join('%19.15fd0' % charg\
                                         for charg in fksborn.real_processes[info['n_me']-1].charges) ))
+                tag_lines.append( \
+                    'DATA (PARTICLE_TAG_D(%d, IPOS), IPOS=1, NEXTERNAL) / %s /' \
+                    % (i + 1, ', '.join(bool_dict[tag] for tag in fksborn.real_processes[info['n_me']-1].particle_tags) ))
                 fks_j_from_i_lines.extend(self.get_fks_j_from_i_lines(fksborn.real_processes[info['n_me']-1],\
                                                 i + 1))
                 split_type_lines.append( \
@@ -2884,6 +2966,7 @@
             pdgs = [l.get('id') for l in bornproc.get('legs')] + [-21]
             colors = [l.get('color') for l in bornproc.get('legs')] + [8]
             charges = [l.get('charge') for l in bornproc.get('legs')] + [0.]
+            tags = [l.get('is_tagged') for l in bornproc.get('legs')] + [False]
 
             fks_i = len(colors)
             # fist look for a colored legs (set j to 1 otherwise)
@@ -2919,6 +3002,8 @@
                             % ', '.join([str(pdg) for pdg in pdgs])]
             charge_lines = ['DATA (PARTICLE_CHARGE_D(1, IPOS), IPOS=1, NEXTERNAL) / %s /' \
                             % ', '.join('%19.15fd0' % charg for charg in charges)]
+            tag_lines = ['DATA (PARTICLE_TAG_D(1, IPOS), IPOS=1, NEXTERNAL) / %s /' \
+                            %  ', '.join(bool_dict[tag] for tag in tags)]
             fks_j_from_i_lines = ['DATA (FKS_J_FROM_I_D(1, %d, JPOS), JPOS = 0, 1)  / 1, %d /' \
                             % (fks_i, fks_j)]
             split_type_lines = [ \
@@ -2930,6 +3015,7 @@
         replace_dict['col_lines'] = '\n'.join(col_lines)
         replace_dict['pdg_lines'] = '\n'.join(pdg_lines)
         replace_dict['charge_lines'] = '\n'.join(charge_lines)
+        replace_dict['tag_lines'] = '\n'.join(tag_lines)
         replace_dict['fks_j_from_i_lines'] = '\n'.join(fks_j_from_i_lines)
         replace_dict['split_type_lines'] = '\n'.join(split_type_lines)
 
 
=== modified file 'madgraph/iolibs/template_files/fks_info.inc'
--- madgraph/iolibs/template_files/fks_info.inc	2014-09-19 17:34:43 +0000
+++ madgraph/iolibs/template_files/fks_info.inc	2021-10-15 11:42:16 +0000
@@ -3,6 +3,7 @@
       INTEGER EXTRA_CNT_D(%(nconfs)d), ISPLITORDER_BORN_D(%(nconfs)d), ISPLITORDER_CNT_D(%(nconfs)d)
       INTEGER FKS_J_FROM_I_D(%(nconfs)d, NEXTERNAL, 0:NEXTERNAL)
       INTEGER PARTICLE_TYPE_D(%(nconfs)d, NEXTERNAL), PDG_TYPE_D(%(nconfs)d, NEXTERNAL)
+      LOGICAL PARTICLE_TAG_D(%(nconfs)d, NEXTERNAL)
       REAL*8 PARTICLE_CHARGE_D(%(nconfs)d, NEXTERNAL)
       LOGICAL SPLIT_TYPE_D(%(nconfs)d, %(nsplitorders)d)
       LOGICAL NEED_COLOR_LINKS_D(%(nconfs)d), NEED_CHARGE_LINKS_D(%(nconfs)d)
@@ -40,3 +41,9 @@
 C     Particle charge:
 C     charge is set 0. with QCD corrections, which is irrelevant
 %(charge_lines)s
+
+C     
+C     Tagged particles:
+C     
+%(tag_lines)s
+
 
=== added file 'madgraph/iolibs/template_files/rescale_alpha_tagged.inc'
--- madgraph/iolibs/template_files/rescale_alpha_tagged.inc	1970-01-01 00:00:00 +0000
+++ madgraph/iolibs/template_files/rescale_alpha_tagged.inc	2021-10-15 11:42:16 +0000
@@ -0,0 +1,40 @@
+double precision function get_rescale_alpha_factor(ntag,qed_pow)
+C returns the power of ratios of alpha needed to compensate 
+C for the presence of external tagged photons
+C It is automatically written, and it detects whether the
+C starting model is in the Gmu or alpha0 scheme.
+C ntag is the number of tagged photons, qed_pow the power of gweak
+C of the contribution considered
+implicit none
+integer ntag, qed_pow
+INCLUDE '../../Source/MODEL/input.inc'
+
+get_rescale_alpha_factor = 1d0
+if (ntag.eq.0) return
+
+get_rescale_alpha_factor = %(rescale_fact)s
+
+return
+end
+
+
+double precision function get_virtual_a0Gmu_conv(qed_pow, ntagph, ivirt, born_wgt)
+C returns the piece to compensate the a0<>Gmu change of scheme for the
+C virtuals (single pole if ivirt = 1, finite if ivirt = 0
+implicit none
+integer qed_pow, ntagph, ivirt
+double precision born_wgt
+INCLUDE '../../Source/MODEL/input.inc'
+
+if (ivirt.eq.1) then
+! single 
+get_virtual_a0Gmu_conv = %(virtual_a0Gmu_conv_pole)s
+else if (ivirt.eq.0) then
+! finite part
+get_virtual_a0Gmu_conv = %(virtual_a0Gmu_conv_finite)s
+else
+write(*,*) 'Error get_virtual_a0Gmu_conv: Invalid ivirt', ivirt
+endif
+
+return
+end
 
=== modified file 'madgraph/various/banner.py'
--- madgraph/various/banner.py	2021-10-14 15:00:59 +0000
+++ madgraph/various/banner.py	2021-10-15 11:42:16 +0000
@@ -4528,6 +4528,7 @@
         """Rules
           e+ e- beam -> lpp:0 ebeam:500  
           p p beam -> set maxjetflavor automatically
+          process with tagged photons -> gamma_is_j = false
         """
 
         # check for beam_id
@@ -4552,13 +4553,23 @@
         if proc_characteristic['ninitial'] == 1:
             #remove all cut
             self.remove_all_cut()
+
+        # check for tagged photons
+        tagged_particles = set()
             
         # Check if need matching
         min_particle = 99
         max_particle = 0
         for proc in proc_def:
+            for leg in proc['legs']:
+                if leg['is_tagged']:
+                    tagged_particles.add(leg['id'])
             min_particle = min(len(proc['legs']), min_particle)
             max_particle = max(len(proc['legs']), max_particle)
+
+        if 22 in tagged_particles:
+            self['gamma_is_j'] = False
+
         matching = False
         if min_particle != max_particle:
             #take one of the process with min_particle
 
=== modified file 'tests/acceptance_tests/test_cmd_amcatnlo.py'
--- tests/acceptance_tests/test_cmd_amcatnlo.py	2020-11-27 13:41:46 +0000
+++ tests/acceptance_tests/test_cmd_amcatnlo.py	2021-10-15 11:42:16 +0000
@@ -700,3 +700,38 @@
         
         result = save_load_object.load_from_file('%s/HTML/results.pkl' % self.path)
         return result[run_name]
+
+
+    def test_generate_taggedph_nloew(self):
+        """test the param_card created is correct"""
+        
+
+        text = """
+        import model loop_qcd_qed_sm_Gmu-a0
+        generate u u~ > !a! !a! [QED]
+        output %s
+        launch NLO
+        set lepphreco False
+        set quarkphreco False
+        """ % (self.path)
+        
+        interface = MGCmd.MasterCmd()
+        interface.no_notification()
+        
+        open(pjoin(self.tmpdir,'cmd'),'w').write(text)
+        
+        
+        os.system('rm -rf %s/RunWeb' % self.path)
+        os.system('rm -rf %s/Events/run_*' % self.path)
+
+        interface.exec_cmd('import command %s' % pjoin(self.tmpdir, 'cmd'))
+        # test the lhe event file exists
+        self.assertTrue(os.path.exists('%s/Events/run_01/summary.txt' % self.path))
+        self.assertTrue(os.path.exists('%s/Events/run_01/run_01_tag_1_banner.txt' % self.path))
+        self.assertTrue(os.path.exists('%s/Events/run_01/res_0.txt' % self.path))
+        self.assertTrue(os.path.exists('%s/Events/run_01/res_1.txt' % self.path))
+        self.assertTrue(os.path.exists('%s/Events/run_01/alllogs_0.html' % self.path))
+        self.assertTrue(os.path.exists('%s/Events/run_01/alllogs_1.html' % self.path))
+
+        check_html_page(self, pjoin(self.path, 'crossx.html'))
+        check_html_page(self, pjoin(self.path, 'HTML', 'run_01', 'results.html'))
 
=== added file 'tests/unit_tests/fks/test_fks_taggedphotons.py'
--- tests/unit_tests/fks/test_fks_taggedphotons.py	1970-01-01 00:00:00 +0000
+++ tests/unit_tests/fks/test_fks_taggedphotons.py	2021-10-15 11:42:16 +0000
@@ -0,0 +1,80 @@
+##############################################################################
+#
+# Copyright (c) 2010 The MadGraph5_aMC@NLO Development team and Contributors
+#
+# This file is a part of the MadGraph5_aMC@NLO project, an application which 
+# automatically generates Feynman diagrams and matrix elements for arbitrary
+# high-energy processes in the Standard Model and beyond.
+#
+# It is subject to the MadGraph5_aMC@NLO license which should accompany this 
+# distribution.
+#
+# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
+#
+################################################################################
+from __future__ import absolute_import
+from __future__ import print_function
+from cmd import Cmd
+from six.moves import zip
+""" Basic test of the command interface """
+
+import unittest
+import madgraph
+import madgraph.interface.master_interface as mgcmd
+import madgraph.interface.extended_cmd as ext_cmd
+import madgraph.interface.amcatnlo_interface as amcatnlocmd
+import os
+import madgraph.fks.fks_helas_objects as fks_helas
+import copy
+import madgraph.iolibs.save_load_object as save_load_object
+import madgraph
+import madgraph.various.misc as misc
+root_path = os.path.split(os.path.dirname(os.path.realpath( __file__ )))[0]
+root_path = os.path.dirname(root_path)
+# root_path is ./tests
+pjoin = os.path.join
+
+class TestAMCatNLOEWTagPh(unittest.TestCase):
+    """ a suite of extra tests for the ew stuff """
+    
+    def setUp(self):
+        self.interface = mgcmd.MasterCmd()
+
+    def test_generate_fks_tagph(self):
+        """check that the generate command works as expected.
+        In particular the correct number of real-emission processes
+        and FKS partners"""
+        cmd_list = [
+            'u u~ > a a [real=QED]',
+            'u u~ > !a! !a! [real=QED]',
+            'u u~ > !2a!  [real=QED]',
+            'u u~ > 2!a! [real=QED]'] 
+
+        nrealproc_list = [9, 3, 3, 3]
+        fks_j_list = [[1,2,3], [1,2], [1,2], [1,2]] # 4 is not there in the first set for symmetry
+
+
+        for cmd, nrealproc, fks_j in \
+                zip(cmd_list, nrealproc_list, fks_j_list):
+            self.interface.do_generate(cmd)
+
+            fksprocess = self.interface._fks_multi_proc['born_processes'][0]
+
+            self.assertEqual(len(fksprocess.real_amps), nrealproc)
+            self.assertEqual(set([r.fks_infos[0]['j'] for r in fksprocess.real_amps]), set(fks_j))
+            
+    def test_invalid_syntax_tag(self):
+
+        cmd = "u u~ > !a! !a!"
+        self.assertRaises(madgraph.InvalidCmd,  self.interface.do_generate, cmd)
+
+        cmd = "u u~ > !z! a [real=QED]"
+        self.assertRaises(madgraph.InvalidCmd,  self.interface.do_generate, cmd)
+
+        cmd = "u u~ > !t! t~ a [real=QED]"
+        self.assertRaises(madgraph.InvalidCmd,  self.interface.do_generate, cmd)
+
+        cmd = "u u~ > !a! a [real=QCD]" # is this valid ? I guess NOT [note this is QCD]
+        self.assertRaises(madgraph.InvalidCmd,  self.interface.do_generate, cmd)
+
+