Merge lp:~maddevelopers/mg5amcnlo/2.6.5_polarization into lp:~mg5core1/mg5amcnlo/2.6.8

[Olivier Mattelaer]

1) support for polarization @LO
2) better identification of zero result at LO

=== modified file 'HELAS/v2xxxx.f'
--- HELAS/v2xxxx.f	2010-10-30 03:26:37 +0000
+++ HELAS/v2xxxx.f	2020-01-09 23:16:36 +0000
@@ -30,11 +30,11 @@
             vc(2) = p(1)/p(0)
             vc(3) = p(2)/p(0)
             vc(4) = p(3)/p(0)
-         else
+      else
             vc(1) = p(0)/vmass
             vc(2) = p(1)/vmass
             vc(3) = p(2)/vmass
             vc(4) = p(3)/vmass
-         endif
+      endif
       return
       end

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2019-11-16 11:59:49 +0000
+++ MadSpin/decay.py	2020-01-09 23:16:36 +0000
@@ -2245,7 +2245,7 @@
             nb_mc_masses=len(indices_for_mc_masses)
 
             p, p_str=self.curr_event.give_momenta(event_map)
-            stdin_text=' %s %s %s %s \n' % ('2', self.options['BW_cut'], self.Ecollider, decay_me['max_weight'])
+            stdin_text=' %s %s %s %s %s \n' % ('2', self.options['BW_cut'], self.Ecollider, decay_me['max_weight'], self.options['frame_id'])
             stdin_text+=p_str
             # here I also need to specify the Monte Carlo Masses
             stdin_text+=" %s \n" % nb_mc_masses
@@ -2373,7 +2373,7 @@
         #no decays for this production mode, run in passthrough mode, only adding the helicities to the events
         nb_mc_masses=0
         p, p_str=self.curr_event.give_momenta(event_map)
-        stdin_text=' %s %s %s %s \n' % ('2', self.options['BW_cut'], self.Ecollider, 1.0)
+        stdin_text=' %s %s %s %s \n' % ('2', self.options['BW_cut'], self.Ecollider, 1.0, self.options['frame_id'])
         stdin_text+=p_str
         # here I also need to specify the Monte Carlo Masses
         stdin_text+=" %s \n" % nb_mc_masses
@@ -3278,8 +3278,7 @@
         """return the max. weight associated with me decay['path']"""
 
         p, p_str=self.curr_event.give_momenta(event_map)
-
-        std_in=" %s  %s %s %s  \n" % ("1",BWcut, self.Ecollider, nbpoints)
+        std_in=" %s  %s %s %s %s \n" % ("1",BWcut, self.Ecollider, nbpoints, self.options['frame_id'])
         std_in+=p_str
         max_weight = self.loadfortran('maxweight',
                                path, std_in)
@@ -3348,7 +3347,7 @@
                 if nb < cut:
                     if key[0]=='full':
                         path=key[1]
-                        end_signal="5 0 0 0 \n"  # before closing, write down the seed 
+                        end_signal="5 0 0 0 0\n"  # before closing, write down the seed 
                         external.stdin.write(end_signal)
                         ranmar_state=external.stdout.readline()
                         ranmar_file=pjoin(path,'ranmar_state.dat')
@@ -3963,7 +3962,7 @@
                     external.terminate()
                     del external
                 elif mode=='full':
-                    stdin_text="5 0 0 0 \n"  # before closing, write down the seed 
+                    stdin_text="5 0 0 0 0\n"  # before closing, write down the seed 
                     external = self.calculator[('full',path)]
                     try:
                         external.stdin.write(stdin_text)
@@ -3990,7 +3989,7 @@
             except Exception:
                 pass
             else:
-                stdin_text="5 0 0 0"
+                stdin_text="5 0 0 0 0"
                 external.stdin.write(stdin_text)
                 external.stdin.close()
                 external.stdout.close()

=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2019-10-16 10:56:33 +0000
+++ MadSpin/interface_madspin.py	2020-01-09 23:16:36 +0000
@@ -75,7 +75,7 @@
         self.add_param('cross_section', {'__type__':0.}, comment="forcing normalization of cross-section after MS (for none/onshell)" )
         self.add_param('new_wgt', 'cross-section' ,allowed=['cross-section', 'BR'], comment="if not consistent number of particles, choose what to do for the weight. (BR: means local according to number of part, cross use the force cross-section")
         self.add_param('input_format', 'auto', allowed=['auto','lhe', 'hepmc', 'lhe_no_banner'])
-        
+        self.add_param('frame_id', 6)
         
     ############################################################################
     ##  Special post-processing of the options                                ## 
@@ -211,7 +211,7 @@
         # Read the banner of the inputfile
         self.events_file = open(os.path.realpath(inputfile))
         self.banner = banner.Banner(self.events_file)
-        
+
 
         # Check the validity of the banner:
         if 'slha' not in self.banner:
@@ -226,14 +226,21 @@
             raise self.InvalidCmd('This event file was already decayed by MS. This is not possible to add to it a second decay')
         
         if 'mgruncard' in self.banner:
+            run_card = self.banner.charge_card('run_card')
             if not self.options['Nevents_for_max_weigth']:
-                nevents = int(self.banner.get_detail('run_card', 'nevents'))
+                nevents = run_card['nevents']
                 N_weight = max([75, int(3*nevents**(1/3))])
                 self.options['Nevents_for_max_weigth'] = N_weight
                 N_sigma = max(4.5, math.log(nevents,7.7))
                 self.options['nb_sigma'] = N_sigma
             if self.options['BW_cut'] == -1:
                 self.options['BW_cut'] = float(self.banner.get_detail('run_card', 'bwcutoff'))
+            
+            if isinstance(run_card, banner.RunCardLO):
+                run_card.update_system_parameter_for_include()
+                self.options['frame_id'] = run_card['frame_id']
+            else:
+                self.options['frame_id'] = 6
         else:
             if not self.options['Nevents_for_max_weigth']:
                 self.options['Nevents_for_max_weigth'] = 75
@@ -313,7 +320,6 @@
                 if model_name in line:
                     final_model = True
                     
