Merge lp:~maddevelopers/mg5amcnlo/aloha_new_C_routine into lp:~maddevelopers/mg5amcnlo/1.4.3

[Olivier Mattelaer]

I propose a specific review on those changes (independently of those in 1.4.3)
since they are potentially quite dangerous (at least for MSSM).

They are three changes with these branches.
First one 'aesthetical' modification (which simplifies one of the bug fixing and is require for spoin3/2)
- a flip of the argument for the C routine. (this is in fact a "technical"
flip which has no impact at MG5 level, but which will be helpfull for the
implementation of spoin3/2)
so before
def FFVC1_0(F1 , F2, V, Coup):
    return F1[4]*F2[3]*V[1] *Coup
pass to
def FFVC1_0(F2 , F1, V, Coup):
    return F2[4]*F1[3]*V[1] *Coup
therefore this doesn't require any change in the way to call the ALOHA function.

Then two bug corrections (both discover due to spin3/2):
- the definition of Epsilon (mu, nu ,rho,sigma) which now is coherent with the FR definition
- The momentum (inside the Lorentz definition) were wrongly define in previous case in case of C routine. Partially due to flipping of the argument that were not present, and due to the fact that momentum convention are different for incoming and outcoming fermion.

Note that FFVC1_1, FFVC1_2 are not modified.
The case of 4 fermion should be dealt correctly (thanks to unittest)

[Johan Alwall (johan-alwall)]

Hi Olivier,

These are bugs and fixes that I can't really assess. Clearly, the first one is harmless, but for the bug fixes, I can't say. How did you find the bugs? And how do you plan to verify them?

I would suggest to create a couple of models which differ only in the order of Majorana particles to make sure we get the same result. This might be a bit too extensive to include in 1.4.3 though if you want to release that within a day or two.

Cheers,
Johan

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

Both are bugs discover due to problem with either goldstino/spin 3/2
(in a new model made by benjamin for spin3/2).
and discover those thanks to the "check" command.

The first one, is not that worrying the Epsilon(1,2,3,4) was never
use, and the definition in the UFO/ALOHA paper didn't specify the
convention for the sign.
The only problem linked to that bug fixing might be if someone
implemented a model by hand (not very probable).

The second one is the most tricky one.
In fact the first modification is not that harmless since this flip
the definition of the momenta. (which is part of the bug)

For the Epsilon, I have just edit the unit-test by flipping the
expecting sign. This is by far enough.
For the flip +sign: I've add a couple of unit-test associated to those
situations.
I think that this is enough for the moment, when the spin3/2
implementation will be part of the official release, then we can add
some
acceptance test using the "check" command associated to this model.

But at current stage, in addition to all the tests, I've checked that
all the "check" were passing for the two model associated with the
goldstino.

Cheers,

Olivier

On 23-févr.-12, at 20:25, Johan Alwall wrote:

> Hi Olivier,
>
> These are bugs and fixes that I can't really assess. Clearly, the
> first one is harmless, but for the bug fixes, I can't say. How did
> you find the bugs? And how do you plan to verify them?
>
> I would suggest to create a couple of models which differ only in
> the order of Majorana particles to make sure we get the same result.
> This might be a bit too extensive to include in 1.4.3 though if you
> want to release that within a day or two.
>
> Cheers,
> Johan
> --
> https://code.launchpad.net/~maddevelopers/madgraph5/
> aloha_new_C_routine/+merge/94406
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/
> madgraph5/1.4.3.

[lp:~maddevelopers/mg5amcnlo/aloha_new_C_routine]

211. By Johan Alwall on 2012-02-24

Add 7545 command for forbidden s-channels, and keep $ for forbidden onshell s-channels

212. By Johan Alwall on 2012-02-24

Ensure that ptj and mmjj < xqcut for xqcut > 0. If set larger, they are set to 0 unless auto_ptj_and_mjj is T.

213. By mattelaer-olivier on 2012-02-24

fix interface issue, allow cluster to be 5 time usresponsive for the submit/control function
add seed to the banner

214. By mattelaer-olivier on 2012-02-24

fix plot issue

[Johan Alwall (johan-alwall)]

Hello Olivier,

Great work. As I said, I can't really assess this, but if all parallel tests pass and the new version passes all checks (and gives the right result compared to analytical calculations for the models in question), then I'm happy to approve.

Cheers,
Johan

review: Approve

[lp:~maddevelopers/mg5amcnlo/aloha_new_C_routine]

215. By Johan Alwall on 2012-02-25

Fixed acceptance tests (once and for all), fixed symfact.dat writing independent of number of diagrams (works also for > 10000)

216. By Olivier Mattelaer on 2012-02-25

new implementation of sqve command

217. By Olivier Mattelaer on 2012-02-26

merge with the new cmd format.

218. By Olivier Mattelaer on 2012-02-26

merge with the aloha fixes

219. By Olivier Mattelaer on 2012-02-27

more secure way for the html double link trick.
fix parralel test for cross-section

220. By mattelaer-olivier on 2012-02-27

Change html output for the details of the cross-sections. (Merge all results on a single page)
Take care correctly of chrome restricting local call to ajax.

221. By Olivier Mattelaer on 2012-02-28

fix bug #942376

222. By Olivier Mattelaer on 2012-02-29

fix a multiple try problem (cluster only)
fix the ./bin/madevent script to close it properly

223. By mattelaer-olivier on 2012-02-29

fix the delphes problem on the web, simple optimisation, adding a new acceptance test

224. By mattelaer-olivier on 2012-03-01

New C routine for wavefunctions (not working for Cpp + not tested for combine routine)

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

I' re-using this branch for the last (?) aloha problem pointed by the spin 3/2.

so now for FF
_1 returns a wavefunction associates to the incoming fermion
_2 returns a wavefunction associates to the outcoming fermion

This has off course no effect on SM. (since this is in IO anyway)
For the MSSM, this is probably just equivalent to a renaming of the argumments in ALOHA.
For some NP, this might be important if momentum are present in the lorentz structure.
And off course this is crucial for spin3/2.

This might also explains the failure of Johan to solve the problem associated with the pid.

Note that this is not yet ready for review. Since
1) the Cpp output is not yet fixed
2) I didn't test the combine aloha output with this additional flipping
3) I didn't run any acceptance/parralel test. (but unit test passes, but one linked to Cpp output)

review: Needs Fixing

[lp:~maddevelopers/mg5amcnlo/aloha_new_C_routine]

225. By Olivier Mattelaer on 2012-03-02

fix Cpp writter

226. By Olivier Mattelaer on 2012-03-02

added test to ensure combine routines in presence of C routines

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

The two first point are fixed.
and the acceptance test and the short test are all working.

So I'm updating the spin 3/2 branch to see if this fix the associated problem.
And will run tomorrow parallel tests

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

The spoin3/2 requires some additional fixes compare to that versions,
but nothing that claim that the fix is conceptually wrong.
anyway it might be better to merge this branch in spin 3/2 directly.
(and merge 1.4.3 independently of this.)

