Merge lp:~maddevelopers/mg5amcnlo/NLO into lp:mg5amcnlo/lts

[Valentin Hirschi]

=============================
|| MadLoop 5 version -1.0 ||
=============================

Major upgrade and modifications with respect to the existing trunk.
This branch basically represents MadLoop5 in its current stage of development and it will take part of the 2.0 release of MadGraph 5 along with MadFKS5 which will ultimately allow to have aMC@NLO incorporated in MadGraph 5 2.0.

This version of MadLoop 5 is ***stable*** and implements:

I) Generation and computation of virtual amplitudes for any process in the full QCD sector of the SM.
  > The model file to be imported for this is 'loop_SM_QCD-full'. (default restriction is not working for loop models yet.)
  > The syntax for process generation is now initial_state > final_state 'born_orders' [pert_orders] 'squared_orders' 'process_options' (Notice that the 'add' functionality still work at NLO)
  > There is only one output mode which is fortran standalone, that you can access as usual with the command 'output standalone'
  > Here is an example of test you could perform with this version:
      a) Compile CutTools once and for all by doing:
          'cd loop_material/CutTools/ && make && cd ../..'
      b) Then run the MG5 shell:
          './bin/mg5'
      c) Run some process, like:
          'import model loop_SM_QCD-full'
          'generate g g > t t~ [QCD]'
          'output standalone Trial'
      d) Then compile the test program 'check_sa.f' in your generated process:
          'cd Trial/SubProcesses/P0_gg_ttx/'
          'make'
          './check'
         You should see some trial PS points being evaluated.
  > Some comments are in order:
      a) The interface needs a lot of polishing and acceptance tests must be written.
      b) Decay chains are supported at NLO in principle (so the structure is thought to be compliant with it), but never tested yet.
      c) Multiple particle labels cannot be used for NLO. (MadFKS will take care of that)
      d) Most massive QCD processes have been crosschecked against MadLoop4, NGLuon and soon GoSam. But the proper systematic parallel tests are still missing.
      e) Mixed order perturbation and loop-induced processes (i.e. no born) are supported for generation only right now, and not for output yet.
      f) Fermion-flows are in principle dealt with, but never tested in practice (never tried majorana yet)
      g) Once CutTools is compiled, all tests should be passing except those related to the loop drawing.

II) What is missing?
  > Loop drawings are there in principle but bugged. Olivier will take care of that as soon as he has some time.
  > Interface upgrade with all NLO process features and the corresponding acceptance tests.
  > The parallel tests for the process I have crosschecked already, like gg>ng, n=2,3,4 and q q~ > q q~, g g > q q~, g g > t t~, etc...
  > The UFO model for the full SM (but notice that it could in principle already be handled by this version, if existing)
  > [...]

===================
|| FINAL COMMENT ||
===================

This is not at all a good-looking finished and polished product but rather a raw sketch of a first version of MadLoop 5 which is stable and complete enough to be brought to light.
I could have waited more for the above mentioned missing items to be implemented but I thought that since the portion of new code in MadLoop5 with respect to the MG5 trunk is already quite big, it is crucial that you guys give me some feedback and remarks on the actual implementation.
This will ease the merging of MadFKS5 and the future upgrade of the existing MG5 trunk.

[lp:~maddevelopers/mg5amcnlo/NLO]

213. By Valentin Hirschi <email address hidden> on 2012-01-05

1. Fixed a bug in the loop model restriction so that it now correctly handles the coupling present in the particle counterterms.
2. Implemented parallel tests for loop processes. (Try it in NLO/tests/parallel_tests/loop_sample_script.py)

214. By Olivier Mattelaer on 2012-02-01

merge with 1.4 version

215. By Olivier Mattelaer on 2012-02-02

correct some of the unit_test failure

216. By Olivier Mattelaer on 2012-02-10

fix a couple of unit tests

217. By mattelaer-olivier on 2012-02-15

try to merge

218. By mattelaer-olivier on 2012-02-15

acceptance tests working

219. By mattelaer-olivier on 2012-02-16

only 3 tests didn't pass

220. By Olivier Mattelaer on 2012-02-16

merge with NLO branch

221. By Valentin Hirschi <email address hidden> on 2012-02-16

1. Upgraded to CutTools v1.6.9

222. By Valentin Hirschi <email address hidden> on 2012-02-21

1. Added two tests for the odd behavior of ALOHA with its new 'tag' setup.

223. By Valentin Hirschi <email address hidden> on 2012-02-21

1. Fixed an error in a modification of the loop_exporters test in last revision.

224. By Valentin Hirschi <email address hidden> on 2012-02-22

1. Fixed aloha problem and clean up for checks.

225. By Valentin Hirschi <email address hidden> on 2012-02-22

1. Merged with Olivier's fix for the aloha problem.

226. By Valentin Hirschi <email address hidden> on 2012-02-27

1. Corrected loop_sample_script.py

227. By Valentin Hirschi <email address hidden> on 2012-02-27

1. Merged with the new interface setup of MG5 2.0

Unmerged revisions

=== modified file 'aloha/aloha_writers.py'
--- aloha/aloha_writers.py	2012-01-31 16:28:22 +0000
+++ aloha/aloha_writers.py	2012-02-27 13:42:21 +0000
@@ -2,7 +2,8 @@
     import madgraph.iolibs.file_writers as writers 
 except:
     import aloha.file_writers as writers
-    
+
+
 import os
 import re 
 from numbers import Number
@@ -10,6 +11,7 @@
 class WriteALOHA: 
     """ Generic writing functions """ 
     
+    LOOP_MODE = False
     power_symbol = '**'
     change_var_format = str
     change_number_format = str
@@ -17,6 +19,7 @@
     type_to_variable = {2:'F',3:'V',5:'T',1:'S'}
     type_to_size = {'S':3, 'T':18, 'V':6, 'F':6}
     
+    
     def __init__(self, abstract_routine, dirpath):
 
 
@@ -37,6 +40,7 @@
         self.offshell = abstract_routine.outgoing 
         self.symmetries = abstract_routine.symmetries
         self.tag = abstract_routine.tag
+        self.loop_routine = 'L' in self.tag
 
         #prepare the necessary object
         self.collect_variables() # Look for the different variables
@@ -310,7 +314,10 @@
         if len(OverM) > 0: 
             str_out += 'double complex ' + ','.join(OverM) + '\n'
         if len(Momenta) > 0:
-            str_out += 'double precision ' + '(0:3),'.join(Momenta) + '(0:3)\n'
+            if self.loop_routine:
+                str_out += 'double complex ' + '(0:3),'.join(Momenta) + '(0:3)\n'
+            else:
+                str_out += 'double precision ' + '(0:3),'.join(Momenta) + '(0:3)\n'
 
         # Add entry for symmetry
         #str_out += '\n'
@@ -344,9 +351,13 @@
             offshelltype = self.particles[self.offshell -1]
             offshell_size = self.type_to_size[offshelltype]            
             #Implement the conservation of Energy Impulsion
-            for i in range(-1,1):
+            if self.LOOP_MODE:
+                max = 3 # each component is one momentum (since it's complex)
+            else:
+                max = 1 # one componenet correspond to 2 component
+            for i in range(-1, max):
                 str_out += '%s%d(%d)= ' % (offshelltype, self.offshell, \
-                                                              offshell_size + i)
+                                                          offshell_size + i)
                 
                 pat=re.compile(r'^[-+]?(?P<spin>\w)')
                 for elem in momentum_conservation:
@@ -363,19 +374,25 @@
             sign = ''
             if self.offshell == index and type in ['V','S']:
                 sign = '-'
-                            
-            str_out += '%s(0) = %s dble(%s%d(%d))\n' % (mom, sign, type, index, energy_pos)
-            str_out += '%s(1) = %s dble(%s%d(%d))\n' % (mom, sign, type, index, energy_pos + 1)
-            str_out += '%s(2) = %s dimag(%s%d(%d))\n' % (mom, sign, type, index, energy_pos + 1)
-            str_out += '%s(3) = %s dimag(%s%d(%d))\n' % (mom, sign, type, index, energy_pos)            
+
+            if self.LOOP_MODE:
+                str_out += '%s(0) = %s %s%d(%d)\n' % (mom, sign, type, index, energy_pos)
+                str_out += '%s(1) = %s %s%d(%d)\n' % (mom, sign, type, index, energy_pos + 1)
+                str_out += '%s(2) = %s %s%d(%d)\n' % (mom, sign, type, index, energy_pos + 2)
+                str_out += '%s(3) = %s %s%d(%d)\n' % (mom, sign, type, index, energy_pos + 3)            
+            else:  
+                str_out += '%s(0) = %s dble(%s%d(%d))\n' % (mom, sign, type, index, energy_pos)
+                str_out += '%s(1) = %s dble(%s%d(%d))\n' % (mom, sign, type, index, energy_pos + 1)
+                str_out += '%s(2) = %s dimag(%s%d(%d))\n' % (mom, sign, type, index, energy_pos + 1)
+                str_out += '%s(3) = %s dimag(%s%d(%d))\n' % (mom, sign, type, index, energy_pos)            
             
                    
-        # Definition for the One Over Mass**2 terms
-        for elem in overm:
-            #Mom is in format OMX with X the number of the particle
-            index = int(elem[2:])
-            str_out += 'OM%d = 0d0\n' % (index)
-            str_out += 'if (M%d .ne. 0d0) OM%d' % (index, index) + '=1d0/M%d**2\n' % (index) 
+                # Definition for the One Over Mass**2 terms
+                for elem in overm:
+                    #Mom is in format OMX with X the number of the particle
+                    index = int(elem[2:])
+                    str_out += 'OM%d = 0d0\n' % (index)
+                    str_out += 'if (M%d .ne. 0d0) OM%d' % (index, index) + '=1d0/M%d**2\n' % (index) 
         
         # Returning result
         return str_out
@@ -429,10 +446,13 @@
             denominator = self.obj.denominator
             for ind in denominator.listindices():
                 denom = self.write_obj(denominator.get_rep(ind))
-            string = 'denom =' + '1d0/(' + denom + ')'
-            string = string.replace('+-', '-')
-            string = re.sub('\((?P<num>[+-]*[0-9])\+(?P<num2>[+-][0-9])[Jj]\)\.', '(\g<num>d0,\g<num2>d0)', string)
-            string = re.sub('(?P<num>[0-9])[Jj]\.', '\g<num>*(0d0,1d0)', string)
+            if self.loop_routine:
+                string = 'denom = 1d0'
+            else:
+                string = 'denom =' + '1d0/(' + denom + ')'
+                string = string.replace('+-', '-')
+                string = re.sub('\((?P<num>[+-]*[0-9])\+(?P<num2>[+-][0-9])[Jj]\)\.', '(\g<num>d0,\g<num2>d0)', string)
+                string = re.sub('(?P<num>[0-9])[Jj]\.', '\g<num>*(0d0,1d0)', string)
             OutString = OutString + string + '\n'
             for ind in numerator.listindices():
                 string = '%s(%d)= COUP*denom*' % (OffShellParticle, self.pass_to_HELAS(ind, start=1))
@@ -509,23 +529,15 @@
         if not offshell:
             text += ' double complex TMP\n'
         else:
-            spin = self.particles[offshell -1] 
-            text += ' double complex TMP(%s)\n integer i' % self.type_to_size[spin]         
+            spin = self.particles[offshell -1]
+            if self.LOOP_MODE:
+                text += ' double complex TMP(%s)\n integer i' % (self.type_to_size[spin]+2)
+            else:
+                text += ' double complex TMP(%s)\n integer i' % self.type_to_size[spin]                   
         
         # Define which part of the routine should be called
         addon = ''.join(self.tag) + '_%s' % self.offshell
         
-#        if 'C' in self.namestring:
-#            short_name, addon = name.split('C',1)
-#            if addon.split('_')[0].isdigit():
-#                addon = 'C' +self.namestring.split('C',1)[1]
-#            elif all([n.isdigit() for n in addon.split('_')[0].split('C')]):
-#                addon = 'C' +self.namestring.split('C',1)[1]
-#            else:
-#                addon = '_%s' % self.offshell
-#        else:
-#            addon = '_%s' % self.offshell
-
         # how to call the routine
         if not offshell:
             main = 'vertex'
@@ -598,6 +610,7 @@
 
 def combine_name(name, other_names, outgoing, tag=None):
     """ build the name for combined aloha function """
+
     
 
     # Two possible scheme FFV1C1_2_X or FFV1__FFV2C1_X

=== modified file 'aloha/create_aloha.py'
--- aloha/create_aloha.py	2012-02-01 20:26:10 +0000
+++ aloha/create_aloha.py	2012-02-27 13:42:21 +0000
@@ -40,6 +40,7 @@
 
 _conjugate_gap = 50
 _spin2_mult = 1000
+LOOP_MODE = False
 
 class ALOHAERROR(Exception): pass
 
@@ -57,8 +58,8 @@
         self.infostr = infostr
         self.symmetries = []
         self.combined = []
+        #self.loop = False # Check if we need the denominator
         self.tag = []
-
         
     def add_symmetry(self, outgoing):
         """ add an outgoing """
@@ -436,9 +437,13 @@
                                                        (lorentzname, outgoing) )
             return None
     
-    def set(self, lorentzname, outgoing, abstract_routine):
+    def set(self, lorentzname, outgoing, abstract_routine, loop=False):
         """ add in the dictionary """
     
+        if loop and not 'L' in abstract_routine.tag:
+            abstract_routine = copy.copy(abstract_routine)
+            abstract_routine.tag = abstract_routine.tag + ['L']
+
         self[(lorentzname, outgoing)] = abstract_routine
     
     def compute_all(self, save=True, wanted_lorentz = []):
@@ -512,15 +517,34 @@
         # reorganize the data (in order to use optimization for a given lorentz
         #structure
         request = {}
+        
+        #Check Loop status
+        aloha_writers.WriteALOHA.LOOP_MODE = any(['L' in tag for l,tag,out in data])
+        # Add loop attribut for those which are not defined
+
         for list_l_name, tag, outgoing in data:
+            conjugate = tuple([int(c[1:]) for c in tag if c.startswith('C')])
+            loop = ('L' in tag)
             for l_name in list_l_name:
                 try:
+<<<<<<< TREE
                     request[l_name][tag].append(outgoing)
+=======
+                    request[l_name][conjugate].append((outgoing, loop))
+>>>>>>> MERGE-SOURCE
                 except:
                     try:
+<<<<<<< TREE
                         request[l_name][tag] = [outgoing]
+=======
+                        request[l_name][conjugate] = [(outgoing, loop)]
+>>>>>>> MERGE-SOURCE
                     except:
+<<<<<<< TREE
                         request[l_name] = {tag: [outgoing]}
+=======
+                        request[l_name] = {conjugate: [(outgoing, loop)]}
+>>>>>>> MERGE-SOURCE
                         
         # Loop on the structure to build exactly what is request
         for l_name in request:
@@ -537,37 +561,43 @@
             
             for conjg in request[l_name]:
                 #ensure that routines are in rising order (for symetries)
-                routines = sorted(request[l_name][conjg])
+                outgoing = sorted([a[0] for a in request[l_name][conjg] if not a[1]])
+                outgoing_loop = sorted([a[0] for a in request[l_name][conjg] if a[1]])
                 if not conjg:
                     # No need to conjugate -> compute directly
-                    self.compute_aloha(builder, routines=routines)
+                    self.compute_aloha(builder, routines=outgoing, 
+                                                routines_loop=outgoing_loop)
                 else:
                     # Define the high level conjugate routine
                     conjg_builder = builder.define_conjugate_builder(conjg)
                     # Compute routines
                     self.compute_aloha(conjg_builder, symmetry=lorentz.name,
-                                        routines=routines)
+                                        routines=outgoing, 
+                                        routines_loop=outgoing_loop )
         
         # Build mutiple lorentz call
-        for list_l_name, conjugate, outgoing in data:
+        for list_l_name, tag, outgoing in data:
             if len(list_l_name) >1:
-                lorentzname = list_l_name[0]
-                for c in conjugate:
-                    lorentzname += 'C%s' % c
+                lorentzname = list_l_name[0] + ''.join(tag)
                 self[(lorentzname, outgoing)].add_combine(list_l_name[1:])
                         
-    def compute_aloha(self, builder, symmetry=None, routines=None):
+    def compute_aloha(self, builder, symmetry=None, routines=None, routines_loop=None):
         """ define all the AbstractRoutine linked to a given lorentz structure
         symmetry authorizes to use the symmetry of anoter lorentz structure.
-        routines to define only a subset of the routines."""
+        routines to define only a subset of the routines.
+        routines loop defines those which should also be written in loop output
+        """
 
         name = builder.name
         if not symmetry:
             symmetry = name
         if not routines:
             routines = range(len(builder.spins) + 1)
-
+        if not routines_loop:
+            routines_loop = []
+  
         # Create the routines
+        # First create the tree-level routines
         for outgoing in routines:
             symmetric = self.has_symmetries(symmetry, outgoing, valid_output=routines)
             if symmetric:
@@ -577,6 +607,19 @@
                 #Store the information
                 realname = name + ''.join(wavefunction.tag)
                 self.set(realname, outgoing, wavefunction)
+        
+        # Creates the Loop routines
+        for outgoing in routines_loop:
+            symmetric = self.has_symmetries(symmetry, outgoing, valid_output=routines_loop)
+            if symmetric:
+                self.get(name+'L', symmetric).add_symmetry(outgoing)
+            elif outgoing in routines:
+                abstract = self.get(name, outgoing)
+                self.set(name+'L', outgoing, abstract, loop=True)
+            else:
+                wavefunction = builder.compute_routine(outgoing)
+                #Store the information
+                self.set(name+'L', outgoing, wavefunction, loop=True)
 
     def compute_aloha_without_kernel(self, builder, symmetry=None, routines=None):
         """define all the AbstractRoutine linked to a given lorentz structure