-            
                 
     def import_model(self, args):
         """syntax: import model NAME CARD_PATH
@@ -400,6 +406,12 @@
         #if self.model and not self.model['case_sensitive']:
         #    decaybranch = decaybranch.lower()
 
+        if self.options['spinmode'] != 'full' and '{' in decaybranch:
+            if self.options['spinmode'] == 'none':
+                logger.warning("polarization option used with spinmode=none. The polarization definition will be done according to the rest-frame of the decaying particles (which is likely not what you expect).")
+            else:
+                logger.warning("polarization option used with spinmode=onshell. This combination is not validated and is by construction using sub-optimal method which can likely lead to bias in some situation. Use at your own risk.")
+
         decay_process, init_part = self.decay.reorder_branch(decaybranch)
         if not self.list_branches.has_key(init_part):
             self.list_branches[init_part] = []
@@ -585,14 +597,16 @@
         self.check_launch(args)
         for part in self.list_branches.keys():
             if part in self.mg5cmd._multiparticles:
+            
                 if any(pid in self.final_state for pid in self.mg5cmd._multiparticles[part]):
                     break
-            pid = self.mg5cmd._curr_model.get('name2pdg')[part]
-            if pid in self.final_state:
-                break
-#        else:
-#            logger.info("Nothing to decay ...")
-#            return
+            else:
+                pid = self.mg5cmd._curr_model.get('name2pdg')[part]
+                if pid in self.final_state:
+                    break
+        else:
+            logger.info("Nothing to decay ...")
+            return
         
 
         model_line = self.banner.get('proc_card', 'full_model_line')
@@ -1443,8 +1457,10 @@
         for curr_event,production in enumerate(orig_lhe):
             if self.options['fixed_order']:
                 production, counterevt= production[0], production[1:]
-            if curr_event and curr_event % 1000 == 0 and float(str(curr_event)[1:]) ==0:
+            if curr_event and self.efficiency:# and curr_event % 10 == 0 and float(str(curr_event)[1:]) ==0:
                 logger.info("decaying event number %s. Efficiency: %s [%s s]" % (curr_event, 1/self.efficiency, time.time()-start))
+            else:
+                logger.info("next event [%s]", time.time()-start)
             while 1:
                 nb_try +=1
                 decays = self.get_decay_from_file(production, evt_decayfile, nb_event-curr_event)
@@ -1579,6 +1595,7 @@
             if len(all_maxwgt) < i:
                 break
             ave_weight, std_weight = decay_tools.get_mean_sd(all_maxwgt[:i])
+            #misc.sprint(ave_weight, std_weight)
             base_max_weight = max(base_max_weight, 1.05 * (ave_weight+self.options['nb_sigma']*std_weight))
                 
             if all_maxwgt[1] > base_max_weight:

=== modified file 'MadSpin/src/driver.f'
--- MadSpin/src/driver.f	2018-09-19 09:41:33 +0000
+++ MadSpin/src/driver.f	2020-01-09 23:16:36 +0000
@@ -18,7 +18,7 @@
       INTEGER I,J,K
       INTEGER HELSET(NEXTERNAL)
       REAL*8 P(0:3,NEXTERNAL_PROD) 
-      REAL*8 PFULL(0:3,NEXTERNAL), Ptrial(0:3,NEXTERNAL) 
+      REAL*8 PFULL(0:3,NEXTERNAL), Ptrial(0:3,NEXTERNAL), P2(0:3,NEXTERNAL) 
       double precision x(36), Ecollider
       CHARACTER*120 BUFF(NEXTERNAL_PROD)
       integer iforest(2,-nexternal:-1,N_MAX_CG)
@@ -39,6 +39,7 @@
       double precision mean, variance, maxweight,weight,std
       double precision temp
       double precision Pprod(0:3,nexternal_prod)
+
       integer nb_mc_masses, indices_mc_masses(nexternal)
       double precision values_mc_masses(nexternal)
 
@@ -54,6 +55,9 @@
        double precision BWcut, maxBW
        common /to_BWcut/BWcut, maxBW
 
+       integer frame_id
+       common /to_me_frame/frame_id
+
 c Conflicting BW stuff
       integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
       double precision cBW_mass(-nexternal:-1,-1:1),
@@ -94,7 +98,8 @@
        read(56,*) P_seed
        close(56)
        iseed = iseed + P_seed
-
+c       open(unit=999,file="/data/DOC/sujays/divisors")
+c       open(unit=999,file="/home/sujay/Downloads/momentaboost")
 cccccccccccccccccccccccccccccccccccccccccccccccccccc
 c   I. read momenta for the production events
 c
@@ -112,7 +117,7 @@
  
 1     continue
       maxBW=0d0
-      read(*,*) mode,  BWcut, Ecollider, temp
+      read(*,*) mode,  BWcut, Ecollider, temp, frame_id
  
 
       if (mode.eq.1) then    ! calculate the maximum weight
@@ -258,10 +263,32 @@
            !do j=1,nexternal
            !   write (*,*) (pfull(k,j), k=0,3)  
            !enddo
-           call SMATRIX(pfull,M_full)
-           call SMATRIX_PROD(pprod,M_prod)
-c           write(*,*) 'M_full ', M_full
-c           write(*,*) 'jac',jac
+
+c=======================Modified part=================================
+
+          
+           call  boost_to_frame(pfull, frame_id, P2)
+           call SMATRIX(P2,M_full)
+
+
+c           call SMATRIX_PROD(pprod,M_prod)
+c           write(999,*) (M_prod)
+c           write(999,*) "---------after-boost------"
+           call  boost_to_frame_prod(pprod, frame_id,nexternal_prod, P2)
+           call SMATRIX_PROD(P2,M_prod)
+
+
+
+c          do j=1,nexternal
+c              write (999,*) (P2(k,j), k=0,3)
+c           enddo
+c              write(999,*) (M_prod)
+c           do j=1,nexternal_prod
+c              write (999,*) (pprod(k,j), k=0,3)
+c           enddo
+c           close(999)
+
+c==========================================================================
 
            weight=M_full*jac/M_prod
            if (weight.gt.maxweight) then
@@ -335,7 +362,6 @@
            do i = 1, 3*(nexternal-nexternal_prod)+1
               call  ntuple(x(i),0d0,1d0,i,1)
            enddo
-
            call generate_momenta_conf(jac,x,itree,qmass,qwidth,pfull,pprod,map_external2res) 
            if (jac.lt.0d0) then
              counter2=counter2+1 
@@ -366,9 +392,32 @@
 
              cycle
            endif
-           call SMATRIX(pfull,M_full)
-           call SMATRIX_PROD(pprod,M_prod)
-
+
+c=======================Modified part=================================
+
+
+           call  boost_to_frame(pfull, frame_id, P2)
+           call SMATRIX(P2,M_full)
+
+
+c           call SMATRIX_PROD(pprod,M_prod)
+c           write(999,*) (M_prod)
+c           write(999,*) "---------after-boost------"
+           call  boost_to_frame_prod(pprod, frame_id,nexternal_prod, P2)
+           call SMATRIX_PROD(P2,M_prod)
+
+
+
+c          do j=1,nexternal
+c              write (999,*) (P2(k,j), k=0,3)
+c           enddo
+c              write(999,*) (M_prod)
+c           do j=1,nexternal_prod
+c              write (999,*) (pprod(k,j), k=0,3)
+c           enddo
+c           close(999)
+
+c==========================================================================
            weight=M_full*jac/M_prod
 
            if (weight.gt.x(3*(nexternal-nexternal_prod)+1)*maxweight) notpass=.false.
@@ -1892,3 +1941,162 @@
       END
 
 
+
+
+c======================================================================
+C     Subroutine to return momenta in a dedicated frame
+C     frame_id is the tag of the particle to put at rest
+C     frame_id follow the convention of cluster.f (sum 2**(N-1))
+C     -----------------------------------------
+
+      subroutine boost_to_frame(P1, frame_id, P2)
+
+      implicit none
+
+      include 'nexternal.inc'
+
+      DOUBLE PRECISION P1(0:3,NEXTERNAL)
+      DOUBLE PRECISION P2(0:3,NEXTERNAL)
+      DOUBLE PRECISION PBOOST(0:3)
+      integer frame_id
+
+      integer ids(nexternal)
+      integer i,j
+      Pboost(0)= 0d0
+      Pboost(1)= 0d0
+      Pboost(2)= 0d0
+      Pboost(3)= 0d0
+
+c     uncompress
+      call mapid(frame_id, ids)
+c     find the boost momenta --sum of particles--
+      do i=1,nexternal
+       if (ids(i).eq.1)then
+c             write (999,*) (Pboost(j), j=0,3)
+            do j=0,3
+                   Pboost(j) = Pboost(j) + P1(j,i)
+            enddo
+         endif
+      enddo
+      do j=1,3
+          Pboost(j) = -1 * Pboost(j)
+      enddo
+       
+c              write (999,*) (Pboost(j), j=0,3)
+      do i=1, nexternal
+         call boostx(p1(0,i), pboost, p2(0,i))
+      enddo
+      return
+      end
+
+
+
+
+c========================================================================
+
+
+      subroutine mapid(id,ids)
+c**************************************************************************
+c     input:
+c            id        compressed particle id
+c            ids       array of particle ids
+c**************************************************************************
+      implicit none
+      include 'nexternal.inc'
+      integer i, icd, id, ids(nexternal)
+
+      icd=id
+      do i=1,nexternal
+         ids(i)=0
+         if (btest(id,i)) then
+            ids(i)=1
+         endif
+c         write(*,*) 'cluster.f: uncompressed code ',i,' is ',ids(i)
+      enddo
+
+      return
+      end
+
+
+
+
+C     -----------------------------------------
+
+      subroutine boost_to_frame_prod(P1, frame_id, nexternal_prod, P2)
+
+      implicit none
+
+      include 'nexternal.inc'
+
+      DOUBLE PRECISION P1(0:3,*)
+      DOUBLE PRECISION P2(0:3,*)
+      DOUBLE PRECISION PBOOST(0:3)
+      integer frame_id, nexternal_prod
+
+      integer ids(nexternal_prod)
+      integer i,j
+      Pboost(0)= 0d0
+      Pboost(1)= 0d0
+      Pboost(2)= 0d0
+      Pboost(3)= 0d0
+
+c     uncompress
+      call mapid_prod(frame_id, ids,nexternal_prod)
+c     find the boost momenta --sum of particles--
+      do i=1,nexternal_prod
+       if (ids(i).eq.1)then
+c             write (999,*) (Pboost(j), j=0,3)
+            do j=0,3
+                   Pboost(j) = Pboost(j) + P1(j,i)
+            enddo
+         endif
+      enddo
+      do j=1,3
+          Pboost(j) = -1 * Pboost(j)
+      enddo
+
+c              write (999,*) (Pboost(j), j=0,3)
+      do i=1, nexternal_prod
+         call boostx(p1(0,i), pboost, p2(0,i))
+      enddo
+      return
+      end
+
+
+
+
+
+
+
+
+      subroutine mapid_prod(id,ids,nexternal_prod)
+c**************************************************************************
+c     input:
+c            id        compressed particle id
+c            ids       array of particle ids
+c**************************************************************************
+      implicit none
+      include 'nexternal.inc'
+      integer i, icd, id, nexternal_prod, ids(nexternal_prod)
+
+      icd=id
+      do i=1,nexternal_prod
+         ids(i)=0
+         if (btest(id,i)) then
+            ids(i)=1
+         endif
+c         write(*,*) 'cluster.f: uncompressed code ',i,' is ',ids(i)
+      enddo
+
+      return
+      end
+
+
+
+
+
+
+
+
+
+

=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2019-05-01 19:53:40 +0000
+++ Template/LO/Cards/run_card.dat	2020-01-09 23:16:36 +0000
@@ -39,6 +39,7 @@
      %(ebeam1)s     = ebeam1  ! beam 1 total energy in GeV
      %(ebeam2)s     = ebeam2  ! beam 2 total energy in GeV
 $beam_pol
+$frame
 #*********************************************************************
 # PDF CHOICE: this automatically fixes also alpha_s and its evol.    *
 #*********************************************************************

=== modified file 'Template/LO/Source/dsample.f'
--- Template/LO/Source/dsample.f	2019-01-25 10:01:29 +0000
+++ Template/LO/Source/dsample.f	2020-01-09 23:16:36 +0000
@@ -134,6 +134,7 @@
             ievent=ievent+1
             call x_to_f_arg(ndim,ipole,mincfig,maxcfig,ninvar,wgt,x,p)
             if (pass_point(p)) then
+               nb_pass_cuts = nb_pass_cuts + 1
                fx = dsig(p,wgt,0) !Evaluate function
                wgt = wgt*fx
                if (wgt .ne. 0d0) call graph_point(p,wgt) !Update graphs
@@ -316,6 +317,7 @@
                fx =0d0
                wgt=0d0
             endif
+            
             if (nzoom .le. 0) then
                call sample_put_point(wgt,x(1),iter,ipole,itmin) !Store result
             else
@@ -705,6 +707,7 @@
       tmean = 0d0
       trmean = 0d0
       tsigma = 0d0
+      nb_pass_cuts = 0
       kn = 0
       cur_it = 1
       do j=1,ng
@@ -1756,6 +1759,11 @@
          if (kn .ge. max_events .and. non_zero .le. 5) then
             call none_pass(max_events)
          endif
+         if (iteration.eq.1) then
+           if (nb_pass_cuts.eq.1000 .and. non_zero.eq.0) then
+              call none_pass(1000)
+           endif
+         endif
          if (non_zero .ge. events .or. (kn .gt. 200*events .and.
      $        non_zero .gt. 5)) then
 
@@ -2267,8 +2275,15 @@
 c----
 c  Begin Code
 c----
-      write(*,*) 'No points passed cuts!'
-      write(*,*) 'Loosen cuts or increase max_events',max_events
+      if (1000.eq.max_events) then
+         write(*,*) nb_pass_cuts, 
+     &    ' points passed the cut but all returned zero'
+         write(*,*) 'therefore considering this contribution as zero'
+      else
+         write(*,*) 'No points passed cuts!'
+         write(*,*) 'Loosen cuts or increase max_events',max_events
+      endif
+
 c      open(unit=22,file=result_file,status='old',access='append',
 c     &           err=23)
 c      write(22,222) 'Iteration',0,'Mean: ',0d0,
@@ -2540,6 +2555,7 @@
 
       write(*,*) "RESET CUMULATIVE VARIABLE"
       non_zero = 0
+      nb_pass_cuts = 0
       do j=1,maxinvar
          do i=1,ng -1
             inon_zero = 0

=== modified file 'Template/LO/Source/genps.inc'
--- Template/LO/Source/genps.inc	2016-04-02 14:23:21 +0000
+++ Template/LO/Source/genps.inc	2020-01-09 23:16:36 +0000
@@ -38,3 +38,8 @@
       double precision hel_jacobian
 	  INTEGER HEL_PICKED
 	  COMMON/HEL_PICKED/HEL_PICKED,hel_jacobian
+c**********************************************************************************
+c       parameters to control is zero-cross section is due to PS or matrix-element
+c**********************************************************************************
+      integer nb_pass_cuts
+      common/cut_efficiency/nb_pass_cuts
\ No newline at end of file

=== modified file 'Template/LO/Source/run.inc'
--- Template/LO/Source/run.inc	2018-10-07 20:06:38 +0000
+++ Template/LO/Source/run.inc	2020-01-09 23:16:36 +0000
@@ -16,6 +16,9 @@
       integer                                        lpp(2)
       double precision    ebeam(2), xbk(2),q2fact(2)
       common/to_collider/ ebeam   , xbk   ,q2fact,   lpp
+
+      integer frame_id
+      common/to_frame_me/frame_id
 c
 c     Number of events
 c

=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f	2019-06-17 14:27:09 +0000
+++ Template/LO/SubProcesses/genps.f	2020-01-09 23:16:36 +0000
@@ -1264,5 +1264,47 @@
       END
       
 
+C     -----------------------------------------
+C     Subroutine to return momenta in a dedicated frame
+C     frame_id is the tag of the particle to put at rest
+C     frame_id follow the convention of cluster.f (sum 2**(N-1))
+C     -----------------------------------------
+
+      subroutine boost_to_frame(P1, frame_id, P2)
+
+      implicit none
+
+      include 'nexternal.inc'
+
+      DOUBLE PRECISION P1(0:3,NEXTERNAL)
+      DOUBLE PRECISION P2(0:3,NEXTERNAL)
+      DOUBLE PRECISION PBOOST(0:3)
+      integer frame_id
+
+      integer ids(nexternal)
+      integer i,j
+
+c     uncompress
+      call mapid(frame_id, ids)
+      pboost(:) = 0d0
+      p2(:,:) = 0d0
+c     find the boost momenta --sum of particles--
+      do i=1,nexternal
+       if (ids(i).eq.1)then
+            do j=0,3
+	           Pboost(j) = Pboost(j) + P1(j,i)
+            enddo
+         endif
+      enddo
+      do j=1,3	
+          Pboost(j) = -1 * Pboost(j)
+      enddo	    
+      do i=1, nexternal
+         call boostx(p1(0,i), pboost, p2(0,i))
+      enddo   
+      return
+      end
+
+
 
 

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2020-01-09 09:35:56 +0000
+++ UpdateNotes.txt	2020-01-09 23:16:36 +0000
@@ -1,8 +1,16 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
+
+
 2.7.0()
      OM: Allow for a new syntax:
          generate p p > 3j replaces p p > j j j 
+      OM: Allow syntax for (fully) polarized particles at LO:
+           ex: p p > w+{0} j, e+{L} e-{R} > mu+ mu- Z{T}
+               p p > w+{T} w-{0} j j, w+ > e+ ve, w- > j j
+      OM: distinguish in the code if zero contribution are related to no point passing cuts or if
+          they are related to vanishing matrix-element. In the later case, allow for a lower threshold.
+          (This allow to fasten the computation of such zero contribution)
 
 
 2.6.7(16/10/19)
@@ -30,7 +38,6 @@
       OM: Energy cut (at LO) are now hidden by default for LHC type of run but visible for lepton collider ones.
       OM: Same for beam polarization
 
-
 2.6.5 (03/02/19)
       OM: Fix some speed issue with the generated gridpack --speed issue introduced in 2.6.1--
       OM: Fix a bug in the computation of systematics when running with Python 2.6.

=== modified file 'aloha/aloha_object.py'
--- aloha/aloha_object.py	2018-04-26 18:04:25 +0000
+++ aloha/aloha_object.py	2020-01-09 23:16:36 +0000
@@ -69,8 +69,140 @@
     def get_unique_name(self, lorentz1, particle):
         return '_P^%s_%s' % (particle, lorentz1)
 
-
-
+#===============================================================================
+# P (Momenta)
+#===============================================================================
+class L_PVec(aloha_lib.LorentzObject):
+    """ Helas Object for an Impulsion """
+    
+    contract_first = 1
+    
+    def __init__(self, name, lorentz1, particle):
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self, name,[lorentz1], [],['PVec%s'%particle])
+        aloha_lib.KERNEL.add_tag((name,))
+    
+    def create_representation(self):
+        #self.sub0 = aloha_lib.DVariable('P%s_0' % self.particle)
+        self.sub1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        self.sub2 = aloha_lib.DVariable('P%s_2' % self.particle)
+        self.sub3 = aloha_lib.DVariable('P%s_3' % self.particle)
+
+        self.representation= aloha_lib.LorentzObjectRepresentation(
+                                    {(0,): 0, (1,): self.sub1, \
+                                     (2,): self.sub2, (3,): self.sub3},                              
+                                    self.lorentz_ind, [])
+
+
+class PVec(aloha_lib.FactoryLorentz):
+    """ Helas Object for an Impulsion """
+    
+    object_class = L_PVec
+    
+    #def __init__(self, lorentz1, particle):
+    @classmethod
+    def get_unique_name(self, lorentz1, particle):
+        return '_Pvec^%s_%s' % (particle, lorentz1)
+
+
+    
+#===============================================================================
+# Three Norm
+#===============================================================================
+class L_Tnorm(aloha_lib.LorentzObject):
+    """ Helas Object for DSQRT(Px^2+Py^2+pz^2) """
+      
+    def __init__(self, name, particle):
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self, name, [], [], tags=['|p%s|' % particle])
+    
+    def create_representation(self):
+        mass = aloha_lib.DVariable('Tnorm%s' % self.particle)
+
+        self.representation = aloha_lib.LorentzObjectRepresentation(
+                                mass, self.lorentz_ind, self.spin_ind)
+
+class Tnorm(aloha_lib.FactoryLorentz):
+    
+    object_class = L_Tnorm
+    
+    @classmethod
+    def get_unique_name(self, particle):
+        return '_tnorm_%s' % particle
+
+#===============================================================================
+# |p| + pz
+#===============================================================================
+class L_TnormZ(aloha_lib.LorentzObject):
+    """ Helas Object for DSQRT(Px^2+Py^2+pz^2) + PZ """
+      
+    def __init__(self, name, particle):
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self, name, [], [], tags=['|p%s|_z' % particle])
+    
+    def create_representation(self):
+        mass = aloha_lib.DVariable('TnormZ%s' % self.particle)
+
+        self.representation = aloha_lib.LorentzObjectRepresentation(
+                                mass, self.lorentz_ind, self.spin_ind)
+
+class TnormZ(aloha_lib.FactoryLorentz):
+    
+    object_class = L_TnormZ
+    
+    @classmethod
+    def get_unique_name(self, particle):
+        return '_tnormz_%s' % particle
+    
+#===============================================================================
+# DSQRT(E+- |p|)
+#===============================================================================
+class L_FermionWP(aloha_lib.LorentzObject):
+    """ Helas Object for DSQRT(Px^2+Py^2+pz^2) """
+      
+    def __init__(self, name, particle):
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self, name, [], [], tags=['FWP_%s' % particle])
+    
+    def create_representation(self):
+        mass = aloha_lib.DVariable('FWP%s' % self.particle)
+
+        self.representation = aloha_lib.LorentzObjectRepresentation(
+                                mass, self.lorentz_ind, self.spin_ind)
+
+class FermionWP(aloha_lib.FactoryLorentz):
+    
+    object_class = L_FermionWP
+    
+    @classmethod
+    def get_unique_name(self, particle):
+        return '_FermionWP_%s' % particle
+    
+#===============================================================================
+# DSQRT(E+- |p|)
+#===============================================================================
+class L_FermionWM(aloha_lib.LorentzObject):
+    """ Helas Object for DSQRT(Px^2+Py^2+pz^2) """
+      
+    def __init__(self, name, particle):
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self, name, [], [], tags=['FWM_%s' % particle])
+    
+    def create_representation(self):
+        mass = aloha_lib.DVariable('FWM%s' % self.particle)
+
+        self.representation = aloha_lib.LorentzObjectRepresentation(
+                                mass, self.lorentz_ind, self.spin_ind)
+
+class FermionWM(aloha_lib.FactoryLorentz):
+    
+    object_class = L_FermionWM
+    
+    @classmethod
+    def get_unique_name(self, particle):
+        return '_FermionWM_%s' % particle
+
+    
 #===============================================================================
 # Pslash
 #===============================================================================
@@ -219,6 +351,45 @@
     def get_unique_name(self, particle):
         return '_OM2_%s' % particle
 
+
+#===============================================================================
+# P (Momenta)
+#===============================================================================
+class L_PT(aloha_lib.LorentzObject):
+    """ Helas Object for an Impulsion """
+    
+    contract_first = 1
+    
+    def __init__(self, name, lorentz1, particle):
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self, name,[lorentz1], [],['PT%s'%particle])
+        aloha_lib.KERNEL.add_tag((name,))
+    
+    def create_representation(self):
+        #self.sub0 = aloha_lib.DVariable('P%s_0' % self.particle)
+        self.sub1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        self.sub2 = aloha_lib.DVariable('P%s_2' % self.particle)
+        #self.sub3 = aloha_lib.DVariable('P%s_3' % self.particle)
+
+        self.representation= aloha_lib.LorentzObjectRepresentation(
+                                    {(0,): 0, (1,): self.sub1, \
+                                     (2,): self.sub2, (3,): 0},                              
+                                    self.lorentz_ind, [])
+
+
+class PT(aloha_lib.FactoryLorentz):
+    """ Helas Object for an Impulsion """
+    
+    object_class = L_PT
+    
+    #def __init__(self, lorentz1, particle):
+    @classmethod
+    def get_unique_name(self, lorentz1, particle):
+        return '_PT^%s_%s' % (particle, lorentz1)
+
+
+
+
 #===============================================================================
 # Width
 #===============================================================================
@@ -928,7 +1099,513 @@
     def get_unique_name(self, spin1, spin2):
         
         return 'PROJP_%s_%s' % (spin1, spin2)
-        
+
+#===============================================================================
+# EPSL (longitudinal part of the epsilon normmaized to avoid denominator)
+#===============================================================================
+class L_EPSL(aloha_lib.LorentzObject):
+    """ eps^mu_L (polarization vector) 
+        (kx**2+ky**2+kz**2, E*kz, E*ky, E*kz)
+        the normalization 1/m *|\vec{k} is NOT include
+    
+    """
+    
+    def __init__(self, name, lor1, particle):
+        
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self,name, [lor1], [])
+    
+    def create_representation(self):
+        """create representation"""
+        E = aloha_lib.DVariable('P%s_0' % self.particle)
+        m = aloha_lib.DVariable('M%s' % self.particle)
+        p1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        p2 = aloha_lib.DVariable('P%s_2' % self.particle)
+        p3 = aloha_lib.DVariable('P%s_3' % self.particle)    
+
+    
+        norm = (p1**2+p2**2+p3**2)
+        eps = {
+            (0,): norm,
+            (1,): p1*E,
+            (2,): p2*E ,
+            (3,): p3*E ,
+            }
+#        eps = {
+#            (0,): norm-E*E,
+#            (1,): 0,
+#            (2,): 0 ,
+#            (3,): 0 ,
+#            }
+    
+        self.representation = aloha_lib.LorentzObjectRepresentation(eps,
+                                self.lorentz_ind,self.spin_ind)
+
+class EPSL(aloha_lib.FactoryLorentz):
+
+    object_class = L_EPSL
+    
+    @classmethod
+    def get_unique_name(self, lor1, particle):
+        return '_EPSL%s_%s_' % (particle, lor1)
+
+#===============================================================================
+# EPSTR (transverse part of the epsilon normmalized to avoid denominator -- one of the eigenstate)
+#        Real part
+#===============================================================================
+class L_EPST1(aloha_lib.LorentzObject):
+    """ eps^mu_1 (polarization vector) 
+        (0, kx kz  , ky kz,  -k_T^2)
+        the normalization 1/k_T / norm(k) is NOT include
+    
+    """
+    
+    def __init__(self, name, lor1, particle):
+        
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self,name, [lor1], [])
+    
+    def create_representation(self):
+        """create representation"""
+#        E = aloha_lib.DVariable('P%s_0' % self.particle)
+        p1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        p2 = aloha_lib.DVariable('P%s_2' % self.particle)
+        p3 = aloha_lib.DVariable('P%s_3' % self.particle)    
+
+    
+        kT = (p1**2+p2**2)
+
+        eps = {
+            (0,): 0,
+            (1,): (p1)*(p3),
+            (2,): p2*(p3),
+            (3,):  -kT ,
+            }
+#        eps = {
+#            (0,): norm-E*E,
+#            (1,): 0,
+#            (2,): 0 ,
+#            (3,): 0 ,
+#            }
+    
+        self.representation = aloha_lib.LorentzObjectRepresentation(eps,
+                                self.lorentz_ind,self.spin_ind)
+
+class EPST1(aloha_lib.FactoryLorentz):
+
+    object_class = L_EPST1
+    
+    @classmethod
+    def get_unique_name(self, lor1, particle):
+        return '_EPST1%s_%s_' % (particle, lor1)
+
+#===============================================================================
+# EPSTI (transverse part of the epsilon normmalized to avoid denominator -- one of the eigenstate)
+#        Imaginary part
+#===============================================================================
+
+class L_EPST2(aloha_lib.LorentzObject):
+    """ eps^mu_+ (polarization vector) 
+        (0,  ky k_T ,  -i kx k_T, 0)
+ 
+        the normalization 1/(kT*norm(k)) is NOT include here
+        as well as the k_T term (to avoid the square-root)
+    
+    """
+    
+    def __init__(self, name, lor1, particle):
+        
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self,name, [lor1], [])
+    
+    def create_representation(self):
+        """create representation"""
+ #       E = aloha_lib.DVariable('P%s_0' % self.particle)
+ #       m = aloha_lib.DVariable('M%s' % self.particle)
+        p1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        p2 = aloha_lib.DVariable('P%s_2' % self.particle)
+ #       p3 = aloha_lib.DVariable('P%s_3' % self.particle)    
+
+
+        eps = {
+            (0,): 0,
+            (1,): -p2,
+            (2,): (p1),
+            (3,):  0 ,
+            }
+#        eps = {
+#            (0,): norm-E*E,
+#            (1,): 0,
+#            (2,): 0 ,
+#            (3,): 0 ,
+#            }
+    
+        self.representation = aloha_lib.LorentzObjectRepresentation(eps,
+                                self.lorentz_ind,self.spin_ind)
+        
+class EPST2(aloha_lib.FactoryLorentz):
+
+    object_class = L_EPST2
+    
+    @classmethod
+    def get_unique_name(self, lor1, particle):
+        return '_EPST2%s_%s_' % (particle, lor1)
+#===============================================================================
+# UFP U fermion plus component (first two component)
+#       the denominator is include in the denominator of the propagator
+#===============================================================================
+fsign = -1
+    
+class L_UFP(aloha_lib.LorentzObject):
+    """ 
+
+    """
+    
+    def __init__(self, name, spin1, particle):
+        
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self,name, [], [spin1],
+                                         tags=['|p%s|_z' % particle,
+                                               'FWP%s' % particle,
+                                               'FWM%s' % particle])
+    
+    def create_representation(self):
+        """create representation"""
+        #E = aloha_lib.DVariable('P%s_0' % self.particle)
+        p1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        p2 = aloha_lib.DVariable('P%s_2' % self.particle)
+        #p3 = aloha_lib.DVariable('P%s_3' % self.particle)    
+        normz = aloha_lib.DVariable('TnormZ%s' % self.particle)
+        fwp = aloha_lib.DVariable('FWP%s' % self.particle)
+        fwm = aloha_lib.DVariable('FWM%s' % self.particle)
+    
+
+        u = {
+            (0,): fwm * normz,
+            (1,): fwm * (fsign * p1 + complex(0,fsign)*p2),
+            (2,): fwp * normz,
+            (3,): fwp * (fsign* p1 + complex(0,fsign)*p2),
+            }
+
+    
+        self.representation = aloha_lib.LorentzObjectRepresentation(u,
+                                self.lorentz_ind,self.spin_ind)
+
+
+class UFP(aloha_lib.FactoryLorentz):
+
+    object_class = L_UFP
+    
+    @classmethod
+    def get_unique_name(self, spin1, particle):
+        return '_UFP%s_%s_' % (particle, spin1)
+
+class L_UFM(aloha_lib.LorentzObject):
+    """ 
+
+    """
+    
+    def __init__(self, name, spin1, particle):
+        
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self,name, [], [spin1],
+                                         tags=['|p%s|_z' % particle,
+                                               'FWP%s' % particle,
+                                               'FWM%s' % particle])
+    
+    def create_representation(self):
+        """create representation"""
+        #E = aloha_lib.DVariable('P%s_0' % self.particle)
+        p1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        p2 = aloha_lib.DVariable('P%s_2' % self.particle)
+        #p3 = aloha_lib.DVariable('P%s_3' % self.particle)    
+        normz = aloha_lib.DVariable('TnormZ%s' % self.particle)
+        fwp = aloha_lib.DVariable('FWP%s' % self.particle)
+        fwm = aloha_lib.DVariable('FWM%s' % self.particle)
+    
+
+        u = {
+            (0,): fwp * (-1*fsign*p1 + complex(0,fsign)*p2),
+            (1,): fwp * normz,
+            (2,): fwm * (-1*fsign*p1 + complex(0,fsign)*p2),
+            (3,): fwm * normz,
+            }
+
+    
+        self.representation = aloha_lib.LorentzObjectRepresentation(u,
+                                self.lorentz_ind,self.spin_ind)
+        
+class UFM(aloha_lib.FactoryLorentz):
+
+    object_class = L_UFM
+    
+    @classmethod
+    def get_unique_name(self, spin1, particle):
+        return '_UFM%s_%s_' % (particle, spin1)
+   
+class L_UFPC(aloha_lib.LorentzObject):
+    """ 
+
+    """
+    
+    def __init__(self, name, spin1, particle):
+        
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self,name, [], [spin1],
+                                         tags=['|p%s|_z' % particle,
+                                               'FWP%s' % particle,
+                                               'FWM%s' % particle])
+    
+    def create_representation(self):
+        """create representation"""
+        #E = aloha_lib.DVariable('P%s_0' % self.particle)
+        p1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        p2 = aloha_lib.DVariable('P%s_2' % self.particle)
+        #p3 = aloha_lib.DVariable('P%s_3' % self.particle)    
+        normz = aloha_lib.DVariable('TnormZ%s' % self.particle)
+        fwp = aloha_lib.DVariable('FWP%s' % self.particle)
+        fwm = aloha_lib.DVariable('FWM%s' % self.particle)
+    
+
+        u = {
+            (0,): fwp * normz,
+            (1,): fwp * (fsign* p1 + complex(0,-fsign)*p2),
+            (2,): fwm * normz,
+            (3,): fwm * (fsign* p1 + complex(0,-fsign)*p2),
+            }
+
+    
+        self.representation = aloha_lib.LorentzObjectRepresentation(u,
+                                self.lorentz_ind,self.spin_ind)
+
+class UFPC(aloha_lib.FactoryLorentz):
+
+    object_class = L_UFPC
+    
+    @classmethod
+    def get_unique_name(self, spin1, particle):
+        return '_UFPC%s_%s_' % (particle, spin1)
+
+class L_UFMC(aloha_lib.LorentzObject):
+    """ 
+
+    """
+    
+    def __init__(self, name, spin1, particle):
+        
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self,name, [], [spin1],
+                                         tags=['|p%s|_z' % particle,
+                                               'FWP%s' % particle,
+                                               'FWM%s' % particle])
+    
+    def create_representation(self):
+        """create representation"""
+        #E = aloha_lib.DVariable('P%s_0' % self.particle)
+        p1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        p2 = aloha_lib.DVariable('P%s_2' % self.particle)
+        #p3 = aloha_lib.DVariable('P%s_3' % self.particle)    
+        normz = aloha_lib.DVariable('TnormZ%s' % self.particle)
+        fwp = aloha_lib.DVariable('FWP%s' % self.particle)
+        fwm = aloha_lib.DVariable('FWM%s' % self.particle)
+    
+
+        u = {
+            (0,): fwm * (-fsign*p1 + complex(0,-fsign)*p2),
+            (1,): fwm * normz,
+            (2,): fwp * (-fsign*p1 + complex(0,-fsign)*p2),
+            (3,): fwp * normz,
+            }
+
+    
+        self.representation = aloha_lib.LorentzObjectRepresentation(u,
+                                self.lorentz_ind,self.spin_ind)
+
+class UFMC(aloha_lib.FactoryLorentz):
+
+    object_class = L_UFMC
+    
+    @classmethod
+    def get_unique_name(self, spin1, particle):
+        return '_UFMC%s_%s_' % (particle, spin1)
+
+class L_VFP(aloha_lib.LorentzObject):
+    """ 
+
+    """
+    
+    def __init__(self, name, spin1, particle):
+        
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self,name, [], [spin1],
+                                         tags=['|p%s|_z' % particle,
+                                               'FWP%s' % particle,
+                                               'FWM%s' % particle])
+    
+    def create_representation(self):
+        """create representation"""
+        #E = aloha_lib.DVariable('P%s_0' % self.particle)
+        p1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        p2 = aloha_lib.DVariable('P%s_2' % self.particle)
+        #p3 = aloha_lib.DVariable('P%s_3' % self.particle)    
+        normz = aloha_lib.DVariable('TnormZ%s' % self.particle)
+        fwp = aloha_lib.DVariable('FWP%s' % self.particle)
+        fwm = aloha_lib.DVariable('FWM%s' % self.particle)
+    
+
+        u = {
+            (0,): -1*fwp * (-1 *fsign * p1 + complex(0,fsign)*p2),
+            (1,): -1*fwp * normz,
+            (2,): 1*fwm * (-1 *fsign * p1 + complex(0,fsign)*p2),
+            (3,): 1*fwm * normz,
+            }
+
+    
+        self.representation = aloha_lib.LorentzObjectRepresentation(u,
+                                self.lorentz_ind,self.spin_ind)
+
+
+class VFP(aloha_lib.FactoryLorentz):
+
+    object_class = L_VFP
+    
+    @classmethod
+    def get_unique_name(self, spin1, particle):
+        return '_VFP%s_%s_' % (particle, spin1)
+
+class L_VFPC(aloha_lib.LorentzObject):
+    """ 
+
+    """
+    
+    def __init__(self, name, spin1, particle):
+        
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self,name, [], [spin1],
+                                         tags=['|p%s|_z' % particle,
+                                               'FWP%s' % particle,
+                                               'FWM%s' % particle])
+    
+    def create_representation(self):
+        """create representation"""
+        #E = aloha_lib.DVariable('P%s_0' % self.particle)
+        p1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        p2 = aloha_lib.DVariable('P%s_2' % self.particle)
+        #p3 = aloha_lib.DVariable('P%s_3' % self.particle)    
+        normz = aloha_lib.DVariable('TnormZ%s' % self.particle)
+        fwp = aloha_lib.DVariable('FWP%s' % self.particle)
+        fwm = aloha_lib.DVariable('FWM%s' % self.particle)
+    
+
+        u = {
+            (0,): fwm * (-1 *fsign * p1 + complex(0,-1*fsign)*p2),
+            (1,): fwm * normz,
+            (2,): -1*fwp * (-1 *fsign * p1 + complex(0,-1*fsign)*p2),
+            (3,): -1*fwp * normz,
+            }
+
+    
+        self.representation = aloha_lib.LorentzObjectRepresentation(u,
+                                self.lorentz_ind,self.spin_ind)
+
+
+class VFPC(aloha_lib.FactoryLorentz):
+
+    object_class = L_VFPC
+    
+    @classmethod
+    def get_unique_name(self, spin1, particle):
+        return '_VFPC%s_%s_' % (particle, spin1)
+
+
+class L_VFM(aloha_lib.LorentzObject):
+    """ 
+
+    """
+    
+    def __init__(self, name, spin1, particle):
+        
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self,name, [], [spin1],
+                                         tags=['|p%s|_z' % particle,
+                                               'FWP%s' % particle,
+                                               'FWM%s' % particle])
+    
+    def create_representation(self):
+        """create representation"""
+        #E = aloha_lib.DVariable('P%s_0' % self.particle)
+        p1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        p2 = aloha_lib.DVariable('P%s_2' % self.particle)
+        #p3 = aloha_lib.DVariable('P%s_3' % self.particle)    
+        normz = aloha_lib.DVariable('TnormZ%s' % self.particle)
+        fwp = aloha_lib.DVariable('FWP%s' % self.particle)
+        fwm = aloha_lib.DVariable('FWM%s' % self.particle)
+    
+
+        u = {
+            (0,): fwm * normz,
+            (1,): fwm * (fsign * p1 + complex(0,fsign)*p2),
+            (2,): -1*fwp * normz,
+            (3,): -1*fwp * (fsign * p1 + complex(0,fsign)*p2),
+            }
+
+    
+        self.representation = aloha_lib.LorentzObjectRepresentation(u,
+                                self.lorentz_ind,self.spin_ind)
+
+
+class VFM(aloha_lib.FactoryLorentz):
+
+    object_class = L_VFM
+    
+    @classmethod
+    def get_unique_name(self, spin1, particle):
+        return '_VFM%s_%s_' % (particle, spin1)
+
+class L_VFMC(aloha_lib.LorentzObject):
+    """ 
+
+    """
+    
+    def __init__(self, name, spin1, particle):
+        
+        self.particle = particle
+        aloha_lib.LorentzObject.__init__(self,name, [], [spin1],
+                                         tags=['|p%s|_z' % particle,
+                                               'FWP%s' % particle,
+                                               'FWM%s' % particle])
+    
+    def create_representation(self):
+        """create representation"""
+        #E = aloha_lib.DVariable('P%s_0' % self.particle)
+        p1 = aloha_lib.DVariable('P%s_1' % self.particle)
+        p2 = aloha_lib.DVariable('P%s_2' % self.particle)
+        #p3 = aloha_lib.DVariable('P%s_3' % self.particle)    
+        normz = aloha_lib.DVariable('TnormZ%s' % self.particle)
+        fwp = aloha_lib.DVariable('FWP%s' % self.particle)
+        fwm = aloha_lib.DVariable('FWM%s' % self.particle)
+    
+
+        u = {
+            (0,): -1*fwp * normz,
+            (1,): -1*fwp * (fsign * p1 + complex(0,-1*fsign)*p2),
+            (2,): fwm * normz,
+            (3,): fwm * (fsign * p1 + complex(0,-1*fsign)*p2),
+
+            }
+
+    
+        self.representation = aloha_lib.LorentzObjectRepresentation(u,
+                                self.lorentz_ind,self.spin_ind)
+
+
+class VFMC(aloha_lib.FactoryLorentz):
+
+    object_class = L_VFMC
+    
+    @classmethod
+    def get_unique_name(self, spin1, particle):
+        return '_VFMC%s_%s_' % (particle, spin1)
 
 #===============================================================================
 # Denominator Propagator 
@@ -979,11 +1656,12 @@
 SpinorPropagatorin = lambda spin1, spin2, particle: (Gamma('mu', spin1, spin2) * \
                     P('mu', particle) + Mass(particle) * Identity(spin1, spin2))
 
+
 VectorPropagator = lambda l1, l2, part: complex(0,1)*(-1 * Metric(l1, l2) + OverMass2(part) * \
                                     Metric(l1,'I3')* P('I3', part) * P(l2, part))
 
 VectorPropagatorMassless= lambda l1, l2, part: complex(0,-1) * Metric(l1, l2)
-
+    
 
 Spin3halfPropagatorin =  lambda mu, nu, s1, s2, part: (\
     -1/3 * (Gamma(mu,s1,-2) + Identity(s1, -2) *  P(mu, part) * Mass(part) * OverMass2(part))* \

=== modified file 'aloha/aloha_parsers.py'
--- aloha/aloha_parsers.py	2019-11-14 08:33:54 +0000
+++ aloha/aloha_parsers.py	2020-01-09 23:16:36 +0000
@@ -206,10 +206,13 @@
 
 class ALOHAExpressionParser(UFOExpressionParser):
 
-    aloha_object = ['P','Gamma','Gamma5','Sigma','Mass','PSlash',
+    aloha_object = ['P', 'PVec','Gamma','Gamma5','Sigma','Mass','PSlash',
                     'OverMass2','Width','Scalar','Spinor','Vector',
                     'Spin2','Spin32','C','Epsilon','Metric','Identity',
-                    'ProjM','ProjP','Coup']
+                    'ProjM','ProjP','Coup','Norm', 'EPSL', 'EPST1', 'EPST2', 'PT',
+                    'UFP', 'UFM', 'UFPC', 'UFMC',
+                    'VFP', 'VFM', 'VFPC', 'VFMC',
+                     'Tnorm', 'TnormZ']
 
     def p_expression_pi(self, p):
         '''expression : PI'''

=== modified file 'aloha/aloha_writers.py'
--- aloha/aloha_writers.py	2018-09-19 09:41:33 +0000
+++ aloha/aloha_writers.py	2020-01-09 23:16:36 +0000
@@ -594,12 +594,12 @@
         p = [] # a list for keeping track how to write the momentum
         
         signs = self.get_momentum_conservation_sign()
-        
         for i,type in enumerate(self.particles):
             if self.declaration.is_used('OM%s' % (i+1)):
                 out.write("    OM{0} = {1}\n    if (M{0}.ne.{1}) OM{0}={2}/M{0}**2\n".format( 
                          i+1, self.change_number_format(0), self.change_number_format(1)))
-            
+
+
             if i+1 == self.outgoing:
                 out_type = type
                 out_size = self.type_to_size[type] 
@@ -612,6 +612,8 @@
                 
         # define the resulting momenta
         if self.offshell:
+            
+            
             energy_pos = out_size -2
             type = self.particles[self.outgoing-1]
             
@@ -622,6 +624,27 @@
             if self.declaration.is_used('P%s' % self.outgoing):
                 self.get_one_momenta_def(self.outgoing, out)
 
+            if "P1T" in self.tag or "P1L" in self.tag:
+                for i in range(1,4):
+                    P = "P%s" % (self.outgoing)
+                    value = ["1d-30", "0d0", "1d-15"]
+                    out.write("  IF (DABS(%(P)s(0))*1e-10.gt.DABS(%(P)s(%(i)s))) %(P)s(%(i)s)=%(val)s\n"
+                              % {"P": P, "i":i, 'val':value[i-1]})
+            i = self.outgoing -1
+            if self.declaration.is_used('Tnorm%s' % (i+1)):
+                out.write("    TNORM{0} = DSQRT(P{0}(1)*P{0}(1)+P{0}(2)*P{0}(2)+P{0}(3)*P{0}(3))\n".format(
+                        i+1))
+            if self.declaration.is_used('TnormZ%s' % (i+1)):
+                out.write("    TNORMZ{0} =  TNORM{0} - P{0}(3)\n".format(
+                        i+1))
+
+            if self.declaration.is_used('FWP%s' % (i+1)):
+                out.write("     FWP{0} = DSQRT(-P{0}(0) + TNORM{0})\n"\
+                          .format(i+1))
+            if self.declaration.is_used('FWM%s' % (i+1)):
+                out.write("     FWM{0} = DSQRT(-P{0}(0) - TNORM{0})\n"\
+                          .format(i+1))
+                #out.write("     FWM{0} = M{0}/FWP{0}\n".format(i+1))
         
         # Returning result
         return out.getvalue()
@@ -715,7 +738,10 @@
             tmp = Fraction(str(number))
             tmp = tmp.limit_denominator(100)
             if not abs(tmp - number) / abs(tmp + number) < 1e-8:
-                out = '%s%s' % (number, self.format)
+                if 'e' in str(number):
+                    out = str(number).replace('e','d')
+                else:
+                    out = '%s%s' % (number, self.format)
             else:
                 out = '%s%s/%s%s' % (tmp.numerator, self.format, tmp.denominator, self.format)
         return out
@@ -745,7 +771,6 @@
         keys.sort(sort_fct)
         for name in keys:
             fct, objs = self.routine.fct[name]
-
             format = ' %s = %s\n' % (name, self.get_fct_format(fct))
             try:
                 text = format % ','.join([self.write_obj(obj) for obj in objs])
@@ -774,7 +799,14 @@
         else:
             OffShellParticle = '%s%d' % (self.particles[self.offshell-1],\
                                                                   self.offshell)
-            if 'L' not in self.tag:
+            is_loop = False
+            if 'L' in self.tag:
+                if self.tag.count('L') == 1 and 'PL' in self.tag:
+                    is_loop = False
+                else:
+                    is_loop = True
+                    
+            if not is_loop:
                 coeff = 'denom*'    
                 if not aloha.complex_mass:
                     if self.routine.denominator:
@@ -917,6 +949,7 @@
             writer.writelines(text)
         return text
 
+
 class QP(object): 
     """routines for writing out Fortran"""
     

=== modified file 'aloha/create_aloha.py'
--- aloha/create_aloha.py	2018-06-13 07:15:48 +0000
+++ aloha/create_aloha.py	2020-01-09 23:16:36 +0000
@@ -238,6 +238,7 @@
                                           for name in aloha_lib.KERNEL.use_tag
                                           if name.startswith('TMP')])
         
+
         output.fct = dict([(name, aloha_lib.KERNEL.reduced_expr2[name])
                                           for name in aloha_lib.KERNEL.use_tag
                                           if name.startswith('FCT')])
@@ -255,7 +256,7 @@
         
         if need_P_sign:
             expr = re.sub(r'\b(P|PSlash)\(', r'-\1(', expr)
-        
+
         calc = aloha_parsers.ALOHAExpressionParser()
         lorentz_expr = calc.parse(expr)
         return lorentz_expr
@@ -431,9 +432,43 @@
     def get_custom_propa(self, propa, spin, id):
         """Return the ALOHA object associated to the user define propagator"""
 
-        propagator = getattr(self.model.propagators, propa)
-        numerator = propagator.numerator
-        denominator = propagator.denominator      
+        if not propa.startswith('1'):
+            propagator = getattr(self.model.propagators, propa)
+            numerator = propagator.numerator
+            denominator = propagator.denominator      
+        elif propa == "1L":
+            numerator = "EPSL(1,id) * EPSL(2,id)"
+            denominator = "-1*PVec(-2,id)*PVec(-2,id)*Mass(id)**2 * (P(-1,id)**2 - Mass(id) * Mass(id) + complex(0,1) * Mass(id) * Width(id))"
+        elif propa == "1T":
+            numerator = "-1*PVec(-2,id)*PVec(-2,id) * EPST2(1,id)*EPST2(2,id) + EPST1(1,id)*EPST1(2,id)"
+            denominator = "PVec(-2,id)*PVec(-2,id) * PT(-3,id)*PT(-3,id) * (P(-1,id)**2 - Mass(id) * Mass(id) + complex(0,1) * Mass(id) * Width(id))"
+        elif propa == "1A":
+            numerator = "(P(-2,id)**2 - Mass(id)**2) * P(1,id) * P(2,id)"
+            denominator = "P(-2,id)**2 * Mass(id)**2 * (P(-1,id)**2 - Mass(id) * Mass(id) + complex(0,1) * Mass(id) * Width(id))"
+        elif propa in ["1P"]:
+            # shift and flip the tag if we multiply by C matrices
+            spin_id = id
+            if (id + 1) // 2 in self.conjg:
+                spin_id += _conjugate_gap + id % 2 - (id +1) % 2
+            if (spin_id % 2):
+                numerator =  "UFP(1,id)*UFPC(2,id)"
+            else:
+                numerator =  "VFP(1,id)*VFPC(2,id)"
+               
+            denominator = "(2*Tnorm(id)*TnormZ(id))*(P(-1,id)*P(-1,id) - Mass(id) * Mass(id) + complex(0,1) * Mass(id) * Width(id))"
+        
+        elif propa == "1M":
+            # shift and flip the tag if we multiply by C matrices
+            spin_id = id
+            if (id + 1) // 2 in self.conjg:
+                spin_id += _conjugate_gap + id % 2 - (id +1) % 2
+            if (spin_id % 2):
+                numerator =  "UFM(1,id)*UFMC(2,id)"
+            else:
+                numerator =  "VFM(1,id)*VFMC(2,id)"
+            denominator = "(2*Tnorm(id)*TnormZ(id))*(P(-1,id)*P(-1,id) - Mass(id) * Mass(id) + complex(0,1) * Mass(id) * Width(id))"
+        else:
+            raise Exception
 
         # Find how to make the replacement for the various tag in the propagator expression
         needPflipping = False
@@ -475,15 +510,16 @@
         
         numerator = self.mod_propagator_expression(tag, numerator)
         if denominator:
-            denominator = self.mod_propagator_expression(tag, denominator)      
+            denominator = self.mod_propagator_expression(tag, denominator)  
+                
         numerator = self.parse_expression(numerator, needPflipping)
-        
+      
         if denominator:
             self.denominator = self.parse_expression(denominator, needPflipping)
             self.denominator = eval(self.denominator)
             if not isinstance(self.denominator, numbers.Number):
                 self.denominator = self.denominator.simplify().expand().simplify().get((0,))
-                
+        needPflipping = False
         if spin ==4:
             return eval(numerator) * propaR
         else:

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2019-09-12 07:01:20 +0000
+++ madgraph/core/base_objects.py	2020-01-09 23:16:36 +0000
@@ -24,6 +24,7 @@
 import re
 import StringIO
 import madgraph.core.color_algebra as color
+import collections
 from madgraph import MadGraph5Error, MG5DIR, InvalidCmd
 import madgraph.various.misc as misc 
 
@@ -1968,6 +1969,8 @@
         self['from_group'] = True
         # onshell: decaying leg (True), forbidden s-channel (False), none (None)
         self['onshell'] = None
+        # filter on the helicty
+        self['polarization'] = []
 
     def filter(self, name, value):
         """Filter for valid leg property values."""
@@ -2000,12 +2003,24 @@
                 raise self.PhysicsObjectError, \
                         "%s is not a valid boolean for leg flag onshell" % \
                                                                     str(value)
+                                                                    
+        
+        if name == 'polarization':
+            if not isinstance(value, list):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid list" % str(value)
+            for i in value:
+                if i not in [-1, 1, 2,-2, 3,-3, 0, 99]:
+                    raise self.PhysicsObjectError, \
+                          "%s is not a valid polarization" % str(value)
+                                                                    
+        
         return True
 
     def get_sorted_keys(self):
         """Return particle property names as a nicely sorted list."""
 
-        return ['id', 'number', 'state', 'from_group', 'loop_line', 'onshell']
+        return ['id', 'number', 'state', 'from_group', 'loop_line', 'onshell', 'polarization']
 
     def is_fermion(self, model):
         """Returns True if the particle corresponding to the leg is a