Spin 3/2 is not yet ready for the merge since they are still some problem for some processes,
like a n1 > grv. But this rather points as a external routines problem (that's my current best guess).

Cheers,

Olivier

[Olivier Mattelaer (olivier-mattelaer)]

I don't know why it passed to merged here.
let put in rejected (since this is now include in spin 3/2 branch)

review: Disapprove

Unmerged revisions

226. By Olivier Mattelaer on 2012-03-02

added test to ensure combine routines in presence of C routines

225. By Olivier Mattelaer on 2012-03-02

fix Cpp writter

224. By mattelaer-olivier on 2012-03-01

New C routine for wavefunctions (not working for Cpp + not tested for combine routine)

223. By mattelaer-olivier on 2012-02-29

fix the delphes problem on the web, simple optimisation, adding a new acceptance test

222. By Olivier Mattelaer on 2012-02-29

fix a multiple try problem (cluster only)
fix the ./bin/madevent script to close it properly

221. By Olivier Mattelaer on 2012-02-28

fix bug #942376

220. By mattelaer-olivier on 2012-02-27

Change html output for the details of the cross-sections. (Merge all results on a single page)
Take care correctly of chrome restricting local call to ajax.

219. By Olivier Mattelaer on 2012-02-27

more secure way for the html double link trick.
fix parralel test for cross-section

218. By Olivier Mattelaer on 2012-02-26

merge with the aloha fixes

217. By Olivier Mattelaer on 2012-02-26

merge with the new cmd format.

=== modified file 'aloha/aloha_writers.py'
--- aloha/aloha_writers.py	2012-02-23 16:07:14 +0000
+++ aloha/aloha_writers.py	2012-03-02 06:37:23 +0000
@@ -34,10 +34,20 @@
         self.namestring = name
         self.abstractname = abstract_routine.name
         self.comment = abstract_routine.infostr
-        self.offshell = abstract_routine.outgoing 
+        self.offshell = abstract_routine.outgoing
+         
         self.symmetries = abstract_routine.symmetries
         self.tag = abstract_routine.tag
 
+        self.outgoing = self.offshell
+        if 'C%s' %((self.outgoing + 1) // 2) in self.tag:
+            #flip the outgoing tag if in conjugate
+            self.outgoing = self.outgoing + self.outgoing % 2 - (self.outgoing +1) % 2
+            
+
+
+
+
         #prepare the necessary object
         self.collect_variables() # Look for the different variables
         self.make_all_lists()    # Compute the expression for the call ordering
@@ -182,12 +192,11 @@
             
             if index in conjugate:
                 index2, spin2 = index+1, self.particles[index+1]
-                if index2 + 1 != outgoing: 
-                    call_arg.append('%s%d' % (spin2, index2 +1)) 
-                call_arg.append('%s%d' % (spin, index +1)) 
+                call_arg.append('%s%d' % (spin2, index2 +1)) 
+                #call_arg.append('%s%d' % (spin, index +1)) 
             elif index-1 in conjugate:
-                if index == outgoing:
-                    call_arg.append('%s%d' % (spin, index +1)) 
+                index2, spin2 = index-1, self.particles[index]
+                call_arg.append('%s%d' % (spin2, index2 +1)) 
             else:
                 call_arg.append('%s%d' % (spin, index +1)) 
                 
@@ -226,14 +235,14 @@
             global_sign = -1
         
         flipped = [2*(int(c[1:])-1) for c in self.tag if c.startswith('C')]
-        if self.offshell % 2:
-             not_flip = self.offshell - 1
-             if not_flip in flipped:
-                 flipped.remove(not_flip)
-        else:
-             not_flip = self.offshell - 2
-             if not_flip in flipped:
-                 flipped.remove(not_flip)           
+#        if self.offshell % 2:
+#             not_flip = self.offshell - 1
+#             if not_flip in flipped:
+#                 flipped.remove(not_flip)
+#        else:
+#             not_flip = self.offshell - 2
+#             if not_flip in flipped:
+#                 flipped.remove(not_flip)           
         
         for index, spin in enumerate(self.particles): 
             assert(spin in ['S','F','V','T'])  
@@ -252,7 +261,7 @@
                 if index-1 in flipped:
                     sign *= -1
             # No need to include the outgoing particles in the definitions
-            if index == self.offshell -1:
+            if index == self.outgoing -1:
                 continue 
             
             # write the
@@ -320,11 +329,11 @@
             if not alredy_update:
                 declare_list.append('double complex denom')
                 declare_list.append('double precision M%(id)d, W%(id)d' % 
-                                                          {'id': self.offshell})
+                                                          {'id': self.outgoing})
             call_arg = '%(args)s, COUP, M%(id)d, W%(id)d, %(spin)s%(id)d' % \
                     {'args': ', '.join(CallList), 
                      'spin': self.particles[self.offshell -1],
-                     'id': self.offshell}
+                     'id': self.outgoing}
             str_out = ' subroutine %s(%s)\n' % (name, call_arg) 
 
         # Forcing implicit None
@@ -371,7 +380,7 @@
             offshell_size = self.type_to_size[offshelltype]            
             #Implement the conservation of Energy Impulsion
             for i in range(-1,1):
-                str_out += '%s%d(%d)= ' % (offshelltype, self.offshell, \
+                str_out += '%s%d(%d)= ' % (offshelltype, self.outgoing, \
                                                               offshell_size + i)
                 
                 pat=re.compile(r'^[-+]?(?P<spin>\w)')
@@ -390,11 +399,11 @@
             if self.offshell == index and type in ['V','S']:
                 sign = -1
             if 'C%s' % ((index +1) // 2)  in self.tag: 
-                if index == self.offshell:
-                    pass
-                elif index % 2 and index -1 != self.offshell:
-                    pass
-                elif index %2 == 1 and index + 1  != self.offshell:
+                if index == self.outgoing:
+                    pass
+                elif index % 2 and index -1 != self.outgoing:
+                    pass
+                elif index % 2 == 1 and index + 1  != self.outgoing:
                     pass
                 else:
                     sign *= -1
@@ -464,7 +473,7 @@
                 OutString = OutString + string + '\n'
         else:
             OffShellParticle = '%s%d' % (self.particles[self.offshell-1],\
-                                                                  self.offshell)
+                                                                  self.outgoing)
             numerator = self.obj.numerator
             denominator = self.obj.denominator
             for ind in denominator.listindices():
@@ -577,11 +586,11 @@
         else:
             main = '%(spin)s%(id)d' % \
                           {'spin': self.particles[offshell -1],
-                           'id': self.offshell}
+                           'id': self.outgoing}
             call_arg = '%(args)s, %(COUP)s, M%(id)d, W%(id)d, %(LAST)s' % \
                     {'args': ', '.join(self.calllist['CallList']), 
                      'COUP':'COUP%d',
-                     'id': self.offshell,
+                     'id': self.outgoing,
                      'LAST': '%s'}
 
         # make the first call
@@ -723,8 +732,8 @@
             str_out = 'void %(name)s(%(args)s, double M%(number)d, double W%(number)d, complex<double>%(out)s%(number)d[])' % \
               {'name': name,
                'args': ','.join(CallList+ ['complex<double> COUP']),
-               'out': self.particles[OffShellParticle],
-               'number': OffShellParticle + 1 
+               'out': self.particles[self.outgoing - 1],
+               'number': self.outgoing 
                }
 
         h_string = str_out + ";\n\n"
@@ -755,7 +764,7 @@
             offshell_size = self.type_to_size[offshelltype]            
             #Implement the conservation of Energy Impulsion
             for i in range(-2,0):
-                str_out += '%s%d[%d]= ' % (offshelltype, self.offshell,
+                str_out += '%s%d[%d]= ' % (offshelltype, self.outgoing,
                                            offshell_size + i)
                 
                 pat=re.compile(r'^[-+]?(?P<spin>\w)')
@@ -775,11 +784,11 @@
             if self.offshell == index and type in ['V', 'S']:
                 sign = -1
             if 'C%s' % ((index +1) // 2)  in self.tag: 
-                if index == self.offshell:
-                    pass
-                elif index % 2 and index -1 != self.offshell:
-                    pass
-                elif index %2 == 1 and index + 1  != self.offshell:
+                if index == self.outgoing:
+                    pass
+                elif index % 2 and index -1 != self.outgoing:
+                    pass
+                elif index %2 == 1 and index + 1  != self.outgoing:
                     pass
                 else:
                     sign *= -1
@@ -845,7 +854,7 @@
                 string = string.replace('+-', '-')
                 OutString = OutString + string + ';\n'
         else:
-            OffShellParticle = self.particles[self.offshell-1]+'%s'%(self.offshell)
+            OffShellParticle = self.particles[self.offshell-1]+'%s'%(self.outgoing)
             numerator = self.obj.numerator
             denominator = self.obj.denominator
             for ind in denominator.listindices():
@@ -1115,7 +1124,7 @@
         
         WriteALOHA.__init__(self, abstract_routine, dirpath)
         self.outname = '%s%s' % (self.particles[self.offshell -1], \
-                                                               self.offshell)
+                                                               self.outgoing)
     
         
     def change_var_format(self, name): 
@@ -1192,7 +1201,7 @@
             str_out += 'def %(name)s(%(args)s, COUP, M%(id)d, W%(id)d):\n' % \
                 {'name': name,
                  'args': ', '.join(CallList), 
-                     'id': self.offshell}            
+                     'id': self.outgoing}            
         return str_out     
 
     def make_declaration_list(self):
@@ -1238,16 +1247,17 @@
             sign = 1
             if self.offshell == index and type in ['V','S']:
                 sign = -1