=== added file 'aloha/template_files/aloha_functions_loop.f'
--- aloha/template_files/aloha_functions_loop.f	1970-01-01 00:00:00 +0000
+++ aloha/template_files/aloha_functions_loop.f	2012-02-27 13:42:21 +0000
@@ -0,0 +1,1690 @@
+C###############################################################################
+C
+C Copyright (c) 2010 The ALOHA Development team and Contributors
+C
+C This file is a part of the MadGraph 5 project, an application which
+C automatically generates Feynman diagrams and matrix elements for arbitrary
+C high-energy processes in the Standard Model and beyond.
+C
+C It is subject to the ALOHA license which should accompany this
+C distribution.
+C
+C###############################################################################
+      subroutine ixxxxx(p, fmass, nhel, nsf ,fi)
+c
+c This subroutine computes a fermion wavefunction with the flowing-IN
+c fermion number.
+c
+c input:
+c       real    p(0:3)         : four-momentum of fermion
+c       real    fmass          : mass          of fermion
+c       integer nhel = -1 or 1 : helicity      of fermion
+c       integer nsf  = -1 or 1 : +1 for particle, -1 for anti-particle
+c
+c output:
+c       complex fi(6)          : fermion wavefunction               |fi>
+c
+      implicit none
+      double complex fi(8),chi(2)
+      double precision p(0:3),sf(2),sfomeg(2),omega(2),fmass,
+     &     pp,pp3,sqp0p3,sqm(0:1)
+      integer nhel,nsf,ip,im,nh
+
+      double precision rZero, rHalf, rTwo
+      parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
+
+c#ifdef HELAS_CHECK
+c      double precision p2
+c      double precision epsi
+c      parameter( epsi = 2.0d-5 )
+c      integer stdo
+c      parameter( stdo = 6 )
+c#endif
+c
+c#ifdef HELAS_CHECK
+c      pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
+c      if ( abs(p(0))+pp.eq.rZero ) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in ixxxxx is zero momentum'
+c      endif
+c      if ( p(0).le.rZero ) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in ixxxxx has non-positive energy'
+c         write(stdo,*)
+c     &        '             : p(0) = ',p(0)
+c      endif
+c      p2 = (p(0)-pp)*(p(0)+pp)
+c      if ( abs(p2-fmass**2).gt.p(0)**2*epsi ) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in ixxxxx has inappropriate mass'
+c         write(stdo,*)
+c     &        '             : p**2 = ',p2,' : fmass**2 = ',fmass**2
+c      endif
+c      if (abs(nhel).ne.1) then
+c         write(stdo,*) ' helas-error : nhel in ixxxxx is not -1,1'
+c         write(stdo,*) '             : nhel = ',nhel
+c      endif
+c      if (abs(nsf).ne.1) then
+c         write(stdo,*) ' helas-error : nsf in ixxxxx is not -1,1'
+c         write(stdo,*) '             : nsf = ',nsf
+c      endif
+c#endif
+
+c     Convention for trees
+c      fi(5) = dcmplx(p(0),p(3))*nsf
+c      fi(6) = dcmplx(p(1),p(2))*nsf
+
+c     Convention for loop computations
+      fi(5) = dcmplx(p(0),0.D0)*nsf
+      fi(6) = dcmplx(p(1),0.D0)*nsf
+      fi(7) = dcmplx(p(2),0.D0)*nsf
+      fi(8) = dcmplx(p(3),0.D0)*nsf
+
+      nh = nhel*nsf
+
+      if ( fmass.ne.rZero ) then
+
+         pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
+
+         if ( pp.eq.rZero ) then
+
+            sqm(0) = dsqrt(abs(fmass)) ! possibility of negative fermion masses
+            sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
+            ip = (1+nh)/2
+            im = (1-nh)/2
+
+            fi(1) = ip     * sqm(ip)
+            fi(2) = im*nsf * sqm(ip)
+            fi(3) = ip*nsf * sqm(im)
+            fi(4) = im     * sqm(im)
+
+         else
+
+            sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
+            sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
+            omega(1) = dsqrt(p(0)+pp)
+            omega(2) = fmass/omega(1)
+            ip = (3+nh)/2
+            im = (3-nh)/2
+            sfomeg(1) = sf(1)*omega(ip)
+            sfomeg(2) = sf(2)*omega(im)
+            pp3 = max(pp+p(3),rZero)
+            chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
+            if ( pp3.eq.rZero ) then
+               chi(2) = dcmplx(-nh )
+            else
+               chi(2) = dcmplx( nh*p(1) , p(2) )/dsqrt(rTwo*pp*pp3)
+            endif
+
+            fi(1) = sfomeg(1)*chi(im)
+            fi(2) = sfomeg(1)*chi(ip)
+            fi(3) = sfomeg(2)*chi(im)
+            fi(4) = sfomeg(2)*chi(ip)
+
+         endif
+
+      else
+
+         if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
+            sqp0p3 = 0d0
+         else
+            sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsf
+         end if
+         chi(1) = dcmplx( sqp0p3 )
+         if ( sqp0p3.eq.rZero ) then
+            chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
+         else
+            chi(2) = dcmplx( nh*p(1), p(2) )/sqp0p3
+         endif
+         if ( nh.eq.1 ) then
+            fi(1) = dcmplx( rZero )
+            fi(2) = dcmplx( rZero )
+            fi(3) = chi(1)
+            fi(4) = chi(2)
+         else
+            fi(1) = chi(2)
+            fi(2) = chi(1)
+            fi(3) = dcmplx( rZero )
+            fi(4) = dcmplx( rZero )
+         endif
+      endif
+c
+      return
+      end
+
+
+      subroutine oxxxxx(p,fmass,nhel,nsf , fo)
+c
+c This subroutine computes a fermion wavefunction with the flowing-OUT
+c fermion number.
+c
+c input:
+c       real    p(0:3)         : four-momentum of fermion
+c       real    fmass          : mass          of fermion
+c       integer nhel = -1 or 1 : helicity      of fermion
+c       integer nsf  = -1 or 1 : +1 for particle, -1 for anti-particle
+c
+c output:
+c       complex fo(6)          : fermion wavefunction               <fo|
+c
+      implicit none
+      double complex fo(8),chi(2)
+      double precision p(0:3),sf(2),sfomeg(2),omega(2),fmass,
+     &     pp,pp3,sqp0p3,sqm(0:1)
+      integer nhel,nsf,nh,ip,im
+
+      double precision rZero, rHalf, rTwo
+      parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
+
+c#ifdef HELAS_CHECK
+c      double precision p2
+c      double precision epsi
+c      parameter( epsi = 2.0d-5 )
+c      integer stdo
+c      parameter( stdo = 6 )
+c#endif
+c
+c#ifdef HELAS_CHECK
+c      pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
+c      if ( abs(p(0))+pp.eq.rZero ) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in oxxxxx is zero momentum'
+c      endif
+c      if ( p(0).le.rZero ) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in oxxxxx has non-positive energy'
+c         write(stdo,*)
+c     &        '         : p(0) = ',p(0)
+c      endif
+c      p2 = (p(0)-pp)*(p(0)+pp)
+c      if ( abs(p2-fmass**2).gt.p(0)**2*epsi ) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in oxxxxx has inappropriate mass'
+c         write(stdo,*)
+c     &        '             : p**2 = ',p2,' : fmass**2 = ',fmass**2
+c      endif
+c      if ( abs(nhel).ne.1 ) then
+c         write(stdo,*) ' helas-error : nhel in oxxxxx is not -1,1'
+c         write(stdo,*) '             : nhel = ',nhel
+c      endif
+c      if ( abs(nsf).ne.1 ) then
+c         write(stdo,*) ' helas-error : nsf in oxxxxx is not -1,1'
+c         write(stdo,*) '             : nsf = ',nsf
+c      endif
+c#endif
+
+c     Convention for trees
+c      fo(5) = dcmplx(p(0),p(3))*nsf
+c      fo(6) = dcmplx(p(1),p(2))*nsf
+
+c     Convention for loop computations
+      fo(5) = dcmplx(p(0),0.D0)*nsf
+      fo(6) = dcmplx(p(1),0.D0)*nsf
+      fo(7) = dcmplx(p(2),0.D0)*nsf
+      fo(8) = dcmplx(p(3),0.D0)*nsf
+
+
+      nh = nhel*nsf
+
+      if ( fmass.ne.rZero ) then
+
+         pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
+
+         if ( pp.eq.rZero ) then
+
+            sqm(0) = dsqrt(abs(fmass)) ! possibility of negative fermion masses
+            sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
+            ip = -((1+nh)/2)
+            im =  (1-nh)/2
+
+            fo(1) = im     * sqm(im)
+            fo(2) = ip*nsf * sqm(im)
+            fo(3) = im*nsf * sqm(-ip)
+            fo(4) = ip     * sqm(-ip)
+
+         else
+
+            pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
+            sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
+            sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
+            omega(1) = dsqrt(p(0)+pp)
+            omega(2) = fmass/omega(1)
+            ip = (3+nh)/2
+            im = (3-nh)/2
+            sfomeg(1) = sf(1)*omega(ip)
+            sfomeg(2) = sf(2)*omega(im)
+            pp3 = max(pp+p(3),rZero)
+            chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
+            if ( pp3.eq.rZero ) then
+               chi(2) = dcmplx(-nh )
+            else
+               chi(2) = dcmplx( nh*p(1) , -p(2) )/dsqrt(rTwo*pp*pp3)
+            endif
+
+            fo(1) = sfomeg(2)*chi(im)
+            fo(2) = sfomeg(2)*chi(ip)
+            fo(3) = sfomeg(1)*chi(im)
+            fo(4) = sfomeg(1)*chi(ip)
+
+         endif
+
+      else
+
+         if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
+            sqp0p3 = 0d0
+         else
+            sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsf
+         end if
+         chi(1) = dcmplx( sqp0p3 )
+         if ( sqp0p3.eq.rZero ) then
+            chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
+         else
+            chi(2) = dcmplx( nh*p(1), -p(2) )/sqp0p3
+         endif
+         if ( nh.eq.1 ) then
+            fo(1) = chi(1)
+            fo(2) = chi(2)
+            fo(3) = dcmplx( rZero )
+            fo(4) = dcmplx( rZero )
+         else
+            fo(1) = dcmplx( rZero )
+            fo(2) = dcmplx( rZero )
+            fo(3) = chi(2)
+            fo(4) = chi(1)
+         endif
+
+      endif
+c
+      return
+      end
+
+      subroutine pxxxxx(p,tmass,nhel,nst , tc)
+
+c    CP3 2009.NOV
+
+c This subroutine computes a PSEUDOR wavefunction.
+c
+c input:
+c       real    p(0:3)         : four-momentum of tensor boson
+c       real    tmass          : mass          of tensor boson
+c       integer nhel           : helicity      of tensor boson
+c                = -2,-1,0,1,2 : (0 is forbidden if tmass=0.0)
+c       integer nst  = -1 or 1 : +1 for final, -1 for initial
+c
+c output:
+c       complex tc(18)         : PSEUDOR  wavefunction    epsilon^mu^nu(t)
+c
+      implicit none
+      double precision p(0:3), tmass
+      integer nhel, nst
+      double complex tc(20)
+
+      double complex ft(6,4), ep(4), em(4), e0(4)
+      double precision pt, pt2, pp, pzpt, emp, sqh, sqs
+      integer i, j
+
+      double precision rZero, rHalf, rOne, rTwo
+      parameter( rZero = 0.0d0, rHalf = 0.5d0 )
+      parameter( rOne = 1.0d0, rTwo = 2.0d0 )
+
+
+      tc(1)=NHEL
+
+c     Convention for trees
+c      tc(17) = dcmplx(p(0),p(3))*nst
+c      tc(18) = dcmplx(p(1),p(2))*nst
+
+c     Convention for loop computations
+      tc(17) = dcmplx(p(0),0.D0)*nst
+      tc(18) = dcmplx(p(1),0.D0)*nst
+      tc(19) = dcmplx(p(2),0.D0)*nst
+      tc(20) = dcmplx(p(3),0.D0)*nst
+
+      return
+      end
+
+      subroutine sxxxxx(p,nss , sc)
+c
+c This subroutine computes a complex SCALAR wavefunction.
+c
+c input:
+c       real    p(0:3)         : four-momentum of scalar boson
+c       integer nss  = -1 or 1 : +1 for final, -1 for initial
+c
+c output:
+c       complex sc(3)          : scalar wavefunction                   s
+c
+      implicit none
+      double complex sc(5)
+      double precision p(0:3)
+      integer nss
+
+      double precision rOne
+      parameter( rOne = 1.0d0 )
+
+c#ifdef HELAS_CHECK
+c      double precision p2
+c      double precision epsi
+c      parameter( epsi = 2.0d-5 )
+c      double precision rZero
+c      parameter( rZero = 0.0d0 )
+c      integer stdo
+c      parameter( stdo = 6 )
+c#endif
+c
+c#ifdef HELAS_CHECK
+c      if ( abs(p(0))+abs(p(1))+abs(p(2))+abs(p(3)).eq.rZero ) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in sxxxxx is zero momentum'
+c      endif
+c      if ( p(0).le.rZero ) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in sxxxxx has non-positive energy'
+c         write(stdo,*)
+c     &        '             : p(0) = ',p(0)
+c      endif
+c      p2 = p(0)**2-(p(1)**2+p(2)**2+p(3)**2)
+c      if ( p2.lt.-p(0)**2*epsi ) then
+c         write(stdo,*) ' helas-error : p(0:3) in sxxxxx is spacelike'
+c         write(stdo,*) '             : p**2 = ',p2
+c      endif
+c      if ( abs(nss).ne.1 ) then
+c         write(stdo,*) ' helas-error : nss in sxxxxx is not -1,1'
+c         write(stdo,*) '             : nss = ',nss
+c      endif
+c#endif
+
+      sc(1) = dcmplx( rOne )
+
+c     Convention for trees
+c      sc(2) = dcmplx(p(0),p(3))*nss
+c      sc(3) = dcmplx(p(1),p(2))*nss
+
+c     Convention for loop computations
+      sc(2) = dcmplx(p(0),0.D0)*nss
+      sc(3) = dcmplx(p(1),0.D0)*nss
+      sc(4) = dcmplx(p(2),0.D0)*nss
+      sc(5) = dcmplx(p(3),0.D0)*nss
+c
+      return
+      end
+
+      subroutine txxxxx(p,tmass,nhel,nst , tc)
+c
+c This subroutine computes a TENSOR wavefunction.
+c
+c input:
+c       real    p(0:3)         : four-momentum of tensor boson
+c       real    tmass          : mass          of tensor boson
+c       integer nhel           : helicity      of tensor boson
+c                = -2,-1,0,1,2 : (0 is forbidden if tmass=0.0)
+c       integer nst  = -1 or 1 : +1 for final, -1 for initial
+c
+c output:
+c       complex tc(18)         : tensor wavefunction    epsilon^mu^nu(t)
+c
+      implicit none
+      double precision p(0:3), tmass
+      integer nhel, nst
+      double complex tc(20)
+
+      double complex ft(6,4), ep(4), em(4), e0(4)
+      double precision pt, pt2, pp, pzpt, emp, sqh, sqs
+      integer i, j
+
+      double precision rZero, rHalf, rOne, rTwo
+      parameter( rZero = 0.0d0, rHalf = 0.5d0 )
+      parameter( rOne = 1.0d0, rTwo = 2.0d0 )
+
+      integer stdo
+      parameter( stdo = 6 )
+
+      sqh = sqrt(rHalf)
+      sqs = sqrt(rHalf/3.d0)
+
+      pt2 = p(1)**2 + p(2)**2
+      pp = min(p(0),sqrt(pt2+p(3)**2))
+      pt = min(pp,sqrt(pt2))
+
+c     Convention for trees
+c      ft(5,1) = dcmplx(p(0),p(3))*nst
+c      ft(6,1) = dcmplx(p(1),p(2))*nst
+
+c     Convention for loop computations
+      ft(5,1) = dcmplx(p(0),0.D0)*nst
+      ft(6,1) = dcmplx(p(1),0.D0)*nst
+      ft(7,1) = dcmplx(p(2),0.D0)*nst
+      ft(8,1) = dcmplx(p(3),0.D0)*nst
+
+      if ( nhel.ge.0 ) then
+c construct eps+
+         if ( pp.eq.rZero ) then
+            ep(1) = dcmplx( rZero )
+            ep(2) = dcmplx( -sqh )
+            ep(3) = dcmplx( rZero , nst*sqh )
+            ep(4) = dcmplx( rZero )
+         else
+            ep(1) = dcmplx( rZero )
+            ep(4) = dcmplx( pt/pp*sqh )
+            if ( pt.ne.rZero ) then
+               pzpt = p(3)/(pp*pt)*sqh
+               ep(2) = dcmplx( -p(1)*pzpt , -nst*p(2)/pt*sqh )
+               ep(3) = dcmplx( -p(2)*pzpt ,  nst*p(1)/pt*sqh )
+            else
+               ep(2) = dcmplx( -sqh )
+               ep(3) = dcmplx( rZero , nst*sign(sqh,p(3)) )
+            endif
+         endif
+      end if
+
+      if ( nhel.le.0 ) then
+c construct eps-
+         if ( pp.eq.rZero ) then
+            em(1) = dcmplx( rZero )
+            em(2) = dcmplx( sqh )
+            em(3) = dcmplx( rZero , nst*sqh )
+            em(4) = dcmplx( rZero )
+         else
+            em(1) = dcmplx( rZero )
+            em(4) = dcmplx( -pt/pp*sqh )
+            if ( pt.ne.rZero ) then
+               pzpt = -p(3)/(pp*pt)*sqh
+               em(2) = dcmplx( -p(1)*pzpt , -nst*p(2)/pt*sqh )
+               em(3) = dcmplx( -p(2)*pzpt ,  nst*p(1)/pt*sqh )
+            else
+               em(2) = dcmplx( sqh )
+               em(3) = dcmplx( rZero , nst*sign(sqh,p(3)) )
+            endif
+         endif
+      end if
+
+      if ( abs(nhel).le.1 ) then
+c construct eps0
+         if ( pp.eq.rZero ) then
+            e0(1) = dcmplx( rZero )
+            e0(2) = dcmplx( rZero )
+            e0(3) = dcmplx( rZero )
+            e0(4) = dcmplx( rOne )
+         else
+            emp = p(0)/(tmass*pp)
+            e0(1) = dcmplx( pp/tmass )
+            e0(4) = dcmplx( p(3)*emp )
+            if ( pt.ne.rZero ) then
+               e0(2) = dcmplx( p(1)*emp )
+               e0(3) = dcmplx( p(2)*emp )
+            else
+               e0(2) = dcmplx( rZero )
+               e0(3) = dcmplx( rZero )
+            endif
+         end if
+      end if
+
+      if ( nhel.eq.2 ) then
+         do j = 1,4
+            do i = 1,4
+               ft(i,j) = ep(i)*ep(j)
+            end do
+         end do
+      else if ( nhel.eq.-2 ) then
+         do j = 1,4
+            do i = 1,4
+               ft(i,j) = em(i)*em(j)
+            end do
+         end do
+      else if (tmass.eq.0) then
+         do j = 1,4
+            do i = 1,4
+               ft(i,j) = 0
+            end do
+         end do
+      else if (tmass.ne.0) then
+        if  ( nhel.eq.1 ) then
+           do j = 1,4
+              do i = 1,4
+                 ft(i,j) = sqh*( ep(i)*e0(j) + e0(i)*ep(j) )
+              end do
+           end do
+        else if ( nhel.eq.0 ) then
+           do j = 1,4
+              do i = 1,4
+                 ft(i,j) = sqs*( ep(i)*em(j) + em(i)*ep(j)
+     &                                + rTwo*e0(i)*e0(j) )
+              end do
+           end do
+        else if ( nhel.eq.-1 ) then
+           do j = 1,4
+              do i = 1,4
+                 ft(i,j) = sqh*( em(i)*e0(j) + e0(i)*em(j) )
+              end do
+           end do
+        else
+           write(stdo,*) 'invalid helicity in TXXXXX'
+           stop
+        end if
+      end if
+
+      tc(1) = ft(1,1)
+      tc(2) = ft(1,2)
+      tc(3) = ft(1,3)
+      tc(4) = ft(1,4)
+      tc(5) = ft(2,1)
+      tc(6) = ft(2,2)
+      tc(7) = ft(2,3)
+      tc(8) = ft(2,4)
+      tc(9) = ft(3,1)
+      tc(10) = ft(3,2)
+      tc(11) = ft(3,3)
+      tc(12) = ft(3,4)
+      tc(13) = ft(4,1)
+      tc(14) = ft(4,2)
+      tc(15) = ft(4,3)
+      tc(16) = ft(4,4)
+      tc(17) = ft(5,1)
+      tc(18) = ft(6,1)
+      tc(19) = ft(5,1)
+      tc(20) = ft(6,1)
+
+      return
+      end
+
+
+      subroutine vxxxxx(p,vmass,nhel,nsv , vc)
+c
+c This subroutine computes a VECTOR wavefunction.
+c
+c input:
+c       real    p(0:3)         : four-momentum of vector boson
+c       real    vmass          : mass          of vector boson
+c       integer nhel = -1, 0, 1: helicity      of vector boson
+c                                (0 is forbidden if vmass=0.0)
+c       integer nsv  = -1 or 1 : +1 for final, -1 for initial
+c
+c output:
+c       complex vc(6)          : vector wavefunction       epsilon^mu(v)
+c
+      implicit none
+      double complex vc(8)
+      double precision p(0:3),vmass,hel,hel0,pt,pt2,pp,pzpt,emp,sqh
+      integer nhel,nsv,nsvahl
+
+      double precision rZero, rHalf, rOne, rTwo
+      parameter( rZero = 0.0d0, rHalf = 0.5d0 )
+      parameter( rOne = 1.0d0, rTwo = 2.0d0 )
+
+c#ifdef HELAS_CHECK
+c      double precision p2
+c      double precision epsi
+c      parameter( epsi = 2.0d-5 )
+c      integer stdo
+c      parameter( stdo = 6 )
+c#endif
+c
+c#ifdef HELAS_CHECK
+c      pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
+c      if ( abs(p(0))+pp.eq.rZero ) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in vxxxxx is zero momentum'
+c      endif
+c      if ( p(0).le.rZero ) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in vxxxxx has non-positive energy'
+c         write(stdo,*)
+c     &        '             : p(0) = ',p(0)
+c      endif
+c      p2 = (p(0)+pp)*(p(0)-pp)
+c      if ( abs(p2-vmass**2).gt.p(0)**2*2.e-5 ) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in vxxxxx has inappropriate mass'
+c         write(stdo,*)
+c     &        '             : p**2 = ',p2,' : vmass**2 = ',vmass**2
+c      endif
+c      if ( vmass.ne.rZero ) then
+c         if ( abs(nhel).gt.1 ) then
+c            write(stdo,*) ' helas-error : nhel in vxxxxx is not -1,0,1'
+c            write(stdo,*) '             : nhel = ',nhel
+c         endif
+c      else
+c         if ( abs(nhel).ne.1 ) then
+c            write(stdo,*) ' helas-error : nhel in vxxxxx is not -1,1'
+c            write(stdo,*) '             : nhel = ',nhel
+c         endif
+c      endif
+c      if ( abs(nsv).ne.1 ) then
+c         write(stdo,*) ' helas-error : nsv in vmxxxx is not -1,1'
+c         write(stdo,*) '             : nsv = ',nsv
+c      endif
+c#endif
+
+      sqh = dsqrt(rHalf)
+      hel = dble(nhel)
+      nsvahl = nsv*dabs(hel)
+      pt2 = p(1)**2+p(2)**2
+      pp = min(p(0),dsqrt(pt2+p(3)**2))
+      pt = min(pp,dsqrt(pt2))
+
+c     Convention for trees
+c      vc(5) = dcmplx(p(0),p(3))*nsv
+c      vc(6) = dcmplx(p(1),p(2))*nsv
+
+c     Convention for loop computations
+      vc(5) = dcmplx(p(0),0.D0)*nsv
+      vc(6) = dcmplx(p(1),0.D0)*nsv
+      vc(7) = dcmplx(p(2),0.D0)*nsv
+      vc(8) = dcmplx(p(3),0.D0)*nsv
+
+c#ifdef HELAS_CHECK
+c nhel=4 option for scalar polarization
+c      if( nhel.eq.4 ) then
+c         if( vmass.eq.rZero ) then
+c            vc(1) = rOne
+c            vc(2) = p(1)/p(0)
+c            vc(3) = p(2)/p(0)
+c            vc(4) = p(3)/p(0)
+c         else
+c            vc(1) = p(0)/vmass
+c            vc(2) = p(1)/vmass
+c            vc(3) = p(2)/vmass
+c            vc(4) = p(3)/vmass
+c         endif
+c         return
+c      endif
+c#endif
+
+      if ( vmass.ne.rZero ) then
+
+         hel0 = rOne-dabs(hel)
+
+         if ( pp.eq.rZero ) then
+
+            vc(1) = dcmplx( rZero )
+            vc(2) = dcmplx(-hel*sqh )
+            vc(3) = dcmplx( rZero , nsvahl*sqh )
+            vc(4) = dcmplx( hel0 )
+
+         else
+
+            emp = p(0)/(vmass*pp)
+            vc(1) = dcmplx( hel0*pp/vmass )
+            vc(4) = dcmplx( hel0*p(3)*emp+hel*pt/pp*sqh )
+            if ( pt.ne.rZero ) then
+               pzpt = p(3)/(pp*pt)*sqh*hel
+               vc(2) = dcmplx( hel0*p(1)*emp-p(1)*pzpt ,
+     &                         -nsvahl*p(2)/pt*sqh       )
+               vc(3) = dcmplx( hel0*p(2)*emp-p(2)*pzpt ,
+     &                          nsvahl*p(1)/pt*sqh       )
+            else
+               vc(2) = dcmplx( -hel*sqh )
+               vc(3) = dcmplx( rZero , nsvahl*sign(sqh,p(3)) )
+            endif
+
+         endif
+
+      else
+
+         pp = p(0)
+         pt = sqrt(p(1)**2+p(2)**2)
+         vc(1) = dcmplx( rZero )
+         vc(4) = dcmplx( hel*pt/pp*sqh )
+         if ( pt.ne.rZero ) then
+            pzpt = p(3)/(pp*pt)*sqh*hel
+            vc(2) = dcmplx( -p(1)*pzpt , -nsv*p(2)/pt*sqh )
+            vc(3) = dcmplx( -p(2)*pzpt ,  nsv*p(1)/pt*sqh )
+         else
+            vc(2) = dcmplx( -hel*sqh )
+            vc(3) = dcmplx( rZero , nsv*sign(sqh,p(3)) )
+         endif
+
+      endif
+c
+      return
+      end
+
+      subroutine boostx(p,q , pboost)
+c
+c This subroutine performs the Lorentz boost of a four-momentum.  The
+c momentum p is assumed to be given in the rest frame of q.  pboost is
+c the momentum p boosted to the frame in which q is given.  q must be a
+c timelike momentum.
+c
+c input:
+c       real    p(0:3)         : four-momentum p in the q rest  frame
+c       real    q(0:3)         : four-momentum q in the boosted frame
+c
+c output:
+c       real    pboost(0:3)    : four-momentum p in the boosted frame
+c
+      implicit none
+      double precision p(0:3),q(0:3),pboost(0:3),pq,qq,m,lf
+
+      double precision rZero
+      parameter( rZero = 0.0d0 )
+
+c#ifdef HELAS_CHECK
+c      integer stdo
+c      parameter( stdo = 6 )
+c      double precision pp
+c#endif
+c
+      qq = q(1)**2+q(2)**2+q(3)**2
+
+c#ifdef HELAS_CHECK
+c      if (abs(p(0))+abs(p(1))+abs(p(2))+abs(p(3)).eq.rZero) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in boostx is zero momentum'
+c      endif
+c      if (abs(q(0))+qq.eq.rZero) then
+c         write(stdo,*)
+c     &        ' helas-error : q(0:3) in boostx is zero momentum'
+c      endif
+c      if (p(0).le.rZero) then
+c         write(stdo,*)
+c     &        ' helas-warn  : p(0:3) in boostx has not positive energy'
+c         write(stdo,*)
+c     &        '             : p(0) = ',p(0)
+c      endif
+c      if (q(0).le.rZero) then
+c         write(stdo,*)
+c     &        ' helas-error : q(0:3) in boostx has not positive energy'
+c         write(stdo,*)
+c     &        '             : q(0) = ',q(0)
+c      endif
+c      pp=p(0)**2-p(1)**2-p(2)**2-p(3)**2
+c      if (pp.lt.rZero) then
+c         write(stdo,*)
+c     &        ' helas-warn  : p(0:3) in boostx is spacelike'
+c         write(stdo,*)
+c     &        '             : p**2 = ',pp
+c      endif
+c      if (q(0)**2-qq.le.rZero) then
+c         write(stdo,*)
+c     &        ' helas-error : q(0:3) in boostx is not timelike'
+c         write(stdo,*)
+c     &        '             : q**2 = ',q(0)**2-qq
+c      endif
+c      if (qq.eq.rZero) then
+c         write(stdo,*)
+c     &   ' helas-warn  : q(0:3) in boostx has zero spacial components'
+c      endif
+c#endif
+
+      if ( qq.ne.rZero ) then
+         pq = p(1)*q(1)+p(2)*q(2)+p(3)*q(3)
+         m = sqrt(q(0)**2-qq)
+         lf = ((q(0)-m)*pq/qq+p(0))/m
+         pboost(0) = (p(0)*q(0)+pq)/m
+         pboost(1) =  p(1)+q(1)*lf
+         pboost(2) =  p(2)+q(2)*lf
+         pboost(3) =  p(3)+q(3)*lf
+      else
+         pboost(0) = p(0)
+         pboost(1) = p(1)
+         pboost(2) = p(2)
+         pboost(3) = p(3)
+      endif
+c
+      return
+      end
+
+      subroutine momntx(energy,mass,costh,phi , p)
+c
+c This subroutine sets up a four-momentum from the four inputs.
+c
+c input:
+c       real    energy         : energy
+c       real    mass           : mass
+c       real    costh          : cos(theta)
+c       real    phi            : azimuthal angle
+c
+c output:
+c       real    p(0:3)         : four-momentum
+c
+      implicit none
+      double precision p(0:3),energy,mass,costh,phi,pp,sinth
+
+      double precision rZero, rOne
+      parameter( rZero = 0.0d0, rOne = 1.0d0 )
+
+c#ifdef HELAS_CHECK
+c      double precision rPi, rTwo
+c      parameter( rPi = 3.14159265358979323846d0, rTwo = 2.d0 )
+c      integer stdo
+c      parameter( stdo = 6 )
+c#endif
+c
+c#ifdef HELAS_CHECK
+c      if (energy.lt.mass) then
+c         write(stdo,*)
+c     &        ' helas-error : energy in momntx is less than mass'
+c         write(stdo,*)
+c     &        '             : energy = ',energy,' : mass = ',mass
+c      endif
+c      if (mass.lt.rZero) then
+c         write(stdo,*) ' helas-error : mass in momntx is negative'
+c         write(stdo,*) '             : mass = ',mass
+c      endif
+c      if (abs(costh).gt.rOne) then
+c         write(stdo,*) ' helas-error : costh in momntx is improper'
+c         write(stdo,*) '             : costh = ',costh
+c      endif
+c      if (phi.lt.rZero .or. phi.gt.rTwo*rPi) then
+c         write(stdo,*)
+c     &   ' helas-warn  : phi in momntx does not lie on 0.0 thru 2.0*pi'
+c         write(stdo,*)
+c     &   '             : phi = ',phi
+c      endif
+c#endif
+
+      p(0) = energy
+
+      if ( energy.eq.mass ) then
+
+         p(1) = rZero
+         p(2) = rZero
+         p(3) = rZero
+
+      else
+
+         pp = sqrt((energy-mass)*(energy+mass))
+         sinth = sqrt((rOne-costh)*(rOne+costh))
+         p(3) = pp*costh
+         if ( phi.eq.rZero ) then
+            p(1) = pp*sinth
+            p(2) = rZero
+         else
+            p(1) = pp*sinth*cos(phi)
+            p(2) = pp*sinth*sin(phi)
+         endif
+
+      endif
+c
+      return
+      end
+      subroutine rotxxx(p,q , prot)
+c
+c This subroutine performs the spacial rotation of a four-momentum.
+c the momentum p is assumed to be given in the frame where the spacial
+c component of q points the positive z-axis.  prot is the momentum p
+c rotated to the frame where q is given.
+c
+c input:
+c       real    p(0:3)         : four-momentum p in q(1)=q(2)=0 frame
+c       real    q(0:3)         : four-momentum q in the rotated frame
+c
+c output:
+c       real    prot(0:3)      : four-momentum p in the rotated frame
+c
+      implicit none
+      double precision p(0:3),q(0:3),prot(0:3),qt2,qt,psgn,qq,p1
+
+      double precision rZero, rOne
+      parameter( rZero = 0.0d0, rOne = 1.0d0 )
+
+c#ifdef HELAS_CHECK
+c      integer stdo
+c      parameter( stdo = 6 )
+c#endif
+c
+      prot(0) = p(0)
+
+      qt2 = q(1)**2 + q(2)**2
+
+c#ifdef HELAS_CHECK
+c      if (abs(p(0))+abs(p(1))+abs(p(2))+abs(p(3)).eq.rZero) then
+c         write(stdo,*)
+c     &        ' helas-error : p(0:3) in rotxxx is zero momentum'
+c      endif
+c      if (abs(q(0))+abs(q(3))+qt2.eq.rZero) then
+c         write(stdo,*)
+c     &        ' helas-error : q(0:3) in rotxxx is zero momentum'
+c      endif
+c      if (qt2+abs(q(3)).eq.rZero) then
+c         write(stdo,*)
+c     &   ' helas-warn  : q(0:3) in rotxxx has zero spacial momentum'
+c      endif
+c#endif
+
+      if ( qt2.eq.rZero ) then
+          if ( q(3).eq.rZero ) then
+             prot(1) = p(1)
+             prot(2) = p(2)
+             prot(3) = p(3)
+          else
+             psgn = dsign(rOne,q(3))
+             prot(1) = p(1)*psgn
+             prot(2) = p(2)*psgn
+             prot(3) = p(3)*psgn
+          endif
+      else
+          qq = sqrt(qt2+q(3)**2)
+          qt = sqrt(qt2)
+          p1 = p(1)
+          prot(1) = q(1)*q(3)/qq/qt*p1 -q(2)/qt*p(2) +q(1)/qq*p(3)
+          prot(2) = q(2)*q(3)/qq/qt*p1 +q(1)/qt*p(2) +q(2)/qq*p(3)
+          prot(3) =          -qt/qq*p1               +q(3)/qq*p(3)
+      endif
+c
+      return
+      end
+
+      subroutine mom2cx(esum,mass1,mass2,costh1,phi1 , p1,p2)
+c
+c This subroutine sets up two four-momenta in the two particle rest
+c frame.
+c
+c input:
+c       real    esum           : energy sum of particle 1 and 2
+c       real    mass1          : mass            of particle 1
+c       real    mass2          : mass            of particle 2
+c       real    costh1         : cos(theta)      of particle 1
+c       real    phi1           : azimuthal angle of particle 1
+c
+c output:
+c       real    p1(0:3)        : four-momentum of particle 1
+c       real    p2(0:3)        : four-momentum of particle 2
+c     
+      implicit none
+      double precision p1(0:3),p2(0:3),
+     &     esum,mass1,mass2,costh1,phi1,md2,ed,pp,sinth1
+
+      double precision rZero, rHalf, rOne, rTwo
+      parameter( rZero = 0.0d0, rHalf = 0.5d0 )
+      parameter( rOne = 1.0d0, rTwo = 2.0d0 )
+
+c#ifdef HELAS_CHECK
+c      double precision rPi
+c      parameter( rPi = 3.14159265358979323846d0 )
+c      integer stdo
+c      parameter( stdo = 6 )
+c#endif
+cc
+c#ifdef HELAS_CHECK
+c      if (esum.lt.mass1+mass2) then
+c         write(stdo,*)
+c     &        ' helas-error : esum in mom2cx is less than mass1+mass2'
+c         write(stdo,*)
+c     &        '             : esum = ',esum,
+c     &        ' : mass1+mass2 = ',mass1,mass2
+c      endif
+c      if (mass1.lt.rZero) then
+c         write(stdo,*) ' helas-error : mass1 in mom2cx is negative'
+c         write(stdo,*) '             : mass1 = ',mass1
+c      endif
+c      if (mass2.lt.rZero) then
+c         write(stdo,*) ' helas-error : mass2 in mom2cx is negative'
+c         write(stdo,*) '             : mass2 = ',mass2
+c      endif
+c      if (abs(costh1).gt.1.) then
+c         write(stdo,*) ' helas-error : costh1 in mom2cx is improper'
+c         write(stdo,*) '             : costh1 = ',costh1
+c      endif
+c      if (phi1.lt.rZero .or. phi1.gt.rTwo*rPi) then
+c         write(stdo,*)
+c     &   ' helas-warn  : phi1 in mom2cx does not lie on 0.0 thru 2.0*pi'
+c         write(stdo,*) 
+c     &   '             : phi1 = ',phi1
+c      endif
+c#endif
+
+      md2 = (mass1-mass2)*(mass1+mass2)
+      ed = md2/esum
+      if ( mass1*mass2.eq.rZero ) then
+         pp = (esum-abs(ed))*rHalf
+      else
+         pp = sqrt((md2/esum)**2-rTwo*(mass1**2+mass2**2)+esum**2)*rHalf
+      endif
+      sinth1 = sqrt((rOne-costh1)*(rOne+costh1))
+
+      p1(0) = max((esum+ed)*rHalf,rZero)
+      p1(1) = pp*sinth1*cos(phi1)
+      p1(2) = pp*sinth1*sin(phi1)
+      p1(3) = pp*costh1
+
+      p2(0) = max((esum-ed)*rHalf,rZero)
+      p2(1) = -p1(1)
+      p2(2) = -p1(2)
+      p2(3) = -p1(3)
+c
+      return
+      end
+
+C###############################################################################
+C LOOP related universal subroutines
+C###############################################################################
+
+C===============================================================================
+C Subroutines to close the lorentz traces of loops 
+C===============================================================================
+
+      SUBROUTINE CLOSE_V(Q,M,AMPS,RES)
+      
+      COMPLEX*16 Q(0:3)
+      COMPLEX*16 RES
+      COMPLEX*16 AMPS(4)
+
+      IF (M.NE.0.D0) THEN
+        STOP 'Massive vector L-cut particle not supported'
+      ENDIF
+      RES=AMPS(1)-AMPS(2)-AMPS(3)-AMPS(4)
+
+      END
+
+c This subroutine is to recreate the fermion propagator with 4 helicities
+c only. This has problems with certain configuration of the imaginary
+c momentum q, so it is not implemented yet.
+
+      SUBROUTINE CLOSE_F4HEL(Q,M,AMPS,RES)      
+      
+      COMPLEX*16 Q(0:3)
+      COMPLEX*16 RES, QNORM
+      REAL*8 M
+      COMPLEX*16 AMPS(4) 
+      COMPLEX*16 PMM, PPM
+
+      PPM=AMPS(1)+AMPS(2)
+      PMM=AMPS(3)+AMPS(4)
+      write(*,*) 'PPM=',PPM
+      write(*,*) 'PMM=',PMM      
+      IF (M.NE.0.D0) THEN
+        QNORM=SQRT(Q(0)**2-Q(1)**2-Q(2)**2-Q(3)**2)
+        write(*,*) 'Q=',Q        
+        write(*,*) 'QNORM=',QNORM
+        write(*,*) 'M=',M
+        RES=(0.D0,0.5D0)*((PPM+PMM)+(PPM-PMM)*(M/QNORM))
+        write(*,*) 'RES1=',RES
+      ELSE
+        RES=(0.D0,0.5D0)*(PPM+PMM)
+        write(*,*) 'RES2=',RES
+      ENDIF
+
+      END
+
+
+      SUBROUTINE CLOSE_F(Q,M,AMPS,RES)      
+      
+      COMPLEX*16 Q(0:3)
+      COMPLEX*16 RES
+      REAL*8 M
+      COMPLEX*16 AMPS(12) 
+
+      RES=(0.D0,0.D0)
+      RES=(Q(0)-Q(3))*AMPS(1)+
+     &    (-Q(1)+(0.0d0,1.0d0)*Q(2))*AMPS(2)+
+     &    (-Q(1)-(0.0d0,1.0d0)*Q(2))*AMPS(3)+
+     &    (Q(0)+Q(3))*AMPS(4)+
+     &    (Q(0)+Q(3))*AMPS(5)+
+     &    (Q(1)-(0.0d0,1.0d0)*Q(2))*AMPS(6)+
+     &    (Q(1)+(0.0d0,1.0d0)*Q(2))*AMPS(7)+
+     &    (Q(0)-Q(3))*AMPS(8)
+      IF (M.NE.0.D0) THEN
+        RES=RES-M*(AMPS(9)+AMPS(10)+AMPS(11)+AMPS(12))
+      ENDIF
+
+      END
+
+      SUBROUTINE CLOSE_S(Q,AMP,RES)
+
+      COMPLEX*16 Q(0:3)
+      COMPLEX*16 RES
+      COMPLEX*16 AMP
+
+      RES=(1.D0,0.D0)*AMP
+      
+      END
+
+C===============================================================================
+C Subroutines to create the external wavefunctions of the L-cut particles
+C===============================================================================
+
+c This subroutine is to recreate the fermion propagator with 4 helicities
+c only. This has problems with certain configuration of the imaginary
+c momentum q, so it is not implemented yet.
+
+      SUBROUTINE LCUT_F4HEL(Q,M,CFIG,SCD,W)
+
+      COMPLEX*16 Q(0:3)
+      INTEGER CFIG,J
+      LOGICAL SCD
+      REAL*8 M
+      COMPLEX*16 W(20)
+
+      IF (CFIG.EQ.1) THEN
+        IF (SCD) THEN
+C         UBAR, HEL=-1        
+          CALL ILXXXX(Q(0),M,-1,1,W(1))
+        ELSE
+C         U, HEL=-1
+          CALL ILXXXX(Q(0),M,-1,-1,W(1))
+          do J=1,4
+          write(*,*) 'Wcf(',j,',1)=',W(1)   
+          enddo
+        ENDIF
+      ELSEIF (CFIG.EQ.2) THEN
+        IF (SCD) THEN
+C         UBAR, HEL=1        
+          CALL ILXXXX(Q(0),M,1,1,W(1))
+        ELSE
+C         U, HEL=1        
+          CALL ILXXXX(Q(0),M,1,-1,W(1))
+        ENDIF
+      ELSEIF (CFIG.EQ.3) THEN
+        IF (SCD) THEN
+C         VBAR, HEL=-1,        
+          CALL OLXXXX(Q(0),M,-1,-1,W(1))
+        ELSE
+C         V, HEL=-1        
+          CALL OLXXXX(Q(0),M,-1,1,W(1))
+        ENDIF
+      ELSEIF (CFIG.EQ.4) THEN
+        IF (SCD) THEN
+C         VBAR, HEL=1        
+          CALL OLXXXX(Q(0),M,1,-1,W(1))
+        ELSE
+C         V, HEL=1    
+          CALL OLXXXX(Q(0),M,1,1,W(1))
+        ENDIF
+      ENDIF
+C     REVERSE THE MOMENTUM IN THE WF FOR THE SECOND L-CUT SPINORS      
+      IF (SCD) THEN
+        W(5)=-Q(0)
+        W(6)=-Q(1)
+        W(7)=-Q(2)
+        W(8)=-Q(3)
+      ENDIF
+
+      END
+
+      SUBROUTINE LCUT_F(Q,M,CFIG,SCD,W)
+
+      COMPLEX*16 Q(0:3)
+      INTEGER CFIG
+      LOGICAL SCD
+      REAL*8 M
+      COMPLEX*16 W(20)
+
+      W(1)=(0.d0,0.d0)
+      W(2)=(0.d0,0.d0)
+      W(3)=(0.d0,0.d0)
+      W(4)=(0.d0,0.d0)
+
+      IF (CFIG.eq.1) then
+        IF (SCD) then
+          W(3)=(1.d0,0.d0)
+        ELSE
+          W(1)=(1.d0,0.d0)
+        ENDIF
+      ELSEIF (CFIG.eq.2) then
+        IF (SCD) then
+          W(4)=(1.d0,0.d0)
+        ELSE
+          W(1)=(1.d0,0.d0)
+        ENDIF
+      ELSEIF (CFIG.eq.3) then
+        IF (SCD) then
+          W(3)=(1.d0,0.d0)
+        ELSE
+          W(2)=(1.d0,0.d0)
+        ENDIF  
+      ELSEIF (CFIG.eq.4) then
+        IF (SCD) then
+          W(4)=(1.d0,0.d0)
+        ELSE
+          W(2)=(1.d0,0.d0)
+        ENDIF
+      ELSEIF (CFIG.eq.5) then
+        IF (SCD) then
+          W(1)=(1.d0,0.d0)
+        ELSE
+          W(3)=(1.d0,0.d0)
+        ENDIF
+      ELSEIF (CFIG.eq.6) then
+        IF (SCD) then
+          W(2)=(1.d0,0.d0)
+        ELSE
+          W(3)=(1.d0,0.d0)
+        ENDIF
+      ELSEIF (CFIG.eq.7) then
+        IF (SCD) then
+          W(1)=(1.d0,0.d0)
+        ELSE
+          W(4)=(1.d0,0.d0)
+        ENDIF
+      ELSEIF (CFIG.eq.8) then
+        IF (SCD) then
+          W(2)=(1.d0,0.d0)
+        ELSE
+          W(4)=(1.d0,0.d0)
+        ENDIF
+      ELSEIF (CFIG.eq.9) then
+        IF (SCD) then
+          W(1)=(1.d0,0.d0)
+        ELSE
+          W(1)=(1.d0,0.d0)
+        ENDIF
+      ELSEIF (CFIG.eq.10) then
+        IF (SCD) then
+          W(2)=(1.d0,0.d0)
+        ELSE
+          W(2)=(1.d0,0.d0)
+        ENDIF
+      ELSEIF (CFIG.eq.11) then
+        IF (SCD) then
+          W(3)=(1.d0,0.d0)
+        ELSE
+          W(3)=(1.d0,0.d0)
+        ENDIF
+      ELSEIF (CFIG.eq.12) then
+        IF (SCD) then
+          W(4)=(1.d0,0.d0)
+        ELSE
+          W(4)=(1.d0,0.d0)
+        ENDIF
+      ENDIF
+C     REVERSE THE MOMENTUM IN THE WF FOR THE SECOND L-CUT SPINORS      
+      IF (SCD) THEN
+        W(5)=-Q(0)
+        W(6)=-Q(1)
+        W(7)=-Q(2)
+        W(8)=-Q(3)
+      ELSE
+        W(5)=Q(0)
+        W(6)=Q(1)
+        W(7)=Q(2)
+        W(8)=Q(3)
+      ENDIF
+
+      END
+
+      SUBROUTINE LCUT_CF(Q,M,CFIG,SCD,W)
+
+      COMPLEX*16 Q(0:3)
+      INTEGER CFIG
+      LOGICAL SCD
+      REAL*8 M
+      COMPLEX*16 W(20)
+
+      IF (CFIG.EQ.1) THEN
+        IF (SCD) THEN
+C         UBAR, HEL=-1        
+          CALL ICLXXX(Q(0),M,-1,1,W(1))
+        ELSE
+C         U, HEL=-1        
+          CALL ICLXXX(Q(0),M,-1,-1,W(1))
+        ENDIF
+      ELSEIF (CFIG.EQ.2) THEN
+        IF (SCD) THEN
+C         UBAR, HEL=1        
+          CALL ICLXXX(Q(0),M,1,1,W(1))
+        ELSE
+C         U, HEL=1        
+          CALL ICLXXX(Q(0),M,1,-1,W(1))
+        ENDIF
+      ELSEIF (CFIG.EQ.3) THEN
+        IF (SCD) THEN
+C         VBAR, HEL=-1,        
+          CALL OCLXXX(Q(0),M,-1,-1,W(1))
+        ELSE
+C         V, HEL=-1        
+          CALL OCLXXX(Q(0),M,-1,1,W(1))
+        ENDIF
+      ELSEIF (CFIG.EQ.4) THEN
+        IF (SCD) THEN
+C         VBAR, HEL=1        
+          CALL OCLXXX(Q(0),M,1,-1,W(1))
+        ELSE
+C         V, HEL=1    
+          CALL OCLXXX(Q(0),M,1,1,W(1))
+        ENDIF
+      ENDIF
+C     REVERSE THE MOMENTUM IN THE WF FOR THE SECOND L-CUT SPINORS      
+      IF (SCD) THEN
+        W(5)=-Q(0)
+        W(6)=-Q(1)
+        W(7)=-Q(2)
+        W(8)=-Q(3)
+      ENDIF
+
+      END
+
+      SUBROUTINE LCUT_V(Q,M,CFIG,SCD,W)
+
+      COMPLEX*16 Q(0:3)
+      INTEGER CFIG
+      LOGICAL SCD      
+      REAL*8 M
+      COMPLEX*16 W(20)
+
+      IF (M.NE.0.D0) THEN
+        STOP 'Massive vector L-cut particle not supported'
+      ENDIF
+      W(1)=(0.d0,0.d0)
+      W(2)=(0.d0,0.d0)
+      W(3)=(0.d0,0.d0)
+      W(4)=(0.d0,0.d0)
+      IF (CFIG.EQ.1) THEN
+        W(1)=(1.d0,0.d0)
+      ELSEIF (CFIG.EQ.2) THEN
+        W(2)=(1.d0,0.d0)
+      ELSEIF (CFIG.EQ.3) THEN
+        W(3)=(1.d0,0.d0)
+      ELSEIF (CFIG.EQ.4) THEN
+        W(4)=(1.d0,0.d0)
+      ENDIF
+C     REVERSE THE MOMENTUM IN THE WF FOR THE SECOND L-CUT VECTORS      
+      IF (SCD) THEN
+        W(5)=-Q(0)
+        W(6)=-Q(1)
+        W(7)=-Q(2)
+        W(8)=-Q(3)
+      ELSE
+        W(5)=Q(0)
+        W(6)=Q(1)
+        W(7)=Q(2)
+        W(8)=Q(3)
+      ENDIF
+
+      END
+
+      SUBROUTINE LCUT_S(Q,M,CFIG,SCD,W)
+
+      COMPLEX*16 Q(0:3)
+      REAL*8 M
+      INTEGER CFIG
+      LOGICAL SCD      
+      COMPLEX*16 W(20)
+
+      W(1)=(1.D0,0.D0)
+
+C     REVERSE THE MOMENTUM IN THE WF FOR THE SECOND SCALAR      
+      IF (SCD) THEN
+        W(2)=-Q(0)
+        W(3)=-Q(1)
+        W(4)=-Q(2)
+        W(5)=-Q(3)
+      ELSE
+        W(2)=Q(0)
+        W(3)=Q(1)
+        W(4)=Q(2)
+        W(5)=Q(3)
+      ENDIF
+
+      END
+
+C===============================================================================
+C Complex-momentum version of the subroutine to create on-shell
+C wavefunctions of particles with different spins.
+C===============================================================================
+
+      subroutine ilxxxx(p,ffmass,nhel,nsf,fi)
+c
+c This subroutine computes a fermion wavefunction with the flowing-IN
+c fermion number and defined with complex ONSHELL momentium.
+c
+c input:
+c       complex    p(0:3)         : four-momentum of fermion
+c       real    fmass          : mass          of fermion
+c       integer nhel = -1 or 1 : helicity      of fermion
+c       integer nsf  = -1 or 1 : +1 for particle, -1 for anti-particle
+c
+c output:
+c       complex fi(8)          : fermion wavefunction               |fi>
+c       Note: There are 4 components for the spinor and four for the
+c             momentum. 
+      implicit none
+      double complex fi(8),chi(2), fmass
+c      double precision p(0:3),sf(2),sfomeg(2),omega(2),fmass,
+c     &     pp,pp3,sqp0p3,sqm(0:1)
+      double complex sqm(0:1)
+      double precision sf(2),ffmass
+      double complex p(0:3), sfomeg(2),omega(2),
+     &     pp,pp3,sqp0p3 
+      integer nhel,nsf,ip,im,nh
+
+      double precision rZero, rHalf, rTwo
+      parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
+      
+
+      
+c      fi(5) = dcmplx(p(0),p(3))*nsf
+c      fi(6) = dcmplx(p(1),p(2))*nsf
+      fi(5) = p(0)*nsf
+      fi(6) = p(1)*nsf
+      fi(7) = p(2)*nsf
+      fi(8) = p(3)*nsf
+
+      nh = nhel*nsf
+
+      fmass = sqrt(p(0)**2-p(1)**2-p(2)**2-p(3)**2)
+
+      if ( ffmass.ne.rZero ) then
+c special treatment for massless particles.
+c         pp = min(p(0),sqrt(p(1)**2+p(2)**2+p(3)**2))
+         pp=sqrt(p(1)**2+p(2)**2+p(3)**2)
+c for time-like four-momenta we can always think of it as the p_vec^2          
+         if ( abs(pp).eq.rZero ) then
+c particle at rest.            
+            sqm(0) = sqrt(fmass) ! possibility of negative fermion masses
+            sqm(1) = sqm(0) ! possibility of negative fermion masses
+            ip = (1+nh)/2
+            im = (1-nh)/2
+            
+            fi(1) = ip     * sqm(ip)
+            fi(2) = im*nsf * sqm(ip)
+            fi(3) = ip*nsf * sqm(im)
+            fi(4) = im     * sqm(im)
+
+         else
+c standard spinor
+
+            pp=sqrt(p(1)**2+p(2)**2+p(3)**2)
+            write(*,*) 'ppre=',pp            
+c            if( (dble(p(0)) .lt. 0 .and. dble(pp) .gt. 0) .or.
+c     &          (dble(p(0)) .lt. 0 .and. dble(pp) .gt. 0) ) then
+c            pp=-pp
+c            endif 
+            sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
+c fermion spin using HELAS conventions.            
+            sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
+            omega(1) = sqrt(p(0)+pp)
+c the omega of the definition.
+c            omega(2) = fmass/omega(1)
+            omega(2) = sqrt(p(0)-pp)
+c the prefactor
+            ip = (3+nh)/2
+            im = (3-nh)/2
+            sfomeg(1) = sf(1)*omega(ip)
+            sfomeg(2) = sf(2)*omega(im)
+c            pp3 = max(pp+p(3),rZero)
+            pp3=pp+p(3)
+            chi(1) = sqrt(pp3*rHalf/pp)
+            if ( abs(pp3).eq.rZero ) then
+               chi(2) = dcmplx(-nh )
+            else
+               chi(2) = ( (nh*p(1)) + ((0d0,1d0)*p(2)) )/
+     .sqrt(rTwo*pp*pp3)
+            endif
+            
+
+c         Write(*,*) 'Chi=',Chi(1),' and ',Chi(2)
+
+            fi(1) = sfomeg(1)*chi(im)
+            fi(2) = sfomeg(1)*chi(ip)
+c         Write(*,*) 'fi(2)=',fi(2)     
+            fi(3) = sfomeg(2)*chi(im)
+c         Write(*,*) 'fi(3)=',fi(3)
+            fi(4) = sfomeg(2)*chi(ip)
+            
+         endif
+         
+      else
+         
+c         if(zabs(p(1)).eq.0d0.and.zabs(p(2)).eq.0d0.and.
+c     .zabs(p(3)).lt.0d0) then
+c            sqp0p3 = 0d0
+c         else
+            sqp0p3 = sqrt(p(0)+p(3))*nsf
+c        end if
+         chi(1) = sqp0p3
+         if ( abs(sqp0p3).eq.rZero ) then
+            chi(2) = dcmplx(-nhel )*sqrt(rTwo*p(0))
+         else
+            chi(2) = ( nh*p(1) + ((0d0,1d0)*p(2) ) )/sqp0p3
+         endif
+         if ( nh.eq.1 ) then
+            fi(1) = dcmplx( rZero )
+            fi(2) = dcmplx( rZero )
+            fi(3) = chi(1)
+            fi(4) = chi(2)
+         else
+            fi(1) = chi(2)
+            fi(2) = chi(1)
+            fi(3) = dcmplx( rZero )
+            fi(4) = dcmplx( rZero )
+         endif
+      endif     
+
+      return
+      end
+
+      subroutine olxxxx(p,ffmass,nhel,nsf,fo)
+c
+c This subroutine computes a fermion wavefunction with the flowing-OUT
+c fermion number and defined with complex ONSHELL  momentum.
+c
+c input:
+c       complex    p(0:3)         : four-momentum of fermion
+c       real    fmass          : mass          of fermion
+c       integer nhel = -1 or 1 : helicity      of fermion
+c       integer nsf  = -1 or 1 : +1 for particle, -1 for anti-particle
+c
+c output:
+c       complex fo(8)          : fermion wavefunction               <fo|
+c     
+      implicit none
+      double complex fo(8),chi(2), fmass
+
+      double complex sqm(0:1)
+      double precision sf(2), ffmass
+
+      double complex p(0:3),sfomeg(2),omega(2),
+     &     pp,pp3,sqp0p3
+      integer nhel,nsf,nh,ip,im
+
+      double precision rZero, rHalf, rTwo
+      parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
+
+
+c      fo(5) = dcmplx(p(0),p(3))*nsf
+c      fo(6) = dcmplx(p(1),p(2))*nsf
+      fo(5) = p(0)*(nsf)
+      fo(6) = p(1)*(nsf)
+      fo(7) = p(2)*(nsf)
+      fo(8) = p(3)*(nsf)      
+
+      nh = nhel*nsf
+
+      fmass=sqrt(p(0)**2-p(1)**2-p(2)**2-p(3)**2)
+
+      if ( ffmass.ne.rZero ) then
+
+c         pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
+         pp=sqrt(p(1)**2+p(2)**2+p(3)**2)
+
+         if ( abs(pp).eq.rZero ) then
+            
+            sqm(0) = sqrt(fmass) ! possibility of negative fermion masses
+            sqm(1) = sqm(0) ! possibility of negative fermion masses
+            ip = -((1+nh)/2)
+            im =  (1-nh)/2
+            
+            fo(1) = im     * sqm(im)
+            fo(2) = ip*nsf * sqm(im)
+            fo(3) = im*nsf * sqm(-ip)
+            fo(4) = ip     * sqm(-ip)
+            
+         else
+            
+c            pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
+            pp=sqrt(p(1)**2+p(2)**2+p(3)**2) !repetition
+c            if( (dble(p(0)) .lt. 0 .and. dble(pp) .gt. 0) .or.
+c     &          (dble(p(0)) .lt. 0 .and. dble(pp) .gt. 0) ) then
+c            pp=-pp
+c            endif
+            sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
+            sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
+            omega(1) = sqrt(p(0)+pp)
+c            omega(2) = fmass/omega(1)
+            omega(2) = sqrt(p(0)-pp)
+            ip = (3+nh)/2
+            im = (3-nh)/2
+            sfomeg(1) = sf(1)*omega(ip)
+            sfomeg(2) = sf(2)*omega(im)
+            pp3 = pp+p(3)
+            chi(1) = sqrt(pp3*rHalf/pp)
+            if ( abs(pp3).eq.rZero ) then
+               chi(2) = dcmplx(-nh )
+            else
+               chi(2) = ( nh*p(1) - ((0d0,1d0)*p(2)) )/
+     .sqrt(rTwo*pp*pp3)
+            endif
+            
+            fo(1) = sfomeg(2)*chi(im)
+            fo(2) = sfomeg(2)*chi(ip)
+            fo(3) = sfomeg(1)*chi(im)
+            fo(4) = sfomeg(1)*chi(ip)
+
+         endif
+         
+      else
+         
+c         if(zabs(p(1)).eq.0d0.and.zabs(p(2)).eq.0d0.and.zabs(p(3)).lt.0d0) then
+c            sqp0p3 = 0d0
+c         else
+            sqp0p3 = sqrt(p(0)+p(3))*nsf
+c         end if
+         chi(1) = sqp0p3
+         if ( abs(sqp0p3).eq.rZero ) then
+            chi(2) = dcmplx(-nhel )*sqrt(rTwo*p(0))
+         else
+            chi(2) = ( nh*p(1)-((0d0,1d0)*p(2)) )/sqp0p3
+         endif
+         if ( nh.eq.1 ) then
+            fo(1) = chi(1)
+            fo(2) = chi(2)
+            fo(3) = dcmplx( rZero )
+            fo(4) = dcmplx( rZero )
+         else
+            fo(1) = dcmplx( rZero )
+            fo(2) = dcmplx( rZero )
+            fo(3) = chi(2)
+            fo(4) = chi(1)
+         endif
+         
+      endif
+c
+      return
+      end
+
+   
+C     The subroutine with charge conjugation are not yet implemented
+
+      subroutine oclxxx(p,ffmass,nhel,nsf,fo)
+
+      implicit none
+      double complex fo(8),p(0:3)
+      double precision ffmass
+      integer nhel,nsf
+
+      CALL olxxxx(p,ffmass,nhel,nsf,fo)
+
+      end
+
+      subroutine iclxxx(p,ffmass,nhel,nsf,fi)
+
+      implicit none
+      double complex fi(8),p(0:3)
+      double precision ffmass
+      integer nhel,nsf
+
+      CALL ilxxxx(p,ffmass,nhel,nsf,fi)
+
+      end