@@ -2165,6 +2180,7 @@
 
         self['ids'] = []
         self['state'] = True
+        self['polarization'] = []
 
     def filter(self, name, value):
         """Filter for valid multileg property values."""
@@ -2178,6 +2194,15 @@
                     raise self.PhysicsObjectError, \
                           "%s is not a valid list of integers" % str(value)
 
+        if name == 'polarization':
+            if not isinstance(value, list):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid list" % str(value)
+            for i in value:
+                if i not in [-1, 1,  2, -2, 3, -3, 0, 99]:
+                    raise self.PhysicsObjectError, \
+                          "%s is not a valid polarization" % str(value)
+
         if name == 'state':
             if not isinstance(value, bool):
                 raise self.PhysicsObjectError, \
@@ -2189,7 +2214,7 @@
     def get_sorted_keys(self):
         """Return particle property names as a nicely sorted list."""
 
-        return ['ids', 'state']
+        return ['ids', 'state','polarization']
 
 #===============================================================================
 # LegList
@@ -3512,13 +3537,11 @@
         final_legs = filter(lambda leg: leg.get('state') == True, \
                               self.get_legs_with_decays())
 
-        identical_indices = {}
+        identical_indices = collections.defaultdict(int)
         for leg in final_legs:
-            if leg.get('id') in identical_indices:
-                identical_indices[leg.get('id')] = \
-                                    identical_indices[leg.get('id')] + 1
-            else:
-                identical_indices[leg.get('id')] = 1
+            key = (leg.get('id'), tuple(leg.get('polarization')))
+            identical_indices[key] += 1
+
         return reduce(lambda x, y: x * y, [ math.factorial(val) for val in \
                         identical_indices.values() ], 1)
 
@@ -3867,18 +3890,19 @@
         
         # First make sure that the desired particle ids belong to those defined
         # in this process definition.