+
             if 'C%s' % ((index +1) // 2)  in self.tag: 
-                if index == self.offshell:
-                    pass
-                elif index % 2 and index -1 != self.offshell:
-                    pass
-                elif index %2 == 1 and index + 1  != self.offshell:
+                if index == self.outgoing:
+                    pass
+                elif index % 2 and index -1 != self.outgoing:
+                    pass
+                elif index %2 == 1 and index + 1  != self.outgoing:
                     pass
                 else:
                     sign *= -1
-            
+
             if sign == -1 :
                 sign = '- '
             else:
@@ -1340,11 +1350,11 @@
         else:
             main = '%(spin)s%(id)d' % \
                           {'spin': self.particles[self.offshell -1],
-                           'id': self.offshell}
+                           'id': self.outgoing}
             call_arg = '%(args)s, %(COUP)s, M%(id)d, W%(id)d' % \
                     {'args': ', '.join(self.calllist['CallList']), 
                      'COUP':'COUP%d',
-                     'id': self.offshell}
+                     'id': self.outgoing}
 
         # make the first call
         line = "    %s = %s%s("+call_arg+")\n"

=== modified file 'aloha/create_aloha.py'
--- aloha/create_aloha.py	2012-02-21 18:10:48 +0000
+++ aloha/create_aloha.py	2012-03-02 06:37:23 +0000
@@ -199,6 +199,9 @@
         #multiply by the wave functions
         nb_spinor = 0
         outgoing = self.outgoing
+        if (outgoing + 1) // 2 in self.conjg:
+            #flip the outgoing tag if in conjugate
+            outgoing = outgoing + outgoing % 2 - (outgoing +1) % 2
         
         if not self.routine_kernel:
             AbstractRoutineBuilder.counter += 1
@@ -215,21 +218,25 @@
         else:
             lorentz = self.routine_kernel
             aloha_lib.USE_TAG = set(self.kernel_tag) 
-        for (i, spin ) in enumerate(self.spins):
+            
+            
+        for (i, spin ) in enumerate(self.spins):   
             id = i + 1
+                     
             #Check if this is the outgoing particle
             if id == outgoing:
                 if spin == 1: 
                     lorentz *= complex(0,1)
                 elif spin == 2:
                     # shift and flip the tag if we multiply by C matrices
-                    if (id+1) // 2 in self.conjg:
-                        id += _conjugate_gap
-                    nb_spinor += 1
-                    if nb_spinor %2:
-                        lorentz *= SpinorPropagator(id, 'I2', self.outgoing)
+                    if (id + 1) // 2 in self.conjg:
+                        id += _conjugate_gap + id % 2 - (id +1) % 2
+                    if id % 2:
+                        #propagator outcoming
+                        lorentz *= SpinorPropagator(id, 'I2', outgoing)
                     else:
-                        lorentz *= SpinorPropagator('I2', id, self.outgoing) 
+                        #propagator incoming
+                        lorentz *= SpinorPropagator('I2', id, outgoing) 
                 elif spin == 3 :
                     lorentz *= VectorPropagator(id, 'I2', id)
                 elif spin == 5 :
@@ -250,15 +257,9 @@
                 elif spin == 2:
                     # shift the tag if we multiply by C matrices
                     if (id+1) // 2 in self.conjg:
-                        spin_id = id + _conjugate_gap
-                        if id % 2:
-                            if outgoing not in [id, id+1]:
-                                spin_id += 1
-                        elif outgoing not in [id, id-1]:
-                            spin_id -= 1
+                        spin_id = id + _conjugate_gap + id % 2 - (id +1) % 2
                     else:
                         spin_id = id
-                    nb_spinor += 1
                     lorentz *= Spinor(spin_id, id)
                 elif spin == 3:        
                     lorentz *= Vector(id, id)