=== added directory 'loop_material'
=== added directory 'loop_material/CutTools'
=== added directory 'loop_material/CutTools/DOC'
=== added file 'loop_material/CutTools/DOC/cuttools.pdf'
Binary files loop_material/CutTools/DOC/cuttools.pdf	1970-01-01 00:00:00 +0000 and loop_material/CutTools/DOC/cuttools.pdf	2012-02-27 13:42:21 +0000 differ
=== added file 'loop_material/CutTools/DOC/cuttools.ps.gz'
Binary files loop_material/CutTools/DOC/cuttools.ps.gz	1970-01-01 00:00:00 +0000 and loop_material/CutTools/DOC/cuttools.ps.gz	2012-02-27 13:42:21 +0000 differ
=== added file 'loop_material/CutTools/DOC/mpfun90.ps.gz'
Binary files loop_material/CutTools/DOC/mpfun90.ps.gz	1970-01-01 00:00:00 +0000 and loop_material/CutTools/DOC/mpfun90.ps.gz	2012-02-27 13:42:21 +0000 differ
=== added file 'loop_material/CutTools/DOC/mpfun90.tex.gz'
Binary files loop_material/CutTools/DOC/mpfun90.tex.gz	1970-01-01 00:00:00 +0000 and loop_material/CutTools/DOC/mpfun90.tex.gz	2012-02-27 13:42:21 +0000 differ
=== added file 'loop_material/CutTools/README'
--- loop_material/CutTools/README	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/README	2012-02-27 13:42:21 +0000
@@ -0,0 +1,449 @@
+------------------------------------------------------------------------
+------------------------------------------------------------------------
+ 
+              CutTools V1.6.9 (Dec 16 2011)
+
+              contact author: [e-mail]
+
+ - Authors: Giovanni Ossola, Costas G. Papadopoulos, Roberto Pittau
+
+ - When using this code, please quote the following 4 papers:
+
+   \bibitem{CutTools} 
+     G.~Ossola, C.~G.~Papadopoulos and R.~Pittau,
+     %``CutTools: a program implementing the OPP reduction method
+     %to compute one-loop amplitudes,''   
+     arXiv:0711.3596 [hep-ph].
+     %%CITATION = ARXIV:0711.3596;%%
+
+   \bibitem{Ossola:2006us}
+     G.~Ossola, C.~G.~Papadopoulos and R.~Pittau,
+     %``Reducing full one-loop amplitudes to scalar integrals at the integrand
+     %level,''
+     Nucl.\ Phys.\  B {\bf 763} (2007) 147 [arXiv:hep-ph/0609007].
+     %%CITATION = NUPHA,B763,147;%%
+
+   \bibitem{Ossola:2008xq}
+     G.~Ossola, C.~G.~Papadopoulos and R.~Pittau,
+     %``On the Rational Terms of the one-loop amplitudes,''
+     JHEP {\bf 0805} (2008) 004
+     [arXiv:0802.1876 [hep-ph]].
+     %%CITATION = JHEPA,0805,004;%%
+
+   \bibitem{Ossola:2007bb}
+     G.~Ossola, C.~G.~Papadopoulos and R.~Pittau,
+     %``Numerical Evaluation of Six-Photon Amplitudes,''
+     JHEP {\bf 0707} (2007) 085
+     [arXiv:0704.1271 [hep-ph]].
+     %%CITATION = JHEPA,0707,085;%%
+
+ - When using this code, together with the routines in 
+   avh_olo_s4.f, please quote also the 2 papers:
+
+    A. van Hameren,                                                   
+      Comput.Phys.Commun. 182 (2011) 2427-2438, arXiv:1007.4716       
+    A. van Hameren, C.G. Papadopoulos and R. Pittau,                  
+      JHEP 0909:106,2009, arXiv:0903.4665  
+
+ - When using this code, together with the routines in 
+   qcdloop, please quote also the paper:
+
+   \bibitem{Ellis:2007qk}
+     R.~K.~Ellis and G.~Zanderighi,
+     %``Scalar one-loop integrals for QCD,''
+     JHEP {\bf 0802} (2008) 002
+     [arXiv:0712.1851 [hep-ph]].
+     %%CITATION = ARXIV:0712.1851;%%
+
+ - !!Do not forget that the rational part R2 should be provided externally!!
+   (see, for example,  arXiv:1111.4965 [hep-ph]
+                       arXiv:0910.3130 [hep-ph]
+                       arXiv:0903.0356 [hep-ph])
+
+---------------------------------------------------------------------------
+---------------------------------------------------------------------------
+
+To compile and to test the program: 
+----------------------------------
+edit the first line of the makefile and set the option for the multiprecision 
+computation between PRECISION= MP (internal multiprecision routines) and
+PRECISION= QP (quadruple precision compiler, if supported by your system)
+
+Then type:
+
+  > make
+
+This will produce the directory "includects", containing the *.h files, the
+*.mod files and the library "libcts.a". The whole directory "includects" 
+and the library "libcts.a", should be present (on linked with a soft link)
+in your working directory (called "examples" in the following).
+
+To check your installation: 
+--------------------------
+go to the working directory "examples" and type
+
+  > make
+
+Then, check the implementation of the multiprecision routines by typing
+
+  > make mpcheck
+
+To compare with reference outputs contained in "examples/testdir", type
+
+  > make test  
+
+In directory "examples" will then appear the file "testoutput", containing the
+differences between the generates outputs (out*_ccc) and the reference
+ones (out*_ref) contained in "examples/testdir.
+
+How to use the program:
+-----------------------
+Your own numerator function to be "CutToolized" should be written in the file
+"numerators.f", that is also included in the directory "examples".
+When provided, a multi-precision version of it should be put in  
+"mpnumerators.f90".  
+If you want to work in another directory, the whole directory "includects" 
+and the library "libcts.a", should be present (on linked with a soft link)
+in that directory.
+
+To perform the computation, three subroutines should be called
+
+       call ctsinit       (to initialize cuttools)  
+       call ctsxcut       (to get the amplitude amp and R1)
+       call ctsstatistics (to get statistical information on the run)  
+
+In particular:
+
+                 ctsinit:
+
+!---------------------------------------------------------------------
+! To initialize CutTools call ctsinit(limitvalue,scaloop)
+!
+! INPUT:
+!
+!         real*8  limitvalue -> limit of precision below which
+!                               the PS point is considered stable
+!                               by tha A=A test (see below)      
+!
+!         integer scaloop    -> library used to compute the scalar 
+!                               1-loop functions:  
+!                               scaloop= 1 -> looptools (not implemented) 
+!                               scaloop= 2 -> avh (complex masses)   
+!                               scaloop= 3 -> qcdloop.  
+! OUTPUT: none
+!---------------------------------------------------------------------
+
+                 ctsxcut: 
+
+!---------------------------------------------------------------------
+! To compute the 1-loop amplitude
+!
+!      call ctsxcut(imode,rootsvalue,muscale,number_propagators,test,
+!     & mpfortest,rnk,pp,m2,amp,ampcc,ampr1,stable)
+!
+!
+! INPUT: integer imode              -> the running mode of CutTools according
+!                                      to the following scheme:  
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!                                                                            !
+!! imode:|  actions performed by ctsxcut:                                     !
+!!       |                                                                    !
+!!   0   | (dp_dir,dp_inv)-> dp_Atest -> stable -> (only if stable=.false.) ->!
+!!       | (mp_dir,mp_inv)-> mp_Atest -> stable                               ! 
+!!   1   | (dp_dir)       -> dp_Ntest -> stable                               !
+!!   2   | (dp_inv)       -> dp_Ntest -> stable                               !
+!!   3   | (dp_dir,dp_inv)-> dp_Atest -> stable                               !
+!!   4   | (mp_dir)       -> mp_Ntest -> stable                               ! 
+!!   5   | (mp_inv)       -> mp_Ntest -> stable                               ! 
+!!   6   | (mp_dir,mp_inv)-> mp_Atest -> stable                               !
+!!                                                                            !
+!! Legenda:                                                                   !
+!!                                                                            !
+!! dp_dir    = compute amp in double precision with normal   propagator order !
+!! dp_inv    = compute amp in double precision with reversed propagator order !
+!! mp_dir    = compute amp in multi  precision with normal   propagator order !
+!! mp_inv    = compute amp in multi  precision with reversed propagator order !
+!! dp_Atest  = perform the A=A test in double precision                       !
+!! mp_Atest  = perform the A=A test in multi  precision                       !
+!! dp_Ntest  = perform the N=N test in double precision                       !
+!! mp_Ntest  = perform the N=N test in multi  precision                       !
+!! -> stable = set stable=.true. or stable=.false.                            !
+!!             according to the outcome of the test                           !
+!!                                                                            !
+!! Tests:                                                                     !
+!!                                                                            !
+!! -The N=N test is a test on the reconstructed OPP integrand performed       !
+!!  by comparing original and reconstacted integrands at an arbirtary value   !
+!!  of the integration momentum.                                              !
+!!                                                                            ! 
+!! -The A=A test checks the 2 amplitudes obtained with dir and inv orders     !
+!!               of the propagators given in the input                        !
+!! Notes:                                                                     !
+!!                                                                            ! 
+!! a) imode= 0 is recommended, unless you really know what you are doing.     !
+!!                                                                            !
+!! b) When two determinations of amp are available, that one with more        !
+!!    accurate recounstructed numerator (coming from the N=N test) is used.   !
+!!                                                                            !
+!! c) When running in multi precision with scaloop= 3 (qcdloop), the loop     !
+!!    functions are computed in double precision only. A full multi           !
+!!    precision result can only be obtained with scaloop= 2 (OneLoop).        !
+!!                                                                            ! 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!        real*8  rootsvalue         -> the arbitrary OPP scale 
+!                                      set event by event.
+!
+!        real*8  muscale            -> the scale for the 1-loop integrals
+!                                      set event by event.
+!                                      It has dimension of an energy. 
+!
+!        integer number_propagators -> number of propagators.
+!
+!        external test              -> name of the subroutine
+!                                      computing N(q). 
+!
+!        external mpfortest         -> name of the subroutine
+!                                      computing N(q) in  
+!                                      multi-precision (if absent
+!                                      put dummy).
+!
+!        integer rnk                -> the maximum rank of N(q) (if 
+!                                      unknown put number_propagators).
+!
+!        real*8 pp(0:3,0:number_propagators-1)           
+!                                   -> momenta flowing in the internal 
+!                                      propagators.
+!
+!        complex*16 m2(0:number_propagators-1)           
+!                                   -> masses squared of the internal 
+!                                      propagators. When scaloop supports it,
+!                                      they can be complex. When scaloop does
+!                                      not support complex masses, only 
+!                                      the real part of m2 is used.  
+!   
+!               
+! OUTPUT:  complex*16 amp(0:2)      -> Amplitude (without r2):     
+!                                      amp(0) is the total finite part
+!                                      (INCLUDING R1!!!)   
+!                                      amp(1) is the coeff. of 1/eps   pole
+!                                      amp(2) is the coeff. of 1/eps^2 pole.
+!
+!          complex*16 ampcc         -> the Cut Constructible contribution
+!          complex*16 ampr1         -> the R1 contribution
+!                                      (NOTE that ampcc+ampr1 = amp(0))
+!
+!          logical stable           -> .false. if CutTools detects
+!                                      numerical instabilities.         
+!---------------------------------------------------------------------
+
+                 ctsstatistics:
+
+!---------------------------------------------------------------------
+! To know the statistics of the run call ctsstatistics(discarded)
+!
+! INPUT :  none
+!
+! OUTPUT:  logical discarded  -> .true. if there are discarded PS points
+!
+! Print out of the statistics of the run:
+!
+!          n_mp   ->  n.of points evaluated in multi-precision.
+!          n_disc ->  n.of discarded points.               
+!          n_tot  ->  n.of calls to ctsxcut               
+!---------------------------------------------------------------------
+
+An explicit example is presented in the file "examples/example0.f".
+The file "examples/ooutput0" contains the output you should obtain when
+compiling and running example0.f in the directory "examples"
+(to copile it, just go to that directory and type 
+
+  > make)  
+
+---------------------------------------------------------------------------
+---------------------------------------------------------------------------
+
+ Some history of the current version:
+
+   25 Jan 2007:
+
+   Towards a version with m-denominators
+
+   09 Apr 2008:
+
+- in dp_cutting4: 
+  if (real(cb3).le.0.d0) then changed to
+  if (real(cb3).lt.0.d0) then
+
+- in mp_cutting4: 
+  if (mpreal(cb3).le.0.d0) then changed to
+  if (mpreal(cb3).lt.0.d0) 
+
+- corrected test_rat in multiprecision (d_6 was not properly allocated 
+  when printing out its value)
+
+   24 Jun 2008:
+
+- in dp_cutting4 and mp_cutting4: corrected the roots of the quadratic 
+  equations to always get precise results
+- corrected l3 and l4 computation to always get precise results
+- The N=N test changed using now the m momenta appearing in the m denominators 
+- program restructured to interface with Formcalc/Feyncalc 
+  (call dcut,ccut etc...)
+
+   28 July 2008:
+
+- The N=N test now also contains a small component with a generic q
+
+  02 September 2008:
+
+- new algorithm implemented for calculating Rational Terms with inf=.true.  
+
+  04 September 2008:
+
+- dp_getd amd mp_getd modified to call the num function only 2 times.
+
+  24 September 2008:
+
+- Corrected error in  subroutine getdloop (cts_loopfunctions.f90)
+
+  wrong ->  k2 = den(bn4(np,3,i))%p-den(bn4(np,2,1))%p 
+  right ->  k2 = den(bn4(np,3,i))%p-den(bn4(np,2,i))%p 
+
+  October 30 2008:
+
+- interface with general avh routines 
+- corrected wrong scaless b11 function in cts_loopfunctions.f90: it was not 0.
+
+  November  04 2008:
+
+- interface with qcdloop (still to be activated at the compilation level)
+- changed the basis for the 1-point functions
+
+  December  01 2008:
+
+- changed the basis for the 2-point functions
+
+  January 23 2009:
+
+- interface with qcdloop activated 
+- loptools library disactivated (because of conflicts with qcdloop)
+
+  July 2 2009:
+ 
+- try to make of it the official version V1.1 
+
+- August 20 2009:
+
+- makefile simplified, eliminated configure, eliminated LoopTools  
+- qcdloop version 1.9 replaces qcdloop version 1.8 
+  (problem fixed in qlbox15.f)
+- makefile modified to run on DOS/Windows Platforms as well
+
+- November 19 2009:
+
+- m^2 (and not m!) as an input in acut,...,xcut
+- thrs added as an input in ctsinit
+- new bases in the 3,2 and 1 point sectors
+- N(q) is called the minimum amount of times
+- N = N test only performed with q = -momenta in the denominators
+- flag inf1 added to choose amomg 2 different ways to compute R_1 
+  in the 4 point sector in the qt2 -> limit
+- version of avh_olo with complex masses 
+
+- June 1 2010:
+
+- new test for the numerical accuracy based on a re-fitting
+  the re-constructed numerator at different values of qt2.
+  The test is now trustable.
+  The computation time is larger for simple Numerator functions because
+  the fitted ones takes more time. In real life, when N(q) is complicated
+  that should be OK. 
+- The 4-point part of R_1 is obtained with a new algorithm.
+- Most of the times, computing N(q) in double precision and
+  all the fitting procedure in mprec ADJUSTS the bad points.
+  See mpnumerators.f90 and example1.f to see how to achieve that.
+- Routines acut, bcut, ccut, dcut ... eliminated from
+  cts_cutroutines.f90
+- example2.f eliminated  
+- R_1 separated from the rest of the amplitude (HELAC-NLO convention) 
+
+- February 4 2011:
+
+- eliminated any use of dynamical memory 
+- loose N=N test added to avoid large unstable weights
+
+- February 10 2011:
+
+- computing R_1 by colling Num (not Numrec). This gives more numerical 
+  stability  
+- updated version of mprec implemented
+
+- March  1 2011:
+
+- computing R_1 and amp1 by calling Num (not Numrec) ONLY in multiprecision
+  to comprimize between speed and accuracy
+- declarations of real, complex, mp_real and mp_complex 
+  variables put in include files
+- put parameter maxden in cts_combinatorics.f90 and inserted protection
+  when number_propagators > maxden
+- inserted rootsvalue and muscalein in the input of ctsxcut  
+- introduced 2 new files cts_combinatorics.f90 and cts_combinatorics.f90, to
+  organize better the code
+
+- March 18 2011:
+
+- double-double precision routines implemented (significantly faster)
+- mptest program added in diectory examples to test the implementation
+  of the multiprecision routines
+
+- May 9 2011:
+
+- double-double precision routines my be buggy: old multiprecision 
+  re-implemented
+- maxden= 6 in this version (to be re-fixed soon)
+- option dmr= -1 implemented when 1 denominator is present in the numerator
+  function
+
+- June 5 2011:
+
+- generic choice of maxden reinserted in /src/cts/cts_combinatorics.f90
+- an automatic test of the cuttools implementation on the user's platform
+  is now available  
+- forcemp flag added to force the code to switch in multiprecision
+
+- October 19 2011:
+
+- new test for numerical precision introduced, based on reversing the order 
+  of the denominators at the input level  
+
+- the threshold for testmp has been put to 10^-11 (rnk= 4) in this version, 
+  because the loop functions are computed in double precision only
+
+- November 24 2011:
+
+- inclusion of new fortran95 version of OneLOop in this version
+
+- December 16 2011:
+  
+- mp_OneLOop implemented in this version
+
+- the threshold for testmp has been re-put to 10^-21 and rnk= 6
+
+- input/output of ctsxcut, ctsinit and ctsstatistics restructured
+
+- imode options intruduced in ctsxcut
+
+- mp/qp option introduced in the makefile
+
+- DOS/Windows Platforms optiond eliminated from the makefile
+
+- Jan 17 2012:
+
+- Bug corrected in OneLOop-2.2:In "subroutine box16"
+  sm1 = cmplx(abs(rmu)) -> cmplx(abs(rmu),kind=kindc2)
+
+---------------------------------------------------------------------------
+Last modification of this document: Jan 17 2012
+---------------------------------------------------------------------------