-        my_isids = [leg.get('ids') for leg in self.get('legs') \
-              if not leg.get('state')]
-        my_fsids = [leg.get('ids') for leg in self.get('legs') \
-             if leg.get('state')]            
-        for i, is_id in enumerate(initial_state_ids):
-            assert is_id in my_isids[i]
-        for i, fs_id in enumerate(final_state_ids):
-            assert fs_id in my_fsids[i]
+        if __debug__:
+            my_isids = [leg.get('ids') for leg in self.get('legs') \
+                  if not leg.get('state')]
+            my_fsids = [leg.get('ids') for leg in self.get('legs') \
+                 if leg.get('state')]            
+            for i, is_id in enumerate(initial_state_ids):
+                assert is_id in my_isids[i]
+            for i, fs_id in enumerate(final_state_ids):
+                assert fs_id in my_fsids[i]
         
         return self.get_process_with_legs(LegList(\
-               [Leg({'id': id, 'state':False}) for id in initial_state_ids] + \
-               [Leg({'id': id, 'state':True}) for id in final_state_ids]))
+               [Leg({'id': id, 'state':False, 'polarization':[]}) for id in initial_state_ids] + \
+               [Leg({'id': id, 'state':True, 'polarization':[]}) for id in final_state_ids]))
 
     def __eq__(self, other):
         """Overloading the equality operator, so that only comparison

=== modified file 'madgraph/core/diagram_generation.py'
--- madgraph/core/diagram_generation.py	2019-03-15 15:43:38 +0000
+++ madgraph/core/diagram_generation.py	2020-01-09 23:16:36 +0000
@@ -627,7 +627,6 @@
         leglist = self.copy_leglist(process.get('legs'))
 
         for leg in leglist:
-            
             # For the first step, ensure the tag from_group 
             # is true for all legs
             leg.set('from_group', True)
@@ -1690,36 +1689,42 @@
         # identifying identical matrix elements already at this stage.
         model = process_definition['model']
         
+        islegs = [leg for leg in process_definition['legs'] \
+                 if leg['state'] == False]
+        fslegs = [leg for leg in process_definition['legs'] \
+                 if leg['state'] == True]        
+        
         isids = [leg['ids'] for leg in process_definition['legs'] \
                  if leg['state'] == False]
-        fsids = [leg['ids'] for leg in process_definition['legs'] \
+        fsids = [leg['ids']  for leg in process_definition['legs'] \
+                 if leg['state'] == True]
+        polids = [tuple(leg['polarization'])  for leg in process_definition['legs'] \
                  if leg['state'] == True]
         # Generate all combinations for the initial state
-        
         for prod in itertools.product(*isids):
             islegs = [\
-                    base_objects.Leg({'id':id, 'state': False}) \
-                    for id in prod]
+                    base_objects.Leg({'id':id, 'state': False, 
+                                      'polarization': islegs[i]['polarization']})
+                    for i,id in enumerate(prod)]
 
             # Generate all combinations for the final state, and make
             # sure to remove double counting
 
-            red_fsidlist = []
+            red_fsidlist = set()
 
             for prod in itertools.product(*fsids):
-
+                tag = zip(prod, polids)
+                tag = sorted(tag)
                 # Remove double counting between final states
-                if tuple(sorted(prod)) in red_fsidlist:
+                if tuple(tag) in red_fsidlist:
                     continue
                 
-                red_fsidlist.append(tuple(sorted(prod)));
-                
+                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}) \
-                        for id in prod])
+                        base_objects.Leg({'id':id, 'state': True, 'polarization': fslegs[i]['polarization']}) \
+                        for i,id  in enumerate(prod)])
                 
                 legs = base_objects.LegList(leg_list)
 

=== modified file 'madgraph/core/helas_objects.py'
--- madgraph/core/helas_objects.py	2019-08-02 21:35:19 +0000
+++ madgraph/core/helas_objects.py	2020-01-09 23:16:36 +0000
@@ -638,6 +638,9 @@
         # conjugate_indices is a list [1,2,...] with fermion lines
         # that need conjugates. Default is "None"
         self['conjugate_indices'] = None
+        #
+        #
+        self['polarization'] = []
 
     # Customized constructor
     def __init__(self, *arguments):
@@ -652,6 +655,7 @@
                 leg = arguments[0]
                 interaction_id = arguments[1]
                 model = arguments[2]
+                
                 # decay_ids is the pdg codes for particles with decay
                 # chains defined
                 decay_ids = []
@@ -672,11 +676,17 @@
                 # the decay is applied
                 if self['state'] == 'final' and self.get('pdg_code') in decay_ids:
                     self.set('decay', True)
-
+                else:
+                    if 99 in leg.get('polarization'):
+                        raise Exception("polarization A only valid for propagator.")
                 # Set fermion flow state. Initial particle and final
                 # antiparticle are incoming, and vice versa for
                 # outgoing
                 if self.is_fermion():
+                    if leg.get('polarization'):
+                        pol = list(leg.get('polarization'))
+                        self.set('polarization', pol) 
+                    
                     if leg.get('state') == False and \
                            self.get('is_part') or \
                            leg.get('state') == True and \
@@ -684,12 +694,16 @@
                         self.set('state', 'incoming')
                     else:
                         self.set('state', 'outgoing')
+                else:
+                    self.set('polarization', leg.get('polarization'))
                 self.set('interaction_id', interaction_id, model)
         elif arguments:
             super(HelasWavefunction, self).__init__(arguments[0])
         else:
             super(HelasWavefunction, self).__init__()
 
+
+
     def filter(self, name, value):
         """Filter for valid wavefunction property values."""
 
@@ -838,6 +852,15 @@
                 raise self.PhysicsObjectError, \
                         "%s is not a valid int" % str(value) + \
                         " for the lcut_size"
+                        
+        if name == 'polarization':
+            if not isinstance(value, list):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid list" % str(value)
+            for i in value:
+                if i not in [-1, 1, 2, -2, 3, -3, 0, 99]:
+                    raise self.PhysicsObjectError, \
+                      "%s is not a valid polarization" % str(value)
 
         return True
 
@@ -1567,6 +1590,33 @@
         output['propa'] = self.get('particle').get('propagator')
         if output['propa'] not in ['', None]:
             output['propa'] = 'P%s' % output['propa']
+            if self.get('polarization'):
+                raise InvalidCmd, 'particle with custom propagator can not have polarization'
+        elif self.get('polarization'):
+            if self.get('polarization') == [0]:
+                if self.get('spin') != 3:
+                    raise InvalidCmd, 'polarization not handle for decay particle'
+                output['propa'] = 'P1L' 
+            elif self.get('polarization') == [1,-1]:
+                if self.get('spin') != 3:
+                    raise InvalidCmd, 'polarization not handle for decay particle'
+                output['propa'] = 'P1T'
+            elif self.get('polarization') == [99]:
+                if self.get('spin') != 3:
+                    raise InvalidCmd, 'polarization not handle for decay particle'
+                output['propa'] = 'P1A'
+            elif self.get('polarization') == [1]:
+                if self.get('spin') != 2:
+                    raise InvalidCmd, 'polarization not handle for decay particle'
+                output['propa'] = 'P1P'
+            elif self.get('polarization') == [-1]:
+                misc.sprint(self.get('polarization'), self.get('spin'))
+                if self.get('spin') != 2:
+                    raise InvalidCmd, 'Left polarization not handle for decay particle for spin (2s+1=%s) particles' % self.get('spin') 
+                output['propa'] = 'P1M'
+            else:            
+                raise InvalidCmd, 'polarization not handle for decay particle'
+            
         # optimization
         if aloha.complex_mass: 
             if (self.get('width') == 'ZERO' or self.get('mass') == 'ZERO'):
@@ -1641,16 +1691,22 @@
         # Sort according to spin and flow direction
         res.sort()
         res.append(self.get_spin_state_number())
-        res.append(self.find_outgoing_number())
+        outgoing =self.find_outgoing_number()
+        res.append(outgoing)
 
         if self['is_loop']:
             res.append(self.get_loop_index())
             if not self.get('mothers'):
                 res.append(self.get('is_part'))
 
+        res.append(tuple(self.get('polarization')) )
+
         # Check if we need to append a charge conjugation flag
         if self.needs_hermitian_conjugate():
             res.append(self.get('conjugate_indices'))
+            
+
+        
 
         return (tuple(res), tuple(self.get('lorentz')))
 
@@ -1774,6 +1830,19 @@
 
         if self.get('particle').get('propagator') not in ['', None]:
             tags.append('P%s' % str(self.get('particle').get('propagator')))
+        elif self.get('polarization'):
+            if self.get('polarization') == [0]:
+                tags.append('P1L') 
+            elif self.get('polarization') == [1,-1]:
+                tags.append('P1T')
+            elif self.get('polarization') == [99]:
+                tags.append('P1A')
+            elif self.get('polarization') == [1]:
+                tags.append('P1P')
+            elif self.get('polarization') == [-1]:
+                tags.append('P1M')
+            else:
+                raise InvalidCmd, 'polarization not handle for decay particle'
 
         return (tuple(self.get('lorentz')),tuple(tags),self.find_outgoing_number())
 
@@ -3811,7 +3880,7 @@
                           filter(lambda wf: not wf.get('mothers') and \
                                  wf.get('number_external') == number,
                                  self.get_all_wavefunctions())]
-
+        
         # Keep track of wavefunction and amplitude numbers, to ensure
         # unique numbers for all new wfs and amps during manipulations
         numbers = [self.get_all_wavefunctions()[-1].get('number'),
@@ -3971,6 +4040,12 @@
            amplitudes which have this wavefunction as mother.
         """
 
+        #check that decay does not specify polarization
+        wfs = filter(lambda w: w.get('state') == 'initial' , decay.get('diagrams')[0].get('wavefunctions'))
+        if any(wf.get('polarization') for wf in wfs):
+            raise InvalidCmd, 'In decay-chain polarization can only be specified in production not in decay. Please Retry'
+
+        
         len_decay = len(decay.get('diagrams'))
 
         number_external = old_wfs[0].get('number_external')
@@ -4218,6 +4293,11 @@
 
                     old_wf_index = [wf.get('number') for wf in \
                                     diagram_wfs].index(old_wf.get('number'))
+                                    
+                    old_wf_pol = diagram_wfs[old_wf_index].get('polarization')
+                    for w in final_decay_wfs:
+                        w.set('polarization', old_wf_pol)
+                    
 
                     diagram_wfs = diagram_wfs[0:old_wf_index] + \
                                   decay_diag_wfs + diagram_wfs[old_wf_index:]
@@ -4425,6 +4505,7 @@
         """Insert decay chain by simply modifying wavefunction. This
         is possible only if there is only one diagram in the decay."""
 
+
         for key in old_wf.keys():
             old_wf.set(key, new_wf[key])
 
@@ -4546,6 +4627,8 @@
             mothers.append(wf)
 
         return mothers
+    
+
 
     def get_num_configs(self):
         """Get number of diagrams, which is always more than number of
@@ -4648,11 +4731,13 @@
             return None
 
         model = self.get('processes')[0].get('model')
+        hel_per_part = [ len(wf.get('polarization')) if wf.get('polarization') 
+                        else len(model.get('particle_dict')[\
+                                  wf.get('pdg_code')].get_helicity_states())
+            for wf in self.get_external_wavefunctions()]
 
         return reduce(lambda x, y: x * y,
-                      [ len(model.get('particle_dict')[wf.get('pdg_code')].\
-                            get_helicity_states())\
-                        for wf in self.get_external_wavefunctions() ], 1)
+                      hel_per_part)
 
     def get_helicity_matrix(self, allow_reverse=True):
         """Gives the helicity matrix for external wavefunctions"""
@@ -4663,9 +4748,12 @@
         process = self.get('processes')[0]
         model = process.get('model')
 
-        return apply(itertools.product, [ model.get('particle_dict')[\
-                                  wf.get('pdg_code')].get_helicity_states(allow_reverse)\
-                                  for wf in self.get_external_wavefunctions()])
+        hel_per_part = [ wf.get('polarization') if wf.get('polarization') 
+                        else model.get('particle_dict')[\
+                                  wf.get('pdg_code')].get_helicity_states(allow_reverse)
+            for wf in self.get_external_wavefunctions()]
+
+        return apply(itertools.product, hel_per_part)
 
     def get_hel_avg_factor(self):
         """ Calculate the denominator factor due to the average over initial
@@ -4674,11 +4762,28 @@
         model = self.get('processes')[0].get('model')
         initial_legs = filter(lambda leg: leg.get('state') == False, \
                               self.get('processes')[0].get('legs'))