=== modified file 'tests/unit_tests/various/test_aloha.py'
--- tests/unit_tests/various/test_aloha.py	2012-02-26 02:56:25 +0000
+++ tests/unit_tests/various/test_aloha.py	2012-03-02 06:37:23 +0000
@@ -3742,55 +3742,55 @@
         """ test that python writer works """
         
         solution ="""import wavefunctions
-def FFV1C1_1(F2, V3, COUP, M1, W1):
+def FFV1C1_1(F1, V3, COUP, M2, W2):
+    F2 = wavefunctions.WaveFunction(size=6)
+    F2[4] = F1[4]+V3[4]
+    F2[5] = F1[5]+V3[5]
+    P2 = [complex(F2[4]).real, \\
+             complex(F2[5]).real, \\
+             complex(F2[5]).imag, \\
+             complex(F2[4]).imag]
+    denom =1.0/(( (M2*( -M2+1j*W2))+( (P2[0]**2)-(P2[1]**2)-(P2[2]**2)-(P2[3]**2))))
+    F2[0]= COUP*denom*( (F1[1]*( (P2[0]*( 1j*V3[1]-V3[2]))+( (P2[1]*( -1j*V3[0]-1j*V3[3]))+( (P2[2]*( V3[0]+V3[3]))+(P2[3]*( 1j*V3[1]-V3[2]))))))+( (F1[0]*( (P2[0]*( -1j*V3[0]+1j*V3[3]))+( (P2[1]*( 1j*V3[1]+V3[2]))+( (P2[2]*( -V3[1]+1j*V3[2]))+(P2[3]*( -1j*V3[0]+1j*V3[3]))))))+(M2*( (F1[3]*( -1j*V3[1]+V3[2]))+(F1[2]*( -1j*V3[0]-1j*V3[3]))))))
+    F2[1]= COUP*denom*( (F1[1]*( (P2[0]*( -1j*V3[0]-1j*V3[3]))+( (P2[1]*( 1j*V3[1]-V3[2]))+( (P2[2]*( V3[1]+1j*V3[2]))+(P2[3]*( 1j*V3[0]+1j*V3[3]))))))+( (F1[0]*( (P2[0]*( 1j*V3[1]+V3[2]))+( (P2[1]*( -1j*V3[0]+1j*V3[3]))+( (P2[2]*( -V3[0]+V3[3]))+(P2[3]*( -1j*V3[1]-V3[2]))))))+(M2*( (F1[3]*( -1j*V3[0]+1j*V3[3]))+(F1[2]*( -1j*V3[1]-V3[2]))))))
+    F2[2]= COUP*denom*( (F1[3]*( (P2[0]*( -1j*V3[1]+V3[2]))+( (P2[1]*( 1j*V3[0]-1j*V3[3]))+( (P2[2]*( -V3[0]+V3[3]))+(P2[3]*( 1j*V3[1]-V3[2]))))))+( (F1[2]*( (P2[0]*( -1j*V3[0]-1j*V3[3]))+( (P2[1]*( 1j*V3[1]+V3[2]))+( (P2[2]*( -V3[1]+1j*V3[2]))+(P2[3]*( 1j*V3[0]+1j*V3[3]))))))+(M2*( (F1[1]*( 1j*V3[1]-V3[2]))+(F1[0]*( -1j*V3[0]+1j*V3[3]))))))
+    F2[3]= COUP*denom*( (F1[3]*( (P2[0]*( -1j*V3[0]+1j*V3[3]))+( (P2[1]*( 1j*V3[1]-V3[2]))+( (P2[2]*( V3[1]+1j*V3[2]))+(P2[3]*( -1j*V3[0]+1j*V3[3]))))))+( (F1[2]*( (P2[0]*( -1j*V3[1]-V3[2]))+( (P2[1]*( 1j*V3[0]+1j*V3[3]))+( (P2[2]*( V3[0]+V3[3]))+(P2[3]*( -1j*V3[1]-V3[2]))))))+(M2*( (F1[1]*( -1j*V3[0]-1j*V3[3]))+(F1[0]*( 1j*V3[1]+V3[2]))))))
+    return F2
+    
+    
+"""
+        
+        FFV = UFOLorentz(name = 'FFV1',
+                 spins = [ 2, 2, 3 ],
+                 structure = 'Gamma(3,2,1)')        
+        builder = create_aloha.AbstractRoutineBuilder(FFV)
+        builder.apply_conjugation()
+        amp = builder.compute_routine(1)
+        routine = amp.write(output_dir=None, language='Python')
+        
+        split_solution = solution.split('\n')
+        split_routine = routine.split('\n')
+        self.assertEqual(split_solution, split_routine)
+        self.assertEqual(len(split_routine), len(split_solution))
+
+        solution ="""import wavefunctions
+def FFV1C1_2(F2, V3, COUP, M1, W1):
     F1 = wavefunctions.WaveFunction(size=6)
-    F1[4] = F2[4]+V3[4]
-    F1[5] = F2[5]+V3[5]
+    F1[4] = F2[4]-V3[4]
+    F1[5] = F2[5]-V3[5]
     P1 = [complex(F1[4]).real, \\
              complex(F1[5]).real, \\
              complex(F1[5]).imag, \\
              complex(F1[4]).imag]
     denom =1.0/(( (M1*( -M1+1j*W1))+( (P1[0]**2)-(P1[1]**2)-(P1[2]**2)-(P1[3]**2))))
-    F1[0]= COUP*denom*( (F2[1]*( (P1[0]*( 1j*V3[1]-V3[2]))+( (P1[1]*( -1j*V3[0]-1j*V3[3]))+( (P1[2]*( V3[0]+V3[3]))+(P1[3]*( 1j*V3[1]-V3[2]))))))+( (F2[0]*( (P1[0]*( -1j*V3[0]+1j*V3[3]))+( (P1[1]*( 1j*V3[1]+V3[2]))+( (P1[2]*( -V3[1]+1j*V3[2]))+(P1[3]*( -1j*V3[0]+1j*V3[3]))))))+(M1*( (F2[3]*( -1j*V3[1]+V3[2]))+(F2[2]*( -1j*V3[0]-1j*V3[3]))))))
-    F1[1]= COUP*denom*( (F2[1]*( (P1[0]*( -1j*V3[0]-1j*V3[3]))+( (P1[1]*( 1j*V3[1]-V3[2]))+( (P1[2]*( V3[1]+1j*V3[2]))+(P1[3]*( 1j*V3[0]+1j*V3[3]))))))+( (F2[0]*( (P1[0]*( 1j*V3[1]+V3[2]))+( (P1[1]*( -1j*V3[0]+1j*V3[3]))+( (P1[2]*( -V3[0]+V3[3]))+(P1[3]*( -1j*V3[1]-V3[2]))))))+(M1*( (F2[3]*( -1j*V3[0]+1j*V3[3]))+(F2[2]*( -1j*V3[1]-V3[2]))))))
-    F1[2]= COUP*denom*( (F2[3]*( (P1[0]*( -1j*V3[1]+V3[2]))+( (P1[1]*( 1j*V3[0]-1j*V3[3]))+( (P1[2]*( -V3[0]+V3[3]))+(P1[3]*( 1j*V3[1]-V3[2]))))))+( (F2[2]*( (P1[0]*( -1j*V3[0]-1j*V3[3]))+( (P1[1]*( 1j*V3[1]+V3[2]))+( (P1[2]*( -V3[1]+1j*V3[2]))+(P1[3]*( 1j*V3[0]+1j*V3[3]))))))+(M1*( (F2[1]*( 1j*V3[1]-V3[2]))+(F2[0]*( -1j*V3[0]+1j*V3[3]))))))
-    F1[3]= COUP*denom*( (F2[3]*( (P1[0]*( -1j*V3[0]+1j*V3[3]))+( (P1[1]*( 1j*V3[1]-V3[2]))+( (P1[2]*( V3[1]+1j*V3[2]))+(P1[3]*( -1j*V3[0]+1j*V3[3]))))))+( (F2[2]*( (P1[0]*( -1j*V3[1]-V3[2]))+( (P1[1]*( 1j*V3[0]+1j*V3[3]))+( (P1[2]*( V3[0]+V3[3]))+(P1[3]*( -1j*V3[1]-V3[2]))))))+(M1*( (F2[1]*( -1j*V3[0]-1j*V3[3]))+(F2[0]*( 1j*V3[1]+V3[2]))))))
+    F1[0]= COUP*denom*( (F2[0]*( (V3[1]*( 1j*P1[1]+P1[2]))+( (V3[2]*( -P1[1]+1j*P1[2]))+( (V3[0]*( -1j*P1[0]+1j*P1[3]))+(V3[3]*( -1j*P1[0]+1j*P1[3]))))))+( (F2[1]*( (V3[0]*( 1j*P1[1]+P1[2]))+( (V3[3]*( -1j*P1[1]-P1[2]))+( (V3[1]*( -1j*P1[0]+1j*P1[3]))+(V3[2]*( -P1[0]+P1[3]))))))+(M1*( (F2[2]*( -1j*V3[0]+1j*V3[3]))+(F2[3]*( 1j*V3[1]+V3[2]))))))
+    F1[1]= COUP*denom*( (F2[0]*( (V3[1]*( -1j*P1[0]-1j*P1[3]))+( (V3[2]*( P1[0]+P1[3]))+( (V3[0]*( 1j*P1[1]-P1[2]))+(V3[3]*( 1j*P1[1]-P1[2]))))))+( (F2[1]*( (V3[0]*( -1j*P1[0]-1j*P1[3]))+( (V3[3]*( 1j*P1[0]+1j*P1[3]))+( (V3[1]*( 1j*P1[1]-P1[2]))+(V3[2]*( P1[1]+1j*P1[2]))))))+(M1*( (F2[2]*( 1j*V3[1]-V3[2]))+(F2[3]*( -1j*V3[0]-1j*V3[3]))))))
+    F1[2]= COUP*denom*( (F2[2]*( (V3[1]*( 1j*P1[1]+P1[2]))+( (V3[2]*( -P1[1]+1j*P1[2]))+( (V3[0]*( -1j*P1[0]-1j*P1[3]))+(V3[3]*( 1j*P1[0]+1j*P1[3]))))))+( (F2[3]*( (V3[0]*( -1j*P1[1]-P1[2]))+( (V3[3]*( -1j*P1[1]-P1[2]))+( (V3[1]*( 1j*P1[0]+1j*P1[3]))+(V3[2]*( P1[0]+P1[3]))))))+(M1*( (F2[0]*( -1j*V3[0]-1j*V3[3]))+(F2[1]*( -1j*V3[1]-V3[2]))))))
+    F1[3]= COUP*denom*( (F2[2]*( (V3[1]*( 1j*P1[0]-1j*P1[3]))+( (V3[2]*( -P1[0]+P1[3]))+( (V3[0]*( -1j*P1[1]+P1[2]))+(V3[3]*( 1j*P1[1]-P1[2]))))))+( (F2[3]*( (V3[0]*( -1j*P1[0]+1j*P1[3]))+( (V3[3]*( -1j*P1[0]+1j*P1[3]))+( (V3[1]*( 1j*P1[1]-P1[2]))+(V3[2]*( P1[1]+1j*P1[2]))))))+(M1*( (F2[0]*( -1j*V3[1]+V3[2]))+(F2[1]*( -1j*V3[0]+1j*V3[3]))))))
     return F1
     
     
 """