=== added directory 'loop_material/CutTools/examples'
=== added file 'loop_material/CutTools/examples/cts_mpc.h'
--- loop_material/CutTools/examples/cts_mpc.h	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/examples/cts_mpc.h	2012-02-27 13:42:21 +0000
@@ -0,0 +1,1 @@
+ type(mp_complex)&

=== added file 'loop_material/CutTools/examples/cts_mpr.h'
--- loop_material/CutTools/examples/cts_mpr.h	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/examples/cts_mpr.h	2012-02-27 13:42:21 +0000
@@ -0,0 +1,1 @@
+ type(mp_real)&

=== added file 'loop_material/CutTools/examples/cts_mprec.h'
--- loop_material/CutTools/examples/cts_mprec.h	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/examples/cts_mprec.h	2012-02-27 13:42:21 +0000
@@ -0,0 +1,1 @@
+  use mpmodule

=== added file 'loop_material/CutTools/examples/example0.f'
--- loop_material/CutTools/examples/example0.f	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/examples/example0.f	2012-02-27 13:42:21 +0000
@@ -0,0 +1,263 @@
+!---------------------------------------------------------------------
+!  EXAMPLE OF MAIN PROGRAM
+!
+!  Pourpose of the program:  
+!  -----------------------
+! 
+!
+!                     /                 N(q)
+!  Computing    amp=  | d^n q   ------------------------  
+!                     /         D_0 D_1 D_2 D_3 D_4 D_5
+!
+!
+!  General Structure of the program:
+!  --------------------------------
+!
+!  To perform the computation, three subroutines should be called
+!
+!       call ctsinit       (to initialize cuttools)  
+!       call ctsxcut       (to get the amplitude amp and R1)
+!       call ctsstatistics (to get statistical information on the run)  
+!             
+!  as detailed below.
+!---------------------------------------------------------------------
+!
+      program example0
+      implicit none
+!                     !----------------------------------!
+      external test   ! Name of the subroutine computing !
+!                     ! the numerator function N(q).     ! 
+!                     ! Location: numerators.f.          ! 
+!                     !----------------------------------! 
+!
+!                     !---------------------------------------!
+      external mptest ! Name of the mpsubroutine (if present) !        
+!                     ! computing the numerator function N(q).!
+!                     ! in multiprecision.                    !
+!                     ! Location: mpnumerators.f90.           ! 
+!                     ! If absent put 'external dummy'        !
+!                     !---------------------------------------!
+!
+!---------------------------------------------------------------------
+
+      common/rango/rango ! only used by the toy numerators test and mptest
+      complex*16 amp(0:2),ampcc,ampr1
+      real*8 rootsvalue,limitvalue,muscale
+      real*8 pp(0:3,0:5)           
+      complex*16 m2(0:5)            
+      integer number_propagators
+      integer rnk,rango,imode
+      integer scaloop
+      logical stable,discarded
+!
+      rootsvalue= 50.d0
+      limitvalue= 1.d-2
+      imode  = 0     
+      scaloop= 2 
+      muscale= 1.d0 
+!---------------------------------------------------------------------
+! To initialize CutTools call ctsinit(limitvalue,scaloop)
+!
+! INPUT:
+!
+!         real*8  limitvalue -> limit of precision below which
+!                               the PS point is considered stable
+!                               by tha A=A test.      
+!
+!         integer scaloop    -> library used to compute the scalar 
+!                               1-loop functions:  
+!                               scaloop= 1 -> looptools (not implemented) 
+!                               scaloop= 2 -> avh (complex masses)   
+!                               scaloop= 3 -> qcdloop.  
+! OUTPUT: none
+!---------------------------------------------------------------------
+      call ctsinit(limitvalue,scaloop)
+
+      number_propagators= 6 
+      rango= 6
+      rnk= rango 
+!
+!     momenta flowing in the 6 internal propagators: 
+!
+!     0 is the energy component
+!     1 is the x component
+!     2 is the y component
+!     3 is the z component
+!
+      pp(0,0)=    0.d0
+      pp(1,0)=    0.d0
+      pp(2,0)=    0.d0
+      pp(3,0)=    0.d0
+!
+      pp(0,1)=    25.d0
+      pp(1,1)=    0.d0
+      pp(2,1)=    0.d0
+      pp(3,1)=   -25.d0
+!
+      pp(0,2)=    12.2944131682730d0
+      pp(1,2)=    1.03940959319740d0
+      pp(2,2)=   -6.52053527152409d0
+      pp(3,2)=   -35.8551455176305d0
+!
+      pp(0,3)=   -3.11940819171487d0
+      pp(1,3)=    16.2625830020165d0
+      pp(2,3)=   -7.82868841887933d0
+      pp(3,3)=   -33.8229999456535d0
+!
+      pp(0,4)=   -4.47137880909124d0
+      pp(1,4)=    15.9241558606968d0
+      pp(2,4)=   -8.11309412871339d0
+      pp(3,4)=   -35.1006560075423d0
+!
+      pp(0,5)=   -25.d0
+      pp(1,5)=     0.d0
+      pp(2,5)=     0.d0
+      pp(3,5)=   -25.d0
+!
+!     masses of the 6 internal propagators: 
+!
+      m2(0)= 0.d0          
+      m2(1)= 0.d0            
+      m2(2)= 0.d0            
+      m2(3)= 0.d0            
+      m2(4)= 0.d0            
+      m2(5)= 0.d0            
+!---------------------------------------------------------------------
+! To compute the 1-loop amplitude
+!
+!      call ctsxcut(imode,rootsvalue,muscale,number_propagators,test,
+!     & mpfortest,rnk,pp,m2,amp,ampcc,ampr1,stable)
+!
+!
+! INPUT: integer imode              -> the running mode of CutTools according
+!                                      to the following scheme:  
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!                                                                            !
+!! imode:|  actions performed by ctsxcut:                                     !
+!!       |                                                                    !
+!!   0   | (dp_dir,dp_inv)-> dp_Atest -> stable -> (only if stable=.false.) ->!
+!!       | (mp_dir,mp_inv)-> mp_Atest -> stable                               ! 
+!!   1   | (dp_dir)       -> dp_Ntest -> stable                               !
+!!   2   | (dp_inv)       -> dp_Ntest -> stable                               !
+!!   3   | (dp_dir,dp_inv)-> dp_Atest -> stable                               !
+!!   4   | (mp_dir)       -> mp_Ntest -> stable                               ! 
+!!   5   | (mp_inv)       -> mp_Ntest -> stable                               ! 
+!!   6   | (mp_dir,mp_inv)-> mp_Atest -> stable                               !
+!!                                                                            !
+!! Legenda:                                                                   !
+!!                                                                            !
+!! dp_dir    = compute amp in double precision with normal   propagator order !
+!! dp_inv    = compute amp in double precision with reversed propagator order !
+!! mp_dir    = compute amp in multi  precision with normal   propagator order !
+!! mp_inv    = compute amp in multi  precision with reversed propagator order !
+!! dp_Atest  = perform the A=A test in double precision                       !
+!! mp_Atest  = perform the A=A test in multi  precision                       !
+!! dp_Ntest  = perform the N=N test in double precision                       !
+!! mp_Ntest  = perform the N=N test in multi  precision                       !
+!! -> stable = set stable=.true. or stable=.false.                            !
+!!             according to the outcome of the test                           !
+!!                                                                            !
+!! Tests:                                                                     !
+!!                                                                            !
+!! -The N=N test is a test on the reconstructed OPP integrand performed       !
+!!  by comparing original and reconstacted integrands at an arbirtary value   !
+!!  of the integration momentum.                                              !
+!!                                                                            ! 
+!! -The A=A test checks the 2 amplitudes obtained with dir and inv orders     !
+!!               of the propagators given in the input                        !
+!! Notes:                                                                     !
+!!                                                                            ! 
+!! a) imode= 0 is recommended, unless you really know what you are doing.     !
+!!                                                                            !
+!! b) When two determinations of amp are available, that one with more        !
+!!    accurate recounstructed numerator (coming from the N=N test) is used.   !
+!!                                                                            !
+!! c) When running in multi precision with scaloop= 3 (qcdloop), the loop     !
+!!    functions are computed in double precision only. A full multi           !
+!!    precision result can only be obtained with scaloop= 2 (OneLoop).        !
+!!                                                                            ! 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!        real*8  rootsvalue         -> the arbitrary OPP scale 
+!                                      set event by event.
+!
+!        real*8  muscale            -> the scale for the 1-loop integrals
+!                                      set event by event.
+!                                      It has dimension of an energy. 
+!
+!        integer number_propagators -> number of propagators.
+!
+!        external test              -> name of the subroutine
+!                                      computing N(q). 
+!
+!        external mpfortest         -> name of the subroutine
+!                                      computing N(q) in  
+!                                      multi-precision (if absent
+!                                      put dummy).
+!
+!        integer rnk                -> the maximum rank of N(q) (if 
+!                                      unknown put number_propagators).
+!
+!        real*8 pp(0:3,0:number_propagators-1)           
+!                                   -> momenta flowing in the internal 
+!                                      propagators.
+!
+!        complex*16 m2(0:number_propagators-1)           
+!                                   -> masses squared of the internal 
+!                                      propagators. When scaloop supports it,
+!                                      they can be complex. When scaloop does
+!                                      not support complex masses, only 
+!                                      the real part of m2 is used.  
+!   
+!               
+! OUTPUT:  complex*16 amp(0:2)      -> Amplitude (without r2):     
+!                                      amp(0) is the total finite part
+!                                      (INCLUDING R1!!!)   
+!                                      amp(1) is the coeff. of 1/eps   pole
+!                                      amp(2) is the coeff. of 1/eps^2 pole.
+!
+!          complex*16 ampcc         -> the Cut Constructible contribution
+!          complex*16 ampr1         -> the R1 contribution
+!                                      (NOTE that ampcc+ampr1 = amp(0))
+!
+!          logical stable           -> .false. if CutTools detects
+!                                      numerical instabilities.         
+!---------------------------------------------------------------------
+
+      call ctsxcut(imode,rootsvalue,muscale,number_propagators,test,
+     & mptest,rnk,pp,m2,amp,ampcc,ampr1,stable)
+      write(*,*)'               '
+      write(*,*)' Complete Amplitude (without r2):     '
+      write(*,*)'               '
+      write(*,*)'               '
+      write(*,*)' finite part           amp(0)=',amp(0)
+      write(*,*)' coeff of 1/eps   pole amp(1)=',amp(1)
+      write(*,*)' coeff of 1/eps^2 pole amp(2)=',amp(2)
+      write(*,*)'                        ampcc=',ampcc
+      write(*,*)'                           R1=',ampr1
+      write(*,*)'                       stable=',stable  
+      write(*,*)'               '
+!---------------------------------------------------------------------
+! To know the statistics of the run call ctsstatistics(discarded)
+!
+! INPUT :  none
+!
+! OUTPUT:  logical discarded  -> .true. if there are discarded PS points
+!
+! Print out of the statistics of the run:
+!
+!          n_mp   ->  n.of points evaluated in multi-precision.
+!          n_disc ->  n.of discarded points.               
+!          n_tot  ->  n.of calls to ctsxcut               
+!---------------------------------------------------------------------
+      call ctsstatistics(discarded)
+      end program example0
+
+
+ 
+
+
+
+
+

