MadGraph5_aMC@NLO

Merge lp:~maddevelopers/mg5amcnlo/ALOHA_spin2 into lp:~madteam/mg5amcnlo/trunk

  1. ALOHA_spin2
  2. Merge into trunk
Proposed by Olivier Mattelaer
Status: Rejected
Rejected by: Olivier Mattelaer
Proposed branch: lp:~maddevelopers/mg5amcnlo/ALOHA_spin2
Merge into: lp:~madteam/mg5amcnlo/trunk
Diff against target: 1524 lines (+975/-267)
11 files modified
aloha/aloha_lib.py (+4/-1)
aloha/aloha_writers.py (+21/-11)
aloha/create_aloha.py (+40/-28)
aloha/template_files/wavefunctions.py (+341/-210)
madgraph/core/base_objects.py (+10/-4)
madgraph/interface/cmd_interface.py (+22/-2)
madgraph/iolibs/export_v4.py (+1/-0)
madgraph/various/process_checks.py (+1/-1)
tests/unit_tests/core/test_base_objects.py (+4/-4)
tests/unit_tests/various/test_aloha.py (+44/-2)
tests/unit_tests/various/test_rambo.py (+487/-4)
To merge this branch: bzr merge lp:~maddevelopers/mg5amcnlo/ALOHA_spin2
Reviewer Review Type Date Requested Status
Johan Alwall Pending
Review via email: mp+46484@code.launchpad.net

Description of the change

Since a bug was discovered in the ALOHA way to make analytical computation and even if it seems quite hard to face this problem. I suggest to merge this branch in the trunk as soon as possible.

This branch update ALOHA in order to fully support massive spin2.

To post a comment you must log in.