-        
-        FFV = UFOLorentz(name = 'FFV1',
-                 spins = [ 2, 2, 3 ],
-                 structure = 'Gamma(3,2,1)')        
-        builder = create_aloha.AbstractRoutineBuilder(FFV)
-        builder.apply_conjugation()
-        amp = builder.compute_routine(1)
-        routine = amp.write(output_dir=None, language='Python')
-        
-        split_solution = solution.split('\n')
-        split_routine = routine.split('\n')
-        self.assertEqual(split_solution, split_routine)
-        self.assertEqual(len(split_routine), len(split_solution))
-
-        solution ="""import wavefunctions
-def FFV1C1_2(F1, V3, COUP, M2, W2):
-    F2 = wavefunctions.WaveFunction(size=6)
-    F2[4] = F1[4]-V3[4]
-    F2[5] = F1[5]-V3[5]
-    P2 = [complex(F2[4]).real, \\
-             complex(F2[5]).real, \\
-             complex(F2[5]).imag, \\
-             complex(F2[4]).imag]
-    denom =1.0/(( (M2*( -M2+1j*W2))+( (P2[0]**2)-(P2[1]**2)-(P2[2]**2)-(P2[3]**2))))
-    F2[0]= COUP*denom*( (F1[0]*( (V3[1]*( 1j*P2[1]+P2[2]))+( (V3[2]*( -P2[1]+1j*P2[2]))+( (V3[0]*( -1j*P2[0]+1j*P2[3]))+(V3[3]*( -1j*P2[0]+1j*P2[3]))))))+( (F1[1]*( (V3[0]*( 1j*P2[1]+P2[2]))+( (V3[3]*( -1j*P2[1]-P2[2]))+( (V3[1]*( -1j*P2[0]+1j*P2[3]))+(V3[2]*( -P2[0]+P2[3]))))))+(M2*( (F1[2]*( -1j*V3[0]+1j*V3[3]))+(F1[3]*( 1j*V3[1]+V3[2]))))))
-    F2[1]= COUP*denom*( (F1[0]*( (V3[1]*( -1j*P2[0]-1j*P2[3]))+( (V3[2]*( P2[0]+P2[3]))+( (V3[0]*( 1j*P2[1]-P2[2]))+(V3[3]*( 1j*P2[1]-P2[2]))))))+( (F1[1]*( (V3[0]*( -1j*P2[0]-1j*P2[3]))+( (V3[3]*( 1j*P2[0]+1j*P2[3]))+( (V3[1]*( 1j*P2[1]-P2[2]))+(V3[2]*( P2[1]+1j*P2[2]))))))+(M2*( (F1[2]*( 1j*V3[1]-V3[2]))+(F1[3]*( -1j*V3[0]-1j*V3[3]))))))
-    F2[2]= COUP*denom*( (F1[2]*( (V3[1]*( 1j*P2[1]+P2[2]))+( (V3[2]*( -P2[1]+1j*P2[2]))+( (V3[0]*( -1j*P2[0]-1j*P2[3]))+(V3[3]*( 1j*P2[0]+1j*P2[3]))))))+( (F1[3]*( (V3[0]*( -1j*P2[1]-P2[2]))+( (V3[3]*( -1j*P2[1]-P2[2]))+( (V3[1]*( 1j*P2[0]+1j*P2[3]))+(V3[2]*( P2[0]+P2[3]))))))+(M2*( (F1[0]*( -1j*V3[0]-1j*V3[3]))+(F1[1]*( -1j*V3[1]-V3[2]))))))
-    F2[3]= COUP*denom*( (F1[2]*( (V3[1]*( 1j*P2[0]-1j*P2[3]))+( (V3[2]*( -P2[0]+P2[3]))+( (V3[0]*( -1j*P2[1]+P2[2]))+(V3[3]*( 1j*P2[1]-P2[2]))))))+( (F1[3]*( (V3[0]*( -1j*P2[0]+1j*P2[3]))+( (V3[3]*( -1j*P2[0]+1j*P2[3]))+( (V3[1]*( 1j*P2[1]-P2[2]))+(V3[2]*( P2[1]+1j*P2[2]))))))+(M2*( (F1[0]*( -1j*V3[1]+V3[2]))+(F1[1]*( -1j*V3[0]+1j*V3[3]))))))
-    return F2
-    
-    
-"""
 
         amp = builder.compute_routine(2)
         
@@ -3800,6 +3800,52 @@
         split_routine = routine.split('\n')
         self.assertEqual(split_solution, split_routine)
         self.assertEqual(len(split_routine), len(split_solution))
+
+        solution = """import wavefunctions
+def FFV1C1_1(F1, V3, COUP, M2, W2):
+    F2 = wavefunctions.WaveFunction(size=6)
+    F2[4] = F1[4]+V3[4]
+    F2[5] = F1[5]+V3[5]
+    P2 = [complex(F2[4]).real, \\
+             complex(F2[5]).real, \\
+             complex(F2[5]).imag, \\
+             complex(F2[4]).imag]
+    denom =1.0/(( (M2*( -M2+1j*W2))+( (P2[0]**2)-(P2[1]**2)-(P2[2]**2)-(P2[3]**2))))
+    F2[0]= COUP*denom*( (F1[1]*( (P2[0]*( 1j*V3[1]-V3[2]))+( (P2[1]*( -1j*V3[0]-1j*V3[3]))+( (P2[2]*( V3[0]+V3[3]))+(P2[3]*( 1j*V3[1]-V3[2]))))))+( (F1[0]*( (P2[0]*( -1j*V3[0]+1j*V3[3]))+( (P2[1]*( 1j*V3[1]+V3[2]))+( (P2[2]*( -V3[1]+1j*V3[2]))+(P2[3]*( -1j*V3[0]+1j*V3[3]))))))+(M2*( (F1[3]*( -1j*V3[1]+V3[2]))+(F1[2]*( -1j*V3[0]-1j*V3[3]))))))
+    F2[1]= COUP*denom*( (F1[1]*( (P2[0]*( -1j*V3[0]-1j*V3[3]))+( (P2[1]*( 1j*V3[1]-V3[2]))+( (P2[2]*( V3[1]+1j*V3[2]))+(P2[3]*( 1j*V3[0]+1j*V3[3]))))))+( (F1[0]*( (P2[0]*( 1j*V3[1]+V3[2]))+( (P2[1]*( -1j*V3[0]+1j*V3[3]))+( (P2[2]*( -V3[0]+V3[3]))+(P2[3]*( -1j*V3[1]-V3[2]))))))+(M2*( (F1[3]*( -1j*V3[0]+1j*V3[3]))+(F1[2]*( -1j*V3[1]-V3[2]))))))
+    F2[2]= COUP*denom*( (F1[3]*( (P2[0]*( -1j*V3[1]+V3[2]))+( (P2[1]*( 1j*V3[0]-1j*V3[3]))+( (P2[2]*( -V3[0]+V3[3]))+(P2[3]*( 1j*V3[1]-V3[2]))))))+( (F1[2]*( (P2[0]*( -1j*V3[0]-1j*V3[3]))+( (P2[1]*( 1j*V3[1]+V3[2]))+( (P2[2]*( -V3[1]+1j*V3[2]))+(P2[3]*( 1j*V3[0]+1j*V3[3]))))))+(M2*( (F1[1]*( 1j*V3[1]-V3[2]))+(F1[0]*( -1j*V3[0]+1j*V3[3]))))))
+    F2[3]= COUP*denom*( (F1[3]*( (P2[0]*( -1j*V3[0]+1j*V3[3]))+( (P2[1]*( 1j*V3[1]-V3[2]))+( (P2[2]*( V3[1]+1j*V3[2]))+(P2[3]*( -1j*V3[0]+1j*V3[3]))))))+( (F1[2]*( (P2[0]*( -1j*V3[1]-V3[2]))+( (P2[1]*( 1j*V3[0]+1j*V3[3]))+( (P2[2]*( V3[0]+V3[3]))+(P2[3]*( -1j*V3[1]-V3[2]))))))+(M2*( (F1[1]*( -1j*V3[0]-1j*V3[3]))+(F1[0]*( 1j*V3[1]+V3[2]))))))
+    return F2
+    
+    
+def FFV1C1_2_1(F1, V3, COUP1,COUP2, M2, W2):
+
+
+    F2 = FFV1C1_1(F1, V3, COUP1, M2, W2)
+    tmp = FFV2C1_1(F1, V3, COUP2, M2, W2)
+    for i in range(4):
+        F2[i] += tmp[i]
+    return F2
+
+"""
+
+
+
+        FFV = UFOLorentz(name = 'FFV1',
+                 spins = [ 2, 2, 3 ],
+                 structure = 'Gamma(3,2,1)')        
+        builder = create_aloha.AbstractRoutineBuilder(FFV) 
+        builder.apply_conjugation()
+        amp = builder.compute_routine(1)
+        amp.add_combine(['FFV2'])
+        routine = amp.write(output_dir=None, language='Python')
+
+
+        split_solution = solution.split('\n')
+        split_routine = routine.split('\n')
+        self.assertEqual(split_solution, split_routine)
+        self.assertEqual(len(split_routine), len(split_solution))
+
        
     def test_pythonwriter_4_fermion(self):
         """ test that python writer works """