=== added file 'loop_material/CutTools/examples/example1.f'
--- loop_material/CutTools/examples/example1.f	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/examples/example1.f	2012-02-27 13:42:21 +0000
@@ -0,0 +1,289 @@
+!---------------------------------------------------------------------
+!  EXAMPLE OF MAIN PROGRAM
+!---------------------------------------------------------------------
+      program example1
+      implicit none
+!                     !----------------------------------!
+      external test   ! Name of the subroutine computing !
+!                     ! the numerator function.          ! 
+!                     ! Location numerators.f.           ! 
+!                     !----------------------------------! 
+!
+!                     !--------------------------------------!
+      external mptest ! Name of the mpsubroutine (if present)!        
+!                     ! computing the numerator function.    !
+!                     ! Location mpnumerators.f90.           ! 
+!                     ! If absent put 'external dummy'       !
+!                     !--------------------------------------!
+!---------------------------------------------------------------------
+      integer maxden
+      parameter (maxden= 10)
+      complex*16 m2(0:maxden-1)
+      real*8 pp(0:3,0:maxden-1)           
+!---------------------------------------------------------------------
+! Auxiliary variables:
+!---------------------------------------------------------------------
+      common/rango/rango  ! only used by the toy numerator 
+      real*8 xm(1:maxden-2),p(4,maxden-2),k(0:3,maxden)           
+      real*8 rootsvalue,limitvalue,roots,muscale
+      complex*16 amp(0:2),ampcc,ampr1
+      integer number_propagators
+      integer rnk,i,j,l,iter,rango
+      integer scaloop,imode
+      integer npoints
+      logical stable,discarded
+!---------------------------------------------------------------------
+! Read number of MC points, the number of propagators and the rank:
+!--------------------------------------------------------------------- 
+      print*,'enter npoints,number_propagators,rank,scaloop,muscale'
+      print*,'    '
+      print*,'scaloop= 1 -> looptools 1-loop '
+      print*,'scaloop= 2 -> avh 1-loop (massive with complex masses)'
+      print*,'scaloop= 3 -> qcdloop   1-loop (Ellis and Zanderighi)'
+      print*,'muscale (dimension of energy) is the scale' 
+      print*,'for the 1-loop integrals' 
+      print*,'    '
+      read*,npoints,number_propagators,rnk,scaloop,muscale
+      rango= rnk                ! only used by the toy numerators 
+                                ! located in numerators.f and  
+                                ! mpnumerators.f90
+!                               !
+      if (number_propagators.gt.maxden) then
+         stop 'increase maxden in example1.f90'
+      endif
+      roots     = 50.d0         ! value of sqrt(s)
+      rootsvalue= roots
+      limitvalue= 1.d-2 
+      imode     = 0
+!---------------------------------------------------------------------
+! Input momenta and masses in each denominator:
+!---------------------------------------------------------------------
+      k(0,1)= roots/2.d0
+      k(1,1)= 0.d0
+      k(2,1)= 0.d0
+      k(3,1)= roots/2.d0
+      k(0,2)= roots/2.d0
+      k(1,2)= 0.d0
+      k(2,2)= 0.d0
+      k(3,2)=-roots/2.d0
+!
+      do i= 1,number_propagators-2
+         xm(i) = 1.d0
+      enddo
+!---------------------------------------------------------------------
+! To initialize CutTools call ctsinit(limitvalue,scaloop)
+!
+! INPUT:
+!
+!         real*8  limitvalue -> limit of precision below which
+!                               the PS point is considered stable
+!                               by tha A=A test.      
+!
+!         integer scaloop    -> library used to compute the scalar 
+!                               1-loop functions:  
+!                               scaloop= 1 -> looptools (not implemented) 
+!                               scaloop= 2 -> avh (complex masses)   
+!                               scaloop= 3 -> qcdloop.  
+! OUTPUT: none
+!---------------------------------------------------------------------
+      call ctsinit(limitvalue,scaloop)
+      do iter= 1,npoints   ! do-loop over events 
+        if     (number_propagators.ge.4) then
+           call rambo(0,number_propagators-2,roots,xm,p)
+        elseif (number_propagators.eq.3) then
+           p(4,1)= k(0,1)+k(0,2) 
+           p(1,1)= k(1,1)+k(1,2)
+           p(2,1)= k(2,1)+k(2,2)
+           p(3,1)= k(3,1)+k(3,2)
+        elseif (number_propagators.eq.2) then
+        elseif (number_propagators.eq.1) then
+        else
+          print*,'number_propagators=',number_propagators,'not allowed!'
+          stop
+        endif 
+          do j= 0,3
+           if (j.ne.0) then
+              do l= 1,number_propagators-2
+                 k(j,l+2)= p(j,l) 
+              enddo
+           else
+            do l= 1,number_propagators-2
+               k(j,l+2)= p(4,l) 
+            enddo
+         endif
+        enddo
+!
+        do l= 0,3 
+          pp(l,0)= 0.d0
+        enddo
+        do i= 1,number_propagators-1
+           do l= 0,3 
+             pp(l,i)=  k(l,2)
+           enddo
+           do j= 3,i+1 
+              do l= 0,3 
+                pp(l,i)= pp(l,i)-k(l,j)
+              enddo 
+           enddo
+        enddo
+!
+!        do i= 0,number_propagators-1
+!           m2(i)= 0.d0            
+!        enddo
+!
+        m2(0)= 1.d0 
+        m2(1)= 2.d0 
+        m2(2)= 3.d0 
+        m2(3)= 4.d0 
+        m2(4)= 0.d0 
+        m2(5)= 0.d0 
+!
+!comment
+!        do i= 1,number_propagators-1
+!           do l= 0,3
+!             print*,'l,i,pp(l,i)=',l,i,pp(l,i)
+!           enddo
+!        enddo
+!comment
+!---------------------------------------------------------------------
+! To compute the 1-loop amplitude
+!
+!      call ctsxcut(imode,rootsvalue,muscale,number_propagators,test,
+!     & mpfortest,rnk,pp,m2,amp,ampcc,ampr1,stable)
+!
+!
+! INPUT: integer imode              -> the running mode of CutTools according
+!                                      to the following scheme:  
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!                                                                            !
+!! imode:|  actions performed by ctsxcut:                                     !
+!!       |                                                                    !
+!!   0   | (dp_dir,dp_inv)-> dp_Atest -> stable -> (only if stable=.false.) ->!
+!!       | (mp_dir,mp_inv)-> mp_Atest -> stable                               ! 
+!!   1   | (dp_dir)       -> dp_Ntest -> stable                               !
+!!   2   | (dp_inv)       -> dp_Ntest -> stable                               !
+!!   3   | (dp_dir,dp_inv)-> dp_Atest -> stable                               !
+!!   4   | (mp_dir)       -> mp_Ntest -> stable                               ! 
+!!   5   | (mp_inv)       -> mp_Ntest -> stable                               ! 
+!!   6   | (mp_dir,mp_inv)-> mp_Atest -> stable                               !
+!!                                                                            !
+!! Legenda:                                                                   !
+!!                                                                            !
+!! dp_dir    = compute amp in double precision with normal   propagator order !
+!! dp_inv    = compute amp in double precision with reversed propagator order !
+!! mp_dir    = compute amp in multi  precision with normal   propagator order !
+!! mp_inv    = compute amp in multi  precision with reversed propagator order !
+!! dp_Atest  = perform the A=A test in double precision                       !
+!! mp_Atest  = perform the A=A test in multi  precision                       !
+!! dp_Ntest  = perform the N=N test in double precision                       !
+!! mp_Ntest  = perform the N=N test in multi  precision                       !
+!! -> stable = set stable=.true. or stable=.false.                            !
+!!             according to the outcome of the test                           !
+!!                                                                            !
+!! Tests:                                                                     !
+!!                                                                            !
+!! -The N=N test is a test on the reconstructed OPP integrand performed       !
+!!  by comparing original and reconstacted integrands at an arbirtary value   !
+!!  of the integration momentum.                                              !
+!!                                                                            ! 
+!! -The A=A test checks the 2 amplitudes obtained with dir and inv orders     !
+!!               of the propagators given in the input                        !
+!! Notes:                                                                     !
+!!                                                                            ! 
+!! a) imode= 0 is recommended, unless you really know what you are doing.     !
+!!                                                                            !
+!! b) When two determinations of amp are available, that one with more        !
+!!    accurate recounstructed numerator (coming from the N=N test) is used.   !
+!!                                                                            !
+!! c) When running in multi precision with scaloop= 3 (qcdloop), the loop     !
+!!    functions are computed in double precision only. A full multi           !
+!!    precision result can only be obtained with scaloop= 2 (OneLoop).        !
+!!                                                                            ! 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!        real*8  rootsvalue         -> the arbitrary OPP scale 
+!                                      set event by event.
+!
+!        real*8  muscale            -> the scale for the 1-loop integrals
+!                                      set event by event.
+!                                      It has dimension of an energy. 
+!
+!        integer number_propagators -> number of propagators.
+!
+!        external test              -> name of the subroutine
+!                                      computing N(q). 
+!
+!        external mpfortest         -> name of the subroutine
+!                                      computing N(q) in  
+!                                      multi-precision (if absent
+!                                      put dummy).
+!
+!        integer rnk                -> the maximum rank of N(q) (if 
+!                                      unknown put number_propagators).
+!
+!        real*8 pp(0:3,0:number_propagators-1)           
+!                                   -> momenta flowing in the internal 
+!                                      propagators.
+!
+!        complex*16 m2(0:number_propagators-1)           
+!                                   -> masses squared of the internal 
+!                                      propagators. When scaloop supports it,
+!                                      they can be complex. When scaloop does
+!                                      not support complex masses, only 
+!                                      the real part of m2 is used.  
+!   
+!               
+! OUTPUT:  complex*16 amp(0:2)      -> Amplitude (without r2):     
+!                                      amp(0) is the total finite part
+!                                      (INCLUDING R1!!!)   
+!                                      amp(1) is the coeff. of 1/eps   pole
+!                                      amp(2) is the coeff. of 1/eps^2 pole.
+!
+!          complex*16 ampcc         -> the Cut Constructible contribution
+!          complex*16 ampr1         -> the R1 contribution
+!                                      (NOTE that ampcc+ampr1 = amp(0))
+!
+!          logical stable           -> .false. if CutTools detects
+!                                      numerical instabilities.         
+!---------------------------------------------------------------------
+        call ctsxcut(imode,rootsvalue,muscale,number_propagators,test,
+     &   mptest,rnk,pp,m2,amp,ampcc,ampr1,stable)
+        write(*,*)'               '
+        write(*,*)'  iter= ',iter
+        write(*,*)'               '
+        write(*,*)'               '
+        write(*,*)' Complete Amplitude (without r2):     '
+        write(*,*)'               '
+        write(*,*)'               '
+        write(*,*)' finite part           amp(0)=',amp(0)
+        write(*,*)' coeff of 1/eps   pole amp(1)=',amp(1)
+        write(*,*)' coeff of 1/eps^2 pole amp(2)=',amp(2)
+        write(*,*)'                        ampcc=',ampcc
+        write(*,*)'                           R1=',ampr1
+        write(*,*)'                       stable=',stable  
+        write(*,*)'               '
+      enddo
+!---------------------------------------------------------------------
+! To know the statistics of the run call ctsstatistics(discarded)
+!
+! INPUT :  none
+!
+! OUTPUT:  logical discarded  -> .true. if there are discarded PS points
+!
+! Print out of the statistics of the run:
+!
+!          n_mp   ->  n.of points evaluated in multi-precision.
+!          n_disc ->  n.of discarded points.               
+!          n_tot  ->  n.of calls to ctsxcut               
+!---------------------------------------------------------------------
+      call ctsstatistics(discarded)
+      end program example1
+
+
+ 
+
+
+
+
+