Unmerged revisions

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1=== modified file 'aloha/aloha_lib.py'
2--- aloha/aloha_lib.py 2010-11-18 14:21:11 +0000
3+++ aloha/aloha_lib.py 2011-01-17 15:18:53 +0000
4@@ -591,6 +591,7 @@
5 self.prefactor *= obj.prefactor
6 if obj in self:
7 index = self.index(obj)
8+ self[index] = self[index].copy()
9 self[index].power += 1
10 else:
11 obj.prefactor = 1
12@@ -1292,6 +1293,7 @@
13 def __mul__(self, obj):
14 """multiplication performing directly the einstein/spin sommation.
15 """
16+
17 if not hasattr(obj, 'vartype') or not self.vartype:
18 out = LorentzObjectRepresentation({}, self.lorentz_ind, self.spin_ind)
19 for ind in out.listindices():
20@@ -1388,10 +1390,11 @@
21 factor = self.get_rep(self_ind)
22 factor *= obj.get_rep(obj_ind)
23
24+
25 if factor:
26 #compute the prefactor due to the lorentz contraction
27 factor *= (-1) ** (len(l_value) - l_value.count(0))
28- out += factor
29+ out += factor
30 return out
31
32 def combine_indices(self, l_dict, s_dict):
33
34=== modified file 'aloha/aloha_writers.py'
35--- aloha/aloha_writers.py 2010-09-12 01:46:50 +0000
36+++ aloha/aloha_writers.py 2011-01-17 15:18:53 +0000
37@@ -34,8 +34,23 @@
38 #prepare the necessary object
39 self.collect_variables() # Look for the different variables
40 self.make_all_lists() # Compute the expression for the call ordering
41- #the definition of objects,...
42-
43+ #the definition of objects,..
44+
45+ def pass_to_HELAS(self, indices, start=0):
46+ """find the Fortran HELAS position for the list of index"""
47+
48+
49+ if len(indices) == 1:
50+ return indices[0] + start
51+
52+ ind_name = self.obj.numerator.lorentz_ind
53+ if ind_name == ['I3', 'I2']:
54+ return 4 * indices[1] + indices[0] + start
55+ elif len(indices) == 2:
56+ return 4 * indices[0] + indices[1] + start
57+ else:
58+ raise Exception
59+
60 def collect_variables(self):
61 """Collects Momenta,Mass,Width into lists"""
62
63@@ -444,6 +459,7 @@
64 out = '%.9f' % number
65 out = self.zero_pattern.sub('', out)
66 return out
67+
68
69 def define_expression(self):
70 OutString = ''
71@@ -466,15 +482,13 @@
72 string = re.sub('\((?P<num>[+-]*[0-9])\+(?P<num2>[+-][0-9])[Jj]\)\.', '(\g<num>d0,\g<num2>d0)', string)
73 string = re.sub('(?P<num>[0-9])[Jj]\.', '\g<num>*(0d0,1d0)', string)
74 OutString = OutString + string + '\n'
75- counter = 1
76 for ind in numerator.listindices():
77- string = '%s(%d)= C*denom*' % (OffShellParticle, counter)
78+ string = '%s(%d)= C*denom*' % (OffShellParticle, self.pass_to_HELAS(ind, start=1))
79 string += self.write_obj(numerator.get_rep(ind))
80 string = string.replace('+-', '-')
81 string = re.sub('\((?P<num>[+-][0-9])\+(?P<num2>[+-][0-9])[Jj]\)\.', '(\g<num>d0,\g<num2>d0)', string)
82 string = re.sub('(?P<num>[0-9])[Jj]\.', '\g<num>*(0d0,1d0)', string)
83 OutString = OutString + string + '\n'
84- counter += 1
85 return OutString
86
87 def define_symmetry(self, new_nb):
88@@ -674,13 +688,11 @@
89 string = 'denom =' + '1./(' + denom + ')'
90 string = string.replace('+-', '-')
91 OutString = OutString + string + ';\n'
92- counter = 0
93 for ind in numerator.listindices():
94- string = '%s[%d]= C*denom*' % (OffShellParticle, counter)
95+ string = '%s[%d]= C*denom*' % (OffShellParticle, self.pass_to_HELAS(ind))
96 string += self.write_obj(numerator.get_rep(ind))
97 string = string.replace('+-', '-')
98 OutString = OutString + string + ';\n'
99- counter += 1
100 OutString = re.sub('(?P<variable>[A-Za-z]+[0-9]\[*[0-9]*\]*)\*\*(?P<num>[0-9])','pow(\g<variable>,\g<num>)',OutString)
101 return OutString
102
103@@ -824,13 +836,11 @@
104 string = 'denom =' + '1.0/(' + denom + ')'
105 string = string.replace('+-', '-')
106 OutString += string + '\n'
107- counter = 0
108 for ind in numerator.listindices():
109- string = '%s[%d]= C*denom*' % (self.outname, counter)
110+ string = '%s[%d]= C*denom*' % (self.outname, self.pass_to_HELAS(ind))
111 string += self.write_obj(numerator.get_rep(ind))
112 string = string.replace('+-', '-')
113 OutString += string + '\n'
114- counter += 1
115 return OutString
116
117 def define_foot(self):
118
119=== modified file 'aloha/create_aloha.py'
120--- aloha/create_aloha.py 2010-10-11 20:32:01 +0000
121+++ aloha/create_aloha.py 2011-01-17 15:18:53 +0000
122@@ -27,6 +27,7 @@
123 sys.path.append(root_path)
124 from aloha.aloha_object import *
125 import aloha.aloha_writers as aloha_writers
126+import aloha.aloha_lib as aloha_lib
127 try:
128 import madgraph.iolibs.files as files
129 except:
130@@ -36,6 +37,7 @@
131 logger = logging.getLogger('ALOHA')
132
133 _conjugate_gap = 50
134+_spin2_mult = 1000
135
136 class ALOHAERROR(Exception): pass
137
138@@ -91,10 +93,10 @@
139 self.routine_kernel = None
140
141
142- def compute_routine(self, mode):
143+ def compute_routine(self, mode, factorize=True):
144 """compute the expression and return it"""
145 self.outgoing = mode
146- self.expr = self.compute_aloha_high_kernel(mode)
147+ self.expr = self.compute_aloha_high_kernel(mode, factorize)
148 return self.define_simple_output()
149
150 def define_all_conjugate_builder(self, pair_list):
151@@ -148,7 +150,7 @@
152 return AbstractRoutine(self.expr, self.outgoing, self.spins, self.name, \
153 infostr)
154
155- def compute_aloha_high_kernel(self, mode):
156+ def compute_aloha_high_kernel(self, mode, factorize=True):
157 """compute the abstract routine associate to this mode """
158 # reset tag for particles
159 aloha_lib.USE_TAG=set()
160@@ -182,9 +184,9 @@
161 id += _conjugate_gap
162 nb_spinor += 1
163 if nb_spinor %2:
164- lorentz *= SpinorPropagator(id, 'I2', i + 1)
165+ lorentz *= SpinorPropagator(id, 'I2', self.outgoing)
166 else:
167- lorentz *= SpinorPropagator('I2', id, i + 1)
168+ lorentz *= SpinorPropagator('I2', id, self.outgoing)
169 elif spin == 3 :
170 lorentz *= VectorPropagator(id, 'I2', id)
171 elif spin == 5 :
172@@ -199,13 +201,15 @@
173 elif spin == 2:
174 # shift the tag if we multiply by C matrices
175 if (id+1) // 2 in self.conjg:
176- id += _conjugate_gap
177+ spin_id = id + _conjugate_gap
178+ else:
179+ spin_id = id
180 nb_spinor += 1
181- lorentz *= Spinor(id, i + 1)
182+ lorentz *= Spinor(spin_id, id)
183 elif spin == 3:
184 lorentz *= Vector(id, id)
185 elif spin == 5:
186- lorentz *= Spin2(10*id+1, 10*id+2, 'I2', 'I3', id)
187+ lorentz *= Spin2(1 * _spin2_mult + id, 2 * _spin2_mult + id, id)
188 else:
189 raise self.AbstractALOHAError(
190 'The spin value %s is not supported yet' % spin)
191@@ -221,15 +225,19 @@
192
193 #lorentz = lorentz.simplify()
194 lorentz = lorentz.expand()
195-
196- if self.spins[self.outgoing-1] == 5:
197+ if self.outgoing and self.spins[self.outgoing-1] == 5:
198 if not self.aloha_lib:
199 AbstractRoutineBuilder.load_library()
200 lorentz *= self.aloha_lib[('Spin2Prop', id)]
201+ aloha_lib.USE_TAG.add('OM%d' % id)
202+ aloha_lib.USE_TAG.add('P%d' % id)
203+
204+
205
206
207 lorentz = lorentz.simplify()
208- lorentz = lorentz.factorize()
209+ if factorize:
210+ lorentz = lorentz.factorize()
211
212 lorentz.tag = set(aloha_lib.USE_TAG)
213 #raise
214@@ -269,8 +277,12 @@
215 @classmethod
216 def load_library(cls):
217 # load the library
218- fsock = open(os.path.join(aloha_path, 'ALOHALib.pkl'), 'r')
219- cls.aloha_lib = cPickle.load(fsock)
220+ try:
221+ fsock = open(os.path.join(aloha_path, 'ALOHALib.pkl'), 'r')
222+ except IOError:
223+ cls.aloha_lib = create_library()
224+ else:
225+ cls.aloha_lib = cPickle.load(fsock)
226
227
228 class AbstractALOHAModel(dict):
229@@ -638,22 +650,22 @@
230 lib = {} # key: (name, part_nb, special) -> object
231 for i in range(1, 10):
232 logger.info('step %s/9' % i)
233- lib[('Scalar', i)] = create( Scalar(i) )
234- lib[('ScalarProp', i)] = complex(0,1)
235- lib[('Denom', i )] = create( DenominatorPropagator(i) )
236- lib[('Spinor', i )] = create( Spinor(i, i) )
237- lib[('SpinorProp', i, 0)] = create( SpinorPropagator(i, 'I2', i) )
238- lib[('SpinorProp', i, 1)] = create( SpinorPropagator('I2', i, i) )
239- lib[('Vector', i)] = create( Vector(i+1, i+1) )
240- lib[('VectorProp', i)] = create( VectorPropagator(i,'I2', i) )
241- lib[('Spin2', i )] = create( Spin2(10*i+1, 10*i+2, i) )
242- lib[('Spin2Prop',i)] = create( Spin2Propagator(10*i+1, \
243- 10*i+2,'I2','I3', i) )
244- logger.info('writing')
245- fsock = open('./ALOHALib.pkl','wb')
246+ #lib[('Scalar', i)] = create( Scalar(i) )
247+ #lib[('ScalarProp', i)] = complex(0,1)
248+ #lib[('Denom', i )] = create( DenominatorPropagator(i) )
249+ #lib[('Spinor', i )] = create( Spinor(i, i) )
250+ #lib[('SpinorProp', i, 0)] = create( SpinorPropagator(i, 'I2', i) )
251+ #lib[('SpinorProp', i, 1)] = create( SpinorPropagator('I2', i, i) )
252+ #lib[('Vector', i)] = create( Vector(i+1, i+1) )
253+ #lib[('VectorProp', i)] = create( VectorPropagator(i,'I2', i) )
254+ #lib[('Spin2', i )] = create( Spin2(10*i+1, 10*i+2, i) )
255+ lib[('Spin2Prop',i)] = create( Spin2Propagator(_spin2_mult + i, \
256+ 2 * _spin2_mult + i,'I2','I3', i) )
257+ logger.info('writing Spin2 lib')
258+ fsock = open(os.path.join(aloha_path, 'ALOHALib.pkl'),'wb')
259 cPickle.dump(lib, fsock, -1)
260- logger.info('done')
261-
262+ return lib
263+
264 if '__main__' == __name__:
265 logging.basicConfig(level=0)
266 #create_library()
267
268=== modified file 'aloha/template_files/wavefunctions.py'
269--- aloha/template_files/wavefunctions.py 2010-09-05 16:14:23 +0000
270+++ aloha/template_files/wavefunctions.py 2011-01-17 15:18:53 +0000
271@@ -1,224 +1,355 @@
272-import cmath
273+from __future__ import division
274 import math
275+from math import sqrt, pow
276
277 class WaveFunction(list):
278- """a objet for a WaveFunction"""
279-
280- spin_to_size={0:1,
281- 1:3,
282- 2:6,
283- 3:6}
284-
285- def __init__(self, spin= None, size=None):
286- """Init the list with zero value"""
287-
288- if spin:
289- size = self.spin_to_size[spin]
290- list.__init__(self, [0]*size)
291-
292+ """a objet for a WaveFunction"""
293+
294+ spin_to_size={0:1,
295+ 1:3,
296+ 2:6,
297+ 3:6,
298+ 5:18}
299+
300+ def __init__(self, spin= None, size=None):
301+ """Init the list with zero value"""
302+
303+ if spin:
304+ size = self.spin_to_size[spin]
305+ list.__init__(self, [0]*size)
306+
307
308 def ixxxxx(p,fmass,nhel,nsf):
309- """Defines an inflow fermion."""
310-
311- fi = WaveFunction(2)
312-
313- fi[4] = complex(p[0]*nsf,p[3]*nsf)
314- fi[5] = complex(p[1]*nsf,p[2]*nsf)
315-
316- nh = nhel*nsf
317-
318- if (fmass != 0.):
319- pp = min(p[0],math.sqrt(math.pow(p[1],2)+math.pow(p[2],2)+math.pow(p[3],2)))
320- if (pp == 0.):
321- sqm = math.sqrt(abs(fmass))
322- ip = (1+nh)/2
323- im = (1-nh)/2
324-
325- fi[0] = ip*sqm
326- fi[1] = im*nsf*sqm
327- fi[2] = ip*nsf*sqm
328- fi[3] = im*sqm
329-
330- else:
331- sf = [(1+nsf+(1-nsf)*nh)*0.5,(1+nsf-(1-nsf)*nh)*0.5]
332- omega = [math.sqrt(p[0]+pp),fmass/(math.sqrt(p[0]+pp))]
333- ip = (1+nh)/2
334- im = (1-nh)/2
335- sfomeg = [sf[0]*omega[ip],sf[1]*omega[im]]
336- pp3 = max(pp+p[3],0.)
337- if (pp3 == 0.):
338- chi1 = complex(-nh,0.)
339- else:
340- chi1 = complex(nh*p[1]/math.sqrt(2.*pp*pp3),\
341- p[2]/math.sqrt(2.*pp*pp3))
342- chi = [complex(math.sqrt(pp3*0.5/pp)),chi1]
343-
344- fi[0] = sfomeg[0]*chi[im]
345- fi[1] = sfomeg[0]*chi[ip]
346- fi[2] = sfomeg[1]*chi[im]
347- fi[3] = sfomeg[1]*chi[ip]
348-
349- else:
350- sqp0p3 = math.sqrt(max(p[0]+p[3],0.))*nsf
351- if (sqp0p3 == 0.):
352- chi1 = complex(-nhel*math.sqrt(2.*p[0]),0.)
353- else:
354- chi1 = complex(nh*p[1]/sqp0p3,p[2]/sqp0p3)
355- chi = [complex(sqp0p3,0.),chi1]
356- if (nh == 1):
357- fi[0] = complex(0.,0.)
358- fi[1] = complex(0.,0.)
359- fi[2] = chi[0]
360- fi[3] = chi[1]
361- else:
362- fi[0] = chi[1]
363- fi[1] = chi[0]
364- fi[2] = complex(0.,0.)
365- fi[3] = complex(0.,0.)
366-
367- return fi
368+ """Defines an inflow fermion."""
369+
370+ fi = WaveFunction(2)
371+
372+ fi[4] = complex(p[0]*nsf,p[3]*nsf)
373+ fi[5] = complex(p[1]*nsf,p[2]*nsf)
374+
375+ nh = nhel*nsf
376+
377+ if (fmass != 0.):
378+ pp = min(p[0],sqrt(p[1]**2 + p[2]**2 + p[3]**2 ))
379+ if (pp == 0.):
380+ sqm = sqrt(abs(fmass))
381+ ip = (1+nh)/2
382+ im = (1-nh)/2
383+
384+ fi[0] = ip*sqm
385+ fi[1] = im*nsf*sqm
386+ fi[2] = ip*nsf*sqm
387+ fi[3] = im*sqm
388+
389+ else:
390+ sf = [(1+nsf+(1-nsf)*nh)*0.5,(1+nsf-(1-nsf)*nh)*0.5]
391+ omega = [sqrt(p[0]+pp),fmass/(sqrt(p[0]+pp))]
392+ ip = (1+nh)//2
393+ im = (1-nh)//2
394+ sfomeg = [sf[0]*omega[ip],sf[1]*omega[im]]
395+ pp3 = max(pp+p[3],0.)
396+ if (pp3 == 0.):
397+ chi1 = complex(-nh,0.)
398+ else:
399+ chi1 = complex(nh*p[1]/sqrt(2.*pp*pp3),\
400+ p[2]/sqrt(2.*pp*pp3))
401+ chi = [complex(sqrt(pp3*0.5/pp)),chi1]
402+
403+ fi[0] = sfomeg[0]*chi[im]
404+ fi[1] = sfomeg[0]*chi[ip]
405+ fi[2] = sfomeg[1]*chi[im]
406+ fi[3] = sfomeg[1]*chi[ip]
407+
408+ else:
409+ sqp0p3 = sqrt(max(p[0]+p[3],0.))*nsf
410+ if (sqp0p3 == 0.):
411+ chi1 = complex(-nhel*sqrt(2.*p[0]),0.)
412+ else:
413+ chi1 = complex(nh*p[1]/sqp0p3,p[2]/sqp0p3)
414+ chi = [complex(sqp0p3,0.),chi1]
415+ if (nh == 1):
416+ fi[0] = complex(0.,0.)
417+ fi[1] = complex(0.,0.)
418+ fi[2] = chi[0]
419+ fi[3] = chi[1]
420+ else:
421+ fi[0] = chi[1]
422+ fi[1] = chi[0]
423+ fi[2] = complex(0.,0.)
424+ fi[3] = complex(0.,0.)
425+
426+ return fi
427
428 def oxxxxx(p,fmass,nhel,nsf):
429- """ initialize an outgoing fermion"""
430-
431- fo = WaveFunction(2)
432-
433- fo[4] = complex(p[0]*nsf,p[3]*nsf)
434- fo[5] = complex(p[1]*nsf,p[2]*nsf)
435-
436- nh = nhel*nsf
437-
438- if (fmass != 0.):
439- pp = min(p[0],math.sqrt(math.pow(p[1],2)+math.pow(p[2],2)+math.pow(p[3],2)))
440- if (pp == 0.):
441- sqm = math.sqrt(abs(fmass))
442- ip = -((1+nh)/2)
443- im = (1-nh)/2
444-
445- fo[0] = im*sqm
446- fo[1] = ip*nsf*sqm
447- fo[2] = im*nsf*sqm
448- fo[3] = ip*sqm
449-
450- else:
451- sf = [(1+nsf+(1-nsf)*nh)*0.5,(1+nsf-(1-nsf)*nh)*0.5]
452- omega = [math.sqrt(p[0]+pp),fmass/(math.sqrt(p[0]+pp))]
453- ip = (1+nh)/2
454- im = (1-nh)/2
455- sfomeg = [sf[0]*omega[ip],sf[1]*omega[im]]
456- pp3 = max(pp+p[3],0.)
457- if (pp3 == 0.):
458- chi1 = complex(-nh,0.)
459- else:
460- chi1 = complex(nh*p[1]/math.sqrt(2.*pp*pp3),\
461- -p[2]/math.sqrt(2.*pp*pp3))
462- chi = [complex(math.sqrt(pp3*0.5/pp)),chi1]
463-
464- fo[0] = sfomeg[1]*chi[im]
465- fo[1] = sfomeg[1]*chi[ip]
466- fo[2] = sfomeg[0]*chi[im]
467- fo[3] = sfomeg[0]*chi[ip]
468-
469- else:
470- sqp0p3 = math.sqrt(max(p[0]+p[3],0.))*nsf
471- if (sqp0p3 == 0.):
472- chi1 = complex(-nhel*math.sqrt(2.*p[0]),0.)
473- else:
474- chi1 = complex(nh*p[1]/sqp0p3,-p[2]/sqp0p3)
475- chi = [complex(sqp0p3,0.),chi1]
476- if (nh == 1):
477- fo[0] = chi[0]
478- fo[1] = chi[1]
479- fo[2] = complex(0.,0.)
480- fo[3] = complex(0.,0.)
481- else:
482- fo[0] = complex(0.,0.)
483- fo[1] = complex(0.,0.)
484- fo[2] = chi[1]
485- fo[3] = chi[0]
486-
487- return fo
488+ """ initialize an outgoing fermion"""
489+
490+ fo = WaveFunction(2)
491+
492+ fo[4] = complex(p[0]*nsf,p[3]*nsf)
493+ fo[5] = complex(p[1]*nsf,p[2]*nsf)
494+
495+ nh = nhel*nsf
496+
497+ if (fmass != 0.):
498+ pp = min(p[0],sqrt(p[1]**2 + p[2]**2 + p[3]**2 ))
499+ if (pp == 0.):
500+ sqm = sqrt(abs(fmass))
501+ ip = -((1+nh)/2)
502+ im = (1-nh)/2
503+
504+ fo[0] = im*sqm
505+ fo[1] = ip*nsf*sqm
506+ fo[2] = im*nsf*sqm
507+ fo[3] = ip*sqm
508+
509+ else:
510+ sf = [(1+nsf+(1-nsf)*nh)*0.5,(1+nsf-(1-nsf)*nh)*0.5]
511+ omega = [sqrt(p[0]+pp),fmass/(sqrt(p[0]+pp))]
512+ ip = (1+nh)//2
513+ im = (1-nh)//2
514+ sfomeg = [sf[0]*omega[ip],sf[1]*omega[im]]
515+ pp3 = max(pp+p[3],0.)
516+ if (pp3 == 0.):
517+ chi1 = complex(-nh,0.)
518+ else:
519+ chi1 = complex(nh*p[1]/sqrt(2.*pp*pp3),\
520+ -p[2]/sqrt(2.*pp*pp3))
521+ chi = [complex(sqrt(pp3*0.5/pp)),chi1]
522+
523+ fo[0] = sfomeg[1]*chi[im]
524+ fo[1] = sfomeg[1]*chi[ip]
525+ fo[2] = sfomeg[0]*chi[im]
526+ fo[3] = sfomeg[0]*chi[ip]
527+
528+ else:
529+ sqp0p3 = sqrt(max(p[0]+p[3],0.))*nsf
530+ if (sqp0p3 == 0.):
531+ chi1 = complex(-nhel*sqrt(2.*p[0]),0.)
532+ else:
533+ chi1 = complex(nh*p[1]/sqp0p3,-p[2]/sqp0p3)
534+ chi = [complex(sqp0p3,0.),chi1]
535+ if (nh == 1):
536+ fo[0] = chi[0]
537+ fo[1] = chi[1]
538+ fo[2] = complex(0.,0.)
539+ fo[3] = complex(0.,0.)
540+ else:
541+ fo[0] = complex(0.,0.)
542+ fo[1] = complex(0.,0.)
543+ fo[2] = chi[1]
544+ fo[3] = chi[0]
545+
546+ return fo
547
548 def vxxxxx(p,vmass,nhel,nsv):
549- """ initialize a vector wavefunction. nhel=4 is for checking BRST"""
550-
551- vc = WaveFunction(3)
552-
553- sqh = math.sqrt(0.5)
554- nsvahl = nsv*abs(nhel)
555- pt2 = math.pow(p[1],2)+math.pow(p[2],2)
556- pp = min(p[0],math.sqrt(pt2+pow(p[3],2)))
557- pt = min(pp,math.sqrt(pt2))
558-
559- vc[4] = complex(p[0]*nsv,p[3]*nsv)
560- vc[5] = complex(p[1]*nsv,p[2]*nsv)
561-
562- if (nhel == 4):
563- if (vmass == 0.):
564- vc[0] = 1.
565- vc[1]=p[1]/p[0]
566- vc[2]=p[2]/p[0]
567- vc[3]=p[3]/p[0]
568- else:
569- vc[0] = p[0]/vmass
570- vc[1] = p[1]/vmass
571- vc[2] = p[2]/vmass
572- vc[3] = p[3]/vmass
573-
574- return vc
575-
576- if (vmass != 0.):
577- hel0 = 1.-abs(nhel)
578-
579- if (pp == 0.):
580- vc[0] = complex(0.,0.)
581- vc[1] = complex(-nhel*sqh,0.)
582- vc[2] = complex(0.,nsvahl*sqh)
583- vc[3] = complex(hel0,0.)
584-
585- else:
586- emp = p[0]/(vmass*pp)
587- vc[0] = complex(hel0*pp/vmass,0.)
588- vc[3] = complex(hel0*p[3]*emp+nhel*pt/pp*sqh)
589- if (pt != 0.):
590- pzpt = p[3]/(pp*pt)*sqh*nhel
591- vc[1] = complex(hel0*p[1]*emp-p[1]*pzpt, \
592- -nsvahl*p[2]/pt*sqh)
593- vc[2] = complex(hel0*p[2]*emp-p[2]*pzpt, \
594- nsvahl*p[1]/pt*sqh)
595- else:
596- vc[1] = complex(-nhel*sqh,0.)
597- vc[2] = complex(0.,nsvahl*sign(sqh,p[3]))
598- else:
599- pp = p[0]
600- pt = math.sqrt(math.pow(p[1],2)+math.pow(p[2],2))
601- vc[0] = complex(0.,0.)
602- vc[3] = complex(nhel*pt/pp*sqh)
603- if (pt != 0.):
604- pzpt = p[3]/(pp*pt)*sqh*nhel
605- vc[1] = complex(-p[1]*pzpt,-nsv*p[2]/pt*sqh)
606- vc[2] = complex(-p[2]*pzpt,nsv*p[1]/pt*sqh)
607- else:
608- vc[1] = complex(-nhel*sqh,0.)
609- vc[2] = complex(0.,nsv*sign(sqh,p[3]))
610-
611- return vc
612+ """ initialize a vector wavefunction. nhel=4 is for checking BRST"""
613+
614+ vc = WaveFunction(3)
615+
616+ sqh = sqrt(0.5)
617+ nsvahl = nsv*abs(nhel)
618+ pt2 = p[1]**2 + p[2]**2
619+ pp = min(p[0],sqrt(pt2 + p[3]**2))
620+ pt = min(pp,sqrt(pt2))
621+
622+ vc[4] = complex(p[0]*nsv,p[3]*nsv)
623+ vc[5] = complex(p[1]*nsv,p[2]*nsv)
624+
625+ if (nhel == 4):
626+ if (vmass == 0.):
627+ vc[0] = 1.
628+ vc[1]=p[1]/p[0]
629+ vc[2]=p[2]/p[0]
630+ vc[3]=p[3]/p[0]
631+ else:
632+ vc[0] = p[0]/vmass
633+ vc[1] = p[1]/vmass
634+ vc[2] = p[2]/vmass
635+ vc[3] = p[3]/vmass
636+
637+ return vc
638+
639+ if (vmass != 0.):
640+ hel0 = 1.-abs(nhel)
641+
642+ if (pp == 0.):
643+ vc[0] = complex(0.,0.)
644+ vc[1] = complex(-nhel*sqh,0.)
645+ vc[2] = complex(0.,nsvahl*sqh)
646+ vc[3] = complex(hel0,0.)
647+
648+ else:
649+ emp = p[0]/(vmass*pp)
650+ vc[0] = complex(hel0*pp/vmass,0.)
651+ vc[3] = complex(hel0*p[3]*emp+nhel*pt/pp*sqh)
652+ if (pt != 0.):
653+ pzpt = p[3]/(pp*pt)*sqh*nhel
654+ vc[1] = complex(hel0*p[1]*emp-p[1]*pzpt, \
655+ -nsvahl*p[2]/pt*sqh)
656+ vc[2] = complex(hel0*p[2]*emp-p[2]*pzpt, \
657+ nsvahl*p[1]/pt*sqh)
658+ else:
659+ vc[1] = complex(-nhel*sqh,0.)
660+ vc[2] = complex(0.,nsvahl*sign(sqh,p[3]))
661+ else:
662+ pp = p[0]
663+ pt = sqrt(p[1]**2 + p[2]**2)
664+ vc[0] = complex(0.,0.)
665+ vc[3] = complex(nhel*pt/pp*sqh)
666+ if (pt != 0.):
667+ pzpt = p[3]/(pp*pt)*sqh*nhel
668+ vc[1] = complex(-p[1]*pzpt,-nsv*p[2]/pt*sqh)
669+ vc[2] = complex(-p[2]*pzpt,nsv*p[1]/pt*sqh)
670+ else:
671+ vc[1] = complex(-nhel*sqh,0.)
672+ vc[2] = complex(0.,nsv*sign(sqh,p[3]))
673+
674+ return vc
675
676 def sign(x,y):
677- """Fortran's sign transfer function"""
678- if (y < 0.):
679- return -abs(x)
680- else:
681- return abs(x)
682-
683+ """Fortran's sign transfer function"""
684+ if (y < 0.):
685+ return -abs(x)
686+ else:
687+ return abs(x)
688+
689
690 def sxxxxx(p,nss):
691- """initialize a scalar wavefunction"""
692-
693- sc = WaveFunction(1)
694-
695- sc[0] = complex(1.,0.)
696- sc[1] = complex(p[0]*nss,p[3]*nss)
697- sc[2] = complex(p[1]*nss,p[2]*nss)
698- return sc
699+ """initialize a scalar wavefunction"""
700+
701+ sc = WaveFunction(1)
702+
703+ sc[0] = complex(1.,0.)
704+ sc[1] = complex(p[0]*nss,p[3]*nss)
705+ sc[2] = complex(p[1]*nss,p[2]*nss)
706+ return sc
707+
708+
709+def txxxxx(p, tmass, nhel, nst):
710+ """ initialize a tensor wavefunction"""
711+
712+ tc = WaveFunction(5)
713+
714+ sqh = sqrt(0.5)
715+ sqs = sqrt(1/6)
716+
717+ pt2 = p[1]**2 + p[2]**2
718+ pp = min(p[0],sqrt(pt2+p[3]**2))
719+ pt = min(pp,sqrt(pt2))
720+
721+ ft = {}
722+ ft[(4,0)] = complex(p[0], p[3]) * nst
723+ ft[(5,0)] = complex(p[1], p[2]) * nst
724+
725+ if ( nhel >= 0 ):
726+ #construct eps+
727+ ep = [0] * 4
728+
729+ if ( pp == 0 ):
730+ #ep[0] = 0
731+ ep[1] = -sqh
732+ ep[2] = complex(0, nst*sqh)
733+ #ep[3] = 0
734+ else:
735+ #ep[0] = 0
736+ ep[3] = pt/pp*sqh
737+ if (pt != 0):
738+ pzpt = p[3]/(pp*pt)*sqh
739+ ep[1] = complex( -p[1]*pzpt , -nst*p[2]/pt*sqh )
740+ ep[2] = complex( -p[2]*pzpt , nst*p[1]/pt*sqh )
741+ else:
742+ ep[1] = -sqh
743+ ep[2] = complex( 0 , nst*sign(sqh,p[3]) )
744+
745+
746+
747+ if ( nhel <= 0 ):
748+ #construct eps-
749+ em = [0] * 4
750+ if ( pp == 0 ):
751+ #em[0] = 0
752+ em[1] = sqh
753+ em[2] = cpmplex( 0 , nst*sqh )
754+ #em[3] = 0
755+ else:
756+ #em[0] = 0
757+ em[3] = -pt/pp*sqh
758+ if pt:
759+ pzpt = -p[3]/(pp*pt)*sqh
760+ em[1] = complex( -p[1]*pzpt , -nst*p[2]/pt*sqh )
761+ em[2] = complex( -p[2]*pzpt , nst*p[1]/pt*sqh )
762+ else:
763+ em[1] = sqh
764+ em[2] = complex( 0 , nst*sign(sqh,p[3]) )
765+
766+
767+ if ( abs(nhel) <= 1 ):
768+ #construct eps0
769+ e0 = [0] * 4
770+ if ( pp == 0 ):
771+ #e0[0] = dcmplx( rZero )
772+ #e0[1] = dcmplx( rZero )
773+ #e0[2] = dcmplx( rZero )
774+ e0[3] = 1
775+ else:
776+ emp = p[0]/(tmass*pp)
777+ e0[0] = pp/tmass
778+ e0[3] = p[3]*emp
779+ if pt:
780+ e0[1] = p[1]*emp
781+ e0[2] = p[2]*emp
782+ #else:
783+ # e0[1] = dcmplx( rZero )
784+ # e0[2] = dcmplx( rZero )
785+
786+ if nhel == 2:
787+ for j in range(4):
788+ for i in range(4):
789+ ft[(i,j)] = ep[i]*ep[j]
790+ elif nhel == -2:
791+ for j in range(4):
792+ for i in range(4):
793+ ft[(i,j)] = em[i]*em[j]
794+ elif tmass == 0:
795+ for j in range(4):
796+ for i in range(4):
797+ ft[(i,j)] = 0
798+ elif nhel == 1:
799+ for j in range(4):
800+ for i in range(4):
801+ ft[(i,j)] = sqh*( ep[i]*e0[j] + e0[i]*ep[j] )
802+ elif nhel == 0:
803+ for j in range(4):
804+ for i in range(4):
805+ ft[(i,j)] = sqs*( ep[i]*em[j] + em[i]*ep[j] + 2 *e0[i]*e0[j] )
806+ elif nhel == -1:
807+ for j in range(4):
808+ for i in range(4):
809+ ft[(i,j)] = sqh*( em[i]*e0[j] + e0[i]*em[j] )
810+
811+ else:
812+ raise Exception, 'invalid helicity TXXXXXX'
813+
814+ tc[0] = ft[(0,0)]
815+ tc[1] = ft[(0,1)]
816+ tc[2] = ft[(0,2)]
817+ tc[3] = ft[(0,3)]
818+ tc[4] = ft[(1,0)]
819+ tc[5] = ft[(1,1)]
820+ tc[6] = ft[(1,2)]
821+ tc[7] = ft[(1,3)]
822+ tc[8] = ft[(2,0)]
823+ tc[9] = ft[(2,1)]
824+ tc[10] = ft[(2,2)]
825+ tc[11] = ft[(2,3)]
826+ tc[12] = ft[(3,0)]
827+ tc[13] = ft[(3,1)]
828+ tc[14] = ft[(3,2)]
829+ tc[15] = ft[(3,3)]
830+ tc[16] = ft[(4,0)]
831+ tc[17] = ft[(5,0)]
832+
833+ return tc
834+
835+
836
837
838=== modified file 'madgraph/core/base_objects.py'
839--- madgraph/core/base_objects.py 2010-12-15 21:03:07 +0000
840+++ madgraph/core/base_objects.py 2011-01-17 15:18:53 +0000
841@@ -341,16 +341,22 @@
842 if spin == 1:
843 # Scalar
844 return [ 0 ]
845- if spin == 2:
846+ elif spin == 2:
847 # Spinor
848 return [ -1, 1 ]
849- if spin == 3 and self.get('mass').lower() == 'zero':
850+ elif spin == 3 and self.get('mass').lower() == 'zero':
851 # Massless vector
852 return [ -1, 1 ]
853- if spin == 3:
854+ elif spin == 3:
855 # Massive vector
856 return [ -1, 0, 1 ]
857-
858+ elif spin == 5 and self.get('mass').lower() == 'zero':
859+ # Massless tensor
860+ return [-2, -1, 1, 2]
861+ elif spin == 5:
862+ # Massive tensor
863+ return [-2, -1, 0, 1, 2]
864+
865 raise self.PhysicsObjectError, \
866 "No helicity state assignment for spin %d particles" % spin
867
868
869=== modified file 'madgraph/interface/cmd_interface.py'
870--- madgraph/interface/cmd_interface.py 2010-12-18 22:39:25 +0000
871+++ madgraph/interface/cmd_interface.py 2011-01-17 15:18:53 +0000
872@@ -1304,7 +1304,8 @@
873
874 # Options and formats available
875 _display_opts = ['particles', 'interactions', 'processes', 'diagrams',
876- 'multiparticles', 'couplings', 'lorentz', 'checks']
877+ 'multiparticles', 'couplings', 'lorentz', 'checks',
878+ 'parameters']
879 _add_opts = ['process']
880 _save_opts = ['model', 'processes']
881 _tutorial_opts = ['start', 'stop']
882@@ -1510,7 +1511,26 @@
883 print "Interactions %s has the following property:" % arg
884 print self._curr_model['interactions'][int(arg)-1]
885
886-
887+ elif args[0] == 'parameters' and len(args) == 1:
888+ text = "Current model contains %i parameters\n" % \
889+ sum([len(part) for part in
890+ self._curr_model['parameters'].values()])
891+
892+ for key, item in self._curr_model['parameters'].items():
893+ text += '\nparameter type: %s\n' % str(key)
894+ for value in item:
895+ if hasattr(value, 'expr'):
896+ if value.value is not None:
897+ text+= ' %s = %s = %s\n' % (value.name, value.expr ,value.value)
898+ else:
899+ text+= ' %s = %s\n' % (value.name, value.expr)
900+ else:
901+ if value.value is not None:
902+ text+= ' %s = %s\n' % (value.name, value.value)
903+ else:
904+ text+= ' %s \n' % (value.name)
905+ pydoc.pager(text)
906+
907 elif args[0] == 'processes':
908 for amp in self._curr_amps:
909 print amp.nice_string_processes()
910
911=== modified file 'madgraph/iolibs/export_v4.py'
912--- madgraph/iolibs/export_v4.py 2011-01-14 18:19:02 +0000
913+++ madgraph/iolibs/export_v4.py 2011-01-17 15:18:53 +0000
914@@ -1616,6 +1616,7 @@
915 one_mass = particle.get('mass')
916 if one_mass.lower() != 'zero':
917 masses.add(one_mass)
918+
919 # find width
920 one_width = particle.get('width')
921 if one_width.lower() != 'zero':
922
923=== modified file 'madgraph/various/process_checks.py'
924--- madgraph/various/process_checks.py 2010-12-15 21:03:07 +0000
925+++ madgraph/various/process_checks.py 2011-01-17 15:18:53 +0000
926@@ -46,7 +46,7 @@
927 import models.model_reader as model_reader
928 import aloha.template_files.wavefunctions as wavefunctions
929 from aloha.template_files.wavefunctions import \
930- ixxxxx, oxxxxx, vxxxxx, sxxxxx
931+ ixxxxx, oxxxxx, vxxxxx, sxxxxx, txxxxx
932
933 #===============================================================================
934 # Logger for process_checks
935
936=== modified file 'tests/unit_tests/core/test_base_objects.py'
937--- tests/unit_tests/core/test_base_objects.py 2010-11-05 02:28:04 +0000
938+++ tests/unit_tests/core/test_base_objects.py 2011-01-17 15:18:53 +0000
939@@ -197,10 +197,10 @@
940 test_part.set('mass', 'Zero')
941 self.assertEqual(test_part.get_helicity_states(), [-1, 1])
942 test_part.set('spin', 5)
943- self.assertRaises(test_part.PhysicsObjectError,
944- test_part.get_helicity_states)
945-
946-
947+ self.assertEqual(test_part.get_helicity_states(), [-2, -1, 1, 2])
948+ test_part.set('mass', 'M')
949+ self.assertEqual(test_part.get_helicity_states(), [-2, -1, 0, 1, 2])
950+
951 def test_particle_list(self):
952 """Test particle list initialization, search and dict generation
953 functions."""
954
955=== modified file 'tests/unit_tests/various/test_aloha.py'
956--- tests/unit_tests/various/test_aloha.py 2010-12-19 01:37:23 +0000
957+++ tests/unit_tests/various/test_aloha.py 2011-01-17 15:18:53 +0000
958@@ -620,8 +620,46 @@
959 else:
960 self.assertFalse(obj1 is term[1])
961
962-
963-
964+ def testdealingwithpower3(self):
965+ """Check that the power is correctly set in a product in the full chain"""
966+
967+ F1_1, F1_2, F1_3, F1_4 = 1,2,3,4
968+
969+ P1_0, P1_1, P1_2, P1_3 = 12, 0, 0, 12
970+ P2_0, P2_1, P2_2, P2_3 = 12, 0, 12, 0
971+ P3_0, P3_1, P3_2, P3_3 = 20, 0, 12, 12
972+ M1, M2, M3 = 0, 0, 100
973+
974+ F2_1, F2_2, F2_3, F2_4 = 5,5,6,7
975+ T3_1, T3_2, T3_3, T3_4 = 8,9,10,11
976+ T3_5, T3_6, T3_7, T3_8 = 8,9,10,11
977+ T3_9, T3_10, T3_11, T3_12 = 8,9,10,11
978+ T3_13, T3_14, T3_15, T3_16 = 8,9,10,11
979+
980+
981+
982+ p1 = aloha_obj.P('mu',2)
983+ gamma1 = aloha_obj.Gamma('mu','a','b')
984+ metric = aloha_obj.Spin2('nu','rho',3)
985+ p2 = aloha_obj.P('rho',2)
986+ gamma2 = aloha_obj.Gamma('nu','b','c')
987+ F1 = aloha_obj.Spinor('c',1)
988+
989+
990+ lor1 = p1 * gamma1 * gamma2 * F1
991+ lor2 = metric * p2
992+ lor1.simplify()
993+ new_lor = lor1.expand()
994+
995+ lor2.simplify()
996+ new_lor2 = lor2.expand()
997+
998+ expr = new_lor * new_lor2
999+
1000+ self.assertEqual((-864+288j), eval(str(expr.get_rep([0]))))
1001+ self.assertEqual((288+864j), eval(str(expr.get_rep([1]))))
1002+ self.assertEqual((2016+288j), eval(str(expr.get_rep([2]))))
1003+ self.assertEqual((-288+2016j), eval(str(expr.get_rep([3]))))
1004
1005
1006 def test_obj_are_not_modified(self):
1007@@ -2373,6 +2411,10 @@
1008 self.assertEqual(complex(0,-1)*ufo_value, v4_value)
1009
1010
1011+
1012+
1013+
1014+
1015
1016
1017
1018
1019=== modified file 'tests/unit_tests/various/test_rambo.py'
1020--- tests/unit_tests/various/test_rambo.py 2010-09-07 23:51:48 +0000
1021+++ tests/unit_tests/various/test_rambo.py 2011-01-17 15:18:53 +0000
1022@@ -14,6 +14,7 @@
1023 ################################################################################
1024
1025 import madgraph.various.rambo as rambo
1026+import aloha.template_files.wavefunctions as wavefunctions
1027 import tests.unit_tests as unittest
1028
1029
1030@@ -56,7 +57,489 @@
1031 for i in range(1,3):
1032 self.assertAlmostEqual(self.m2[i]**2, P[(4,i)]**2 - P[(1,i)]**2 - P[(2,i)]**2 - P[(3,i)]**2)
1033
1034-
1035-
1036-
1037-
1038+
1039+class test_wavefunctions(unittest.TestCase):
1040+ """ check that the wavefunctions TXXXX, IXXXXX, ... are correctly define"""
1041+
1042+ def test_T(self):
1043+ """check T wavefunctions against fortran output"""
1044+
1045+ #input
1046+ P = [624.99999999999989, 83.193213332874592, 333.62309211609107, -149.66469744815910 ]
1047+ M = 500
1048+
1049+ # first
1050+ NHEL = 2
1051+ NST = 1
1052+
1053+ solution = [
1054+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1055+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1056+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1057+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1058+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1059+ complex(-0.46606677614840392 ,-9.36959948731358322E-002),
1060+ complex( 0.13607969464455591 ,-0.17618862729680237 ),
1061+ complex( 4.42705319225060942E-002,-0.44483078948812155 ),
1062+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1063+ complex( 0.13607969464455591 ,-0.17618862729680237 ),
1064+ complex( 4.57095198364004252E-002, 9.36959948731358322E-002),
1065+ complex( 0.17753457473164128 , 0.11092428444383282 ),
1066+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1067+ complex( 4.42705319225060942E-002,-0.44483078948812155 ),
1068+ complex( 0.17753457473164128 , 0.11092428444383282 ),
1069+ complex( 0.42035725631200371 , 0.0000000000000000 ),
1070+ complex( 624.99999999999989 , -149.66469744815910 ),
1071+ complex( 83.193213332874592 , 333.62309211609107 )]
1072+
1073+ values = wavefunctions.txxxxx(P, M, NHEL, NST)
1074+ self.assertEqual(len(values), len(solution))
1075+ for i in range(len(values)):
1076+ self.assertAlmostEqual(values[i], solution[i])
1077+
1078+ # NEXT
1079+ NHEL = 1
1080+ solution = [
1081+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1082+ complex( 3.62119349359586659E-002,-0.36385791873139811 ),
1083+ complex( 0.14521782752280418 , 9.07326566152192315E-002),
1084+ complex( 0.34383932052304755 , 0.0000000000000000 ),
1085+ complex( 3.62119349359586659E-002,-0.36385791873139811 ),
1086+ complex( 2.67785531406523024E-002,-0.26907119516335065 ),
1087+ complex( 0.10738789070969987 ,-0.50596916746687115 ),
1088+ complex( 0.10304644292578606 , 0.24202971253800695 ),
1089+ complex( 0.14521782752280418 , 9.07326566152192315E-002),
1090+ complex( 0.10738789070969987 ,-0.50596916746687115 ),
1091+ complex( 0.43064907243145895 , 0.26907119516335065 ),
1092+ complex( 0.41323891148317965 ,-6.03532248932644386E-002),
1093+ complex( 0.34383932052304755 , 0.0000000000000000 ),
1094+ complex( 0.10304644292578606 , 0.24202971253800695 ),
1095+ complex( 0.41323891148317965 ,-6.03532248932644386E-002),
1096+ complex(-0.45742762557211131 , 0.0000000000000000 ),
1097+ complex( 624.99999999999989 , -149.66469744815910 ),
1098+ complex( 83.193213332874592 , 333.62309211609107 )]
1099+
1100+ values = wavefunctions.txxxxx(P, M, NHEL, NST)
1101+ self.assertEqual(len(values), len(solution))
1102+ for i in range(len(values)):
1103+ self.assertAlmostEqual(values[i], solution[i])
1104+
1105+ #NEXT
1106+ NHEL = 0
1107+ solution = [
1108+ complex( 0.45927932677184580 , 0.0000000000000000 ),
1109+ complex( 0.16981743560670751 , 0.0000000000000000 ),
1110+ complex( 0.68100528507831026 , 0.0000000000000000 ),
1111+ complex(-0.30550178438001113 , 0.0000000000000000 ),
1112+ complex( 0.16981743560670751 , 0.0000000000000000 ),
1113+ complex(-0.32536602932851599 , 0.0000000000000000 ),
1114+ complex( 0.33237610537903167 , 0.0000000000000000 ),
1115+ complex(-0.14910529404613407 , 0.0000000000000000 ),
1116+ complex( 0.68100528507831026 , 0.0000000000000000 ),
1117+ complex( 0.33237610537903167 , 0.0000000000000000 ),
1118+ complex( 0.92465303140679633 , 0.0000000000000000 ),
1119+ complex(-0.59794503971747703 , 0.0000000000000000 ),
1120+ complex(-0.30550178438001113 , 0.0000000000000000 ),
1121+ complex(-0.14910529404613407 , 0.0000000000000000 ),
1122+ complex(-0.59794503971747703 , 0.0000000000000000 ),
1123+ complex(-0.14000767530643479 , 0.0000000000000000 ),
1124+ complex( 624.99999999999989 , -149.66469744815910 ),
1125+ complex( 83.193213332874592 , 333.62309211609107 )]
1126+
1127+ values = wavefunctions.txxxxx(P, M, NHEL, NST)
1128+ self.assertEqual(len(values), len(solution))
1129+ for i in range(len(values)):
1130+ self.assertAlmostEqual(values[i], solution[i])
1131+
1132+ NHEL=-1
1133+ solution = [
1134+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1135+ complex(-3.62119349359586659E-002,-0.36385791873139811 ),
1136+ complex(-0.14521782752280418 , 9.07326566152192315E-002),
1137+ complex(-0.34383932052304755 , 0.0000000000000000 ),
1138+ complex(-3.62119349359586659E-002,-0.36385791873139811 ),
1139+ complex(-2.67785531406523024E-002,-0.26907119516335065 ),
1140+ complex(-0.10738789070969987 ,-0.50596916746687115 ),
1141+ complex(-0.10304644292578606 , 0.24202971253800695 ),
1142+ complex(-0.14521782752280418 , 9.07326566152192315E-002),
1143+ complex(-0.10738789070969987 ,-0.50596916746687115 ),
1144+ complex(-0.43064907243145895 , 0.26907119516335065 ),
1145+ complex(-0.41323891148317965 ,-6.03532248932644386E-002),
1146+ complex(-0.34383932052304755 , 0.0000000000000000 ),
1147+ complex(-0.10304644292578606 , 0.24202971253800695 ),
1148+ complex(-0.41323891148317965 ,-6.03532248932644386E-002),
1149+ complex( 0.45742762557211131 , -0.0000000000000000 ),
1150+ complex( 624.99999999999989 , -149.66469744815910 ),
1151+ complex( 83.193213332874592 , 333.62309211609107 )]
1152+
1153+ values = wavefunctions.txxxxx(P, M, NHEL, NST)
1154+ self.assertEqual(len(values), len(solution))
1155+ for i in range(len(values)):
1156+ self.assertAlmostEqual(values[i], solution[i])
1157+
1158+ NHEL=-2
1159+ solution = [
1160+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1161+ complex( 0.0000000000000000 , -0.0000000000000000 ),
1162+ complex( -0.0000000000000000 , 0.0000000000000000 ),
1163+ complex( -0.0000000000000000 , 0.0000000000000000 ),
1164+ complex( 0.0000000000000000 , -0.0000000000000000 ),
1165+ complex(-0.46606677614840392 , 9.36959948731358322E-002),
1166+ complex( 0.13607969464455591 , 0.17618862729680237 ),
1167+ complex( 4.42705319225060942E-002, 0.44483078948812155 ),
1168+ complex( -0.0000000000000000 , 0.0000000000000000 ),
1169+ complex( 0.13607969464455591 , 0.17618862729680237 ),
1170+ complex( 4.57095198364004252E-002,-9.36959948731358322E-002),
1171+ complex( 0.17753457473164128 ,-0.11092428444383282 ),
1172+ complex( -0.0000000000000000 , 0.0000000000000000 ),
1173+ complex( 4.42705319225060942E-002, 0.44483078948812155 ),
1174+ complex( 0.17753457473164128 ,-0.11092428444383282 ),
1175+ complex( 0.42035725631200371 , -0.0000000000000000 ),
1176+ complex( 624.99999999999989 , -149.66469744815910 ),
1177+ complex( 83.193213332874592 , 333.62309211609107 )]
1178+
1179+ values = wavefunctions.txxxxx(P, M, NHEL, NST)
1180+ self.assertEqual(len(values), len(solution))
1181+ for i in range(len(values)):
1182+ self.assertAlmostEqual(values[i], solution[i])
1183+
1184+
1185+ # first
1186+ NHEL = 2
1187+ NST = -1
1188+ solution = [
1189+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1190+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1191+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1192+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1193+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1194+ complex(-0.46606677614840392 , 9.36959948731358322E-002),
1195+ complex( 0.13607969464455591 , 0.17618862729680237 ),
1196+ complex( 4.42705319225060942E-002, 0.44483078948812155 ),
1197+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1198+ complex( 0.13607969464455591 , 0.17618862729680237 ),
1199+ complex( 4.57095198364004252E-002,-9.36959948731358322E-002),
1200+ complex( 0.17753457473164128 ,-0.11092428444383282 ),
1201+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1202+ complex( 4.42705319225060942E-002, 0.44483078948812155 ),
1203+ complex( 0.17753457473164128 ,-0.11092428444383282 ),
1204+ complex( 0.42035725631200371 , 0.0000000000000000 ),
1205+ complex( -624.99999999999989 , 149.66469744815910 ),
1206+ complex( -83.193213332874592 , -333.62309211609107 )]
1207+
1208+ values = wavefunctions.txxxxx(P, M, NHEL, NST)
1209+ self.assertEqual(len(values), len(solution))
1210+ for i in range(len(values)):
1211+ self.assertAlmostEqual(values[i], solution[i])
1212+
1213+ # first
1214+ NHEL = 1
1215+ solution = [
1216+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1217+ complex( 3.62119349359586659E-002, 0.36385791873139811 ),
1218+ complex( 0.14521782752280418 ,-9.07326566152192315E-002),
1219+ complex( 0.34383932052304755 , 0.0000000000000000 ),
1220+ complex( 3.62119349359586659E-002, 0.36385791873139811 ),
1221+ complex( 2.67785531406523024E-002, 0.26907119516335065 ),
1222+ complex( 0.10738789070969987 , 0.50596916746687115 ),
1223+ complex( 0.10304644292578606 ,-0.24202971253800695 ),
1224+ complex( 0.14521782752280418 ,-9.07326566152192315E-002),
1225+ complex( 0.10738789070969987 , 0.50596916746687115 ),
1226+ complex( 0.43064907243145895 ,-0.26907119516335065 ),
1227+ complex( 0.41323891148317965 , 6.03532248932644386E-002),
1228+ complex( 0.34383932052304755 , 0.0000000000000000 ),
1229+ complex( 0.10304644292578606 ,-0.24202971253800695 ),
1230+ complex( 0.41323891148317965 , 6.03532248932644386E-002),
1231+ complex(-0.45742762557211131 , 0.0000000000000000 ),
1232+ complex( -624.99999999999989 , 149.66469744815910 ),
1233+ complex( -83.193213332874592 , -333.62309211609107 )]
1234+
1235+ values = wavefunctions.txxxxx(P, M, NHEL, NST)
1236+ self.assertEqual(len(values), len(solution))
1237+ for i in range(len(values)):
1238+ self.assertAlmostEqual(values[i], solution[i])
1239+
1240+ NHEL = 0
1241+ solution = [
1242+ complex( 0.45927932677184580 , 0.0000000000000000 ),
1243+ complex( 0.16981743560670751 , 0.0000000000000000 ),
1244+ complex( 0.68100528507831026 , 0.0000000000000000 ),
1245+ complex(-0.30550178438001113 , 0.0000000000000000 ),
1246+ complex( 0.16981743560670751 , 0.0000000000000000 ),
1247+ complex(-0.32536602932851599 , 0.0000000000000000 ),
1248+ complex( 0.33237610537903167 , 0.0000000000000000 ),
1249+ complex(-0.14910529404613407 , 0.0000000000000000 ),
1250+ complex( 0.68100528507831026 , 0.0000000000000000 ),
1251+ complex( 0.33237610537903167 , 0.0000000000000000 ),
1252+ complex( 0.92465303140679633 , 0.0000000000000000 ),
1253+ complex(-0.59794503971747703 , 0.0000000000000000 ),
1254+ complex(-0.30550178438001113 , 0.0000000000000000 ),
1255+ complex(-0.14910529404613407 , 0.0000000000000000 ),
1256+ complex(-0.59794503971747703 , 0.0000000000000000 ),
1257+ complex(-0.14000767530643479 , 0.0000000000000000 ),
1258+ complex( -624.99999999999989 , 149.66469744815910 ),
1259+ complex( -83.193213332874592 , -333.62309211609107 )]
1260+
1261+ values = wavefunctions.txxxxx(P, M, NHEL, NST)
1262+ self.assertEqual(len(values), len(solution))
1263+ for i in range(len(values)):
1264+ self.assertAlmostEqual(values[i], solution[i])
1265+
1266+ NHEL = -1
1267+ solution = [
1268+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1269+ complex(-3.62119349359586659E-002, 0.36385791873139811 ),
1270+ complex(-0.14521782752280418 ,-9.07326566152192315E-002),
1271+ complex(-0.34383932052304755 , 0.0000000000000000 ),
1272+ complex(-3.62119349359586659E-002, 0.36385791873139811 ),
1273+ complex(-2.67785531406523024E-002, 0.26907119516335065 ),
1274+ complex(-0.10738789070969987 , 0.50596916746687115 ),
1275+ complex(-0.10304644292578606 ,-0.24202971253800695 ),
1276+ complex(-0.14521782752280418 ,-9.07326566152192315E-002),
1277+ complex(-0.10738789070969987 , 0.50596916746687115 ),
1278+ complex(-0.43064907243145895 ,-0.26907119516335065 ),
1279+ complex(-0.41323891148317965 , 6.03532248932644386E-002),
1280+ complex(-0.34383932052304755 , 0.0000000000000000 ),
1281+ complex(-0.10304644292578606 ,-0.24202971253800695 ),
1282+ complex(-0.41323891148317965 , 6.03532248932644386E-002),
1283+ complex( 0.45742762557211131 , -0.0000000000000000 ),
1284+ complex( -624.99999999999989 , 149.66469744815910 ),
1285+ complex( -83.193213332874592 , -333.62309211609107 )]
1286+
1287+
1288+ values = wavefunctions.txxxxx(P, M, NHEL, NST)
1289+ self.assertEqual(len(values), len(solution))
1290+ for i in range(len(values)):
1291+ self.assertAlmostEqual(values[i], solution[i])
1292+
1293+ NHEL = -2
1294+ solution = [
1295+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1296+ complex( -0.0000000000000000 , 0.0000000000000000 ),
1297+ complex( 0.0000000000000000 , -0.0000000000000000 ),
1298+ complex( -0.0000000000000000 , 0.0000000000000000 ),
1299+ complex( -0.0000000000000000 , 0.0000000000000000 ),
1300+ complex(-0.46606677614840392 ,-9.36959948731358322E-002),
1301+ complex( 0.13607969464455591 ,-0.17618862729680237 ),
1302+ complex( 4.42705319225060942E-002,-0.44483078948812155 ),
1303+ complex( 0.0000000000000000 , -0.0000000000000000 ),
1304+ complex( 0.13607969464455591 ,-0.17618862729680237 ),
1305+ complex( 4.57095198364004252E-002, 9.36959948731358322E-002),
1306+ complex( 0.17753457473164128 , 0.11092428444383282 ),
1307+ complex( -0.0000000000000000 , 0.0000000000000000 ),
1308+ complex( 4.42705319225060942E-002,-0.44483078948812155 ),
1309+ complex( 0.17753457473164128 , 0.11092428444383282 ),
1310+ complex( 0.42035725631200371 , -0.0000000000000000 ),
1311+ complex( -624.99999999999989 , 149.66469744815910 ),
1312+ complex(-83.193213332874592 , -333.62309211609107)]
1313+
1314+ values = wavefunctions.txxxxx(P, M, NHEL, NST)
1315+ self.assertEqual(len(values), len(solution))
1316+ for i in range(len(values)):
1317+ self.assertAlmostEqual(values[i], solution[i])
1318+
1319+ def test_I(self):
1320+ """check I wavefunctions against fortran output"""
1321+
1322+ P = [2499.9999999999991, 537.98548101331721, 2157.4401624693646, -967.83657009607577]
1323+ M = 607.71370000000002
1324+
1325+ NHEL=-1
1326+ NRST=-1
1327+ solution = [
1328+ complex( 4.7465641991565066 , 0.0000000000000000 ),
1329+ complex( 1.7524192771692548 , 7.0275869209877353 ),
1330+ complex( -38.466940073898407 , -0.0000000000000000 ),
1331+ complex( -14.201895200573350 , -56.952724878721050 ),
1332+ complex( -2499.9999999999991 , 967.83657009607577 ),
1333+ complex( -537.98548101331721 , -2157.4401624693646 )]
1334+
1335+
1336+ values = wavefunctions.ixxxxx(P, M, NHEL, NRST)
1337+ self.assertEqual(len(values), len(solution))
1338+ for i in range(len(values)):
1339+ self.assertAlmostEqual(values[i], solution[i])
1340+
1341+
1342+ NHEL=1
1343+ NRST=-1
1344+ solution = [
1345+ complex( 14.201895200573350 , -56.952724878721050 ),
1346+ complex( -38.466940073898407 , -0.0000000000000000 ),
1347+ complex( -1.7524192771692548 , 7.0275869209877353 ),
1348+ complex( 4.7465641991565066 , 0.0000000000000000 ),
1349+ complex( -2499.9999999999991 , 967.83657009607577 ),
1350+ complex( -537.98548101331721 , -2157.4401624693646 )]
1351+
1352+ values = wavefunctions.ixxxxx(P, M, NHEL, NRST)
1353+ self.assertEqual(len(values), len(solution))
1354+ for i in range(len(values)):
1355+ self.assertAlmostEqual(values[i], solution[i])
1356+
1357+ NHEL=1
1358+ NRST=1
1359+ solution = [
1360+ complex( 4.7465641991565066 , 0.0000000000000000 ),
1361+ complex( 1.7524192771692548 , 7.0275869209877353 ),
1362+ complex( 38.466940073898407 , 0.0000000000000000 ),
1363+ complex( 14.201895200573350 , 56.952724878721050 ),
1364+ complex( 2499.9999999999991 , -967.83657009607577 ),
1365+ complex( 537.98548101331721 , 2157.4401624693646 )]
1366+
1367+ values = wavefunctions.ixxxxx(P, M, NHEL, NRST)
1368+ self.assertEqual(len(values), len(solution))
1369+ for i in range(len(values)):
1370+ self.assertAlmostEqual(values[i], solution[i])
1371+
1372+ NHEL=-1
1373+ NRST=1
1374+ solution = [
1375+ complex( -14.201895200573350 , 56.952724878721050 ),
1376+ complex( 38.466940073898407 , 0.0000000000000000 ),
1377+ complex( -1.7524192771692548 , 7.0275869209877353 ),
1378+ complex( 4.7465641991565066 , 0.0000000000000000 ),
1379+ complex( 2499.9999999999991 , -967.83657009607577 ),
1380+ complex( 537.98548101331721 , 2157.4401624693646 )]
1381+
1382+
1383+
1384+
1385+ values = wavefunctions.ixxxxx(P, M, NHEL, NRST)
1386+ self.assertEqual(len(values), len(solution))
1387+ for i in range(len(values)):
1388+ self.assertAlmostEqual(values[i], solution[i])
1389+
1390+ def test_O(self):
1391+ """check O wavefunctions against fortran output"""
1392+
1393+
1394+ P = [2499.9999999999991, 537.98548101331721, 2157.4401624693646, -967.83657009607577]
1395+ M = 607.71370000000002
1396+
1397+ NHEL=1
1398+ NRST=1
1399+ solution=[
1400+ complex( 38.466940073898407 , 0.0000000000000000 ),
1401+ complex( 14.201895200573350 , -56.952724878721050 ),
1402+ complex( 4.7465641991565066 , 0.0000000000000000 ),
1403+ complex( 1.7524192771692548 , -7.0275869209877353 ),
1404+ complex( 2499.9999999999991 , -967.83657009607577 ),
1405+ complex( 537.98548101331721 , 2157.4401624693646 )]
1406+
1407+ values = wavefunctions.oxxxxx(P, M, NHEL, NRST)
1408+ self.assertEqual(len(values), len(solution))
1409+ for i in range(len(values)):
1410+ self.assertAlmostEqual(values[i], solution[i])
1411+
1412+
1413+ NHEL=1
1414+ NRST=-1
1415+ solution=[
1416+ complex( -1.7524192771692548 , -7.0275869209877353 ),
1417+ complex( 4.7465641991565066 , 0.0000000000000000 ),
1418+ complex( 14.201895200573350 , 56.952724878721050 ),
1419+ complex( -38.466940073898407 , -0.0000000000000000 ),
1420+ complex( -2499.9999999999991 , 967.83657009607577 ),
1421+ complex( -537.98548101331721 , -2157.4401624693646 )]
1422+
1423+ values = wavefunctions.oxxxxx(P, M, NHEL, NRST)
1424+ self.assertEqual(len(values), len(solution))
1425+ for i in range(len(values)):
1426+ self.assertAlmostEqual(values[i], solution[i])
1427+
1428+ NHEL=-1
1429+ NRST=1
1430+ solution=[
1431+ complex( -1.7524192771692548 , -7.0275869209877353 ),
1432+ complex( 4.7465641991565066 , 0.0000000000000000 ),
1433+ complex( -14.201895200573350 , -56.952724878721050 ),
1434+ complex( 38.466940073898407 , 0.0000000000000000 ),
1435+ complex( 2499.9999999999991 , -967.83657009607577 ),
1436+ complex( 537.98548101331721 , 2157.4401624693646 )]
1437+
1438+ values = wavefunctions.oxxxxx(P, M, NHEL, NRST)
1439+ self.assertEqual(len(values), len(solution))
1440+ for i in range(len(values)):
1441+ self.assertAlmostEqual(values[i], solution[i])
1442+
1443+
1444+ NHEL=-1
1445+ NRST=-1
1446+ solution=[
1447+ complex( -38.466940073898407 , -0.0000000000000000 ),
1448+ complex( -14.201895200573350 , 56.952724878721050 ),
1449+ complex( 4.7465641991565066 , 0.0000000000000000 ),
1450+ complex( 1.7524192771692548 , -7.0275869209877353 ),
1451+ complex( -2499.9999999999991 , 967.83657009607577 ),
1452+ complex( -537.98548101331721 , -2157.4401624693646 )]
1453+
1454+ values = wavefunctions.oxxxxx(P, M, NHEL, NRST)
1455+ self.assertEqual(len(values), len(solution))
1456+ for i in range(len(values)):
1457+ self.assertAlmostEqual(values[i], solution[i])
1458+
1459+ def test_V(self):
1460+ """check V wavefunctions against fortran output"""
1461+
1462+ P = [2500, 0, 0, 2500]
1463+ M = 0
1464+
1465+ NHEL=1
1466+ NRST=1
1467+ solution=[
1468+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1469+ complex(-0.70710678118654757 , 0.0000000000000000 ),
1470+ complex( 0.0000000000000000 , 0.70710678118654757 ),
1471+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1472+ complex( 2500.0000000000000 , 2500.0000000000000 ),
1473+ complex( 0.0000000000000000 , 0.0000000000000000 )]
1474+
1475+ values = wavefunctions.vxxxxx(P, M, NHEL, NRST)
1476+ self.assertEqual(len(values), len(solution))
1477+ for i in range(len(values)):
1478+ self.assertAlmostEqual(values[i], solution[i])
1479+
1480+ NHEL=1
1481+ NRST=-1
1482+ solution=[
1483+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1484+ complex(-0.70710678118654757 , 0.0000000000000000 ),
1485+ complex( 0.0000000000000000 ,-0.70710678118654757 ),
1486+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1487+ complex( -2500.0000000000000 , -2500.0000000000000 ),
1488+ complex( -0.0000000000000000 , -0.0000000000000000 )]
1489+
1490+ values = wavefunctions.vxxxxx(P, M, NHEL, NRST)
1491+ self.assertEqual(len(values), len(solution))
1492+ for i in range(len(values)):
1493+ self.assertAlmostEqual(values[i], solution[i])
1494+
1495+ NHEL=-1
1496+ NRST=-1
1497+ solution=[
1498+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1499+ complex( 0.70710678118654757 , 0.0000000000000000 ),
1500+ complex( 0.0000000000000000 ,-0.70710678118654757 ),
1501+ complex( -0.0000000000000000 , 0.0000000000000000 ),
1502+ complex( -2500.0000000000000 , -2500.0000000000000 ),
1503+ complex( -0.0000000000000000 , -0.0000000000000000 )]
1504+
1505+ values = wavefunctions.vxxxxx(P, M, NHEL, NRST)
1506+ self.assertEqual(len(values), len(solution))
1507+ for i in range(len(values)):
1508+ self.assertAlmostEqual(values[i], solution[i])
1509+
1510+ NHEL=-1
1511+ NRST=1
1512+ solution=[
1513+ complex( 0.0000000000000000 , 0.0000000000000000 ),
1514+ complex( 0.70710678118654757 , 0.0000000000000000 ),
1515+ complex( 0.0000000000000000 , 0.70710678118654757 ),
1516+ complex( -0.0000000000000000 , 0.0000000000000000 ),
1517+ complex( 2500.0000000000000 , 2500.0000000000000 ),
1518+ complex( 0.0000000000000000 , 0.0000000000000000 )]
1519+
1520+ values = wavefunctions.vxxxxx(P, M, NHEL, NRST)
1521+ self.assertEqual(len(values), len(solution))
1522+ for i in range(len(values)):
1523+ self.assertAlmostEqual(values[i], solution[i])
1524\ No newline at end of file

Subscribers

People subscribed via source and target branches