@@ -3838,20 +3884,20 @@
         
         
         solution ="""import wavefunctions
-def FFFF1C1_1(F2, F3, F4, COUP, M1, W1):
-    F1 = wavefunctions.WaveFunction(size=6)
-    F1[4] = F2[4]-F3[4]+F4[4]
-    F1[5] = F2[5]-F3[5]+F4[5]
-    P1 = [complex(F1[4]).real, \\
-             complex(F1[5]).real, \\
-             complex(F1[5]).imag, \\
-             complex(F1[4]).imag]
-    denom =1.0/(( (M1*( -M1+1j*W1))+( (P1[0]**2)-(P1[1]**2)-(P1[2]**2)-(P1[3]**2))))
-    F1[0]= COUP*denom*( (F2[3]*( (P1[1]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))+(P1[2]*( -(F4[0]*F3[0])-(F4[1]*F3[1])-(F4[2]*F3[2])-(F4[3]*F3[3])))))+( (F2[2]*( (P1[0]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))+(P1[3]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))))+(M1*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3]))*F2[0])))
-    F1[1]= COUP*denom*( (F2[2]*( (P1[1]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))+(P1[2]*( (F4[0]*F3[0])+(F4[1]*F3[1])+(F4[2]*F3[2])+(F4[3]*F3[3])))))+( (F2[3]*( (P1[0]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))+(P1[3]*( -1j*(F4[0]*F3[0])-1j*(F4[1]*F3[1])-1j*(F4[2]*F3[2])-1j*(F4[3]*F3[3])))))+(M1*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3]))*F2[1])))
-    F1[2]= COUP*denom*( (F2[1]*( (P1[1]*( -1j*(F4[0]*F3[0])-1j*(F4[1]*F3[1])-1j*(F4[2]*F3[2])-1j*(F4[3]*F3[3])))+(P1[2]*( (F4[0]*F3[0])+(F4[1]*F3[1])+(F4[2]*F3[2])+(F4[3]*F3[3])))))+( (F2[0]*( (P1[0]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))+(P1[3]*( -1j*(F4[0]*F3[0])-1j*(F4[1]*F3[1])-1j*(F4[2]*F3[2])-1j*(F4[3]*F3[3])))))+(M1*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3]))*F2[2])))
-    F1[3]= COUP*denom*( (F2[0]*( (P1[1]*( -1j*(F4[0]*F3[0])-1j*(F4[1]*F3[1])-1j*(F4[2]*F3[2])-1j*(F4[3]*F3[3])))+(P1[2]*( -(F4[0]*F3[0])-(F4[1]*F3[1])-(F4[2]*F3[2])-(F4[3]*F3[3])))))+( (F2[1]*( (P1[0]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))+(P1[3]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))))+(M1*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3]))*F2[3])))
-    return F1
+def FFFF1C1_1(F1, F3, F4, COUP, M2, W2):
+    F2 = wavefunctions.WaveFunction(size=6)
+    F2[4] = F1[4]-F3[4]+F4[4]
+    F2[5] = F1[5]-F3[5]+F4[5]
+    P2 = [complex(F2[4]).real, \\
+             complex(F2[5]).real, \\
+             complex(F2[5]).imag, \\
+             complex(F2[4]).imag]
+    denom =1.0/(( (M2*( -M2+1j*W2))+( (P2[0]**2)-(P2[1]**2)-(P2[2]**2)-(P2[3]**2))))
+    F2[0]= COUP*denom*( (F1[3]*( (P2[1]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))+(P2[2]*( -(F4[0]*F3[0])-(F4[1]*F3[1])-(F4[2]*F3[2])-(F4[3]*F3[3])))))+( (F1[2]*( (P2[0]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))+(P2[3]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))))+(M2*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3]))*F1[0])))
+    F2[1]= COUP*denom*( (F1[2]*( (P2[1]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))+(P2[2]*( (F4[0]*F3[0])+(F4[1]*F3[1])+(F4[2]*F3[2])+(F4[3]*F3[3])))))+( (F1[3]*( (P2[0]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))+(P2[3]*( -1j*(F4[0]*F3[0])-1j*(F4[1]*F3[1])-1j*(F4[2]*F3[2])-1j*(F4[3]*F3[3])))))+(M2*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3]))*F1[1])))
+    F2[2]= COUP*denom*( (F1[1]*( (P2[1]*( -1j*(F4[0]*F3[0])-1j*(F4[1]*F3[1])-1j*(F4[2]*F3[2])-1j*(F4[3]*F3[3])))+(P2[2]*( (F4[0]*F3[0])+(F4[1]*F3[1])+(F4[2]*F3[2])+(F4[3]*F3[3])))))+( (F1[0]*( (P2[0]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))+(P2[3]*( -1j*(F4[0]*F3[0])-1j*(F4[1]*F3[1])-1j*(F4[2]*F3[2])-1j*(F4[3]*F3[3])))))+(M2*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3]))*F1[2])))
+    F2[3]= COUP*denom*( (F1[0]*( (P2[1]*( -1j*(F4[0]*F3[0])-1j*(F4[1]*F3[1])-1j*(F4[2]*F3[2])-1j*(F4[3]*F3[3])))+(P2[2]*( -(F4[0]*F3[0])-(F4[1]*F3[1])-(F4[2]*F3[2])-(F4[3]*F3[3])))))+( (F1[1]*( (P2[0]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))+(P2[3]*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3])))))+(M2*( 1j*(F4[0]*F3[0])+1j*(F4[1]*F3[1])+1j*(F4[2]*F3[2])+1j*(F4[3]*F3[3]))*F1[3])))
+    return F2
     
     
 """
@@ -3941,6 +3987,181 @@
 
 
 