=== added file 'loop_material/CutTools/examples/makefile'
--- loop_material/CutTools/examples/makefile	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/examples/makefile	2012-02-27 13:42:21 +0000
@@ -0,0 +1,38 @@
+all: testmp$(EXE) example0$(EXE) example1$(EXE) 
+
+include ../makefile
+
+NUMS0= numerators.o mpnumerators.o 
+
+NUMS = numerators.o mpnumerators.o rambo.o random.o 
+
+mpcheck:
+	./testmp
+
+test:
+	cd testdir; cp input* ../; cp do ../; cp cdiff ../; cd ../; sh do; \
+ cp out* testdir; rm out*; rm input* ; sh cdiff > testoutput; rm cdiff; rm do
+
+testmp$(EXE): $(NUMS) libcts.a testmp.o
+	$(FC) $(FFLAGS) -Iincludects -o testmp$(EXE) \
+ testmp.f $(NUMS) libcts.a
+
+example0$(EXE): $(NUMS0) libcts.a example0.o
+	$(FC) $(FFLAGS) -Iincludects -o example0$(EXE) \
+ example0.f $(NUMS0) libcts.a
+
+example1$(EXE): $(NUMS) libcts.a example1.o 
+	$(FC) $(FFLAGS) -Iincludects -o example1$(EXE) \
+ example1.f $(NUMS) libcts.a
+
+libcts.a:
+	cd .. ; $(MAKE)
+
+.SUFFIXES: .f90
+
+.f90.o:
+	$(FC) $(FFLAGS) -Iincludects -c $<
+
+clean:
+	$(RM) testmp$(EXE) example0$(EXE) example1$(EXE)  *.o testoutput
+

=== added file 'loop_material/CutTools/examples/mpnumerators.f90'
--- loop_material/CutTools/examples/mpnumerators.f90	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/examples/mpnumerators.f90	2012-02-27 13:42:21 +0000
@@ -0,0 +1,86 @@
+!---------------------------------------------------------------------
+!
+!
+!     In this file (mpnumerators.f90) 
+!     you should place all multi-precision numerator functions
+!
+!
+!---------------------------------------------------------------------
+!
+      subroutine dummy(q,amp)
+!
+!-----------------------------------------------
+!                               
+!     dummy subroutine, to be calld when idig= 0           
+!                               
+!----------------------------------------------
+
+      include 'cts_mprec.h'
+      implicit none
+      include 'cts_mpc.h'
+       , intent(in), dimension(0:3) :: q
+      include 'cts_mpc.h'
+       , intent(out) :: amp
+      amp = 0
+      end subroutine dummy
+!      
+      subroutine mptest(q,amp)
+!
+!-------------------------------
+!                               
+!     "Test" numerator           
+!                               
+!-------------------------------
+!
+      include 'cts_mprec.h'
+      implicit none
+      include 'cts_mpc.h'
+       , intent(in), dimension(0:3) :: q
+      include 'cts_mpc.h'
+       , intent(out) :: amp
+      complex(kind(1.d0)) :: dpamp,dpq(0:3)
+      integer ::  j
+      integer rango
+      common/rango/rango
+!
+! comment: the double precision version of the function is called here 
+!
+!      amp  = 2.d0+(4.d0+q(0)*2.22333d0-q(1)*30.666d0+q(2)*1.55d0 &
+!     &       +q(3)*3.6541122d0)**rango
+
+!      amp= (2.d0+(4.d0+q(0)*2.22333d0-q(1)*30.666d0+q(2)*1.55d0 &
+!                  +q(3)*3.6541122d0)**2                         &
+!           )*(q(0)**2-q(1)**2-q(2)**2-q(3)**2+7.d0)*            &
+!           (q(0)*1.12333d0-q(1)*10.666d0+q(2)*5.55d0            &  
+!                  -q(3)*2.6541122d0)
+
+      do j= 0,3; dpq(j)= q(j); enddo
+      call test(dpq,dpamp)
+      amp= dpamp
+! comment
+      end subroutine mptest
+
+      subroutine mpfortest(q,amp)
+!
+!-------------------------------
+!                               
+!     Numerator for testing the 
+!     implementation of the 
+!     multiprecision routines           
+!                               
+!-------------------------------
+!
+      include 'cts_mprec.h'
+      implicit none
+      include 'cts_mpc.h'
+       , intent(in), dimension(0:3) :: q
+      include 'cts_mpc.h'
+       , intent(out) :: amp
+      complex(kind(1.d0)) :: dpamp,dpq(0:3)
+      integer ::  j
+      integer rango
+      common/rango/rango
+      amp  = 2.d0+(4.d0+q(0)*2.22333d0-q(1)*30.666d0+q(2)*1.55d0 &
+     &       +q(3)*3.6541122d0)**rango
+
+      end subroutine mpfortest

=== added file 'loop_material/CutTools/examples/numerators.f'
--- loop_material/CutTools/examples/numerators.f	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/examples/numerators.f	2012-02-27 13:42:21 +0000
@@ -0,0 +1,31 @@
+!---------------------------------------------------------------------
+!
+!
+!     In this file (numerators.f) 
+!     you should place all numerator functions
+!
+!
+!---------------------------------------------------------------------
+!
+      subroutine test(q,amp)
+!
+!-------------------------------
+!                               
+!     "Test" numerator           
+!                               
+!-------------------------------
+!
+      implicit none 
+      complex*16 q(0:3)
+      complex*16 amp
+      integer rango
+      common/rango/rango
+      amp  = 2.d0+(4.d0+q(0)*2.22333d0-q(1)*30.666d0+q(2)*1.55d0 
+     &       +q(3)*3.6541122d0)**rango
+
+!      amp= (2.d0+(4.d0+q(0)*2.22333d0-q(1)*30.666d0+q(2)*1.55d0 
+!     &            +q(3)*3.6541122d0)**2
+!     &     )*(q(0)**2-q(1)**2-q(2)**2-q(3)**2+7.d0)*
+!     &     (q(0)*1.12333d0-q(1)*10.666d0+q(2)*5.55d0 
+!     &            -q(3)*2.6541122d0)
+      end subroutine test

=== added file 'loop_material/CutTools/examples/rambo.f'
--- loop_material/CutTools/examples/rambo.f	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/examples/rambo.f	2012-02-27 13:42:21 +0000
@@ -0,0 +1,185 @@
+      SUBROUTINE RAMBO(LFLAG,N,ET,XM,P)
+c------------------------------------------------------
+c
+c                       RAMBO
+c
+c    RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)
+c
+c    A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR
+c    AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING
+c    THIS IS VERSION 1.0 -  WRITTEN BY R. KLEISS
+c    (MODIFIED BY R. PITTAU)
+c
+c                INPUT                 OUTPUT
+c
+c    LFLAG= 0:   N, ET, XM             P, (DJ)
+c    LFLAG= 1:   N, ET, XM, P          (DJ)
+c
+c    N  = NUMBER OF PARTICLES (>1, IN THIS VERSION <101)
+c    ET = TOTAL CENTRE-OF-MASS ENERGY
+c    XM = PARTICLE MASSES ( DIM=100 )
+c    P  = PARTICLE MOMENTA ( DIM=(4,100) )
+c    DJ = 1/(WEIGHT OF THE EVENT)
+c
+c------------------------------------------------------
+      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
+      DIMENSION XM(100),P(4,100),Q(4,100),Z(100),R(4),
+     .   B(3),P2(100),XM2(100),E(100),V(100),IWARN(5)
+      SAVE ACC,ITMAX,IBEGIN,IWARN,Z,TWOPI,PO2LOG
+      DATA ACC/1.D-14/,ITMAX/10/,IBEGIN/0/,IWARN/5*0/
+C
+C INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
+      IF(IBEGIN.NE.0) GOTO 103
+      IBEGIN=1
+      TWOPI=8.*DATAN(1.D0)
+      PO2LOG=LOG(TWOPI/4.)
+      Z(2)=PO2LOG
+      DO 101 K=3,100
+  101 Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2))
+      DO 102 K=3,100
+  102 Z(K)=(Z(K)-LOG(DFLOAT(K-1)))
+C
+C CHECK ON THE NUMBER OF PARTICLES
+  103 IF(N.GT.1.AND.N.LT.101) GOTO 104
+      PRINT 1001,N
+      STOP
+C
+C CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
+  104 XMT=0.
+      NM=0
+      DO 105 I=1,N
+      IF(XM(I).NE.0.D0) NM=NM+1
+  105 XMT=XMT+ABS(XM(I))
+      IF(XMT.LE.ET) GOTO 201
+      PRINT 1002,XMT,ET
+      STOP
+
+  201 CONTINUE 
+      if (lflag.eq.1) then
+        w0= exp((2.*N-4.)*LOG(ET)+Z(N))
+        do j= 1,N
+          v(j)= sqrt(p(1,j)**2+p(2,j)**2+p(3,j)**2)
+        enddo
+
+        a1= 0.d0
+        a3= 0.d0
+        a2= 1.d0
+        do j= 1,N
+          a1= a1+v(j)/ET
+          a2= a2*v(j)/p(4,j)
+          a3= a3+v(j)*v(j)/p(4,j)/ET
+        enddo
+        wm= a1**(2*N-3)*a2/a3
+        dj= 1.d0/w0/wm
+        return
+      endif
+C
+C THE PARAMETER VALUES ARE NOW ACCEPTED
+C
+C GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
+
+      DO 202 I=1,N
+      call rans(RAN1)
+      call rans(RAN2)
+      call rans(RAN3)
+      call rans(RAN4)
+      C=2.*RAN1-1.
+      S=SQRT(1.-C*C)
+      F=TWOPI*RAN2
+      Q(4,I)=-LOG(RAN3*RAN4)
+      Q(3,I)=Q(4,I)*C
+      Q(2,I)=Q(4,I)*S*COS(F)
+  202 Q(1,I)=Q(4,I)*S*SIN(F)
+C
+C CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
+      DO 203 I=1,4
+  203 R(I)=0.
+      DO 204 I=1,N
+      DO 204 K=1,4
+  204 R(K)=R(K)+Q(K,I)
+      RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
+      DO 205 K=1,3
+  205 B(K)=-R(K)/RMAS
+      G=R(4)/RMAS
+      A=1./(1.+G)
+      X=ET/RMAS
+C
+C TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
+      DO 207 I=1,N
+      BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
+      DO 206 K=1,3
+  206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
+  207 P(4,I)=X*(G*Q(4,I)+BQ)
+C
+C CALCULATE WEIGHT AND POSSIBLE WARNINGS
+      WT=PO2LOG
+      IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
+      IF(WT.GE.-180.D0) GOTO 208
+      IF(IWARN(1).LE.5) PRINT 1004,WT
+      IWARN(1)=IWARN(1)+1
+  208 IF(WT.LE. 174.D0) GOTO 209
+      IF(IWARN(2).LE.5) PRINT 1005,WT
+      IWARN(2)=IWARN(2)+1
+C
+C RETURN FOR WEIGHTED MASSLESS MOMENTA
+  209 IF(NM.NE.0) GOTO 210
+      WT=EXP(WT)
+      DJ= 1.d0/WT
+      RETURN
+C
+C MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
+  210 XMAX=SQRT(1.-(XMT/ET)**2)
+      DO 301 I=1,N
+      XM2(I)=XM(I)**2
+  301 P2(I)=P(4,I)**2
+      ITER=0
+      X=XMAX
+      ACCU=ET*ACC
+  302 F0=-ET
+      G0=0.
+      X2=X*X
+      DO 303 I=1,N
+      E(I)=SQRT(XM2(I)+X2*P2(I))
+      F0=F0+E(I)
+  303 G0=G0+P2(I)/E(I)
+      IF(ABS(F0).LE.ACCU) GOTO 305
+      ITER=ITER+1
+      IF(ITER.LE.ITMAX) GOTO 304
+      PRINT 1006,ITMAX
+      GOTO 305
+  304 X=X-F0/(X*G0)
+      GOTO 302
+  305 DO 307 I=1,N
+      V(I)=X*P(4,I)
+      DO 306 K=1,3
+  306 P(K,I)=X*P(K,I)
+  307 P(4,I)=E(I)
+C
+C CALCULATE THE MASS-EFFECT WEIGHT FACTOR
+      WT2=1.
+      WT3=0.
+      DO 308 I=1,N
+      WT2=WT2*V(I)/E(I)
+  308 WT3=WT3+V(I)**2/E(I)
+      WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
+C
+C RETURN FOR  WEIGHTED MASSIVE MOMENTA
+      WT=WT+WTM
+      IF(WT.GE.-180.D0) GOTO 309
+      IF(IWARN(3).LE.5) PRINT 1004,WT
+      IWARN(3)=IWARN(3)+1
+  309 IF(WT.LE. 174.D0) GOTO 310
+      IF(IWARN(4).LE.5) PRINT 1005,WT
+      IWARN(4)=IWARN(4)+1
+  310 WT=EXP(WT)
+      DJ= 1.d0/WT
+      RETURN
+C
+ 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
+ 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',
+     . ' SMALLER THAN TOTAL ENERGY =',D15.6)
+ 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
+ 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY  OVERFLOW')
+ 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE',
+     . ' DESIRED ACCURACY =',D15.6)
+      END

=== added file 'loop_material/CutTools/examples/random.f'
--- loop_material/CutTools/examples/random.f	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/examples/random.f	2012-02-27 13:42:21 +0000
@@ -0,0 +1,52 @@
+      subroutine rans(r)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Random number generator                                    c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      real*8 r,rangen
+      r = rangen(1)
+      end
+
+*-- Author :    F. James, modified by Mike Seymour
+C-----------------------------------------------------------------------
+      FUNCTION RANGEN(I)
+C-----------------------------------------------------------------------
+C     MAIN RANDOM NUMBER GENERATOR
+C     USES METHOD OF l'Ecuyer, (VIA F.JAMES, COMP PHYS COMM 60(1990)329)
+C-----------------------------------------------------------------------
+      IMPLICIT NONE
+      DOUBLE PRECISION RANGEN,RANSET,RANGET
+      INTEGER I,ISEED(2),K,IZ,JSEED(2)
+      SAVE ISEED
+      DATA ISEED/12345,67890/
+      K=ISEED(1)/53668
+      ISEED(1)=40014*(ISEED(1)-K*53668)-K*12211
+      IF (ISEED(1).LT.0) ISEED(1)=ISEED(1)[phone]
+      K=ISEED(2)/52774
+      ISEED(2)=40692*(ISEED(2)-K*52774)-K*3791
+      IF (ISEED(2).LT.0) ISEED(2)=ISEED(2)[phone]
+      IZ=ISEED(1)-ISEED(2)
+      IF (IZ.LT.1) IZ=IZ+2147483562
+      RANGEN=DBLE(IZ)*4.656613001013252D-10
+C--->                (4.656613001013252D-10 = 1.D0/2147483589)
+      RETURN
+C-----------------------------------------------------------------------
+      ENTRY RANSET(JSEED)
+C-----------------------------------------------------------------------
+      IF (JSEED(1).EQ.0.OR.JSEED(2).EQ.0) then
+         write(6,*) 'Jseeds=0, wrong settings for RANSET'
+         stop
+      endif
+      ISEED(1)=JSEED(1)
+      ISEED(2)=JSEED(2)
+ 999  RETURN
+C-----------------------------------------------------------------------
+      ENTRY RANGET(JSEED)
+C-----------------------------------------------------------------------
+      JSEED(1)=ISEED(1)
+      JSEED(2)=ISEED(2)
+      RETURN
+      END
+

=== added file 'loop_material/CutTools/examples/testmp.f'
--- loop_material/CutTools/examples/testmp.f	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/examples/testmp.f	2012-02-27 13:42:21 +0000
@@ -0,0 +1,288 @@
+!---------------------------------------------------------------------
+! PROGRAM FOR TESTING THE IMPLEMENTATION OF THE MP ROUTINES
+!---------------------------------------------------------------------
+      program testmp
+      implicit none
+!                     !----------------------------------!
+      external test   ! Name of the subroutine computing !
+!                     ! the numerator function.          ! 
+!                     ! Location numerators.f.           ! 
+!                     !----------------------------------! 
+!
+!                         !--------------------------------------!
+      external mpfortest  ! Name of the mpsubroutine (if present)!        
+!                         ! computing the numerator function.    !
+!                         ! Location mpnumerators.f90.           ! 
+!                         ! If absent put 'external dummy'       !
+!                         !--------------------------------------!
+!---------------------------------------------------------------------
+      integer maxden,imode
+      parameter (maxden= 10)
+      complex*16 m2(0:maxden-1)
+      real*8 pp(0:3,0:maxden-1)           
+!---------------------------------------------------------------------
+! Auxiliary variables:
+!---------------------------------------------------------------------
+      common/rango/rango  ! only used by the toy numerator 
+      real*8 xm(1:maxden-2),p(4,maxden-2),k(0:3,maxden)           
+      real*8 rootsvalue,limitvalue,roots,muscale,thrs
+      complex*16 amp(0:2),ampcc,ampr1
+      integer number_propagators
+      integer rnk,i,j,l,iter,idig,rango
+      integer scaloop
+      integer npoints
+      logical stable,discarded
+!
+      npoints           = 10
+      imode             = 0
+      number_propagators= 6    
+      rnk               = 6    
+      scaloop           = 2   
+      muscale           = 1.d0 
+      limitvalue        = 1.d-21
+!
+      rango= rnk        ! only used by the toy numerators 
+                        ! located in numerators.f and  
+                        ! mpnumerators.f90
+!                       !
+      if (number_propagators.gt.maxden) then
+         stop 'increase maxden in testmp.f'
+      endif
+!
+      roots     = 50.d0  ! value of sqrt(s)
+      rootsvalue= roots  ! used as an internal scale of the OPP algorithm
+!---------------------------------------------------------------------
+! Input momenta and masses in each denominator:
+!---------------------------------------------------------------------
+      k(0,1)= roots/2.d0
+      k(1,1)= 0.d0
+      k(2,1)= 0.d0
+      k(3,1)= roots/2.d0
+      k(0,2)= roots/2.d0
+      k(1,2)= 0.d0
+      k(2,2)= 0.d0
+      k(3,2)=-roots/2.d0
+      do i= 1,number_propagators-2
+         xm(i) = 0.d0
+      enddo
+!---------------------------------------------------------------------
+! To initialize CutTools call ctsinit(limitvalue,scaloop)
+!
+! INPUT:
+!
+!         real*8  limitvalue -> limit of precision below which
+!                               the PS point is considered stable
+!                               by tha A=A test.      
+!
+!         integer scaloop    -> library used to compute the scalar 
+!                               1-loop functions:  
+!                               scaloop= 1 -> looptools (not implemented) 
+!                               scaloop= 2 -> avh (complex masses)   
+!                               scaloop= 3 -> qcdloop.  
+! OUTPUT: none
+!---------------------------------------------------------------------
+      call ctsinit(limitvalue,scaloop)
+      print*,'  '
+      do iter= 1,npoints   ! do-loop over events 
+        print*,'I am running to check event n.',iter,' over 10'
+        if     (number_propagators.ge.4) then
+           call rambo(0,number_propagators-2,roots,xm,p)
+        elseif (number_propagators.eq.3) then
+           p(4,1)= k(0,1)+k(0,2) 
+           p(1,1)= k(1,1)+k(1,2)
+           p(2,1)= k(2,1)+k(2,2)
+           p(3,1)= k(3,1)+k(3,2)
+        elseif (number_propagators.eq.2) then
+        elseif (number_propagators.eq.1) then
+        else
+          print*,'number_propagators=',number_propagators,'not allowed!'
+          stop
+        endif 
+          do j= 0,3
+           if (j.ne.0) then
+              do l= 1,number_propagators-2
+                 k(j,l+2)= p(j,l) 
+              enddo
+           else
+            do l= 1,number_propagators-2
+               k(j,l+2)= p(4,l) 
+            enddo
+         endif
+        enddo
+!
+        do l= 0,3 
+          pp(l,0)= 0.d0
+        enddo
+        do i= 1,number_propagators-1
+           do l= 0,3 
+             pp(l,i)=  k(l,2)
+           enddo
+           do j= 3,i+1 
+              do l= 0,3 
+                pp(l,i)= pp(l,i)-k(l,j)
+              enddo 
+           enddo
+        enddo
+
+        do i= 0,number_propagators-1
+           m2(i)= 0.d0            
+        enddo
+!
+        m2(0)= 1.d0 
+        m2(1)= 2.d0
+        m2(2)= 3.d0
+        m2(3)= 4.d0
+        m2(4)= 0.d0 
+        m2(5)= 0.d0 
+!comment
+!
+!        do i= 1,number_propagators-1
+!           do l= 0,3
+!             print*,'l,i,pp(l,i)=',l,i,pp(l,i)
+!           enddo
+!        enddo
+!
+!comment
+
+!---------------------------------------------------------------------
+! To compute the 1-loop amplitude
+!
+!      call ctsxcut(imode,rootsvalue,muscale,number_propagators,test,
+!     & mpfortest,rnk,pp,m2,amp,ampcc,ampr1,stable)
+!
+!
+! INPUT: integer imode              -> the running mode of CutTools according
+!                                      to the following scheme:  
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!                                                                            !
+!! imode:|  actions performed by ctsxcut:                                     !
+!!       |                                                                    !
+!!   0   | (dp_dir,dp_inv)-> dp_Atest -> stable -> (only if stable=.false.) ->!
+!!       | (mp_dir,mp_inv)-> mp_Atest -> stable                               ! 
+!!   1   | (dp_dir)       -> dp_Ntest -> stable                               !
+!!   2   | (dp_inv)       -> dp_Ntest -> stable                               !
+!!   3   | (dp_dir,dp_inv)-> dp_Atest -> stable                               !
+!!   4   | (mp_dir)       -> mp_Ntest -> stable                               ! 
+!!   5   | (mp_inv)       -> mp_Ntest -> stable                               ! 
+!!   6   | (mp_dir,mp_inv)-> mp_Atest -> stable                               !
+!!                                                                            !
+!! Legenda:                                                                   !
+!!                                                                            !
+!! dp_dir    = compute amp in double precision with normal   propagator order !
+!! dp_inv    = compute amp in double precision with reversed propagator order !
+!! mp_dir    = compute amp in multi  precision with normal   propagator order !
+!! mp_inv    = compute amp in multi  precision with reversed propagator order !
+!! dp_Atest  = perform the A=A test in double precision                       !
+!! mp_Atest  = perform the A=A test in multi  precision                       !
+!! dp_Ntest  = perform the N=N test in double precision                       !
+!! mp_Ntest  = perform the N=N test in multi  precision                       !
+!! -> stable = set stable=.true. or stable=.false.                            !
+!!             according to the outcome of the test                           !
+!!                                                                            !
+!! Tests:                                                                     !
+!!                                                                            !
+!! -The N=N test is a test on the reconstructed OPP integrand performed       !
+!!  by comparing original and reconstacted integrands at an arbirtary value   !
+!!  of the integration momentum.                                              !
+!!                                                                            ! 
+!! -The A=A test checks the 2 amplitudes obtained with dir and inv orders     !
+!!               of the propagators given in the input                        !
+!! Notes:                                                                     !
+!!                                                                            ! 
+!! a) imode= 0 is recommended, unless you really know what you are doing.     !
+!!                                                                            !
+!! b) When two determinations of amp are available, that one with more        !
+!!    accurate recounstructed numerator (coming from the N=N test) is used.   !
+!!                                                                            !
+!! c) When running in multi precision with scaloop= 3 (qcdloop), the loop     !
+!!    functions are computed in double precision only. A full multi           !
+!!    precision result can only be obtained with scaloop= 2 (OneLoop).        !
+!!                                                                            ! 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!        real*8  rootsvalue         -> the arbitrary OPP scale 
+!                                      set event by event.
+!
+!        real*8  muscale            -> the scale for the 1-loop integrals
+!                                      set event by event.
+!                                      It has dimension of an energy. 
+!
+!        integer number_propagators -> number of propagators.
+!
+!        external test              -> name of the subroutine
+!                                      computing N(q). 
+!
+!        external mpfortest         -> name of the subroutine
+!                                      computing N(q) in  
+!                                      multi-precision (if absent
+!                                      put dummy).
+!
+!        integer rnk                -> the maximum rank of N(q) (if 
+!                                      unknown put number_propagators).
+!
+!        real*8 pp(0:3,0:number_propagators-1)           
+!                                   -> momenta flowing in the internal 
+!                                      propagators.
+!
+!        complex*16 m2(0:number_propagators-1)           
+!                                   -> masses squared of the internal 
+!                                      propagators. When scaloop supports it,
+!                                      they can be complex. When scaloop does
+!                                      not support complex masses, only 
+!                                      the real part of m2 is used.  
+!   
+!               
+! OUTPUT:  complex*16 amp(0:2)      -> Amplitude (without r2):     
+!                                      amp(0) is the total finite part
+!                                      (INCLUDING R1!!!)   
+!                                      amp(1) is the coeff. of 1/eps   pole
+!                                      amp(2) is the coeff. of 1/eps^2 pole.
+!
+!          complex*16 ampcc         -> the Cut Constructible contribution
+!          complex*16 ampr1         -> the R1 contribution
+!                                      (NOTE that ampcc+ampr1 = amp(0))
+!
+!          logical stable           -> .false. if CutTools detects
+!                                      numerical instabilities.         
+!---------------------------------------------------------------------
+      call ctsxcut(imode,rootsvalue,muscale,number_propagators,test,
+     & mpfortest,rnk,pp,m2,amp,ampcc,ampr1,stable)
+      enddo
+      print*,'  '
+!---------------------------------------------------------------------
+! To know the statistics of the run call ctsstatistics(discarded)
+!
+! INPUT :  none
+!
+! OUTPUT:  logical discarded  -> .true. if there are discarded PS points
+!
+! Print out of the statistics of the run:
+!
+!          n_mp   ->  n.of points evaluated in multi-precision.
+!          n_disc ->  n.of discarded points.               
+!          n_tot  ->  n.of calls to ctsxcut               
+!---------------------------------------------------------------------
+      call ctsstatistics(discarded)
+      if (.not.discarded) then
+        print*,'  '  
+        print*,'  '  
+        print*,'Test on Multiprecision passed'
+        print*,'  '  
+        print*,'  '  
+      else
+        print*,'  '  
+        print*,'  '  
+        print*,'Test on Multiprecision NOT passed !!!'
+        print*,'  '  
+        print*,'  '  
+      endif
+      end program testmp
+
+
+ 
+
+
+
+
+