-        
-        return reduce(lambda x, y: x * y,
-                      [ len(model.get('particle_dict')[leg.get('id')].\
-                                   get_helicity_states())\
-                        for leg in initial_legs ])
+        hel_per_part = [ len(leg.get('polarization')) if leg.get('polarization') 
+                        else len(model.get('particle_dict')[\
+                                  leg.get('id')].get_helicity_states())
+            for leg in initial_legs]
+        
+        return reduce(lambda x, y: x * y, hel_per_part, 1)
+
+    def get_spin_state_initial(self):
+        """Gives (number of state for each initial particle)"""
+
+        model = self.get('processes')[0].get('model')
+        initial_legs = filter(lambda leg: leg.get('state') == False, \
+                              self.get('processes')[0].get('legs'))
+        hel_per_part = [ len(leg.get('polarization')) if leg.get('polarization') 
+                        else len(model.get('particle_dict')[\
+                                  leg.get('id')].get_helicity_states())
+            for leg in initial_legs]
+        
+        if len(hel_per_part) == 1:
+            hel_per_part.append(0)
+        
+        return hel_per_part
 
     def get_beams_hel_avg_factor(self):
         """ Calculate the denominator factor due to the average over initial
@@ -4714,7 +4819,9 @@
             
         spin_factor = reduce(lambda x, y: x * y,
                              [ len(model.get('particle_dict')[leg.get('id')].\
-                                   get_helicity_states())\
+                                   get_helicity_states()) 
+                              if not leg.get('polarization') else 
+                              len(leg.get('polarization'))
                                for leg in initial_legs ])
 
         return spin_factor * color_factor * self['identical_particle_factor']

=== modified file 'madgraph/fks/fks_common.py'
--- madgraph/fks/fks_common.py	2019-02-26 09:15:55 +0000
+++ madgraph/fks/fks_common.py	2020-01-09 23:16:36 +0000
@@ -667,6 +667,7 @@
          'number': fksleg.get('number'),
          'state': fksleg.get('state'),
          'from_group': fksleg.get('from_group'),
+         'polarization': fksleg.get('polarization')
           })
     return leg
 

=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py	2020-01-09 13:09:41 +0000
+++ madgraph/interface/madgraph_interface.py	2020-01-09 23:16:36 +0000
@@ -1062,6 +1062,49 @@
             if re.search('\D\$', particles):
                 raise self.InvalidCmd(
                 'wrong process format: restriction should be place after the final states')
+                
+        # '{}' should only be used for onshell particle (including initial/final state)
+        # check first that polarization are not include between > >
+        if nbsep == 2:
+            if '{' in particles_parts[1]:
+                raise self.InvalidCmd('Polarization restriction can not be used as required s-channel')
+        split = re.split('\D[$|/]',particles_parts[-1],1)
+        if len(split)==2:
+            if '{' in split[1]:
+                raise self.InvalidCmd('Polarization restriction can not be used in forbidding particles')
+            
+        if '[' in process and '{' in process:
+            valid = False
+            if 'noborn' in process:
+                valid = True
+            else:
+                raise self.InvalidCmd('Polarization restriction can not be used for NLO processes')
+
+            # below are the check when [QCD] will be valid for computation            
+            order = process.split('[')[1].split(']')[0]
+            if '=' in order:
+                order = order.split('=')[1]
+            if order.strip().lower() != 'qcd':
+                raise self.InvalidCmd('Polarization restriction can not be used for generic NLO computations')
+
+
+            for p in particles_parts[1].split():
+                if '{' in p:
+                    part = p.split('{')[0]
+                else:
+                    continue
+                if self._curr_model:
+                    p = self._curr_model.get_particle(part)
+                    if not p:
+                        if part in self._multiparticles:
+                            for part2 in self._multiparticles[part]:
+                                p = self._curr_model.get_particle(part2)
+                                if p.get('color') != 1:
+                                    raise self.InvalidCmd('Polarization restriction can not be used for color charged particles')
+                        continue
+                    if p.get('color') != 1:
+                        raise self.InvalidCmd('Polarization restriction can not be used for color charged particles')
+        
 
 
     def check_tutorial(self, args):
@@ -4645,6 +4688,7 @@
                 else:
                     orders[order]=99
         
+
         if not self._curr_model['case_sensitive']:
             # Particle names lowercase
             line = line.lower()
@@ -4700,6 +4744,82 @@
                 continue
 
             mylegids = []
+            polarization = []
+            if '{' in part_name:
+                part_name, pol = part_name.split('{',1)
+                pol, rest = pol.split('}',1)
+                
+                no_dup_name = part_name
+                while True:
+                    try:
+                        spin = self._curr_model.get_particle(no_dup_name).get('spin')
+                        break
+                    except AttributeError:
+                        if no_dup_name in self._multiparticles:
+                            spins = set([self._curr_model.get_particle(p).get('spin') for p in self._multiparticles[no_dup_name]])
+                            if len(spins) > 1:
+                                raise self.InvalidCmd('Can not use polarised on multi-particles for multi-particles with various spin')
+                            else:
+                                spin = spins.pop()
+                                break
+                        elif no_dup_name[0].isdigit():
+                            no_dup_name = no_dup_name[1:]
+                        else:
+                            raise
+                if rest:
+                    raise self.InvalidCmd('A space is required after the "}" symbol to separate particles')
+                ignore  =False
+                for i,p in enumerate(pol):
+                    if ignore or p==',':
+                        ignore= False
+                        continue
+                    if p in ['t','T']:
+                        if spin == 3:
+                            polarization += [1,-1]
+                        else:
+                            raise self.InvalidCmd('"T" (transverse) polarization are only supported for spin one particle.')
+                    elif p in ['l', 'L']:
+                        if spin == 3:
+                            logger.warning('"L" polarization is interpreted as Left for Longitudinal please use "0".')
+                        polarization += [-1]
+                    elif p in ['R','r']:
+                        polarization += [1]
+                    elif p in ["A",'a']:
+                        if spin == 3:
+                            polarization += [99]
+                        else:
+                            raise self.InvalidCmd('"A" (auxiliary) polarization are only supported for spin one particle.')
+                    elif p in ['+']:
+                        if i +1 < len(pol) and pol[i+1].isdigit():
+                            p = int(pol[i+1])
+                            if abs(p) > 3: 
+                                raise self.InvalidCmd("polarization are between -3 and 3")
+                            polarization.append(p)
+                            ignore = True
+                        else:
+                            polarization += [1]
+                    elif p in ['-']:
+                        if i+1 < len(pol) and pol[i+1].isdigit():
+                            p = int(pol[i+1])
+                            if abs(p) > 3: 
+                                raise self.InvalidCmd("polarization are between -3 and 3")
+                            polarization.append(-p)
+                            ignore = True
+                        else:
+                            polarization += [-1]
+                    elif p in [0,'0']:
+                        if spin in [1,2]:
+                            raise self.InvalidCmd('"0" (longitudinal) polarization are not supported for scalar/fermion.')
+                        else:
+                            polarization += [0]
+                    elif p.isdigit():
+                        p = int(p)
+                        if abs(p) > 3: 
+                            raise self.InvalidCmd("polarization are between -3 and 3")
+                        polarization.append(p)
+                    else:
+                        raise self.InvalidCmd('Invalid Polarization')
+
             duplicate =1
             if part_name in self._multiparticles:
                 if isinstance(self._multiparticles[part_name][0], list):
@@ -4735,7 +4855,8 @@
             if mylegids:
                 for _ in range(duplicate):
                     myleglist.append(base_objects.MultiLeg({'ids':mylegids,
-                                                        'state':state}))
+                                                        'state':state,
+                                                        'polarization': polarization}))
             else:
                 raise self.InvalidCmd, "No particle %s in model" % part_name
 
@@ -5063,6 +5184,8 @@
         final_states = re.search(r'> ([^\/\$\=\@>]*)(\[|\s\S+\=|\$|\/|\@|$)', procline)
         particles = final_states.groups()[0]
         for particle in particles.split():
+            if '{' in particle:
+                particle = particle.split('{')[0]
             if particle in pids:
                 final.add(pids[particle])
             elif particle in self._multiparticles:

=== modified file 'madgraph/interface/reweight_interface.py'
--- madgraph/interface/reweight_interface.py	2019-12-11 10:31:48 +0000
+++ madgraph/interface/reweight_interface.py	2020-01-09 23:16:36 +0000
@@ -99,6 +99,7 @@
         self.dedicated_path = {}
         self.soft_threshold = None
         self.systematics = False # allow to run systematics in ouput2.0 mode
+        self.boost_event = False
         self.mg5cmd = master_interface.MasterCmd()
         if mother:
             self.mg5cmd.options.update(mother.options)
@@ -368,6 +369,8 @@
                 self.second_process.append(" ".join(args[1:-1]))
             else:
                 self.second_process = [" ".join(args[1:])]
+        elif args[0] == "boost":
+            self.boost_event = eval(' '.join(args[1:]))
         elif args[0] in ['virtual_path', 'tree_path']:
             self.dedicated_path[args[0]] = os.path.abspath(args[1])
         elif args[0] == "output":
@@ -1196,7 +1199,27 @@
         # For 2>N pass in the center of mass frame
         #   - required for helicity by helicity re-weighitng
         #   - Speed-up loop computation 
-        if (hasattr(event[1], 'status') and event[1].status == -1) or \
+        if (hypp_id == 0 and ('fram_id' in self.banner.run_card and self.banner.run_card['frame_id'] !=6)):
+            import copy
+            new_event = copy.deepcopy(event)
+            pboost = FourMomenta()
+            to_inc = bin(self.banner.run_card['frame_id'])[2:]
+            to_inc.reverse()
+            nb_ext = 0
+            for p in new_event:
+                if p.status in [-1,1]:
+                    nb_ext += 1
+                    if to_inc[nb_ext]:
+                        pboost += p                    
+            new_event.boost(pboost)
+            p = new_event.get_momenta(orig_order)
+        elif (hypp_id == 1 and self.boost_event):
+            if self.boost_event is not True:
+                import copy
+                new_event = copy.deepcopy(event)
+                new_event.boost(self.boost_event)
+                p = new_event.get_momenta(orig_order)        
+        elif (hasattr(event[1], 'status') and event[1].status == -1) or \
            (event[1].px == event[1].py == 0.):
             pboost = lhe_parser.FourMomentum(p[0]) + lhe_parser.FourMomentum(p[1])
             for i,thisp in enumerate(p):

=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2019-06-26 14:12:13 +0000
+++ madgraph/iolibs/export_v4.py	2020-01-09 23:16:36 +0000
@@ -916,7 +916,7 @@
 
         #copy Helas Template
         cp(MG5DIR + '/aloha/template_files/Makefile_F', write_dir+'/makefile')
-        if any([any(['L' in tag for tag in d[1]]) for d in wanted_lorentz]):
+        if any([any([tag.startswith('L') for tag in d[1]]) for d in wanted_lorentz]):
             cp(MG5DIR + '/aloha/template_files/aloha_functions_loop.f', 
                                                  write_dir+'/aloha_functions.f')
             aloha_model.loop_mode = False
@@ -4136,6 +4136,11 @@
                           'iolibs/template_files/%s' % self.matrix_file)
         replace_dict['template_file2'] = pjoin(_file_path, \
                           'iolibs/template_files/split_orders_helping_functions.inc')      
+        
+        s1,s2 = matrix_element.get_spin_state_initial()
+        replace_dict['nb_spin_state1'] = s1
+        replace_dict['nb_spin_state2'] = s2
+        
         if writer:
             file = open(replace_dict['template_file']).read()
             file = file % replace_dict
@@ -5224,6 +5229,10 @@
 
         ncomb=matrix_elements[0].get_helicity_combinations()
         replace_dict['read_write_good_hel'] = self.read_write_good_hel(ncomb)
+
+        s1,s2 = matrix_elements[0].get_spin_state_initial()
+        replace_dict['nb_spin_state1'] = s1
+        replace_dict['nb_spin_state2'] = s2
         
         if writer:
             file = open(pjoin(_file_path, \

=== modified file 'madgraph/iolibs/helas_call_writers.py'
--- madgraph/iolibs/helas_call_writers.py	2018-11-08 22:24:49 +0000
+++ madgraph/iolibs/helas_call_writers.py	2020-01-09 23:16:36 +0000
@@ -363,7 +363,7 @@
 
         # Gluon 4-vertex division tensor calls ggT for the FR sm and mssm
 
-        key = ((3, 3, 5, 3), ('A',))
+        key = ((3, 3, 5, 3,tuple()), ('A',))
         call = lambda wf: \
                "CALL UVVAXX(W(1,%d),W(1,%d),%s,zero,zero,zero,W(1,%d))" % \
                (FortranHelasCallWriter.sorted_mothers(wf)[0].get('me_id'),
@@ -372,7 +372,7 @@
                 wf.get('me_id'))
         self.add_wavefunction(key, call)
 
-        key = ((3, 5, 3, 1), ('A',))
+        key = ((3, 5, 3, 1,tuple()), ('A',))
         call = lambda wf: \
                "CALL JVTAXX(W(1,%d),W(1,%d),%s,zero,zero,W(1,%d))" % \
                (FortranHelasCallWriter.sorted_mothers(wf)[0].get('me_id'),
@@ -393,7 +393,7 @@
 
         # SM gluon 4-vertex components
 
-        key = ((3, 3, 3, 3, 1), ('gggg3',))
+        key = ((3, 3, 3, 3, 1,tuple()), ('gggg3',))
         call = lambda wf: \
                "CALL JGGGXX(W(1,%d),W(1,%d),W(1,%d),%s,W(1,%d))" % \
                (FortranHelasCallWriter.sorted_mothers(wf)[1].get('me_id'),
@@ -412,7 +412,7 @@
                 amp.get('coupling')[0],
                 amp.get('number'))
         self.add_amplitude(key, call)
-        key = ((3, 3, 3, 3, 1), ('gggg2',))
+        key = ((3, 3, 3, 3, 1 ,tuple()), ('gggg2',))
         call = lambda wf: \
                "CALL JGGGXX(W(1,%d),W(1,%d),W(1,%d),%s,W(1,%d))" % \
                (FortranHelasCallWriter.sorted_mothers(wf)[0].get('me_id'),
@@ -431,7 +431,7 @@
                 amp.get('coupling')[0],
                 amp.get('number'))
         self.add_amplitude(key, call)
-        key = ((3, 3, 3, 3, 1), ('gggg1',))
+        key = ((3, 3, 3, 3, 1,tuple()), ('gggg1',))
         call = lambda wf: \
                "CALL JGGGXX(W(1,%d),W(1,%d),W(1,%d),%s,W(1,%d))" % \
                (FortranHelasCallWriter.sorted_mothers(wf)[2].get('me_id'),
@@ -453,7 +453,7 @@
 
         # HEFT VVVS calls
 
-        key = ((1, 3, 3, 3, 3), ('',))
+        key = ((1, 3, 3, 3, 3,tuple()), ('',))
         call = lambda wf: \
                "CALL JVVSXX(W(1,%d),W(1,%d),W(1,%d),DUM1,%s,%s,%s,W(1,%d))" % \
                (wf.get('mothers')[0].get('me_id'),
@@ -465,7 +465,7 @@
                 wf.get('me_id'))
         self.add_wavefunction(key, call)
 
-        key = ((3, 3, 3, 1, 4), ('',))
+        key = ((3, 3, 3, 1, 4,tuple()), ('',))
         call = lambda wf: \
                "CALL HVVVXX(W(1,%d),W(1,%d),W(1,%d),DUM1,%s,%s,%s,W(1,%d))" % \
                (wf.get('mothers')[0].get('me_id'),
@@ -490,7 +490,7 @@
 
         # HEFT VVVS calls
 
-        key = ((1, 3, 3, 3, 1), ('',))
+        key = ((1, 3, 3, 3, 1,tuple()), ('',))
         call = lambda wf: \
                "CALL JVVSXX(W(1,%d),W(1,%d),W(1,%d),DUM1,%s,%s,%s,W(1,%d))" % \
                (wf.get('mothers')[0].get('me_id'),
@@ -502,7 +502,7 @@
                 wf.get('me_id'))
         self.add_wavefunction(key, call)
 
-        key = ((3, 3, 3, 1, 4), ('',))
+        key = ((3, 3, 3, 1, 4,tuple()), ('',))
         call = lambda wf: \
                "CALL HVVVXX(W(1,%d),W(1,%d),W(1,%d),DUM1,%s,%s,%s,W(1,%d))" % \
                (wf.get('mothers')[0].get('me_id'),
@@ -537,7 +537,7 @@
                 amp.get('number'))
         self.add_amplitude(key, call)
         
-        key = ((-2, 2, 5, 3), ('',))
+        key = ((-2, 2, 5, 3,tuple()), ('',))
         call = lambda wf: \
                "CALL UIOXXX(W(1,%d),W(1,%d),%s,%s,%s,%s,W(1,%d))" % \
                (wf.get('mothers')[0].get('me_id'),
@@ -1145,6 +1145,7 @@
             outgoing = 0
 
 
+
         # Check if we need to append a charge conjugation flag
         l = [str(l) for l in argument.get('lorentz')]
         flag = []
@@ -1223,7 +1224,7 @@
         # Now we have a line correctly formatted
         call_function = lambda wf: call % wf.get_helas_call_dict(\
           OptimizedOutput=False, specifyHel=self.hel_sum)
-
+        
         # Add the constructed function to wavefunction or amplitude dictionary
         if isinstance(argument, helas_objects.HelasWavefunction):
             self.add_wavefunction(argument.get_call_key(), call_function)
@@ -1510,7 +1511,6 @@
 
         # Creating line formatting:
         call = 'CALL %(routine_name)s(%(wf)s%(coup)s%(mass)s%(out)s)'
-
         arg = {'routine_name': aloha_writers.combine_name(\
                                         '%s' % l[0], l[1:], outgoing, flag, True),
                'coup': ("%%(coup%d)s," * len(argument.get('coupling'))) % \

=== modified file 'madgraph/iolibs/template_files/auto_dsig_v4.inc'
--- madgraph/iolibs/template_files/auto_dsig_v4.inc	2016-08-20 06:27:26 +0000
+++ madgraph/iolibs/template_files/auto_dsig_v4.inc	2020-01-09 23:16:36 +0000
@@ -65,6 +65,11 @@
 %(define_subdiag_lines)s
       include 'coupl.inc'
       include 'run.inc'
+c
+c     local
+c
+	   double precision P1(0:3, nexternal)
+
 C  
 C DATA
 C  
@@ -100,7 +105,12 @@
             DSIG%(proc_id)s = PD(0)
             RETURN
          ENDIF
-         CALL SMATRIX%(proc_id)s(PP,DSIGUU)
+	 if(frame_id.ne.6)then
+           call boost_to_frame(PP, frame_id, P1)
+         else
+           P1 = PP
+         endif
+         CALL SMATRIX%(proc_id)s(P1,DSIGUU)
 		 IF (IMODE.EQ.5) THEN
             IF (DSIGUU.LT.1D199) THEN		 
 		      DSIG%(proc_id)s = DSIGUU*CONV

=== modified file 'madgraph/iolibs/template_files/born_fks.inc'
--- madgraph/iolibs/template_files/born_fks.inc	2017-06-02 12:16:50 +0000
+++ madgraph/iolibs/template_files/born_fks.inc	2020-01-09 23:16:36 +0000
@@ -25,6 +25,7 @@
 C  
       REAL*8 P1(0:3,NEXTERNAL-1)
       COMPLEX*16 ANS(2)
+
 C  
 C LOCAL VARIABLES 
 C  
@@ -34,6 +35,8 @@
       INTEGER NTRY(%(nconfs)d)
       DATA NTRY /%(nconfs)d*0/
       INTEGER NHEL(NEXTERNAL-1,NCOMB)
+      real*8 T_SAVE(NCOMB,2)
+      real*8 T_SUM(2)
 %(helicity_lines)s
 %(den_factor_lines)s
 %(ij_lines)s
@@ -105,6 +108,8 @@
       ANS(1) = 0D0
       ANS(2) = 0D0
       hel_fac=1d0
+      T_SUM(1) = 0d0
+      T_SUM(2) = 0d0
       DO IHEL=1,NCOMB
         ! the following lines are to avoid segfaults when glu_ij=0
         cond_ij=skip(nfksprocess).eq.0
@@ -112,17 +117,32 @@
         !if (nhel(glu_ij,ihel).EQ.NHEL(GLU_IJ,1).or.skip(nfksprocess).eq.0) then
         if (cond_ij) then
           IF ((GOODHEL(IHEL,nFKSprocess) .OR. GOODHEL(IHEL+SKIP(nFKSprocess),nFKSprocess) .OR. NTRY(nFKSprocess) .LT. 2) ) THEN
-            ANS(1)=ANS(1)+BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE,borns)
-            ANS(2)=ANS(2)+BORNTILDE
-            if ( borns(1).ne.0d0 .AND. .NOT. GOODHEL(IHEL,nFKSprocess) ) then
+	    ANS(1) = ANS(1) + BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE,borns) 
+	   ANS(2)=ANS(2)+ BORNTILDE 
+            T_SAVE(IHEL,1) = borns(1)
+	    T_SAVE(IHEL,2) = borns(2)
+	    T_SUM(1) = T_SUM(1) + DABS(borns(1))
+	    T_SUM(2) = T_SUM(2) + DABS(borns(2))
+	  ENDIF
+	  else	
+            T_SAVE(IHEL,1) = 0d0
+            T_SAVE(IHEL,2) = 0d0
+        ENDIF
+      ENDDO
+
+      IF(ntry(nFKSprocess).LT.2) then
+	    if (T_SUM(2).lt.1e-10*T_SUM(1)) T_SUM(2) = T_SUM(1)	
+	do IHEL=1,NCOMB
+            if ( DABS(T_SAVE(IHEL,1)).gt.DABS(T_SUM(1))*1d-8/NCOMB)then
               GOODHEL(IHEL,nFKSprocess)=.TRUE.
             endif
-            if ( borns(2).ne.0d0 .AND. .NOT. GOODHEL(IHEL+SKIP(nFKSprocess),nFKSprocess) ) then
+
+            if ( DABS(T_SAVE(IHEL,2)).gt.DABS(T_SUM(2))*1d-8/NCOMB)then
               GOODHEL(IHEL+SKIP(nFKSprocess),nFKSprocess)=.TRUE.
             endif
-	  ENDIF
-        ENDIF
-      ENDDO
+	 enddo	
+       endif
+
       ANS(1)=ANS(1)/DBLE(IDEN)
       ANS(2)=ANS(2)/DBLE(IDEN)
       wgt_me_born=dble(ans(1))

=== modified file 'madgraph/iolibs/template_files/matrix_madevent_group_v4.inc'
--- madgraph/iolibs/template_files/matrix_madevent_group_v4.inc	2018-10-26 08:20:43 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_group_v4.inc	2020-01-09 23:16:36 +0000
@@ -37,6 +37,7 @@
         LOGICAL GOODHEL(NCOMB,2)
 	INTEGER NTRY(2)
         common/BLOCK_GOODHEL/NTRY,GOODHEL
+
 C  
 C LOCAL VARIABLES 
 C  
@@ -52,6 +53,7 @@
         INTEGER THIS_NTRY(2)
         SAVE THIS_NTRY
         DATA THIS_NTRY /0,0/
+c
 C       This is just to temporarily store the reference grid for helicity of the DiscreteSampler so as to obtain its number of entries with ref_helicity_grid%%n_tot_entries
 	    type(SampledDimension) ref_helicity_grid
 C  
@@ -62,6 +64,9 @@
    
     CHARACTER*101         HEL_BUFF
     COMMON/TO_HELICITY/  HEL_BUFF
+
+    INTEGER NB_SPIN_STATE_in(2)
+    common /nb_hel_state/ nb_spin_state_in
     
     INTEGER IMIRROR
     COMMON/TO_MIRROR/ IMIRROR
@@ -120,9 +125,9 @@
             T=MATRIX%(proc_id)s(P ,NHEL(1,I),JC(1))
             DO JJ=1,nincoming
                IF(POL(JJ).NE.1d0.AND.NHEL(JJ,I).EQ.INT(SIGN(1d0,POL(JJ)))) THEN
-                  T=T*ABS(POL(JJ))
+                  T=T*ABS(POL(JJ))*NB_SPIN_STATE_IN(JJ)/2d0 ! NB_SPIN_STATE(JJ)/2d0 is added for polarised beam
                ELSE IF(POL(JJ).NE.1d0)THEN
-                  T=T*(2d0-ABS(POL(JJ)))
+                  T=T*(2d0-ABS(POL(JJ)))*NB_SPIN_STATE_IN(JJ)/2d0
                ENDIF
             ENDDO
 			IF (ISUM_HEL.NE.0.and.DS_get_dim_status('Helicity').eq.0.and.ALLOW_HELICITY_GRID_ENTRIES) then

=== modified file 'madgraph/iolibs/template_files/matrix_madevent_v4.inc'
--- madgraph/iolibs/template_files/matrix_madevent_v4.inc	2019-06-10 21:44:14 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_v4.inc	2020-01-09 23:16:36 +0000
@@ -37,6 +37,9 @@
     LOGICAL GOODHEL(NCOMB)
     INTEGER NTRY
     common/BLOCK_GOODHEL/NTRY,GOODHEL
+    integer nb_spin_state(2)
+    data  nb_spin_state /%(nb_spin_state1)i,%(nb_spin_state2)i/
+    common /nb_hel_state/ nb_spin_state
 C  
 C LOCAL VARIABLES 
 C  
@@ -102,9 +105,9 @@
                T=MATRIX%(proc_id)s(P ,NHEL(1,I),JC(1))            
              DO JJ=1,nincoming
                IF(POL(JJ).NE.1d0.AND.NHEL(JJ,I).EQ.INT(SIGN(1d0,POL(JJ)))) THEN
-                 T=T*ABS(POL(JJ))
+                 T=T*ABS(POL(JJ))*NB_SPIN_STATE(JJ)/2d0 ! NB_SPIN_STATE(JJ)/2d0 is added for polarised beam
                ELSE IF(POL(JJ).NE.1d0)THEN
-                 T=T*(2d0-ABS(POL(JJ)))
+                 T=T*(2d0-ABS(POL(JJ)))*NB_SPIN_STATE(JJ)/2d0
                ENDIF
              ENDDO
 			 IF (ISUM_HEL.NE.0) then

=== modified file 'madgraph/iolibs/template_files/realmatrix_fks.inc'
--- madgraph/iolibs/template_files/realmatrix_fks.inc	2016-09-08 23:15:34 +0000
+++ madgraph/iolibs/template_files/realmatrix_fks.inc	2020-01-09 23:16:36 +0000
@@ -41,8 +41,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C for the first ps-point, check for helicities that give
 C identical matrix elements
             T=MATRIX_%(N_me)d(P ,NHEL(1,IHEL))            
@@ -68,11 +68,17 @@
           ENDIF
 C add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+        IF ( NTRY .LT. 10) THEN
+           DO IHEL=1,NCOMB
+             IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1d-8/NCOMB)) THEN 
+               GOODHEL(IHEL) = .true.
+	     ENDIF
+	   ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       wgt_ME_real=ans
       END

=== modified file 'madgraph/iolibs/template_files/super_auto_dsig_group_v4.inc'
--- madgraph/iolibs/template_files/super_auto_dsig_group_v4.inc	2017-03-12 21:28:46 +0000
+++ madgraph/iolibs/template_files/super_auto_dsig_group_v4.inc	2020-01-09 23:16:36 +0000
@@ -69,6 +69,10 @@
 C  
 C GLOBAL VARIABLES
 C  
+      integer nb_spin_state(2)
+      data  nb_spin_state /%(nb_spin_state1)i,%(nb_spin_state2)i/
+      common /nb_hel_state/ nb_spin_state
+
       include 'coupl.inc'
       include 'run.inc'
 C   ICONFIG has this config number
@@ -132,9 +136,11 @@
           ENDDO
         ENDDO
         WRITE(*,*)'Relative summed weights:'
+	if (SUMPROB.ne.0d0)then
         DO J=1,SYMCONF(0)
           WRITE(*,'(%(nsprocs)dE12.4)')((SUMWGT(K,I,J)/SUMPROB,K=1,2),I=1,MAXSPROC)
         ENDDO
+	endif
         SUMPROB=0D0
         DO J=1,SYMCONF(0)
           DO I=1,MAXSPROC
@@ -144,9 +150,11 @@
           ENDDO
         ENDDO
         WRITE(*,*)'Relative number of events:'
+	if (SUMPROB.ne.0d0)then
         DO J=1,SYMCONF(0)
           WRITE(*,'(%(nsprocs)dE12.4)')((NUMEVTS(K,I,J)/SUMPROB,K=1,2),I=1,MAXSPROC)
         ENDDO
+	ENDIF
         WRITE(*,*)'Events:'
         DO J=1,SYMCONF(0)
           WRITE(*,'(%(nsprocs)dI12)')((NUMEVTS(K,I,J),K=1,2),I=1,MAXSPROC)
@@ -451,7 +459,6 @@
 
       END
 
-
 C     -----------------------------------------
 C     Subroutine to map three positive integers
 C     I, J and K with upper bounds J_bound and

=== modified file 'madgraph/loop/loop_base_objects.py'
--- madgraph/loop/loop_base_objects.py	2017-03-27 08:37:47 +0000
+++ madgraph/loop/loop_base_objects.py	2020-01-09 23:16:36 +0000
@@ -1545,7 +1545,8 @@
     def get_sorted_keys(self):
         """Return process property names as a nicely sorted list."""
 
-        return ['id', 'number', 'state', 'from_group','loop_line','depth']
+        return ['id', 'number', 'state', 'from_group','loop_line','depth', 
+                'polarization']
     
     def convert_to_leg(self):
         """ Converts a DGLoopLeg back to a Leg. Basically removes the extra