+    def test_fortranwriter_C(self):
+        """ test that python writer works """
+
+        solution = """ subroutine FFV1C1_1(F1, V3, COUP, M2, W2, F2)
+implicit none 
+double complex F1(*)
+double complex F2(*)
+double complex V3(*)
+double complex COUP
+double complex denom
+double precision M2, W2
+double precision P2(0:3)
+
+F2(5)= F1(5)+V3(5)
+F2(6)= F1(6)+V3(6)
+P2(0) =  dble(F2(5))
+P2(1) =  dble(F2(6))
+P2(2) =  dimag(F2(6))
+P2(3) =  dimag(F2(5))
+
+denom =1d0/(( (M2*( -M2+(0, 1)*W2))+( (P2(0)**2)-(P2(1)**2)-(P2(2)**2)-(P2(3)**2))))
+F2(1)= COUP*denom*( (F1(2)*( (P2(0)*( (0, 1)*V3(2)-V3(3)))+( (P2(1)*( (0, -1)*V3(1)+(0, -1)*V3(4)))+( (P2(2)*( V3(1)+V3(4)))+(P2(3)*( (0, 1)*V3(2)-V3(3)))))))+( (F1(1)*( (P2(0)*( (0, -1)*V3(1)+(0, 1)*V3(4)))+( (P2(1)*( (0, 1)*V3(2)+V3(3)))+( (P2(2)*( -V3(2)+(0, 1)*V3(3)))+(P2(3)*( (0, -1)*V3(1)+(0, 1)*V3(4)))))))+(M2*( (F1(4)*( (0, -1)*V3(2)+V3(3)))+(F1(3)*( (0, -1)*V3(1)+(0, -1)*V3(4)))))))
+F2(2)= COUP*denom*( (F1(2)*( (P2(0)*( (0, -1)*V3(1)+(0, -1)*V3(4)))+( (P2(1)*( (0, 1)*V3(2)-V3(3)))+( (P2(2)*( V3(2)+(0, 1)*V3(3)))+(P2(3)*( (0, 1)*V3(1)+(0, 1)*V3(4)))))))+( (F1(1)*( (P2(0)*( (0, 1)*V3(2)+V3(3)))+( (P2(1)*( (0, -1)*V3(1)+(0, 1)*V3(4)))+( (P2(2)*( -V3(1)+V3(4)))+(P2(3)*( (0, -1)*V3(2)-V3(3)))))))+(M2*( (F1(4)*( (0, -1)*V3(1)+(0, 1)*V3(4)))+(F1(3)*( (0, -1)*V3(2)-V3(3)))))))
+F2(3)= COUP*denom*( (F1(4)*( (P2(0)*( (0, -1)*V3(2)+V3(3)))+( (P2(1)*( (0, 1)*V3(1)+(0, -1)*V3(4)))+( (P2(2)*( -V3(1)+V3(4)))+(P2(3)*( (0, 1)*V3(2)-V3(3)))))))+( (F1(3)*( (P2(0)*( (0, -1)*V3(1)+(0, -1)*V3(4)))+( (P2(1)*( (0, 1)*V3(2)+V3(3)))+( (P2(2)*( -V3(2)+(0, 1)*V3(3)))+(P2(3)*( (0, 1)*V3(1)+(0, 1)*V3(4)))))))+(M2*( (F1(2)*( (0, 1)*V3(2)-V3(3)))+(F1(1)*( (0, -1)*V3(1)+(0, 1)*V3(4)))))))
+F2(4)= COUP*denom*( (F1(4)*( (P2(0)*( (0, -1)*V3(1)+(0, 1)*V3(4)))+( (P2(1)*( (0, 1)*V3(2)-V3(3)))+( (P2(2)*( V3(2)+(0, 1)*V3(3)))+(P2(3)*( (0, -1)*V3(1)+(0, 1)*V3(4)))))))+( (F1(3)*( (P2(0)*( (0, -1)*V3(2)-V3(3)))+( (P2(1)*( (0, 1)*V3(1)+(0, 1)*V3(4)))+( (P2(2)*( V3(1)+V3(4)))+(P2(3)*( (0, -1)*V3(2)-V3(3)))))))+(M2*( (F1(2)*( (0, -1)*V3(1)+(0, -1)*V3(4)))+(F1(1)*( (0, 1)*V3(2)+V3(3)))))))
+end
+
+
+"""        
+        FFV = UFOLorentz(name = 'FFV1',
+                 spins = [ 2, 2, 3 ],
+                 structure = 'Gamma(3,2,1)')        
+        builder = create_aloha.AbstractRoutineBuilder(FFV)
+        builder.apply_conjugation()
+        amp = builder.compute_routine(1)
+        routine = amp.write(output_dir=None, language='Fortran')
+        
+        split_solution = solution.split('\n')
+        split_routine = routine.split('\n')
+        self.assertEqual(split_solution, split_routine)
+        self.assertEqual(len(split_routine), len(split_solution))
+
+        solution=""" subroutine FFV1C1_2(F2, V3, COUP, M1, W1, F1)
+implicit none 
+double complex F1(*)
+double complex F2(*)
+double complex V3(*)
+double complex COUP
+double complex denom
+double precision M1, W1
+double precision P1(0:3)
+
+F1(5)= F2(5)-V3(5)
+F1(6)= F2(6)-V3(6)
+P1(0) =  dble(F1(5))
+P1(1) =  dble(F1(6))
+P1(2) =  dimag(F1(6))
+P1(3) =  dimag(F1(5))
+
+denom =1d0/(( (M1*( -M1+(0, 1)*W1))+( (P1(0)**2)-(P1(1)**2)-(P1(2)**2)-(P1(3)**2))))
+F1(1)= COUP*denom*( (F2(1)*( (V3(2)*( (0, 1)*P1(1)+P1(2)))+( (V3(3)*( -P1(1)+(0, 1)*P1(2)))+( (V3(1)*( (0, -1)*P1(0)+(0, 1)*P1(3)))+(V3(4)*( (0, -1)*P1(0)+(0, 1)*P1(3)))))))+( (F2(2)*( (V3(1)*( (0, 1)*P1(1)+P1(2)))+( (V3(4)*( (0, -1)*P1(1)-P1(2)))+( (V3(2)*( (0, -1)*P1(0)+(0, 1)*P1(3)))+(V3(3)*( -P1(0)+P1(3)))))))+(M1*( (F2(3)*( (0, -1)*V3(1)+(0, 1)*V3(4)))+(F2(4)*( (0, 1)*V3(2)+V3(3)))))))
+F1(2)= COUP*denom*( (F2(1)*( (V3(2)*( (0, -1)*P1(0)+(0, -1)*P1(3)))+( (V3(3)*( P1(0)+P1(3)))+( (V3(1)*( (0, 1)*P1(1)-P1(2)))+(V3(4)*( (0, 1)*P1(1)-P1(2)))))))+( (F2(2)*( (V3(1)*( (0, -1)*P1(0)+(0, -1)*P1(3)))+( (V3(4)*( (0, 1)*P1(0)+(0, 1)*P1(3)))+( (V3(2)*( (0, 1)*P1(1)-P1(2)))+(V3(3)*( P1(1)+(0, 1)*P1(2)))))))+(M1*( (F2(3)*( (0, 1)*V3(2)-V3(3)))+(F2(4)*( (0, -1)*V3(1)+(0, -1)*V3(4)))))))
+F1(3)= COUP*denom*( (F2(3)*( (V3(2)*( (0, 1)*P1(1)+P1(2)))+( (V3(3)*( -P1(1)+(0, 1)*P1(2)))+( (V3(1)*( (0, -1)*P1(0)+(0, -1)*P1(3)))+(V3(4)*( (0, 1)*P1(0)+(0, 1)*P1(3)))))))+( (F2(4)*( (V3(1)*( (0, -1)*P1(1)-P1(2)))+( (V3(4)*( (0, -1)*P1(1)-P1(2)))+( (V3(2)*( (0, 1)*P1(0)+(0, 1)*P1(3)))+(V3(3)*( P1(0)+P1(3)))))))+(M1*( (F2(1)*( (0, -1)*V3(1)+(0, -1)*V3(4)))+(F2(2)*( (0, -1)*V3(2)-V3(3)))))))
+F1(4)= COUP*denom*( (F2(3)*( (V3(2)*( (0, 1)*P1(0)+(0, -1)*P1(3)))+( (V3(3)*( -P1(0)+P1(3)))+( (V3(1)*( (0, -1)*P1(1)+P1(2)))+(V3(4)*( (0, 1)*P1(1)-P1(2)))))))+( (F2(4)*( (V3(1)*( (0, -1)*P1(0)+(0, 1)*P1(3)))+( (V3(4)*( (0, -1)*P1(0)+(0, 1)*P1(3)))+( (V3(2)*( (0, 1)*P1(1)-P1(2)))+(V3(3)*( P1(1)+(0, 1)*P1(2)))))))+(M1*( (F2(1)*( (0, -1)*V3(2)+V3(3)))+(F2(2)*( (0, -1)*V3(1)+(0, 1)*V3(4)))))))
+end
+
+
+"""
+        amp = builder.compute_routine(2)
+        
+        routine = amp.write(output_dir=None, language='Fortran')
+        
+        split_solution = solution.split('\n')
+        split_routine = routine.split('\n')
+        self.assertEqual(split_solution, split_routine)
+        self.assertEqual(len(split_routine), len(split_solution))
+
+    def test_Cppwriter_C(self):
+        """ test that python writer works """
+
+ 
+        solution_h = """#ifndef FFV1C1_1_guard
+#define FFV1C1_1_guard
+#include <complex>
+using namespace std;
+
+void FFV1C1_1(complex<double> F1[],complex<double> V3[],complex<double> COUP, double M2, double W2, complex<double>F2[]);
+
+#endif
+
+"""
+        solution_c = """#include "FFV1C1_1.h"
+
+void FFV1C1_1(complex<double> F1[],complex<double> V3[],complex<double> COUP, double M2, double W2, complex<double>F2[]){
+complex<double> denom;
+double P2[4];
+F2[4]= F1[4]+V3[4];
+F2[5]= F1[5]+V3[5];
+P2[0] = F2[4].real();
+P2[1] = F2[5].real();
+P2[2] = F2[5].imag();
+P2[3] = F2[4].imag();
+denom =1./(( (M2*( -M2+complex<double>(0., 1.)*W2))+( (pow(P2[0],2))-(pow(P2[1],2))-(pow(P2[2],2))-(pow(P2[3],2)))));
+F2[0]= COUP*denom*( (F1[1]*( (P2[0]*( complex<double>(0., 1.)*V3[1]-V3[2]))+( (P2[1]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., -1.)*V3[3]))+( (P2[2]*( V3[0]+V3[3]))+(P2[3]*( complex<double>(0., 1.)*V3[1]-V3[2]))))))+( (F1[0]*( (P2[0]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., 1.)*V3[3]))+( (P2[1]*( complex<double>(0., 1.)*V3[1]+V3[2]))+( (P2[2]*( -V3[1]+complex<double>(0., 1.)*V3[2]))+(P2[3]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., 1.)*V3[3]))))))+(M2*( (F1[3]*( complex<double>(0., -1.)*V3[1]+V3[2]))+(F1[2]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., -1.)*V3[3]))))));
+F2[1]= COUP*denom*( (F1[1]*( (P2[0]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., -1.)*V3[3]))+( (P2[1]*( complex<double>(0., 1.)*V3[1]-V3[2]))+( (P2[2]*( V3[1]+complex<double>(0., 1.)*V3[2]))+(P2[3]*( complex<double>(0., 1.)*V3[0]+complex<double>(0., 1.)*V3[3]))))))+( (F1[0]*( (P2[0]*( complex<double>(0., 1.)*V3[1]+V3[2]))+( (P2[1]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., 1.)*V3[3]))+( (P2[2]*( -V3[0]+V3[3]))+(P2[3]*( complex<double>(0., -1.)*V3[1]-V3[2]))))))+(M2*( (F1[3]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., 1.)*V3[3]))+(F1[2]*( complex<double>(0., -1.)*V3[1]-V3[2]))))));
+F2[2]= COUP*denom*( (F1[3]*( (P2[0]*( complex<double>(0., -1.)*V3[1]+V3[2]))+( (P2[1]*( complex<double>(0., 1.)*V3[0]+complex<double>(0., -1.)*V3[3]))+( (P2[2]*( -V3[0]+V3[3]))+(P2[3]*( complex<double>(0., 1.)*V3[1]-V3[2]))))))+( (F1[2]*( (P2[0]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., -1.)*V3[3]))+( (P2[1]*( complex<double>(0., 1.)*V3[1]+V3[2]))+( (P2[2]*( -V3[1]+complex<double>(0., 1.)*V3[2]))+(P2[3]*( complex<double>(0., 1.)*V3[0]+complex<double>(0., 1.)*V3[3]))))))+(M2*( (F1[1]*( complex<double>(0., 1.)*V3[1]-V3[2]))+(F1[0]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., 1.)*V3[3]))))));
+F2[3]= COUP*denom*( (F1[3]*( (P2[0]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., 1.)*V3[3]))+( (P2[1]*( complex<double>(0., 1.)*V3[1]-V3[2]))+( (P2[2]*( V3[1]+complex<double>(0., 1.)*V3[2]))+(P2[3]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., 1.)*V3[3]))))))+( (F1[2]*( (P2[0]*( complex<double>(0., -1.)*V3[1]-V3[2]))+( (P2[1]*( complex<double>(0., 1.)*V3[0]+complex<double>(0., 1.)*V3[3]))+( (P2[2]*( V3[0]+V3[3]))+(P2[3]*( complex<double>(0., -1.)*V3[1]-V3[2]))))))+(M2*( (F1[1]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., -1.)*V3[3]))+(F1[0]*( complex<double>(0., 1.)*V3[1]+V3[2]))))));
+}
+
+"""
+        
+        FFV = UFOLorentz(name = 'FFV1',
+                 spins = [ 2, 2, 3 ],
+                 structure = 'Gamma(3,2,1)')        
+        builder = create_aloha.AbstractRoutineBuilder(FFV)
+        builder.apply_conjugation()
+        amp = builder.compute_routine(1)
+        routine = amp.write(output_dir=None, language='CPP')
+        
+        split_solution = solution_h.split('\n')
+        split_routine = routine[0].split('\n')
+        self.assertEqual(split_solution, split_routine)
+        self.assertEqual(len(split_routine), len(split_solution))
+
+        split_solution = solution_c.split('\n')
+        split_routine = routine[1].split('\n')
+        self.assertEqual(split_solution, split_routine)
+        self.assertEqual(len(split_routine), len(split_solution))
+
+        solution_h = """#ifndef FFV1C1_2_guard
+#define FFV1C1_2_guard
+#include <complex>
+using namespace std;
+
+void FFV1C1_2(complex<double> F2[],complex<double> V3[],complex<double> COUP, double M1, double W1, complex<double>F1[]);
+
+#endif
+
+"""
+
+        solution_c = """#include "FFV1C1_2.h"
+
+void FFV1C1_2(complex<double> F2[],complex<double> V3[],complex<double> COUP, double M1, double W1, complex<double>F1[]){
+complex<double> denom;
+double P1[4];
+F1[4]= F2[4]-V3[4];
+F1[5]= F2[5]-V3[5];
+P1[0] = F1[4].real();
+P1[1] = F1[5].real();
+P1[2] = F1[5].imag();
+P1[3] = F1[4].imag();
+denom =1./(( (M1*( -M1+complex<double>(0., 1.)*W1))+( (pow(P1[0],2))-(pow(P1[1],2))-(pow(P1[2],2))-(pow(P1[3],2)))));
+F1[0]= COUP*denom*( (F2[0]*( (V3[1]*( complex<double>(0., 1.)*P1[1]+P1[2]))+( (V3[2]*( -P1[1]+complex<double>(0., 1.)*P1[2]))+( (V3[0]*( complex<double>(0., -1.)*P1[0]+complex<double>(0., 1.)*P1[3]))+(V3[3]*( complex<double>(0., -1.)*P1[0]+complex<double>(0., 1.)*P1[3]))))))+( (F2[1]*( (V3[0]*( complex<double>(0., 1.)*P1[1]+P1[2]))+( (V3[3]*( complex<double>(0., -1.)*P1[1]-P1[2]))+( (V3[1]*( complex<double>(0., -1.)*P1[0]+complex<double>(0., 1.)*P1[3]))+(V3[2]*( -P1[0]+P1[3]))))))+(M1*( (F2[2]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., 1.)*V3[3]))+(F2[3]*( complex<double>(0., 1.)*V3[1]+V3[2]))))));
+F1[1]= COUP*denom*( (F2[0]*( (V3[1]*( complex<double>(0., -1.)*P1[0]+complex<double>(0., -1.)*P1[3]))+( (V3[2]*( P1[0]+P1[3]))+( (V3[0]*( complex<double>(0., 1.)*P1[1]-P1[2]))+(V3[3]*( complex<double>(0., 1.)*P1[1]-P1[2]))))))+( (F2[1]*( (V3[0]*( complex<double>(0., -1.)*P1[0]+complex<double>(0., -1.)*P1[3]))+( (V3[3]*( complex<double>(0., 1.)*P1[0]+complex<double>(0., 1.)*P1[3]))+( (V3[1]*( complex<double>(0., 1.)*P1[1]-P1[2]))+(V3[2]*( P1[1]+complex<double>(0., 1.)*P1[2]))))))+(M1*( (F2[2]*( complex<double>(0., 1.)*V3[1]-V3[2]))+(F2[3]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., -1.)*V3[3]))))));
+F1[2]= COUP*denom*( (F2[2]*( (V3[1]*( complex<double>(0., 1.)*P1[1]+P1[2]))+( (V3[2]*( -P1[1]+complex<double>(0., 1.)*P1[2]))+( (V3[0]*( complex<double>(0., -1.)*P1[0]+complex<double>(0., -1.)*P1[3]))+(V3[3]*( complex<double>(0., 1.)*P1[0]+complex<double>(0., 1.)*P1[3]))))))+( (F2[3]*( (V3[0]*( complex<double>(0., -1.)*P1[1]-P1[2]))+( (V3[3]*( complex<double>(0., -1.)*P1[1]-P1[2]))+( (V3[1]*( complex<double>(0., 1.)*P1[0]+complex<double>(0., 1.)*P1[3]))+(V3[2]*( P1[0]+P1[3]))))))+(M1*( (F2[0]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., -1.)*V3[3]))+(F2[1]*( complex<double>(0., -1.)*V3[1]-V3[2]))))));
+F1[3]= COUP*denom*( (F2[2]*( (V3[1]*( complex<double>(0., 1.)*P1[0]+complex<double>(0., -1.)*P1[3]))+( (V3[2]*( -P1[0]+P1[3]))+( (V3[0]*( complex<double>(0., -1.)*P1[1]+P1[2]))+(V3[3]*( complex<double>(0., 1.)*P1[1]-P1[2]))))))+( (F2[3]*( (V3[0]*( complex<double>(0., -1.)*P1[0]+complex<double>(0., 1.)*P1[3]))+( (V3[3]*( complex<double>(0., -1.)*P1[0]+complex<double>(0., 1.)*P1[3]))+( (V3[1]*( complex<double>(0., 1.)*P1[1]-P1[2]))+(V3[2]*( P1[1]+complex<double>(0., 1.)*P1[2]))))))+(M1*( (F2[0]*( complex<double>(0., -1.)*V3[1]+V3[2]))+(F2[1]*( complex<double>(0., -1.)*V3[0]+complex<double>(0., 1.)*V3[3]))))));
+}
+
+"""
+
+        amp = builder.compute_routine(2)
+        
+        routine = amp.write(output_dir=None, language='CPP')
+        
+        split_solution = solution_h.split('\n')
+        split_routine = routine[0].split('\n')
+        self.assertEqual(split_solution, split_routine)
+        self.assertEqual(len(split_routine), len(split_solution))
+
+        split_solution = solution_c.split('\n')
+        split_routine = routine[1].split('\n')
+        self.assertEqual(split_solution, split_routine)
+        self.assertEqual(len(split_routine), len(split_solution))
+        
+