=== added file 'loop_material/CutTools/makefile'
--- loop_material/CutTools/makefile	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/makefile	2012-02-27 13:42:21 +0000
@@ -0,0 +1,61 @@
+PRECISION= MP
+#PRECISION= QP
+CTS_VERSION = v1.7.0
+SRC = ./src
+EXE = 
+FC  = gfortran
+FFLAGS = -fno-automatic -O2 -funroll-all-loops  
+BLD = includects
+CTS_DIR = cuttools_$(CTS_VERSION)
+CTS_TAR = $(CTS_DIR).tar.gz
+
+ARGS = \
+  EXE="$(EXE)" \
+  FC="$(FC)" \
+  FFLAGS="$(FFLAGS)" \
+
+ifeq ($(PRECISION),MP)
+#
+# For building of the version with internal multiprecision routines:
+#
+mp: cpmp clean$(BLD)
+else
+#
+# For building of the version with quadruple precision compiler (if present):
+#
+qp: cpqp clean$(BLD)
+endif
+
+cpmp:   
+	cp  -p ./src/cts/cts_mpr.in ./src/cts/cts_mpr.h
+	cp  -p ./src/cts/cts_mpc.in ./src/cts/cts_mpc.h
+
+cpqp:   
+	cp  -p ./src/cts/cts_qpr.in ./src/cts/cts_mpr.h
+	cp  -p ./src/cts/cts_qpc.in ./src/cts/cts_mpc.h
+
+clean$(BLD): default
+	rm -fr  $(BLD)/*.f
+	rm -fr  $(BLD)/*.f90
+	rm -fr  $(BLD)/*.o
+	rm -fr  $(BLD)/makefile
+
+default: force
+	cd $(BLD) && $(MAKE) $(ARGS) $@
+
+force: $(BLD)/version.h
+
+$(BLD)/version.h: 
+	-mkdir -p $(BLD)
+	cp  -p ./src/avh/* $(BLD)/ 
+	cp  -p ./src/cts/* $(BLD)/ 
+	cp  -p ./src/mpfun90/* $(BLD)/ 
+	cp  -p ./src/qcdloop/* $(BLD)/ 
+	cp  -p ./src/makefile $(BLD)/ 
+
+tar:
+	tar -czvf $(CTS_TAR) *
+
+clean:
+	rm -fr $(BLD) $(CTS_TAR)
+

=== added directory 'loop_material/CutTools/src'
=== added directory 'loop_material/CutTools/src/avh'
=== added file 'loop_material/CutTools/src/avh/avh_olo.f90'
--- loop_material/CutTools/src/avh/avh_olo.f90	1970-01-01 00:00:00 +0000
+++ loop_material/CutTools/src/avh/avh_olo.f90	2012-02-27 13:42:21 +0000
@@ -0,0 +1,8454 @@
+!
+!  Copyright (C) 2011 Andreas van Hameren. 
+!
+!  This is the package OneLOop-2.2.
+! 
+!  OneLOop-2.2 is free software: you can redistribute it and/or modify
+!  it under the terms of the GNU General Public License as published by
+!  the Free Software Foundation, either version 3 of the License, or
+!  (at your option) any later version.
+!
+!  OneLOop-2.2 is distributed in the hope that it will be useful,
+!  but WITHOUT ANY WARRANTY; without even the implied warranty of
+!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+!  GNU General Public License for more details.
+!
+!  You should have received a copy of the GNU General Public License
+!  along with OneLOop-2.2.  If not, see <http://www.gnu.org/licenses/>.
+!
+
+
+module avh_olo_xkind
+! use !XKIND_MODULE
+!
+  implicit none
+  private
+  public :: olo_xkind
+!
+  integer ,parameter :: olo_xkind=kind(1d0) !XKIND_SETPAR
+!  
+end module
+
+
+module avh_olo_kinds
+  use avh_olo_xkind
+!
+  implicit none
+  private
+  public :: kindr2,kindc2 &
+           ,R0P0,R1P0,R5M1,TWOPI,SQRT2,C0P0,C1P0,CiP0
+!
+  integer ,parameter :: kindr2 = olo_xkind
+  integer ,parameter :: kindc2 = kindr2
+!
+  include 'cts_dpr.h'
+   ,parameter :: R0P0=0._kindr2
+  include 'cts_dpr.h'
+   ,parameter :: R1P0=1._kindr2
+  include 'cts_dpr.h'
+   ,parameter :: R5M1=0.5_kindr2
+!                                  1 2345678901234567890123456789012 
+  include 'cts_dpr.h'
+   ,parameter :: TWOPI=6.2831853071795864769252867665590_kindr2
+  include 'cts_dpr.h'
+   ,parameter :: SQRT2=1.4142135623730950488016887242097_kindr2
+  include 'cts_dpc.h'
+   ,parameter :: C0P0 = (0._kindr2,0._kindr2)
+  include 'cts_dpc.h'
+   ,parameter :: C1P0 = (1._kindr2,0._kindr2)
+  include 'cts_dpc.h'
+   ,parameter :: CiP0 = (0._kindr2,1._kindr2)
+!  
+end module
+
+
+module avh_olo_units
+  implicit none
+  integer :: eunit=6 !PROTECTED
+  integer :: wunit=6 !PROTECTED
+  integer :: munit=6 !PROTECTED
+  integer :: punit=0 !PROTECTED ! print all
+contains
+  subroutine set_unit( message ,val )
+!***********************************************************************
+! message is intended to be one of the following:
+! 'printall', 'message' ,'warning' ,'error'
+!***********************************************************************
+  character(*) ,intent(in) :: message
+  integer      ,intent(in) :: val
+  if (.false.) then
+  elseif (message(1:8).eq.'printall') then ;punit=val
+  elseif (message(1:7).eq.'message' ) then ;munit=val
+  elseif (message(1:7).eq.'warning' ) then ;wunit=val
+  elseif (message(1:5).eq.'error'   ) then ;eunit=val
+  else
+    eunit=val
+    wunit=val
+    munit=val
+    punit=0
+  endif
+  end subroutine
+end module
+
+
+module avh_olo_print
+  use avh_olo_kinds
+  implicit none
+  private
+  public :: myprint,init_print
+!
+  integer ,parameter :: noverh=10 !maximally 6 decimals for exponent
+  integer :: ndigits=19
+  integer :: nefrmt=19+noverh
+!
+  interface myprint
+    module procedure printr,printc,printi
+  end interface
+!
+contains
+!
+  subroutine init_print( ndig )
+  integer ,intent(in) :: ndig
+  ndigits = ndig+ndig/4+1
+  nefrmt  = ndigits+noverh
+  end subroutine
+!
+  function printc( zz ) result(rslt)
+  include 'cts_dpc.h'
+   ,intent(in) :: zz
+  character(nefrmt*2+3) :: rslt
+  rslt = '('//trim(printr(real(zz))) &
+       //','//trim(printr(aimag(zz)           )) &
+       //')'
+  rslt = adjustl(rslt)
+  end function
+!
+  function printr( xx ) result(rslt)
+  include 'cts_dpr.h'
+   ,intent(in) :: xx
+  character(nefrmt  ) :: rslt
+  character(nefrmt+1) :: cc
+  character(10) :: aa,bb
+  write(aa,'(i10)') nefrmt+1  ;aa=adjustl(aa)
+  write(bb,'(i10)') ndigits   ;bb=adjustl(bb)
+  aa = '(e'//trim(aa)//'.'//trim(bb)//')'
+  write(cc,aa) xx  ;cc=adjustl(cc)
+  if (cc(1:2).eq.'-0') then ;rslt = '-'//cc(3:ndigits*2)
+  else                      ;rslt = ' '//cc(2:ndigits*2)
+  endif
+  end function
+!
+  function printi( ii ) result(rslt)
+  integer ,intent(in) :: ii
+  character(ndigits) :: rslt
+  character(ndigits) :: cc
+  character(10) :: aa
+  write(aa,'(i10)') ndigits ;aa=adjustl(aa)
+  aa = '(i'//trim(aa)//')'
+  write(cc,aa) ii ;cc=adjustl(cc)
+  if (cc(1:1).ne.'-') then ;rslt=' '//cc
+  else                     ;rslt=cc 
+  endif
+  end function
+!
+end module
+
+
+module avh_olo_func
+  use avh_olo_kinds
+  use avh_olo_units
+!
+  implicit none
+!
+  type :: qmplx_type
+   include 'cts_dpc.h'
+    :: c
+   integer         :: p
+  end type
+!
+  interface mysqrt
+    module procedure mysqrt_0,mysqrt_r,mysqrt_i
+  end interface
+!
+  interface qonv
+    module procedure qonv_r,qonv_0,qonv_i
+  end interface
+!
+  interface operator (*)
+    module procedure prduct,prduct_r
+  end interface
+  interface operator (/)
+    module procedure ratio,ratio_r
+  end interface
+!
+contains
+!
+!
+   function mysqrt_0(xx) result(rslt)
+!*******************************************************************
+! Returns the square-root of xx .
+! If  Im(xx)  is equal zero and  Re(xx)  is negative, the result is
+! negative imaginary.
+!*******************************************************************
+   include 'cts_dpc.h'
+    ,intent(in) :: xx
+   include 'cts_dpc.h'
+    :: rslt ,zz
+   include 'cts_dpr.h'
+    :: xim,xre
+   xim = aimag(xx)
+   if (xim.eq.R0P0) then
+     xre = real(xx)
+     if (xre.ge.R0P0) then
+       zz = cmplx(sqrt(xre),R0P0,kind=kindc2)
+     else
+       zz = cmplx(R0P0,-sqrt(-xre),kind=kindc2)
+     endif
+   else
+     zz = sqrt(xx)
+   endif
+   rslt = zz
+   end function
+
+   function mysqrt_r(xx,sgn) result(rslt)
+!*******************************************************************
+! Returns the square-root of xx .
+! If  Im(xx)  is equal zero and  Re(xx)  is negative, the result is
+! imaginary and has the same sign as  sgn .
+!*******************************************************************
+   include 'cts_dpc.h'
+    ,intent(in) :: xx
+   include 'cts_dpr.h'
+    ,intent(in) :: sgn
+   include 'cts_dpc.h'
+    :: rslt ,zz
+   include 'cts_dpr.h'
+    :: xim,xre
+   xim = aimag(xx)
+   if (xim.eq.R0P0) then
+     xre = real(xx)
+     if (xre.ge.R0P0) then
+       zz = cmplx(sqrt(xre),R0P0,kind=kindc2)
+     else
+       zz = cmplx(R0P0,sign(sqrt(-xre),sgn),kind=kindc2)
+     endif
+   else
+     zz = sqrt(xx)
+   endif
+   rslt = zz
+   end function
+
+   function mysqrt_i(xx,sgn) result(rslt)
+!*******************************************************************
+! Returns the square-root of xx .
+! If  Im(xx)  is equal zero and  Re(xx)  is negative, the result is
+! imaginary and has the same sign as  sgn .
+!*******************************************************************
+   include 'cts_dpc.h'
+    ,intent(in) :: xx
+   integer         ,intent(in) :: sgn
+   include 'cts_dpc.h'
+    :: rslt ,zz
+   include 'cts_dpr.h'
+    :: xim,xre
+   xim = aimag(xx)
+   if (xim.eq.R0P0) then
+     xre = real(xx)
+     if (xre.ge.R0P0) then
+       zz = cmplx(sqrt(xre),R0P0,KIND=kindc2)
+     else
+       zz = cmplx(R0P0,sign(sqrt(-xre),real(sgn,KIND=kindr2)),KIND=kindc2)
+     endif
+   else
+     zz = sqrt(xx)
+   endif
+   rslt = zz
+   end function
+
+
+   subroutine solabc( x1,x2 ,dd ,aa,bb,cc ,imode )
+!*******************************************************************
+! Returns the solutions  x1,x2  to the equation  aa*x^2+bb*x+cc=0
+! Also returns  dd = aa*(x1-x2)
+! If  imode=/=0  it uses  dd  as input as value of  sqrt(b^2-4*a*c)
+!*******************************************************************
+   include 'cts_dpc.h'
+    ,intent(out)   :: x1,x2
+   include 'cts_dpc.h'
+    ,intent(inout) :: dd
+   include 'cts_dpc.h'
+    ,intent(in) :: aa,bb,cc
+   integer         ,intent(in) :: imode
+   include 'cts_dpc.h'
+    :: qq,hh
+   include 'cts_dpr.h'
+    :: r1,r2
+!
+   if (aa.eq.C0P0) then
+     if (bb.eq.C0P0) then
+       if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop solabc: ' &
+         ,'no solutions, returning 0'
+       x1 = C0P0
+       x2 = C0P0
+       dd = C0P0
+     else
+       x1 = -cc/bb
+       x2 = x1
+       dd = bb
+     endif
+   elseif (cc.eq.C0P0) then
+     dd = -bb
+     x1 = dd/aa
+     x2 = C0P0
+   else
+     if (imode.eq.0) dd = sqrt(bb*bb - 4*aa*cc)
+     qq = -bb+dd
+     hh = -bb-dd
+     r1 = abs(qq)
+     r2 = abs(hh)
+     if (r1.ge.r2) then
+       x1 = qq/(2*aa)
+       x2 = (2*cc)/qq
+     else
+       qq = hh
+       x2 = qq/(2*aa)
+       x1 = (2*cc)/qq
+     endif
+   endif
+   end subroutine
+
+
+   subroutine rfun(rr,dd ,qq)
+!*******************************************************************
+! Returns  rr  such that  qq = rr + 1/rr  and  Im(rr)  has the same
+! sign as  Im(qq) .
+! If  Im(qq)  is zero, then  Im(rr)  is negative or zero.
+! If  Im(rr)  is zero, then  |rr| > 1/|rr| .
+! Also returns  dd = rr - 1/rr .
+!*******************************************************************
+   include 'cts_dpc.h'
+    ,intent(out) :: rr,dd
+   include 'cts_dpc.h'
+    ,intent(in)  :: qq
+   include 'cts_dpc.h'
+    :: r2
+   include 'cts_dpr.h'
+    :: aa,bb
+   integer :: ir,ik
+   include 'cts_dpc.h'
+    ,parameter :: two=2*C1P0,four=4*C1P0
+   dd = sqrt(qq*qq-four)
+   rr = qq+dd
+   r2 = qq-dd
+   aa = abs(rr)
+   bb = abs(r2)
+   if (bb.gt.aa) then
+     rr = r2
+     dd = -dd
+   endif
+   aa = aimag(qq)
+   bb = aimag(rr)
+   if (aa.eq.R0P0) then
+     if (bb.le.R0P0) then
+       rr = rr/two
+     else
+       rr = two/rr
+       dd = -dd
+     endif
+   else
+     ik = int(sign(R1P0,aa))
+     ir = int(sign(R1P0,bb))
+     if (ir.eq.ik) then
+       rr = rr/two
+     else
+       rr = two/rr
+       dd = -dd
+     endif
+   endif
+   end subroutine
+
+   subroutine rfun0(rr ,dd,qq)
+!*******************************************************************
+! Like rfun, but now  dd  is input, which may get a minus sign
+!*******************************************************************
+   include 'cts_dpc.h'
+    ,intent(out)   :: rr
+   include 'cts_dpc.h'
+    ,intent(inout) :: dd
+   include 'cts_dpc.h'
+    ,intent(in)  :: qq
+   include 'cts_dpc.h'
+    :: r2
+   include 'cts_dpr.h'
+    :: aa,bb
+   integer :: ir,ik
+   include 'cts_dpc.h'
+    ,parameter :: two=2*C1P0
+   rr = qq+dd
+   r2 = qq-dd
+   aa = abs(rr)
+   bb = abs(r2)
+   if (bb.gt.aa) then
+     rr = r2
+     dd = -dd
+   endif
+   aa = aimag(qq)
+   bb = aimag(rr)
+   if (aa.eq.R0P0) then
+     if (bb.le.R0P0) then
+       rr = rr/two
+     else
+       rr = two/rr
+       dd = -dd
+     endif
+   else
+     ik = int(sign(R1P0,aa))
+     ir = int(sign(R1P0,bb))
+     if (ir.eq.ik) then
+       rr = rr/two
+     else
+       rr = two/rr
+       dd = -dd
+     endif
+   endif
+   end subroutine
+
+
+   function qonv_r(xx,sgn) result(rslt)
+!*******************************************************************
+! zz=rslt%c ,iz=rslt%p
+! Determine  zz,iz  such that  xx = zz*exp(iz*imag*pi)  and  Re(zz)
+! is positive. If  Im(x)=0  and  Re(x)<0  then  iz  becomes the
+! sign of  sgn .
+!*******************************************************************
+   include 'cts_dpc.h'
+    ,intent(in) :: xx
+   include 'cts_dpr.h'
+    ,intent(in) :: sgn
+   type(qmplx_type) :: rslt
+   include 'cts_dpr.h'
+    :: xre,xim
+   xre = real(xx)
+   if (xre.ge.R0P0) then
+     rslt%c = xx
+     rslt%p = 0
+   else
+     xim = aimag(xx)
+     if (xim.eq.R0P0) then
+       rslt%c = cmplx(-xre,R0P0,kind=kindc2)
+       rslt%p = int(sign(R1P0,sgn))
+     else
+       rslt%c = -xx
+       rslt%p = int(sign(R1P0,xim)) ! xim = -Im(rslt%c)
+     endif
+   endif
+   end function
+!
+   function qonv_i(xx,sgn) result(rslt)
+!*******************************************************************
+! zz=rslt%c ,iz=rslt%p
+! Determine  zz,iz  such that  xx = zz*exp(iz*imag*pi)  and  Re(zz)
+! is positive. If  Im(x)=0  and  Re(x)<0  then  iz  becomes the
+! sign of  sgn .
+!*******************************************************************
+   include 'cts_dpc.h'
+    ,intent(in) :: xx
+   integer         ,intent(in) :: sgn
+   type(qmplx_type) :: rslt
+   include 'cts_dpr.h'
+    :: xre,xim
+   xre = real(xx)
+   if (xre.ge.R0P0) then
+     rslt%c = xx
+     rslt%p = 0
+   else
+     xim = aimag(xx)
+     if (xim.eq.R0P0) then
+       rslt%c = cmplx(-xre,R0P0,kind=kindc2)
+       rslt%p = sign(1,sgn)
+     else
+       rslt%c = -xx
+       rslt%p = int(sign(R1P0,xim)) ! xim = -Im(rslt%c)
+     endif
+   endif
+   end function
+!
+   function qonv_0(xx) result(rslt)
+!*******************************************************************
+! zz=rslt%c ,iz=rslt%p
+! Determine  zz,iz  such that  xx = zz*exp(iz*imag*pi)  and  Re(zz)
+! is positive. If  Im(x)=0  and  Re(x)<0  then  iz=1
+!*******************************************************************
+   include 'cts_dpc.h'
+    ,intent(in) :: xx
+   type(qmplx_type) :: rslt
+   include 'cts_dpr.h'
+    :: xre,xim
+   xre = real(xx)
+   if (xre.ge.R0P0) then
+     rslt%c = xx
+     rslt%p = 0
+   else
+     xim = aimag(xx)
+     if (xim.eq.R0P0) then
+       if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop qonv: ' &
+         ,'negative input with undefined sign for the imaginary part, ' &
+         ,'putting +ieps'
+       rslt%c = cmplx(-xre,R0P0,kind=kindc2)
+       rslt%p = 1
+     else
+       rslt%c = -xx
+       rslt%p = int(sign(R1P0,xim)) ! xim = -Im(rslt%c)
+     endif
+   endif
+   end function
+!
+   function directly(xx,ix) result(rslt)
+!*******************************************************************
+!*******************************************************************
+   include 'cts_dpc.h'
+    ,intent(in) :: xx
+   integer         ,intent(in) :: ix
+   type(qmplx_type) :: rslt
+   rslt%c = xx
+   rslt%p = ix
+   end function
+
+
+   function sheet(xx) result(ii)
+!*******************************************************************
+! Returns the number of the Riemann-sheet (times 2) for the complex
+! number  xx*exp(ix*imag*pi) . The real part of xx is assumed to be
+! positive or zero. Examples:
+! xx=1+imag, ix=-1 -> ii= 0 
+! xx=1+imag, ix= 1 -> ii= 2 
+! xx=1-imag, ix=-1 -> ii=-2 
+! xx=1-imag, ix= 1 -> ii= 0 
+! xx=1     , ix= 1 -> ii= 0  convention that log(-1)=pi on
+! xx=1     , ix=-1 -> ii=-2  the principal Riemann-sheet
+!*******************************************************************
+   type(qmplx_type) ,intent(in) :: xx
+   integer :: ii,jj
+   include 'cts_dpr.h'
+    :: xim
+   ii = xx%p/2*2
+   jj = xx%p-ii
+   xim = aimag(xx%c)
+   if (xim.le.R0P0) then ! also xim=0 <==> log(-1)=pi, not -pi
+     if (jj.eq.-1) ii = ii-2
+   else
+     if (jj.eq. 1) ii = ii+2
+   endif
+   end function
+
+
+   function prduct(yy,xx) result(zz)
+!*******************************************************************
+! Return the product  zz  of  yy  and  xx  
+! keeping track of (the multiple of pi of) the phase %p such that
+! the real part of  zz%c  remains positive 
+!*******************************************************************
+   type(qmplx_type) ,intent(in) :: yy,xx
+   type(qmplx_type) :: zz
+   zz%c = yy%c*xx%c
+   zz%p = yy%p+xx%p
+   if (real(zz%c).lt.R0P0) then
+     zz%p = zz%p + int(sign(R1P0,aimag(xx%c)))
+     zz%c = -zz%c
+   endif
+   end function
+
+   function prduct_r(yy,xx) result(zz)
+!*******************************************************************
+! Return the product  zz  of  yy  and  xx  
+! keeping track of (the multiple of pi of) the phase %p such that
+! the real part of  zz%c  remains positive 
+!*******************************************************************
+   type(qmplx_type) ,intent(in) :: yy
+   include 'cts_dpr.h'
+    ,intent(in) :: xx
+   type(qmplx_type) :: zz
+   zz%c = yy%c*abs(xx)
+   zz%p = yy%p
+   end function
+
+   function ratio(yy,xx) result(zz)
+!*******************************************************************
+! Return the ratio  zz  of  yy  and  xx  
+! keeping track of (the multiple of pi of) the phase %p such that
+! the real part of  zz%c  remains positive 
+!*******************************************************************
+   type(qmplx_type) ,intent(in) :: yy,xx
+   type(qmplx_type) :: zz
+   zz%c = yy%c/xx%c
+   zz%p = yy%p-xx%p
+   if (real(zz%c).lt.R0P0) then
+     zz%p = zz%p - int(sign(R1P0,aimag(xx%c)))
+     zz%c = -zz%c
+   endif
+   end function
+!
+   function ratio_r(yy,xx) result(zz)
+!*******************************************************************
+!*******************************************************************
+   type(qmplx_type) ,intent(in) :: yy
+   include 'cts_dpr.h'
+    ,intent(in) :: xx
+   type(qmplx_type) :: zz
+   zz%c = yy%c/abs(xx)
+   zz%p = yy%p
+   end function
+!
+!
+   function eta3( aa,sa ,bb,sb ,cc,sc ) result(rslt)
+!*******************************************************************
+! 2*pi*imag times the result of
+!     theta(-Im(a))*theta(-Im(b))*theta( Im(c))
+!   - theta( Im(a))*theta( Im(b))*theta(-Im(c))
+! where a,b,c are interpreted as a+i|eps|sa, b+i|eps|sb, c+i|eps|sc
+!*******************************************************************
+   include 'cts_dpc.h'
+    ,intent(in) :: aa,bb,cc
+   include 'cts_dpr.h'
+    ,intent(in) :: sa,sb,sc
+   include 'cts_dpc.h'
+    :: rslt
+   include 'cts_dpr.h'
+    :: ima,imb,imc
+   ima = aimag(aa)
+   imb = aimag(bb)
+   imc = aimag(cc)
+   if (ima.eq.R0P0) ima = sa
+   if (imb.eq.R0P0) imb = sb
+   if (imc.eq.R0P0) imc = sc
+   ima = sign(R1P0,ima)
+   imb = sign(R1P0,imb)
+   imc = sign(R1P0,imc)
+   if (ima.eq.imb.and.ima.ne.imc) then
+     rslt = cmplx(R0P0,imc*TWOPI,kind=kindc2)
+   else
+     rslt = R0P0
+   endif
+   end function
+ 
+   function eta2( aa,sa ,bb,sb ) result(rslt)
+!*******************************************************************
+! The same as  eta3, but with  c=a*b, so that
+!   eta(a,b) = log(a*b) - log(a) - log(b)
+!*******************************************************************
+   include 'cts_dpc.h'
+    ,intent(in) :: aa,bb
+   include 'cts_dpr.h'
+    ,intent(in) :: sa,sb
+   include 'cts_dpr.h'
+    :: rslt ,rea,reb,ima,imb,imab
+   rea = real(aa)  ;ima = aimag(aa)
+   reb = real(bb)  ;imb = aimag(bb)
+   imab = rea*imb + reb*ima
+   if (ima.eq.R0P0) ima = sa
+   if (imb.eq.R0P0) imb = sb
+   if (imab.eq.R0P0) imab = sign(rea,sb) + sign(reb,sa)
+   ima  = sign(R1P0,ima)
+   imb  = sign(R1P0,imb)
+   imab = sign(R1P0,imab)
+   if (ima.eq.imb.and.ima.ne.imab) then
+     rslt = cmplx(R0P0,imab*TWOPI,kind=kindc2)
+   else
+     rslt = R0P0
+   endif
+   end function 
+!
+end module
+
+
+module avh_olo_loga
+!*******************************************************************
+! log( |xx|*exp(imag*pi*iph) ) = log|xx| + imag*pi*iph
+!*******************************************************************
+  use avh_olo_kinds
+  use avh_olo_units
+  use avh_olo_func
+  implicit none
+  private
+  public :: loga
+  include 'cts_dpr.h'
+   ,parameter :: pi=TWOPI/2
+contains
+!
+  function loga(xx,iph) result(rslt)
+  include 'cts_dpr.h'
+   ,intent(in) :: xx
+  integer      ,intent(in) :: iph
+  include 'cts_dpc.h'
+   :: rslt
+  include 'cts_dpr.h'
+   :: rr
+!
+  rr = abs(xx)
+  if (rr.eq.R0P0) then
+    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop loga: ' &
+      ,'|xx|=',rr
+  endif
+  rslt = cmplx(log(rr),iph*pi,kind=kindc2)
+  end function
+!
+end module
+
+
+module avh_olo_bern
+!*******************************************************************
+! the first nn Bernoulli numbers
+!*******************************************************************
+  use avh_olo_kinds
+  implicit none
+  private
+  public :: init_bern,rbern,cbern
+  integer ,parameter :: nn=40
+  include 'cts_dpr.h'
+   :: rbern(nn) !PROTECTED
+  include 'cts_dpc.h'
+   :: cbern(nn) !PROTECTED
+  integer :: ndigits=0
+contains
+!  
+  subroutine init_bern(ndig)
+  integer ,intent(in) :: ndig
+  integer :: jj
+  integer ,parameter :: d=kindr2
+  if (ndigits.eq.ndig) return ;ndigits=ndig
+  rbern(1:nn) = R0P0
+  rbern( 1) = -1._d/2._d
+  rbern( 2) =  1._d/6._d
+  rbern( 4) = -1._d/30._d
+  rbern( 6) =  1._d/42._d
+  rbern( 8) = -1._d/30._d
+  rbern(10) =  5._d/66._d
+  rbern(12) = -691._d/2730._d
+  rbern(14) =  7._d/6._d
+  rbern(16) = -3617._d/510._d
+  rbern(18) =  43867._d/798._d
+  rbern(20) = -174611._d/330._d
+  rbern(22) =  854513._d/138._d
+  rbern(24) = -236364091._d/2730._d
+  rbern(26) =  8553103._d/6._d
+  rbern(28) = -23749461029._d/870._d
+  rbern(30) =  8615841276005._d/14322._d
+  rbern(32) = -7709321041217._d/510._d
+  rbern(34) =  2577687858367._d/6._d
+  rbern(36) = -26315271553053477373._d/1919190._d
+  rbern(38) =  2929993913841559._d/6._d
+  rbern(40) = -261082718496449122051._d/13530._d
+  do jj=1,nn
+    cbern(jj) = cmplx(rbern(jj),kind=kindc2)
+  enddo
+  end subroutine
+!
+end module
+
+
+module avh_olo_li2a
+!*******************************************************************
+!                  /1    ln(1-zz*t)
+! avh_olo_li2a = - |  dt ---------- 
+!                  /0        t
+! with  zz = 1 - |xx|*exp(imag*pi*iph)
+! Examples:
+!   In order to get the dilog of  1.1  use  xx=1.1, iph=0
+!   In order to get the dilog of -1.1  use  xx=1.1, iph=1
+! Add multiples of  2  to  iph  in order to get the result on
+! different Riemann-sheets.
+!*******************************************************************
+  use avh_olo_kinds
+  use avh_olo_func
+  use avh_olo_bern
+  implicit none
+  private
+  public :: init_li2a,li2a
+  include 'cts_dpr.h'
+   ,parameter :: pi=TWOPI/2
+  include 'cts_dpc.h'
+   ,parameter :: pi2o6=C1P0*TWOPI*TWOPI/24
+  integer :: nn=16
+  integer :: ndigits=0
+contains
+!
+  subroutine init_li2a(ndig)
+  integer ,intent(in) :: ndig
+  if (ndigits.eq.ndig) return ;ndigits=ndig
+  call init_bern(ndigits)
+  if     (ndigits.lt.24) then
+    nn = 16
+  else
+    nn = 30
+  endif
+  end subroutine
+
+  function li2a(xx,iph) result(rslt)
+  include 'cts_dpr.h'
+   ,intent(in) :: xx
+  integer      ,intent(in) :: iph
+  include 'cts_dpc.h'
+   :: rslt
+  include 'cts_dpr.h'
+   :: rr,yy,lyy,loy,zz,z2,liox
+  integer :: ii,ntwo,ione
+  logical :: positive , r_gt_1 , y_lt_h
+!
+  rr = abs(xx)
+  ntwo = iph/2*2
+  ione = iph - ntwo
+  positive = (ione.eq.0)
+! 
+  if     (rr.eq.R0P0) then
+    rslt = pi2o6
+  elseif (rr.eq.R1P0.and.positive) then
+    rslt = C0P0
+  else
+    yy  = rr
+    lyy = log(rr)
+    if (.not.positive) yy = -yy
+!
+    r_gt_1 = (rr.gt.R1P0)
+    if (r_gt_1) then
+      yy   = R1P0/yy
+      lyy  = -lyy
+      ntwo = -ntwo
+      ione = -ione
+    endif
+    loy = log(R1P0-yy) ! log(1-yy) is always real
+!
+    y_lt_h = (yy.lt.R5M1)
+    if (y_lt_h) then
+      zz = -loy ! log(1-yy) is real
+    else
+      zz = -lyy ! yy>0.5 => log(yy) is real
+    endif
+!
+    z2 = zz*zz
+    liox = rbern(nn)
+    do ii=nn,4,-2
+      liox = rbern(ii-2) + liox*z2/(ii*(ii+1))
+    enddo
+    liox = rbern(1) + liox*zz/3
+    liox = zz + liox*z2/2
+!
+    rslt = cmplx(liox,kind=kindc2)
+!
+    if (y_lt_h) then
+      rslt = pi2o6 - rslt - cmplx(loy*lyy,loy*pi*ione,kind=kindc2)
+    endif
+!
+    rslt = rslt + cmplx( R0P0 , -loy*pi*ntwo ,kind=kindc2)
+!
+    if (r_gt_1) rslt = -rslt - cmplx(-lyy,iph*pi,kind=kindc2)**2/2
+  endif
+  end function
+!
+end module
+
+
+module avh_olo_loga2
+!*******************************************************************
+! log(xx)/(1-xx)  with  xx = log|xx| + imag*pi*iph
+!*******************************************************************
+  use avh_olo_kinds
+  use avh_olo_units
+  use avh_olo_func
+  implicit none
+  private
+  public :: init_loga2,loga2
+  include 'cts_dpr.h'
+   :: thrs=epsilon(R1P0)
+  integer :: ndigits=0
+contains
+!
+  subroutine init_loga2(ndig)
+  integer ,intent(in) :: ndig
+  if (ndigits.eq.ndig) return ;ndigits=ndig
+  thrs = 10*thrs
+  end subroutine
+!
+  function loga2(xx,iph) result(rslt)
+  use avh_olo_loga ,only : loga
+  include 'cts_dpr.h'
+   ,intent(in) :: xx
+  integer      ,intent(in) :: iph
+  include 'cts_dpc.h'
+   :: rslt
+  include 'cts_dpr.h'
+   :: omx
+!
+  if (mod(iph,2).eq.0) then
+    omx = R1P0-abs(xx)
+  else
+    omx = R1P0+abs(xx)
+  endif
+!
+  if (iph.ne.0) then
+    if (omx.eq.R0P0) then
+      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop loga2: ' &
+        ,'1-xx,iph=',omx,iph
+      rslt = C0P0
+    else
+      rslt = loga(xx,iph)/cmplx(omx,kind=kindc2)
+    endif
+  else
+    if (abs(omx).lt.thrs) then
+      rslt = cmplx(-R1P0-omx/2,kind=kindc2)
+    else
+      rslt = loga(xx,iph)/cmplx(omx,kind=kindc2)
+    endif
+  endif
+  end function
+!
+end module
+
+
+module avh_olo_logc
+!*******************************************************************
+! Returns  log( |Re(xx)| + imag*Im(xx) ) + imag*pi*iph
+!*******************************************************************
+  use avh_olo_kinds
+  use avh_olo_units
+  use avh_olo_func
+  implicit none
+  private
+  public :: logc
+  include 'cts_dpc.h'
+   ,parameter :: ipi=CiP0*TWOPI/2
+contains
+!
+  function logc(xx) result(rslt)
+  type(qmplx_type) ,intent(in) :: xx
+  include 'cts_dpc.h'
+   :: rslt
+  if (xx%c.eq.C0P0) then
+    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop logc: xx%c =',xx%c
+    rslt = C0P0
+  else
+    rslt = log( cmplx(abs(real(xx%c)),aimag(xx%c),kind=kindc2) ) &
+         + ipi*xx%p
+  endif
+  end function
+!
+end module
+
+
+module avh_olo_li2c
+!*******************************************************************
+!                  /1    ln(1-zz*t)
+! avh_olo_li2c = - |  dt ---------- 
+!                  /0        t
+! with  zz = 1 - ( |Re(xx)| + imag*Im(xx) )*exp(imag*pi*iph)
+! Examples:
+!   In order to get the dilog of  1+imag  use  xx=1+imag, iph= 0
+!   In order to get the dilog of  1-imag  use  xx=1-imag, iph= 0
+!   In order to get the dilog of -1+imag  use  xx=1-imag, iph= 1
+!   In order to get the dilog of -1-imag  use  xx=1+imag, iph=-1
+! Add multiples of  2  to  iph  in order to get the result on
+! different Riemann-sheets.
+!*******************************************************************
+  use avh_olo_kinds
+  use avh_olo_func
+  use avh_olo_bern
+  use avh_olo_li2a
+  implicit none
+  private
+  public :: init_li2c,li2c
+  include 'cts_dpc.h'
+   ,parameter :: ipi=CiP0*TWOPI/2
+  include 'cts_dpc.h'
+   ,parameter :: pi2o6=C1P0*TWOPI*TWOPI/24
+  integer :: nn=18
+  integer :: ndigits=0
+contains
+!
+  subroutine init_li2c(ndig)
+  integer ,intent(in) :: ndig
+  if (ndigits.eq.ndig) return ;ndigits=ndig
+  call init_li2a(ndigits)
+  call init_bern(ndigits)
+  if     (ndigits.lt.24) then
+    nn = 18
+  else
+    nn = 36
+  endif
+  end subroutine
+
+  function li2c(xx) result(rslt)
+  type(qmplx_type) :: xx
+  include 'cts_dpc.h'
+   :: rslt ,yy,lyy,loy,zz,z2
+  include 'cts_dpr.h'
+   :: rex,imx
+  integer :: ii,iyy
+  logical :: x_gt_1 , y_lt_h
+!
+  rex = real(xx%c)
+  imx = aimag(xx%c)
+! 
+  if (imx.eq.R0P0) then
+    rslt = li2a(rex,xx%p)
+  else
+    rex = abs(rex)
+!
+    if (mod(xx%p,2).eq.0) then
+      yy = cmplx(rex,imx,kind=kindc2)
+      iyy = xx%p
+    else
+      yy = cmplx(-rex,-imx,kind=kindc2)
+! Notice that  iyy=xx%p/2*2  does not deal correctly with the
+! situation when  xx%p-xx%p/2*2 = sign(Im(xx%c)) . The following does:
+      iyy = xx%p + nint(sign(R1P0,imx))
+    endif
+!
+    x_gt_1 = (abs(xx%c).gt.R1P0)
+    if (x_gt_1) then
+      yy = C1P0/yy
+      iyy = -iyy
+    endif
+    lyy = log(yy)
+    loy = log(C1P0-yy)
+!
+    y_lt_h = (real(yy).lt.R5M1)
+    if (y_lt_h) then
+      zz = -loy
+    else
+      zz = -lyy
+    endif
+!
+    z2 = zz*zz
+    rslt = cbern(nn)
+    do ii=nn,4,-2
+      rslt = cbern(ii-2) + rslt*z2/(ii*(ii+1))
+    enddo
+    rslt = cbern(1) + rslt*zz/3
+    rslt = zz + rslt*z2/2
+!
+    if (y_lt_h) rslt = pi2o6 - rslt - loy*lyy
+!
+    rslt = rslt - loy*ipi*iyy
+!
+    if (x_gt_1) rslt = -rslt - (lyy + ipi*iyy)**2/2
+  endif
+  end function
+!
+end module
+
+
+module avh_olo_logc2
+!*******************************************************************
+! log(xx)/(1-xx)
+! with  log(xx) = log( |Re(xx)| + imag*Im(xx) ) + imag*pi*iph
+!*******************************************************************
+  use avh_olo_kinds
+  use avh_olo_units
+  use avh_olo_func
+  implicit none
+  private
+  public :: init_logc2,logc2
+  include 'cts_dpr.h'
+   :: thrs=epsilon(R1P0)
+  integer :: ndigits=0
+contains
+!
+  subroutine init_logc2(ndig)
+  integer ,intent(in) :: ndig
+  if (ndigits.eq.ndig) return ;ndigits=ndig
+  thrs = 10*thrs
+  end subroutine
+!
+  function logc2(xx) result(rslt)
+  use avh_olo_logc ,only : logc
+  type(qmplx_type) ,intent(in) :: xx
+  include 'cts_dpc.h'
+   :: rslt ,omx
+  if (mod(xx%p,2).eq.0) then
+    omx = cmplx(1d0-abs(real(xx%c)),-aimag(xx%c),kind=kindc2)
+  else
+    omx = cmplx(1d0+abs(real(xx%c)), aimag(xx%c),kind=kindc2)
+  endif
+  if (xx%p.ne.0) then
+    if (omx.eq.C0P0) then
+      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop logc2: ' &
+        ,'1-xx%c,xx%p=',omx,xx%p
+      rslt = C0P0
+    else
+      rslt = logc(xx)/omx
+    endif
+  else
+    if (abs(omx).lt.thrs) then
+      rslt = -C1P0-omx/2
+    else
+      rslt = logc(xx)/omx
+    endif
+  endif
+  end function
+!
+end module
+
+
+module avh_olo_li2c2
+!*******************************************************************
+! avh_olo_li2c2 = ( li2(x1) - li2(x2) )/(x1%c-x2%c)
+!
+!                    /1    ln(1-zz*t)
+! where  li2(x1) = - |  dt ----------
+!                    /0        t
+! with  zz = 1 - ( |Re(x1%c)| + imag*Im(x1%c) )*exp(imag*pi*x1%p)
+! and similarly for li2(x2)
+!*******************************************************************
+  use avh_olo_kinds
+  use avh_olo_units
+  use avh_olo_func
+  use avh_olo_li2c
+  use avh_olo_logc2
+  implicit none
+  private
+  public :: init_li2c2,li2c2
+  include 'cts_dpc.h'
+   ,parameter :: ipi=CiP0*TWOPI/2
+  include 'cts_dpr.h'
+   ,parameter :: thrs1=epsilon(R1P0)
+  include 'cts_dpr.h'
+   :: thrs=0.11_kindr2
+  integer      :: nmax=12
+  integer :: ndigits=0
+contains
+!
+  subroutine init_li2c2(ndig)
+  integer ,intent(in) :: ndig
+  if (ndigits.eq.ndig) return ;ndigits=ndig
+  call init_logc2(ndigits)
+  call init_li2c(ndigits)
+  if     (ndigits.lt.16) then
+    thrs = 0.11_kindr2 ! double precision
+    nmax = 12
+  elseif (ndigits.lt.24) then
+    thrs = 0.02_kindr2 ! guess
+    nmax = 12
+  else
+    thrs = 0.008_kindr2 ! quadruple precision
+    nmax = 12
+  endif
+  end subroutine
+!
+  function li2c2(x1,x2) result(rslt)
+  type(qmplx_type) ,intent(in) :: x1,x2
+  include 'cts_dpc.h'
+   :: rslt
+  include 'cts_dpc.h'
+   :: x1r,x2r,delta,xx,xr,omx,del,hh,ff(0:20),zz
+  integer :: ih,ii
+!
+  if (mod(x1%p,2).eq.0) then
+    x1r = cmplx( abs(real(x1%c)), aimag(x1%c),kind=kindc2)
+  else
+    x1r = cmplx(-abs(real(x1%c)),-aimag(x1%c),kind=kindc2)
+  endif     
+  if (mod(x2%p,2).eq.0) then
+    x2r = cmplx( abs(real(x2%c)), aimag(x2%c),kind=kindc2)
+  else
+    x2r = cmplx(-abs(real(x2%c)),-aimag(x2%c),kind=kindc2)
+  endif
+  delta = x1r-x2r
+!
+  if (x1%p.ne.x2%p) then
+    if (delta.eq.C0P0) then
+      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop li2c2: ' &
+        ,'x1%p,x2%p,delta=',x1%p,x2%p,delta
+      rslt = C0P0
+    else
+      rslt = ( li2c(x1)-li2c(x2) )/delta
+    endif
+  else
+    if (abs(delta/x1%c).gt.thrs) then
+      rslt = ( li2c(x1)-li2c(x2) )/delta
+    else
+      xx  = x1%c
+      xr  = x1r
+      omx = C1P0-xr
+      del = delta
+      hh = C1P0-x2r
+      if (abs(hh).gt.abs(omx)) then