@@ -1555,6 +1556,7 @@
         for key in aleg.get_sorted_keys():
             aleg.set(key,self[key])
 
+
         return aleg
 
 #===============================================================================

=== modified file 'madgraph/various/banner.py'
--- madgraph/various/banner.py	2019-12-04 21:34:25 +0000
+++ madgraph/various/banner.py	2020-01-09 23:16:36 +0000
@@ -2624,7 +2624,7 @@
                 continue
             else:
                 pathinc = incname
-                
+
             fsock = file_writers.FortranWriter(pjoin(output_dir,pathinc))  
             for key in self.includepath[incname]:                
                 #define the fortran name
@@ -2794,8 +2794,20 @@
  %(e_max_pdg)s = e_max_pdg ! E cut for other particles (syntax e.g. {6: 100, 25: 50})
 """, 
     template_off= '#\n# For display option for energy cut in the partonic center of mass frame type \'update ecut\'\n#'),
-    ]    
-    
+
+
+#    Frame for polarization
+    runblock(name='frame', fields=('me_frame'),
+              template_on=\
+"""#*********************************************************************
+# Frame where to evaluate the matrix-element (not the cut!) for polarization   
+#*********************************************************************
+  %(me_frame)s  = me_frame     ! list of particles to sum-up to define the rest-frame
+                               ! in which to evaluate the matrix-element
+                               ! [1,2] means the partonic center of mass 
+""", 
+    template_off= ''),
+    ]        
     
     def default_setup(self):
         """default value for the run_card.dat"""
@@ -2857,6 +2869,8 @@
         self.add_param("clusinfo", True)
         self.add_param("lhe_version", 3.0)
         self.add_param("boost_event", "False", hidden=True, include=False,      comment="allow to boost the full event. The boost put at rest the sume of 4-momenta of the particle selected by the filter defined here. example going to the higgs rest frame: lambda p: p.pid==25")
+        self.add_param("me_frame", [1,2], hidden=True, include=False, comment="choose lorentz frame where to evaluate matrix-element [for non lorentz invariant matrix-element/polarization]:\n  - 0: partonic center of mass\n - 1: Multi boson frame\n - 2 : (multi) scalar frame\n - 3 : user custom")
+        self.add_param('frame_id', 6,  system=True)
         self.add_param("event_norm", "average", allowed=['sum','average', 'unity'],
                         include=False, sys_default='sum')
         #cut
@@ -3119,6 +3133,9 @@
 
     def update_system_parameter_for_include(self):
         
+        # polarization
+        self['frame_id'] = sum(2**(n) for n in self['me_frame'])
+        
         # set the pdg_for_cut fortran parameter
         pdg_to_cut = set(self['pt_min_pdg'].keys() +self['pt_max_pdg'].keys() + 
                          self['e_min_pdg'].keys() +self['e_max_pdg'].keys() +
@@ -3280,6 +3297,24 @@
         if no_systematics:
             self['use_syst'] = False
             self['systematics_program'] = 'none'
+        
+        # if polarization is used, set the choice of the frame in the run_card
+        # But only if polarization is used for massive particles
+        for plist in proc_def:
+            for proc in plist:
+                for l in proc.get('legs') + proc.get('legs_with_decays'):
+                    if l.get('polarization'):
+                        model = proc.get('model')
+                        particle = model.get_particle(l.get('id'))
+                        if particle.get('mass').lower() != 'zero':
+                            self.display_block.append('frame') 
+                            break
+                else:
+                    continue
+                break
+            else:
+                continue
+            break
 
         if 'MLM' in proc_characteristic['limitations']:
             if self['dynamical_scale_choice'] ==  -1:

=== modified file 'madgraph/various/lhe_parser.py'
--- madgraph/various/lhe_parser.py	2019-11-16 12:01:49 +0000
+++ madgraph/various/lhe_parser.py	2020-01-09 23:16:36 +0000
@@ -1734,11 +1734,14 @@
         """modify the current event to boost it according to the current filter"""
         if filter is None:
             filter = lambda p: p.status==-1
-            
-        pboost = FourMomentum()
-        for p in self:
-            if filter(p):
-                pboost += p
+        
+        if not isinstance(filter, FourMomentum):
+            pboost = FourMomentum()
+            for p in self:
+                if filter(p):
+                    pboost += p
+        else:
+            pboost = FourMomentum(pboost)
 
         # change sign of three-component due to helas convention
         pboost.px *=-1

=== modified file 'madgraph/various/misc.py'
--- madgraph/various/misc.py	2019-11-24 19:59:28 +0000
+++ madgraph/various/misc.py	2020-01-09 23:16:36 +0000
@@ -611,6 +611,8 @@
     else:
         return True
 
+
+
 #===============================================================================
 # mute_logger (designed to work as with statement)
 #===============================================================================

=== modified file 'madgraph/various/process_checks.py'
--- madgraph/various/process_checks.py	2016-08-23 00:33:58 +0000
+++ madgraph/various/process_checks.py	2020-01-09 23:16:36 +0000
@@ -734,7 +734,7 @@
 
         # If directory doesn't exist, skip and return 0
         if not shell_name:
-            logging.info("Directory hasn't been created for process %s" %proc)
+            logging.info("Directory hasn't been created for process %s: %s", proc, directories)
             return ((0.0, 0.0, 0.0, 0.0, 0), [])
 
         if verbose: logging.debug("Working on process %s in dir %s" % (proc, shell_name))

=== modified file 'models/import_ufo.py'
--- models/import_ufo.py	2019-06-20 08:57:47 +0000
+++ models/import_ufo.py	2020-01-09 23:16:36 +0000
@@ -812,10 +812,13 @@
                         particle.set('line', value)
                 elif key == 'propagator':
                     if value:
-                        if aloha.unitary_gauge:
-                            particle.set(key, str(value[0]))
-                        else: 
-                            particle.set(key, str(value[1]))
+                        if isinstance(value, (list,dict)):
+                            if aloha.unitary_gauge:
+                                particle.set(key, str(value[0]))
+                            else: 
+                                particle.set(key, str(value[1]))
+                        else:
+                            particle.set(key, str(value))
                     else:
                         particle.set(key, '')
                 else:

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%born.f	2017-06-02 12:16:50 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%born.f	2020-01-09 23:16:36 +0000
@@ -27,6 +27,7 @@
 C     
       REAL*8 P1(0:3,NEXTERNAL-1)
       COMPLEX*16 ANS(2)
+
 C     
 C     LOCAL VARIABLES 
 C     
@@ -36,6 +37,8 @@
       INTEGER NTRY(1)
       DATA NTRY /1*0/
       INTEGER NHEL(NEXTERNAL-1,NCOMB)
+      REAL*8 T_SAVE(NCOMB,2)
+      REAL*8 T_SUM(2)
       DATA (NHEL(I,   1),I=1,4) /-1,-1,-1, 1/
       DATA (NHEL(I,   2),I=1,4) /-1,-1,-1,-1/
       DATA (NHEL(I,   3),I=1,4) /-1,-1, 1, 1/
@@ -125,6 +128,8 @@
       ANS(1) = 0D0
       ANS(2) = 0D0
       HEL_FAC=1D0
+      T_SUM(1) = 0D0
+      T_SUM(2) = 0D0
       DO IHEL=1,NCOMB
           ! the following lines are to avoid segfaults when glu_ij=0
         COND_IJ=SKIP(NFKSPROCESS).EQ.0
@@ -135,19 +140,33 @@
           IF ((GOODHEL(IHEL,NFKSPROCESS) .OR. GOODHEL(IHEL
      $     +SKIP(NFKSPROCESS),NFKSPROCESS) .OR. NTRY(NFKSPROCESS) .LT.
      $      2) ) THEN
-            ANS(1)=ANS(1)+BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE,BORNS)
-            ANS(2)=ANS(2)+BORNTILDE
-            IF ( BORNS(1).NE.0D0 .AND. .NOT. GOODHEL(IHEL,NFKSPROCESS)
-     $        ) THEN
-              GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
-            ENDIF
-            IF ( BORNS(2).NE.0D0 .AND. .NOT. GOODHEL(IHEL
-     $       +SKIP(NFKSPROCESS),NFKSPROCESS) ) THEN
-              GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
-            ENDIF
+            ANS(1) = ANS(1) + BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE
+     $       ,BORNS)
+            ANS(2)=ANS(2)+ BORNTILDE
+            T_SAVE(IHEL,1) = BORNS(1)
+            T_SAVE(IHEL,2) = BORNS(2)
+            T_SUM(1) = T_SUM(1) + DABS(BORNS(1))
+            T_SUM(2) = T_SUM(2) + DABS(BORNS(2))
           ENDIF
+        ELSE
+          T_SAVE(IHEL,1) = 0D0
+          T_SAVE(IHEL,2) = 0D0
         ENDIF
       ENDDO
+
+      IF(NTRY(NFKSPROCESS).LT.2) THEN
+        IF (T_SUM(2).LT.1E-10*T_SUM(1)) T_SUM(2) = T_SUM(1)
+        DO IHEL=1,NCOMB
+          IF ( DABS(T_SAVE(IHEL,1)).GT.DABS(T_SUM(1))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
+          ENDIF
+
+          IF ( DABS(T_SAVE(IHEL,2)).GT.DABS(T_SUM(2))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS(1)=ANS(1)/DBLE(IDEN)
       ANS(2)=ANS(2)/DBLE(IDEN)
       WGT_ME_BORN=DBLE(ANS(1))

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%born.f	2017-06-02 12:16:50 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%born.f	2020-01-09 23:16:36 +0000
@@ -30,6 +30,7 @@
 C     
       REAL*8 P1(0:3,NEXTERNAL-1)
       COMPLEX*16 ANS(2)
+
 C     
 C     LOCAL VARIABLES 
 C     
@@ -39,6 +40,8 @@
       INTEGER NTRY(1)
       DATA NTRY /1*0/
       INTEGER NHEL(NEXTERNAL-1,NCOMB)
+      REAL*8 T_SAVE(NCOMB,2)
+      REAL*8 T_SUM(2)
       DATA (NHEL(I,   1),I=1,4) / 1,-1,-1, 1/
       DATA (NHEL(I,   2),I=1,4) / 1,-1,-1,-1/
       DATA (NHEL(I,   3),I=1,4) / 1,-1, 1, 1/
@@ -128,6 +131,8 @@
       ANS(1) = 0D0
       ANS(2) = 0D0
       HEL_FAC=1D0
+      T_SUM(1) = 0D0
+      T_SUM(2) = 0D0
       DO IHEL=1,NCOMB
           ! the following lines are to avoid segfaults when glu_ij=0
         COND_IJ=SKIP(NFKSPROCESS).EQ.0
@@ -138,19 +143,33 @@
           IF ((GOODHEL(IHEL,NFKSPROCESS) .OR. GOODHEL(IHEL
      $     +SKIP(NFKSPROCESS),NFKSPROCESS) .OR. NTRY(NFKSPROCESS) .LT.
      $      2) ) THEN
-            ANS(1)=ANS(1)+BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE,BORNS)
-            ANS(2)=ANS(2)+BORNTILDE
-            IF ( BORNS(1).NE.0D0 .AND. .NOT. GOODHEL(IHEL,NFKSPROCESS)
-     $        ) THEN
-              GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
-            ENDIF
-            IF ( BORNS(2).NE.0D0 .AND. .NOT. GOODHEL(IHEL
-     $       +SKIP(NFKSPROCESS),NFKSPROCESS) ) THEN
-              GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
-            ENDIF
+            ANS(1) = ANS(1) + BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE
+     $       ,BORNS)
+            ANS(2)=ANS(2)+ BORNTILDE
+            T_SAVE(IHEL,1) = BORNS(1)
+            T_SAVE(IHEL,2) = BORNS(2)
+            T_SUM(1) = T_SUM(1) + DABS(BORNS(1))
+            T_SUM(2) = T_SUM(2) + DABS(BORNS(2))
           ENDIF
+        ELSE
+          T_SAVE(IHEL,1) = 0D0
+          T_SAVE(IHEL,2) = 0D0
         ENDIF
       ENDDO
+
+      IF(NTRY(NFKSPROCESS).LT.2) THEN
+        IF (T_SUM(2).LT.1E-10*T_SUM(1)) T_SUM(2) = T_SUM(1)
+        DO IHEL=1,NCOMB
+          IF ( DABS(T_SAVE(IHEL,1)).GT.DABS(T_SUM(1))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
+          ENDIF
+
+          IF ( DABS(T_SAVE(IHEL,2)).GT.DABS(T_SUM(2))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS(1)=ANS(1)/DBLE(IDEN)
       ANS(2)=ANS(2)/DBLE(IDEN)
       WGT_ME_BORN=DBLE(ANS(1))

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%born.f	2017-06-02 12:16:50 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%born.f	2020-01-09 23:16:36 +0000
@@ -30,6 +30,7 @@
 C     
       REAL*8 P1(0:3,NEXTERNAL-1)
       COMPLEX*16 ANS(2)
+
 C     
 C     LOCAL VARIABLES 
 C     
@@ -39,6 +40,8 @@
       INTEGER NTRY(1)
       DATA NTRY /1*0/
       INTEGER NHEL(NEXTERNAL-1,NCOMB)
+      REAL*8 T_SAVE(NCOMB,2)
+      REAL*8 T_SUM(2)
       DATA (NHEL(I,   1),I=1,4) /-1, 1,-1, 1/
       DATA (NHEL(I,   2),I=1,4) /-1, 1,-1,-1/
       DATA (NHEL(I,   3),I=1,4) /-1, 1, 1, 1/
@@ -128,6 +131,8 @@
       ANS(1) = 0D0
       ANS(2) = 0D0
       HEL_FAC=1D0
+      T_SUM(1) = 0D0
+      T_SUM(2) = 0D0
       DO IHEL=1,NCOMB
           ! the following lines are to avoid segfaults when glu_ij=0
         COND_IJ=SKIP(NFKSPROCESS).EQ.0
@@ -138,19 +143,33 @@
           IF ((GOODHEL(IHEL,NFKSPROCESS) .OR. GOODHEL(IHEL
      $     +SKIP(NFKSPROCESS),NFKSPROCESS) .OR. NTRY(NFKSPROCESS) .LT.
      $      2) ) THEN
-            ANS(1)=ANS(1)+BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE,BORNS)
-            ANS(2)=ANS(2)+BORNTILDE
-            IF ( BORNS(1).NE.0D0 .AND. .NOT. GOODHEL(IHEL,NFKSPROCESS)
-     $        ) THEN
-              GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
-            ENDIF
-            IF ( BORNS(2).NE.0D0 .AND. .NOT. GOODHEL(IHEL
-     $       +SKIP(NFKSPROCESS),NFKSPROCESS) ) THEN
-              GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
-            ENDIF
+            ANS(1) = ANS(1) + BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE
+     $       ,BORNS)
+            ANS(2)=ANS(2)+ BORNTILDE
+            T_SAVE(IHEL,1) = BORNS(1)
+            T_SAVE(IHEL,2) = BORNS(2)
+            T_SUM(1) = T_SUM(1) + DABS(BORNS(1))
+            T_SUM(2) = T_SUM(2) + DABS(BORNS(2))
           ENDIF
+        ELSE
+          T_SAVE(IHEL,1) = 0D0
+          T_SAVE(IHEL,2) = 0D0
         ENDIF
       ENDDO
+
+      IF(NTRY(NFKSPROCESS).LT.2) THEN
+        IF (T_SUM(2).LT.1E-10*T_SUM(1)) T_SUM(2) = T_SUM(1)
+        DO IHEL=1,NCOMB
+          IF ( DABS(T_SAVE(IHEL,1)).GT.DABS(T_SUM(1))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
+          ENDIF
+
+          IF ( DABS(T_SAVE(IHEL,2)).GT.DABS(T_SUM(2))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS(1)=ANS(1)/DBLE(IDEN)
       ANS(2)=ANS(2)/DBLE(IDEN)
       WGT_ME_BORN=DBLE(ANS(1))

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%born.f	2017-06-02 12:16:50 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%born.f	2020-01-09 23:16:36 +0000
@@ -27,6 +27,7 @@
 C     
       REAL*8 P1(0:3,NEXTERNAL-1)
       COMPLEX*16 ANS(2)
+
 C     
 C     LOCAL VARIABLES 
 C     
@@ -36,6 +37,8 @@
       INTEGER NTRY(8)
       DATA NTRY /8*0/
       INTEGER NHEL(NEXTERNAL-1,NCOMB)
+      REAL*8 T_SAVE(NCOMB,2)
+      REAL*8 T_SUM(2)
       DATA (NHEL(I,   1),I=1,4) /-1,-1,-1, 1/
       DATA (NHEL(I,   2),I=1,4) /-1,-1,-1,-1/
       DATA (NHEL(I,   3),I=1,4) /-1,-1, 1, 1/
@@ -125,6 +128,8 @@
       ANS(1) = 0D0
       ANS(2) = 0D0
       HEL_FAC=1D0
+      T_SUM(1) = 0D0
+      T_SUM(2) = 0D0
       DO IHEL=1,NCOMB
           ! the following lines are to avoid segfaults when glu_ij=0
         COND_IJ=SKIP(NFKSPROCESS).EQ.0
@@ -135,19 +140,33 @@
           IF ((GOODHEL(IHEL,NFKSPROCESS) .OR. GOODHEL(IHEL
      $     +SKIP(NFKSPROCESS),NFKSPROCESS) .OR. NTRY(NFKSPROCESS) .LT.
      $      2) ) THEN
-            ANS(1)=ANS(1)+BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE,BORNS)
-            ANS(2)=ANS(2)+BORNTILDE
-            IF ( BORNS(1).NE.0D0 .AND. .NOT. GOODHEL(IHEL,NFKSPROCESS)
-     $        ) THEN
-              GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
-            ENDIF
-            IF ( BORNS(2).NE.0D0 .AND. .NOT. GOODHEL(IHEL
-     $       +SKIP(NFKSPROCESS),NFKSPROCESS) ) THEN
-              GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
-            ENDIF
+            ANS(1) = ANS(1) + BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE
+     $       ,BORNS)
+            ANS(2)=ANS(2)+ BORNTILDE
+            T_SAVE(IHEL,1) = BORNS(1)
+            T_SAVE(IHEL,2) = BORNS(2)
+            T_SUM(1) = T_SUM(1) + DABS(BORNS(1))
+            T_SUM(2) = T_SUM(2) + DABS(BORNS(2))
           ENDIF
+        ELSE
+          T_SAVE(IHEL,1) = 0D0
+          T_SAVE(IHEL,2) = 0D0
         ENDIF
       ENDDO
+
+      IF(NTRY(NFKSPROCESS).LT.2) THEN
+        IF (T_SUM(2).LT.1E-10*T_SUM(1)) T_SUM(2) = T_SUM(1)
+        DO IHEL=1,NCOMB
+          IF ( DABS(T_SAVE(IHEL,1)).GT.DABS(T_SUM(1))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
+          ENDIF
+
+          IF ( DABS(T_SAVE(IHEL,2)).GT.DABS(T_SUM(2))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS(1)=ANS(1)/DBLE(IDEN)
       ANS(2)=ANS(2)/DBLE(IDEN)
       WGT_ME_BORN=DBLE(ANS(1))

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_1.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_1.f	2016-02-27 14:33:01 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_1.f	2020-01-09 23:16:36 +0000
@@ -74,8 +74,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_1(P ,NHEL(1,IHEL))
@@ -101,11 +101,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_2.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_2.f	2016-02-27 14:33:01 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_2.f	2020-01-09 23:16:36 +0000
@@ -77,8 +77,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_2(P ,NHEL(1,IHEL))
@@ -104,11 +104,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_3.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_3.f	2016-02-27 14:33:01 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_3.f	2020-01-09 23:16:36 +0000
@@ -77,8 +77,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_3(P ,NHEL(1,IHEL))
@@ -104,11 +104,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_4.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_4.f	2016-02-27 14:33:01 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_4.f	2020-01-09 23:16:36 +0000
@@ -77,8 +77,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_4(P ,NHEL(1,IHEL))
@@ -104,11 +104,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_5.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_5.f	2016-02-27 14:33:01 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_5.f	2020-01-09 23:16:36 +0000
@@ -77,8 +77,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_5(P ,NHEL(1,IHEL))
@@ -104,11 +104,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%born.f	2017-06-02 12:16:50 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%born.f	2020-01-09 23:16:36 +0000
@@ -30,6 +30,7 @@
 C     
       REAL*8 P1(0:3,NEXTERNAL-1)
       COMPLEX*16 ANS(2)
+
 C     
 C     LOCAL VARIABLES 
 C     
@@ -39,6 +40,8 @@
       INTEGER NTRY(6)
       DATA NTRY /6*0/
       INTEGER NHEL(NEXTERNAL-1,NCOMB)
+      REAL*8 T_SAVE(NCOMB,2)
+      REAL*8 T_SUM(2)
       DATA (NHEL(I,   1),I=1,4) / 1,-1,-1, 1/
       DATA (NHEL(I,   2),I=1,4) / 1,-1,-1,-1/
       DATA (NHEL(I,   3),I=1,4) / 1,-1, 1, 1/
@@ -128,6 +131,8 @@
       ANS(1) = 0D0
       ANS(2) = 0D0
       HEL_FAC=1D0
+      T_SUM(1) = 0D0
+      T_SUM(2) = 0D0
       DO IHEL=1,NCOMB
           ! the following lines are to avoid segfaults when glu_ij=0
         COND_IJ=SKIP(NFKSPROCESS).EQ.0
@@ -138,19 +143,33 @@
           IF ((GOODHEL(IHEL,NFKSPROCESS) .OR. GOODHEL(IHEL
      $     +SKIP(NFKSPROCESS),NFKSPROCESS) .OR. NTRY(NFKSPROCESS) .LT.
      $      2) ) THEN
-            ANS(1)=ANS(1)+BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE,BORNS)
-            ANS(2)=ANS(2)+BORNTILDE
-            IF ( BORNS(1).NE.0D0 .AND. .NOT. GOODHEL(IHEL,NFKSPROCESS)
-     $        ) THEN
-              GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
-            ENDIF
-            IF ( BORNS(2).NE.0D0 .AND. .NOT. GOODHEL(IHEL
-     $       +SKIP(NFKSPROCESS),NFKSPROCESS) ) THEN
-              GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
-            ENDIF
+            ANS(1) = ANS(1) + BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE
+     $       ,BORNS)
+            ANS(2)=ANS(2)+ BORNTILDE
+            T_SAVE(IHEL,1) = BORNS(1)
+            T_SAVE(IHEL,2) = BORNS(2)
+            T_SUM(1) = T_SUM(1) + DABS(BORNS(1))
+            T_SUM(2) = T_SUM(2) + DABS(BORNS(2))
           ENDIF
+        ELSE
+          T_SAVE(IHEL,1) = 0D0
+          T_SAVE(IHEL,2) = 0D0
         ENDIF
       ENDDO
+
+      IF(NTRY(NFKSPROCESS).LT.2) THEN
+        IF (T_SUM(2).LT.1E-10*T_SUM(1)) T_SUM(2) = T_SUM(1)
+        DO IHEL=1,NCOMB
+          IF ( DABS(T_SAVE(IHEL,1)).GT.DABS(T_SUM(1))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
+          ENDIF
+
+          IF ( DABS(T_SAVE(IHEL,2)).GT.DABS(T_SUM(2))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS(1)=ANS(1)/DBLE(IDEN)
       ANS(2)=ANS(2)/DBLE(IDEN)
       WGT_ME_BORN=DBLE(ANS(1))

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%matrix_1.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%matrix_1.f	2016-02-27 14:33:01 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%matrix_1.f	2020-01-09 23:16:36 +0000
@@ -77,8 +77,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_1(P ,NHEL(1,IHEL))
@@ -104,11 +104,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%matrix_2.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%matrix_2.f	2016-02-27 14:33:01 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%matrix_2.f	2020-01-09 23:16:36 +0000
@@ -77,8 +77,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_2(P ,NHEL(1,IHEL))
@@ -104,11 +104,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%matrix_3.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%matrix_3.f	2016-02-27 14:33:01 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%matrix_3.f	2020-01-09 23:16:36 +0000
@@ -77,8 +77,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_3(P ,NHEL(1,IHEL))
@@ -104,11 +104,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%born.f	2017-06-02 12:16:50 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%born.f	2020-01-09 23:16:36 +0000
@@ -30,6 +30,7 @@
 C     
       REAL*8 P1(0:3,NEXTERNAL-1)
       COMPLEX*16 ANS(2)
+
 C     
 C     LOCAL VARIABLES 
 C     
@@ -39,6 +40,8 @@
       INTEGER NTRY(6)
       DATA NTRY /6*0/
       INTEGER NHEL(NEXTERNAL-1,NCOMB)
+      REAL*8 T_SAVE(NCOMB,2)
+      REAL*8 T_SUM(2)
       DATA (NHEL(I,   1),I=1,4) /-1, 1,-1, 1/
       DATA (NHEL(I,   2),I=1,4) /-1, 1,-1,-1/
       DATA (NHEL(I,   3),I=1,4) /-1, 1, 1, 1/
@@ -128,6 +131,8 @@
       ANS(1) = 0D0
       ANS(2) = 0D0
       HEL_FAC=1D0
+      T_SUM(1) = 0D0
+      T_SUM(2) = 0D0
       DO IHEL=1,NCOMB
           ! the following lines are to avoid segfaults when glu_ij=0
         COND_IJ=SKIP(NFKSPROCESS).EQ.0
@@ -138,19 +143,33 @@
           IF ((GOODHEL(IHEL,NFKSPROCESS) .OR. GOODHEL(IHEL
      $     +SKIP(NFKSPROCESS),NFKSPROCESS) .OR. NTRY(NFKSPROCESS) .LT.
      $      2) ) THEN
-            ANS(1)=ANS(1)+BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE,BORNS)
-            ANS(2)=ANS(2)+BORNTILDE
-            IF ( BORNS(1).NE.0D0 .AND. .NOT. GOODHEL(IHEL,NFKSPROCESS)
-     $        ) THEN
-              GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
-            ENDIF
-            IF ( BORNS(2).NE.0D0 .AND. .NOT. GOODHEL(IHEL
-     $       +SKIP(NFKSPROCESS),NFKSPROCESS) ) THEN
-              GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
-            ENDIF
+            ANS(1) = ANS(1) + BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE
+     $       ,BORNS)
+            ANS(2)=ANS(2)+ BORNTILDE
+            T_SAVE(IHEL,1) = BORNS(1)
+            T_SAVE(IHEL,2) = BORNS(2)
+            T_SUM(1) = T_SUM(1) + DABS(BORNS(1))
+            T_SUM(2) = T_SUM(2) + DABS(BORNS(2))
           ENDIF
+        ELSE
+          T_SAVE(IHEL,1) = 0D0
+          T_SAVE(IHEL,2) = 0D0
         ENDIF
       ENDDO
+
+      IF(NTRY(NFKSPROCESS).LT.2) THEN
+        IF (T_SUM(2).LT.1E-10*T_SUM(1)) T_SUM(2) = T_SUM(1)
+        DO IHEL=1,NCOMB
+          IF ( DABS(T_SAVE(IHEL,1)).GT.DABS(T_SUM(1))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
+          ENDIF
+
+          IF ( DABS(T_SAVE(IHEL,2)).GT.DABS(T_SUM(2))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS(1)=ANS(1)/DBLE(IDEN)
       ANS(2)=ANS(2)/DBLE(IDEN)
       WGT_ME_BORN=DBLE(ANS(1))

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%matrix_1.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%matrix_1.f	2016-02-27 14:33:01 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%matrix_1.f	2020-01-09 23:16:36 +0000
@@ -77,8 +77,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_1(P ,NHEL(1,IHEL))
@@ -104,11 +104,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%matrix_2.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%matrix_2.f	2016-02-27 14:33:01 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%matrix_2.f	2020-01-09 23:16:36 +0000
@@ -77,8 +77,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_2(P ,NHEL(1,IHEL))
@@ -104,11 +104,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%matrix_3.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%matrix_3.f	2016-02-27 14:33:01 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%matrix_3.f	2020-01-09 23:16:36 +0000
@@ -77,8 +77,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_3(P ,NHEL(1,IHEL))
@@ -104,11 +104,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%born.f	2017-06-02 12:16:50 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%born.f	2020-01-09 23:16:36 +0000
@@ -28,6 +28,7 @@
 C     
       REAL*8 P1(0:3,NEXTERNAL-1)
       COMPLEX*16 ANS(2)
+
 C     
 C     LOCAL VARIABLES 
 C     
@@ -37,6 +38,8 @@
       INTEGER NTRY(4)
       DATA NTRY /4*0/
       INTEGER NHEL(NEXTERNAL-1,NCOMB)
+      REAL*8 T_SAVE(NCOMB,2)
+      REAL*8 T_SUM(2)
       DATA (NHEL(I,   1),I=1,3) /-1, 1,-1/
       DATA (NHEL(I,   2),I=1,3) /-1, 1, 0/
       DATA (NHEL(I,   3),I=1,3) /-1, 1, 1/
@@ -122,6 +125,8 @@
       ANS(1) = 0D0
       ANS(2) = 0D0
       HEL_FAC=1D0
+      T_SUM(1) = 0D0
+      T_SUM(2) = 0D0
       DO IHEL=1,NCOMB
           ! the following lines are to avoid segfaults when glu_ij=0
         COND_IJ=SKIP(NFKSPROCESS).EQ.0
@@ -132,19 +137,33 @@
           IF ((GOODHEL(IHEL,NFKSPROCESS) .OR. GOODHEL(IHEL
      $     +SKIP(NFKSPROCESS),NFKSPROCESS) .OR. NTRY(NFKSPROCESS) .LT.
      $      2) ) THEN
-            ANS(1)=ANS(1)+BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE,BORNS)
-            ANS(2)=ANS(2)+BORNTILDE
-            IF ( BORNS(1).NE.0D0 .AND. .NOT. GOODHEL(IHEL,NFKSPROCESS)
-     $        ) THEN
-              GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
-            ENDIF
-            IF ( BORNS(2).NE.0D0 .AND. .NOT. GOODHEL(IHEL
-     $       +SKIP(NFKSPROCESS),NFKSPROCESS) ) THEN
-              GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
-            ENDIF
+            ANS(1) = ANS(1) + BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE
+     $       ,BORNS)
+            ANS(2)=ANS(2)+ BORNTILDE
+            T_SAVE(IHEL,1) = BORNS(1)
+            T_SAVE(IHEL,2) = BORNS(2)
+            T_SUM(1) = T_SUM(1) + DABS(BORNS(1))
+            T_SUM(2) = T_SUM(2) + DABS(BORNS(2))
           ENDIF
+        ELSE
+          T_SAVE(IHEL,1) = 0D0
+          T_SAVE(IHEL,2) = 0D0
         ENDIF
       ENDDO
+
+      IF(NTRY(NFKSPROCESS).LT.2) THEN
+        IF (T_SUM(2).LT.1E-10*T_SUM(1)) T_SUM(2) = T_SUM(1)
+        DO IHEL=1,NCOMB
+          IF ( DABS(T_SAVE(IHEL,1)).GT.DABS(T_SUM(1))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
+          ENDIF
+
+          IF ( DABS(T_SAVE(IHEL,2)).GT.DABS(T_SUM(2))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS(1)=ANS(1)/DBLE(IDEN)
       ANS(2)=ANS(2)/DBLE(IDEN)
       WGT_ME_BORN=DBLE(ANS(1))

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%matrix_1.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%matrix_1.f	2016-03-31 10:53:12 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%matrix_1.f	2020-01-09 23:16:36 +0000
@@ -67,8 +67,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_1(P ,NHEL(1,IHEL))
@@ -94,11 +94,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%matrix_2.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%matrix_2.f	2016-03-31 10:53:12 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%matrix_2.f	2020-01-09 23:16:36 +0000
@@ -67,8 +67,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_2(P ,NHEL(1,IHEL))
@@ -94,11 +94,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%matrix_3.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%matrix_3.f	2016-03-31 10:53:12 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%matrix_3.f	2020-01-09 23:16:36 +0000
@@ -67,8 +67,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_3(P ,NHEL(1,IHEL))
@@ -94,11 +94,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%born.f	2017-06-02 12:16:50 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%born.f	2020-01-09 23:16:36 +0000
@@ -28,6 +28,7 @@
 C     
       REAL*8 P1(0:3,NEXTERNAL-1)
       COMPLEX*16 ANS(2)
+
 C     
 C     LOCAL VARIABLES 
 C     
@@ -37,6 +38,8 @@
       INTEGER NTRY(4)
       DATA NTRY /4*0/
       INTEGER NHEL(NEXTERNAL-1,NCOMB)
+      REAL*8 T_SAVE(NCOMB,2)
+      REAL*8 T_SUM(2)
       DATA (NHEL(I,   1),I=1,3) / 1,-1,-1/
       DATA (NHEL(I,   2),I=1,3) / 1,-1, 0/
       DATA (NHEL(I,   3),I=1,3) / 1,-1, 1/
@@ -122,6 +125,8 @@
       ANS(1) = 0D0
       ANS(2) = 0D0
       HEL_FAC=1D0
+      T_SUM(1) = 0D0
+      T_SUM(2) = 0D0
       DO IHEL=1,NCOMB
           ! the following lines are to avoid segfaults when glu_ij=0
         COND_IJ=SKIP(NFKSPROCESS).EQ.0
@@ -132,19 +137,33 @@
           IF ((GOODHEL(IHEL,NFKSPROCESS) .OR. GOODHEL(IHEL
      $     +SKIP(NFKSPROCESS),NFKSPROCESS) .OR. NTRY(NFKSPROCESS) .LT.
      $      2) ) THEN
-            ANS(1)=ANS(1)+BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE,BORNS)
-            ANS(2)=ANS(2)+BORNTILDE
-            IF ( BORNS(1).NE.0D0 .AND. .NOT. GOODHEL(IHEL,NFKSPROCESS)
-     $        ) THEN
-              GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
-            ENDIF
-            IF ( BORNS(2).NE.0D0 .AND. .NOT. GOODHEL(IHEL
-     $       +SKIP(NFKSPROCESS),NFKSPROCESS) ) THEN
-              GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
-            ENDIF
+            ANS(1) = ANS(1) + BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE
+     $       ,BORNS)
+            ANS(2)=ANS(2)+ BORNTILDE
+            T_SAVE(IHEL,1) = BORNS(1)
+            T_SAVE(IHEL,2) = BORNS(2)
+            T_SUM(1) = T_SUM(1) + DABS(BORNS(1))
+            T_SUM(2) = T_SUM(2) + DABS(BORNS(2))
           ENDIF
+        ELSE
+          T_SAVE(IHEL,1) = 0D0
+          T_SAVE(IHEL,2) = 0D0
         ENDIF
       ENDDO
+
+      IF(NTRY(NFKSPROCESS).LT.2) THEN
+        IF (T_SUM(2).LT.1E-10*T_SUM(1)) T_SUM(2) = T_SUM(1)
+        DO IHEL=1,NCOMB
+          IF ( DABS(T_SAVE(IHEL,1)).GT.DABS(T_SUM(1))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
+          ENDIF
+
+          IF ( DABS(T_SAVE(IHEL,2)).GT.DABS(T_SUM(2))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS(1)=ANS(1)/DBLE(IDEN)
       ANS(2)=ANS(2)/DBLE(IDEN)
       WGT_ME_BORN=DBLE(ANS(1))

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%matrix_1.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%matrix_1.f	2016-03-31 10:53:12 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%matrix_1.f	2020-01-09 23:16:36 +0000
@@ -67,8 +67,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_1(P ,NHEL(1,IHEL))
@@ -94,11 +94,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%matrix_2.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%matrix_2.f	2016-03-31 10:53:12 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%matrix_2.f	2020-01-09 23:16:36 +0000
@@ -67,8 +67,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_2(P ,NHEL(1,IHEL))
@@ -94,11 +94,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%matrix_3.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%matrix_3.f	2016-03-31 10:53:12 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%matrix_3.f	2020-01-09 23:16:36 +0000
@@ -67,8 +67,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_3(P ,NHEL(1,IHEL))
@@ -94,11 +94,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_tdecay_fksreal/%SubProcesses%P0_t_budx%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_tdecay_fksreal/%SubProcesses%P0_t_budx%born.f	2017-06-02 12:16:50 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_tdecay_fksreal/%SubProcesses%P0_t_budx%born.f	2020-01-09 23:16:36 +0000
@@ -28,6 +28,7 @@
 C     
       REAL*8 P1(0:3,NEXTERNAL-1)
       COMPLEX*16 ANS(2)
+
 C     
 C     LOCAL VARIABLES 
 C     
@@ -37,6 +38,8 @@
       INTEGER NTRY(3)
       DATA NTRY /3*0/
       INTEGER NHEL(NEXTERNAL-1,NCOMB)
+      REAL*8 T_SAVE(NCOMB,2)
+      REAL*8 T_SUM(2)
       DATA (NHEL(I,   1),I=1,4) / 1,-1,-1, 1/
       DATA (NHEL(I,   2),I=1,4) / 1,-1,-1,-1/
       DATA (NHEL(I,   3),I=1,4) / 1,-1, 1, 1/
@@ -126,6 +129,8 @@
       ANS(1) = 0D0
       ANS(2) = 0D0
       HEL_FAC=1D0
+      T_SUM(1) = 0D0
+      T_SUM(2) = 0D0
       DO IHEL=1,NCOMB
           ! the following lines are to avoid segfaults when glu_ij=0
         COND_IJ=SKIP(NFKSPROCESS).EQ.0
@@ -136,19 +141,33 @@
           IF ((GOODHEL(IHEL,NFKSPROCESS) .OR. GOODHEL(IHEL
      $     +SKIP(NFKSPROCESS),NFKSPROCESS) .OR. NTRY(NFKSPROCESS) .LT.
      $      2) ) THEN
-            ANS(1)=ANS(1)+BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE,BORNS)
-            ANS(2)=ANS(2)+BORNTILDE
-            IF ( BORNS(1).NE.0D0 .AND. .NOT. GOODHEL(IHEL,NFKSPROCESS)
-     $        ) THEN
-              GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
-            ENDIF
-            IF ( BORNS(2).NE.0D0 .AND. .NOT. GOODHEL(IHEL
-     $       +SKIP(NFKSPROCESS),NFKSPROCESS) ) THEN
-              GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
-            ENDIF
+            ANS(1) = ANS(1) + BORN(P1,NHEL(1,IHEL),IHEL,BORNTILDE
+     $       ,BORNS)
+            ANS(2)=ANS(2)+ BORNTILDE
+            T_SAVE(IHEL,1) = BORNS(1)
+            T_SAVE(IHEL,2) = BORNS(2)
+            T_SUM(1) = T_SUM(1) + DABS(BORNS(1))
+            T_SUM(2) = T_SUM(2) + DABS(BORNS(2))
           ENDIF
+        ELSE
+          T_SAVE(IHEL,1) = 0D0
+          T_SAVE(IHEL,2) = 0D0
         ENDIF
       ENDDO
+
+      IF(NTRY(NFKSPROCESS).LT.2) THEN
+        IF (T_SUM(2).LT.1E-10*T_SUM(1)) T_SUM(2) = T_SUM(1)
+        DO IHEL=1,NCOMB
+          IF ( DABS(T_SAVE(IHEL,1)).GT.DABS(T_SUM(1))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
+          ENDIF
+
+          IF ( DABS(T_SAVE(IHEL,2)).GT.DABS(T_SUM(2))*1D-8/NCOMB)THEN
+            GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS(1)=ANS(1)/DBLE(IDEN)
       ANS(2)=ANS(2)/DBLE(IDEN)
       WGT_ME_BORN=DBLE(ANS(1))

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_tdecay_fksreal/%SubProcesses%P0_t_budx%matrix_1.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_tdecay_fksreal/%SubProcesses%P0_t_budx%matrix_1.f	2016-08-07 07:16:07 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_tdecay_fksreal/%SubProcesses%P0_t_budx%matrix_1.f	2020-01-09 23:16:36 +0000
@@ -75,8 +75,8 @@
       NTRY=NTRY+1
       ANS = 0D0
       DO IHEL=1,NCOMB
-        IF (GOODHEL(IHEL) .OR. NTRY .LT. 2) THEN
-          IF (NTRY.LT.2) THEN
+        IF (GOODHEL(IHEL) .OR. NTRY .LT. 10) THEN
+          IF (NTRY.LT.10) THEN
 C           for the first ps-point, check for helicities that give
 C           identical matrix elements
             T=MATRIX_1(P ,NHEL(1,IHEL))
@@ -102,11 +102,18 @@
           ENDIF
 C         add to the sum of helicities
           ANS=ANS+T
-          IF (T .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
-            GOODHEL(IHEL)=.TRUE.
-          ENDIF
         ENDIF
       ENDDO
+
+      IF ( NTRY .LT. 10) THEN
+        DO IHEL=1,NCOMB
+          IF (.NOT.GOODHEL(IHEL) .AND. (DABS(T_SAVE(IHEL)).GT.ANS*1D-8
+     $     /NCOMB)) THEN
+            GOODHEL(IHEL) = .TRUE.
+          ENDIF
+        ENDDO
+      ENDIF
+
       ANS=ANS/DBLE(IDEN)
       WGT_ME_REAL=ANS
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/auto_dsig.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/auto_dsig.f	2018-03-16 12:10:36 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/auto_dsig.f	2020-01-09 23:16:36 +0000
@@ -71,6 +71,11 @@
       INCLUDE 'coupl.inc'
       INCLUDE 'run.inc'
 C     
+C     local
+C     
+      DOUBLE PRECISION P1(0:3, NEXTERNAL)
+
+C     
 C     DATA
 C     
       DATA U1/1*1D0/
@@ -119,7 +124,12 @@
         DSIG1 = PD(0)
         RETURN
       ENDIF
-      CALL SMATRIX1(PP,DSIGUU)
+      IF(FRAME_ID.NE.6)THEN
+        CALL BOOST_TO_FRAME(PP, FRAME_ID, P1)
+      ELSE
+        P1 = PP
+      ENDIF
+      CALL SMATRIX1(P1,DSIGUU)
       IF (IMODE.EQ.5) THEN
         IF (DSIGUU.LT.1D199) THEN
           DSIG1 = DSIGUU*CONV

=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/matrix1.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/matrix1.f	2018-11-08 22:24:49 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/matrix1.f	2020-01-09 23:16:36 +0000
@@ -39,6 +39,7 @@
       LOGICAL GOODHEL(NCOMB,2)
       INTEGER NTRY(2)
       COMMON/BLOCK_GOODHEL/NTRY,GOODHEL
+
 C     
 C     LOCAL VARIABLES 
 C     
@@ -54,6 +55,7 @@
       INTEGER THIS_NTRY(2)
       SAVE THIS_NTRY
       DATA THIS_NTRY /0,0/
+C     
 C     This is just to temporarily store the reference grid for
 C      helicity of the DiscreteSampler so as to obtain its number of
 C      entries with ref_helicity_grid%%n_tot_entries
@@ -67,6 +69,9 @@
       CHARACTER*101         HEL_BUFF
       COMMON/TO_HELICITY/  HEL_BUFF
 
+      INTEGER NB_SPIN_STATE_IN(2)
+      COMMON /NB_HEL_STATE/ NB_SPIN_STATE_IN
+
       INTEGER IMIRROR
       COMMON/TO_MIRROR/ IMIRROR
 
@@ -146,9 +151,9 @@
             DO JJ=1,NINCOMING
               IF(POL(JJ).NE.1D0.AND.NHEL(JJ,I).EQ.INT(SIGN(1D0,POL(JJ))
      $         )) THEN
-                T=T*ABS(POL(JJ))
+                T=T*ABS(POL(JJ))*NB_SPIN_STATE_IN(JJ)/2D0  ! NB_SPIN_STATE(JJ)/2d0 is added for polarised beam
               ELSE IF(POL(JJ).NE.1D0)THEN
-                T=T*(2D0-ABS(POL(JJ)))
+                T=T*(2D0-ABS(POL(JJ)))*NB_SPIN_STATE_IN(JJ)/2D0
               ENDIF
             ENDDO
             IF (ISUM_HEL.NE.0.AND.DS_GET_DIM_STATUS('Helicity')

=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/super_auto_dsig.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/super_auto_dsig.f	2017-05-04 12:07:34 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/super_auto_dsig.f	2020-01-09 23:16:36 +0000
@@ -72,6 +72,10 @@
 C     
 C     GLOBAL VARIABLES
 C     
+      INTEGER NB_SPIN_STATE(2)
+      DATA  NB_SPIN_STATE /2,2/
+      COMMON /NB_HEL_STATE/ NB_SPIN_STATE
+
       INCLUDE 'coupl.inc'
       INCLUDE 'run.inc'
 C     ICONFIG has this config number
@@ -140,10 +144,12 @@
           ENDDO
         ENDDO
         WRITE(*,*)'Relative summed weights:'
-        DO J=1,SYMCONF(0)
-          WRITE(*,'(4E12.4)')((SUMWGT(K,I,J)/SUMPROB,K=1,2),I=1
-     $     ,MAXSPROC)
-        ENDDO
+        IF (SUMPROB.NE.0D0)THEN
+          DO J=1,SYMCONF(0)
+            WRITE(*,'(4E12.4)')((SUMWGT(K,I,J)/SUMPROB,K=1,2),I=1
+     $       ,MAXSPROC)
+          ENDDO
+        ENDIF
         SUMPROB=0D0
         DO J=1,SYMCONF(0)
           DO I=1,MAXSPROC
@@ -153,10 +159,12 @@
           ENDDO
         ENDDO
         WRITE(*,*)'Relative number of events:'
-        DO J=1,SYMCONF(0)
-          WRITE(*,'(4E12.4)')((NUMEVTS(K,I,J)/SUMPROB,K=1,2),I=1
-     $     ,MAXSPROC)
-        ENDDO
+        IF (SUMPROB.NE.0D0)THEN
+          DO J=1,SYMCONF(0)
+            WRITE(*,'(4E12.4)')((NUMEVTS(K,I,J)/SUMPROB,K=1,2),I=1
+     $       ,MAXSPROC)
+          ENDDO
+        ENDIF
         WRITE(*,*)'Events:'
         DO J=1,SYMCONF(0)
           WRITE(*,'(4I12)')((NUMEVTS(K,I,J),K=1,2),I=1,MAXSPROC)
@@ -487,7 +495,6 @@
 
       END
 
-
 C     -----------------------------------------
 C     Subroutine to map three positive integers
 C     I, J and K with upper bounds J_bound and

=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/auto_dsig.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/auto_dsig.f	2018-04-23 20:23:22 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/auto_dsig.f	2020-01-09 23:16:36 +0000
@@ -69,6 +69,11 @@
       INCLUDE 'coupl.inc'
       INCLUDE 'run.inc'
 C     
+C     local
+C     
+      DOUBLE PRECISION P1(0:3, NEXTERNAL)
+
+C     
 C     DATA
 C     
       DATA U1/1*1D0/
@@ -117,7 +122,12 @@
           DSIG = PD(0)
           RETURN
         ENDIF
-        CALL SMATRIX(PP,DSIGUU)
+        IF(FRAME_ID.NE.6)THEN
+          CALL BOOST_TO_FRAME(PP, FRAME_ID, P1)
+        ELSE
+          P1 = PP
+        ENDIF
+        CALL SMATRIX(P1,DSIGUU)
         IF (IMODE.EQ.5) THEN
           IF (DSIGUU.LT.1D199) THEN
             DSIG = DSIGUU*CONV

=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f	2019-06-26 14:12:13 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f	2020-01-09 23:16:36 +0000
@@ -39,6 +39,9 @@
       LOGICAL GOODHEL(NCOMB)
       INTEGER NTRY
       COMMON/BLOCK_GOODHEL/NTRY,GOODHEL
+      INTEGER NB_SPIN_STATE(2)
+      DATA  NB_SPIN_STATE /2,2/
+      COMMON /NB_HEL_STATE/ NB_SPIN_STATE
 C     
 C     LOCAL VARIABLES 
 C     
@@ -122,9 +125,9 @@
             DO JJ=1,NINCOMING
               IF(POL(JJ).NE.1D0.AND.NHEL(JJ,I).EQ.INT(SIGN(1D0,POL(JJ))
      $         )) THEN
-                T=T*ABS(POL(JJ))
+                T=T*ABS(POL(JJ))*NB_SPIN_STATE(JJ)/2D0  ! NB_SPIN_STATE(JJ)/2d0 is added for polarised beam
               ELSE IF(POL(JJ).NE.1D0)THEN
-                T=T*(2D0-ABS(POL(JJ)))
+                T=T*(2D0-ABS(POL(JJ)))*NB_SPIN_STATE(JJ)/2D0
               ENDIF
             ENDDO
             IF (ISUM_HEL.NE.0) THEN

=== modified file 'tests/parallel_tests/compare_with_old_mg5_version.py'
--- tests/parallel_tests/compare_with_old_mg5_version.py	2017-05-24 22:17:29 +0000
+++ tests/parallel_tests/compare_with_old_mg5_version.py	2020-01-09 23:16:36 +0000
@@ -38,7 +38,7 @@
     """A class to compare the value of a old MG5 version and the current one"""
     
     old_mg5 = None # link to the previous version of MG5 (prevent multiple build)
-    reference_number = 263 #2.4.0
+    reference_number = 310 #2.4.0
     nb_test = 0
     
     
@@ -219,12 +219,14 @@
            Note that if you need to redo this, this is potentially due to a change
            in the model. In consequence, you need to change the old MG5 comparison
            point. (Since the default use another model)."""
-
+        
+        
         return # By default no need this
         self.create_short_parallel_sqso()
         self.create_short_paralel_sm()
         self.create_short_paralel_mssm()
         self.create_short_paralel_heft()
+        self.create_short_polarization()
         self.assertTrue(False)
         
 
@@ -303,6 +305,26 @@
                              filename = "short_sqso.log",
                              pickle_file = pickle_file)
 
+    def create_short_polarization(self):
+        """Test a short list of processes with polarization"""
+        # Create a list of processes to check automatically
+        my_proc_list = ['u u~ > w+{T} w-{T0}', 
+                        ' u d~ > w+{T}, w+ > e+ ve',
+                        ' u d~ > w+{0}, w+ > e+ ve',
+                        ' g g > t{L} t~{R}, t > w+ b, t~ > w- b~'
+                        'W+{T} > e+ ve',
+                        'e+{L} e-{R} > mu+ mu-',
+                        'e+ e- > mu+ mu-{L}']
+                        
+        
+        # Store list of non-zero processes and results in file
+        pickle_file = os.path.join(_pickle_path, 
+                                              "mg5_short_paralleltest_pol.pkl")
+        self.compare_processes(my_proc_list,
+                             model='sm',
+                             orders = {},
+                             filename = "short_pol.log",
+                             pickle_file = pickle_file)
 
     ############################################################################    
     #  ROUTINE FOR THE SHORT TEST (USE by the release script)
@@ -331,6 +353,33 @@
             # Do some cleanup
             my_comp.cleanup()
 
+    def test_short_polarization(self): 
+        """Test a short list of processes with squared order constraints"""
+
+        #self.create_short_polarization()
+        comparisons = me_comparator.PickleRunner.find_comparisons(\
+            os.path.join(_pickle_path, "mg5_short_paralleltest_pol.pkl"))
+
+        for stored_runner in comparisons:
+            # Create a MERunner object for MG5
+            my_mg5 = me_comparator.MG5_UFO_Runner()
+            my_mg5.setup(MG5DIR, MG5DIR)
+
+            # Create and setup a comparator
+            my_comp = me_comparator.MEComparator()
+            my_comp.set_me_runners(stored_runner, my_mg5)
+
+            # Run the actual comparison
+            my_comp.run_comparison(stored_runner.proc_list,'sm',
+                                   stored_runner.orders,
+                                   stored_runner.energy)
+
+            my_comp.assert_processes(self)
+
+            # Do some cleanup
+            my_comp.cleanup()
+            
+
     def test_short_sqso(self): 
         """Test a short list of processes with squared order constraints"""
 
@@ -488,6 +537,8 @@
                              filename = "sm_mini.log",
                              energy = 1000)
         
+        
+        
     def test_mg5_minitest_mssm(self):
         """Test a minimal list of sm 2->2 processes, mainly to test the test"""
         # Create a list of processes to check automatically

=== added file 'tests/parallel_tests/input_files/mg5_short_paralleltest_pol.pkl'
Binary files tests/parallel_tests/input_files/mg5_short_paralleltest_pol.pkl	1970-01-01 00:00:00 +0000 and tests/parallel_tests/input_files/mg5_short_paralleltest_pol.pkl	2020-01-09 23:16:36 +0000 differ
=== modified file 'tests/parallel_tests/test_aloha.py'
--- tests/parallel_tests/test_aloha.py	2018-11-28 22:12:28 +0000
+++ tests/parallel_tests/test_aloha.py	2020-01-09 23:16:36 +0000
@@ -2301,8 +2301,9 @@
             if data1:
                 data1 = data1.simplify()
                 self.assertAlmostEqual(eval(str(data1)), eval(str(data2))) 
+                                
                 
-    def testGammaAlgebraDefinition(self):
+    def test_short_GammaAlgebraDefinition(self):
         """Test the coherence between gamma/gamma5/sigma/projector"""
         Gamma = aloha_obj.Gamma
         Gamma5 = aloha_obj.Gamma5

=== modified file 'tests/unit_tests/core/test_base_objects.py'
--- tests/unit_tests/core/test_base_objects.py	2016-08-07 07:16:07 +0000
+++ tests/unit_tests/core/test_base_objects.py	2020-01-09 23:16:36 +0000
@@ -1010,7 +1010,8 @@
                       'state':True,
                       'from_group':False,
                       'onshell':None,                       
-                      'loop_line':False}
+                      'loop_line':False,
+                      'polarization':[]}
 
         self.myleg = base_objects.Leg(self.mydict)
 
@@ -1085,14 +1086,16 @@
     def test_representation(self):
         """Test leg object string representation."""
 
-        goal = "{\n"
-        goal = goal + "    \'id\': 3,\n"
-        goal = goal + "    \'number\': 5,\n"
-        goal = goal + "    \'state\': True,\n"
-        goal = goal + "    \'from_group\': False,\n" 
-        goal = goal + "    \'loop_line\': False,\n"
-        goal = goal + "    \'onshell\': None\n}"
-
+        goal ="""{
+    'id': 3,
+    'number': 5,
+    'state': True,
+    'from_group': False,
+    'loop_line': False,
+    'onshell': None,
+    'polarization': []
+}"""
+        
         self.assertEqual(goal, str(self.myleg))
 
     def test_leg_list(self):
@@ -1155,7 +1158,8 @@
     def setUp(self):
 
         self.mydict = {'ids':[3, 2, 5],
-                      'state':True}
+                      'state':True,
+                      'polarization':[]}
 
         self.my_multi_leg = base_objects.MultiLeg(self.mydict)
 
@@ -1227,10 +1231,12 @@
     def test_representation(self):
         """Test multi_leg object string representation."""
 
-        goal = "{\n"
-        goal = goal + "    \'ids\': [3, 2, 5],\n"
-        goal = goal + "    \'state\': True\n}"
 
+        goal ="""{
+    'ids': [3, 2, 5],
+    'state': True,
+    'polarization': []
+}"""
         self.assertEqual(goal, str(self.my_multi_leg))
 
     def test_multi_leg_list(self):

=== modified file 'tests/unit_tests/interface/test_cmd.py'
--- tests/unit_tests/interface/test_cmd.py	2020-01-09 09:35:56 +0000
+++ tests/unit_tests/interface/test_cmd.py	2020-01-09 23:16:36 +0000
@@ -204,6 +204,7 @@
         """check if generate format are correctly supported"""
     
         cmd = self.cmd
+        cmd.do_import('sm')
         
         # valid syntax
         cmd.check_process_format('e+ e- > e+ e-')
@@ -211,7 +212,13 @@
         cmd.check_process_format('e+ e- > mu+ ta- / x $y @1')
         cmd.check_process_format('e+ e- > mu+ ta- $ x /y @1')
         cmd.check_process_format('e+ e- > mu+ ta- $ x /y, (e+ > e-, e-> ta) @1')
+        cmd.check_process_format('e+ e- > Z{L}, Z > mu+ mu- @1')
+        cmd.check_process_format('e+ e- > Z{0}, Z > mu+ mu- @1')
+        cmd.check_process_format('e+{L} e- > mu+{L} mu-{R} @1')
+        cmd.check_process_format('e+ e- > t{L} t~ Z{L}, t > mu+ mu- @1')
+        cmd.check_process_format('g g > Z Z [ noborn=QCD] @1')
         cmd.check_process_format('u u~ > 2w+ 2j')
+        cmd.check_process_format('u u~ > 2w+{0} 2j')
         
         # unvalid syntax
         self.wrong(cmd.check_process_format, ' e+ e-')
@@ -224,6 +231,13 @@
         self.wrong(cmd.check_process_format, ' e+ > e+, (e+ > e- / z, e- > top')   
         self.wrong(cmd.check_process_format, 'e+ > ')
         self.wrong(cmd.check_process_format, 'e+ >')
+        self.wrong(cmd.check_process_format, 'e+ e- > Z{L} > mu+ mu-')
+        self.wrong(cmd.check_process_format, 'e+ e- > Z > mu+ mu- / W+{L}')
+        self.wrong(cmd.check_process_format, 'e+ e- > Z > mu+ mu- $ W+{L}')
+        self.wrong(cmd.check_process_format, 'u u~ > t{L} t~ [QCD]')
+        self.wrong(cmd.check_process_format, 'u u~ > W+{L} vl [ QED QCD]')
+        self.wrong(cmd.check_process_format,'u u~ > w+{L} [QCD]')
+        self.wrong(cmd.check_process_format,'u u~ > e+{L} vl [QCD]')
         
     @test_aloha.set_global()
     def test_output_default(self):
@@ -243,6 +257,7 @@
         """check that simple syntax goes trough and return expected process"""
            
         cmd = self.cmd
+        self.do('import model sm')
         self.do('generate 2p > 2j')
         self.assertTrue(cmd._curr_amps)
         proc = cmd._curr_amps[0].get('process').get('legs')
@@ -327,4 +342,14 @@
         output = self.cmd.get_final_part('p p > 2Z')
         self.assertEqual(output, set([23]))        
         
+        output = self.cmd.get_final_part('p p > Z{L} j')
+        self.assertEqual(output, set([1, 2, 3, 4, -1, 21, -4, -3, -2, 23]))         
+
+        output = self.cmd.get_final_part('p p > Z{L} j, Z > e+ e-')
+        self.assertEqual(output, set([1, 2, 3, 4, -1, 21, -4, -3, -2, 11, -11])) 
         
+        output = self.cmd.get_final_part('p p > 2Z{L} ')
+        self.assertEqual(output, set([23]))         
+
+        output = self.cmd.get_final_part('p p > 2Z{L} j, Z > e+ e-')
+        self.assertEqual(output, set([1, 2, 3, 4, -1, 21, -4, -3, -2, 11, -11]))         

=== modified file 'tests/unit_tests/iolibs/test_helas_call_writers.py'
--- tests/unit_tests/iolibs/test_helas_call_writers.py	2018-06-01 08:46:21 +0000
+++ tests/unit_tests/iolibs/test_helas_call_writers.py	2020-01-09 23:16:36 +0000
@@ -497,14 +497,14 @@
 
         wavefunctions = {}
         # IXXXXXX.Key: (spin, state)
-        key1 = ((-2, 0), ('',))
+        key1 = ((-2, 0, tuple()), ('',))
         wavefunctions[key1] = \
                           lambda wf: 'CALL IXXXXX(P(0,%d),%s,NHEL(%d),%d*IC(%d),W(1,%d))' % \
                           (wf.get('number_external'), wf.get('mass'),
                            wf.get('number_external'), -(-1) ** wf.get_with_flow('is_part'),
                            wf.get('number_external'), wf.get('number'))
         # OXXXXXX.Key: (spin, state)
-        key2 = ((2, 0), ('',))
+        key2 = ((2, 0, tuple()), ('',))
         wavefunctions[key2] = \
                           lambda wf: 'CALL OXXXXX(P(0,%d),%s,NHEL(%d),%d*IC(%d),W(1,%d))' % \
                           (wf.get('number_external'), wf.get('mass'),