← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_970226/PROC_970226/Source/DHELAS/aloha_functions.f

  • Committer: John Doe
  • Date: 2013-03-25 20:27:02 UTC
  • Revision ID: john.doe@gmail.com-20130325202702-5sk3t1r8h33ca4p4
first clean version

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/mardelcourt/PROC_970226/PROC_970226/Source/DHELAS/aloha_functions.f
 
1
C###############################################################################
 
2
C
 
3
C Copyright (c) 2010 The ALOHA Development team and Contributors
 
4
C
 
5
C This file is a part of the MadGraph 5 project, an application which
 
6
C automatically generates Feynman diagrams and matrix elements for arbitrary
 
7
C high-energy processes in the Standard Model and beyond.
 
8
C
 
9
C It is subject to the ALOHA license which should accompany this
 
10
C distribution.
 
11
C
 
12
C###############################################################################
 
13
      subroutine ixxxxx(p, fmass, nhel, nsf ,fi)
 
14
c
 
15
c This subroutine computes a fermion wavefunction with the flowing-IN
 
16
c fermion number.
 
17
c
 
18
c input:
 
19
c       real    p(0:3)         : four-momentum of fermion
 
20
c       real    fmass          : mass          of fermion
 
21
c       integer nhel = -1 or 1 : helicity      of fermion
 
22
c       integer nsf  = -1 or 1 : +1 for particle, -1 for anti-particle
 
23
c
 
24
c output:
 
25
c       complex fi(6)          : fermion wavefunction               |fi>
 
26
c
 
27
      implicit none
 
28
      double complex fi(6),chi(2)
 
29
      double precision p(0:3),sf(2),sfomeg(2),omega(2),fmass,
 
30
     &     pp,pp3,sqp0p3,sqm(0:1)
 
31
      integer nhel,nsf,ip,im,nh
 
32
 
 
33
      double precision rZero, rHalf, rTwo
 
34
      parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
 
35
 
 
36
c#ifdef HELAS_CHECK
 
37
c      double precision p2
 
38
c      double precision epsi
 
39
c      parameter( epsi = 2.0d-5 )
 
40
c      integer stdo
 
41
c      parameter( stdo = 6 )
 
42
c#endif
 
43
c
 
44
c#ifdef HELAS_CHECK
 
45
c      pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
 
46
c      if ( abs(p(0))+pp.eq.rZero ) then
 
47
c         write(stdo,*)
 
48
c     &        ' helas-error : p(0:3) in ixxxxx is zero momentum'
 
49
c      endif
 
50
c      if ( p(0).le.rZero ) then
 
51
c         write(stdo,*)
 
52
c     &        ' helas-error : p(0:3) in ixxxxx has non-positive energy'
 
53
c         write(stdo,*)
 
54
c     &        '             : p(0) = ',p(0)
 
55
c      endif
 
56
c      p2 = (p(0)-pp)*(p(0)+pp)
 
57
c      if ( abs(p2-fmass**2).gt.p(0)**2*epsi ) then
 
58
c         write(stdo,*)
 
59
c     &        ' helas-error : p(0:3) in ixxxxx has inappropriate mass'
 
60
c         write(stdo,*)
 
61
c     &        '             : p**2 = ',p2,' : fmass**2 = ',fmass**2
 
62
c      endif
 
63
c      if (abs(nhel).ne.1) then
 
64
c         write(stdo,*) ' helas-error : nhel in ixxxxx is not -1,1'
 
65
c         write(stdo,*) '             : nhel = ',nhel
 
66
c      endif
 
67
c      if (abs(nsf).ne.1) then
 
68
c         write(stdo,*) ' helas-error : nsf in ixxxxx is not -1,1'
 
69
c         write(stdo,*) '             : nsf = ',nsf
 
70
c      endif
 
71
c#endif
 
72
 
 
73
      fi(1) = dcmplx(p(0),p(3))*nsf*-1
 
74
      fi(2) = dcmplx(p(1),p(2))*nsf*-1
 
75
 
 
76
      nh = nhel*nsf
 
77
 
 
78
      if ( fmass.ne.rZero ) then
 
79
 
 
80
         pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
81
 
 
82
         if ( pp.eq.rZero ) then
 
83
 
 
84
            sqm(0) = dsqrt(abs(fmass)) ! possibility of negative fermion masses
 
85
            sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
 
86
            ip = (1+nh)/2
 
87
            im = (1-nh)/2
 
88
 
 
89
            fi(3) = ip     * sqm(ip)
 
90
            fi(4) = im*nsf * sqm(ip)
 
91
            fi(5) = ip*nsf * sqm(im)
 
92
            fi(6) = im     * sqm(im)
 
93
 
 
94
         else
 
95
 
 
96
            sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
 
97
            sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
 
98
            omega(1) = dsqrt(p(0)+pp)
 
99
            omega(2) = fmass/omega(1)
 
100
            ip = (3+nh)/2
 
101
            im = (3-nh)/2
 
102
            sfomeg(1) = sf(1)*omega(ip)
 
103
            sfomeg(2) = sf(2)*omega(im)
 
104
            pp3 = max(pp+p(3),rZero)
 
105
            chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
 
106
            if ( pp3.eq.rZero ) then
 
107
               chi(2) = dcmplx(-nh )
 
108
            else
 
109
               chi(2) = dcmplx( nh*p(1) , p(2) )/dsqrt(rTwo*pp*pp3)
 
110
            endif
 
111
 
 
112
            fi(3) = sfomeg(1)*chi(im)
 
113
            fi(4) = sfomeg(1)*chi(ip)
 
114
            fi(5) = sfomeg(2)*chi(im)
 
115
            fi(6) = sfomeg(2)*chi(ip)
 
116
 
 
117
         endif
 
118
 
 
119
      else
 
120
 
 
121
         if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
 
122
            sqp0p3 = 0d0
 
123
         else
 
124
            sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsf
 
125
         end if
 
126
         chi(1) = dcmplx( sqp0p3 )
 
127
         if ( sqp0p3.eq.rZero ) then
 
128
            chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
 
129
         else
 
130
            chi(2) = dcmplx( nh*p(1), p(2) )/sqp0p3
 
131
         endif
 
132
         if ( nh.eq.1 ) then
 
133
            fi(3) = dcmplx( rZero )
 
134
            fi(4) = dcmplx( rZero )
 
135
            fi(5) = chi(1)
 
136
            fi(6) = chi(2)
 
137
         else
 
138
            fi(3) = chi(2)
 
139
            fi(4) = chi(1)
 
140
            fi(5) = dcmplx( rZero )
 
141
            fi(6) = dcmplx( rZero )
 
142
         endif
 
143
      endif
 
144
c
 
145
      return
 
146
      end
 
147
 
 
148
 
 
149
      subroutine ixxxso(p, fmass, nhel, nsf ,fi)
 
150
c Identical to ixxxxx, except that fi returns only the spinor (without the momentum)
 
151
      implicit none
 
152
      double complex fi(4),chi(2)
 
153
      double precision p(0:3),sf(2),sfomeg(2),omega(2),fmass,
 
154
     &     pp,pp3,sqp0p3,sqm(0:1)
 
155
      integer nhel,nsf,ip,im,nh
 
156
 
 
157
      double precision rZero, rHalf, rTwo
 
158
      parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
 
159
 
 
160
c#ifdef HELAS_CHECK
 
161
c      double precision p2
 
162
c      double precision epsi
 
163
c      parameter( epsi = 2.0d-5 )
 
164
c      integer stdo
 
165
c      parameter( stdo = 6 )
 
166
c#endif
 
167
c
 
168
c#ifdef HELAS_CHECK
 
169
c      pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
 
170
c      if ( abs(p(0))+pp.eq.rZero ) then
 
171
c         write(stdo,*)
 
172
c     &        ' helas-error : p(0:3) in ixxxxx is zero momentum'
 
173
c      endif
 
174
c      if ( p(0).le.rZero ) then
 
175
c         write(stdo,*)
 
176
c     &        ' helas-error : p(0:3) in ixxxxx has non-positive energy'
 
177
c         write(stdo,*)
 
178
c     &        '             : p(0) = ',p(0)
 
179
c      endif
 
180
c      p2 = (p(0)-pp)*(p(0)+pp)
 
181
c      if ( abs(p2-fmass**2).gt.p(0)**2*epsi ) then
 
182
c         write(stdo,*)
 
183
c     &        ' helas-error : p(0:3) in ixxxxx has inappropriate mass'
 
184
c         write(stdo,*)
 
185
c     &        '             : p**2 = ',p2,' : fmass**2 = ',fmass**2
 
186
c      endif
 
187
c      if (abs(nhel).ne.1) then
 
188
c         write(stdo,*) ' helas-error : nhel in ixxxxx is not -1,1'
 
189
c         write(stdo,*) '             : nhel = ',nhel
 
190
c      endif
 
191
c      if (abs(nsf).ne.1) then
 
192
c         write(stdo,*) ' helas-error : nsf in ixxxxx is not -1,1'
 
193
c         write(stdo,*) '             : nsf = ',nsf
 
194
c      endif
 
195
c#endif
 
196
 
 
197
c$$$      fi(1) = dcmplx(p(0),p(3))*nsf*-1
 
198
c$$$      fi(2) = dcmplx(p(1),p(2))*nsf*-1
 
199
 
 
200
      nh = nhel*nsf
 
201
 
 
202
      if ( fmass.ne.rZero ) then
 
203
 
 
204
         pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
205
 
 
206
         if ( pp.eq.rZero ) then
 
207
 
 
208
            sqm(0) = dsqrt(abs(fmass)) ! possibility of negative fermion masses
 
209
            sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
 
210
            ip = (1+nh)/2
 
211
            im = (1-nh)/2
 
212
 
 
213
            fi(1) = ip     * sqm(ip)
 
214
            fi(2) = im*nsf * sqm(ip)
 
215
            fi(3) = ip*nsf * sqm(im)
 
216
            fi(4) = im     * sqm(im)
 
217
 
 
218
         else
 
219
 
 
220
            sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
 
221
            sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
 
222
            omega(1) = dsqrt(p(0)+pp)
 
223
            omega(2) = fmass/omega(1)
 
224
            ip = (3+nh)/2
 
225
            im = (3-nh)/2
 
226
            sfomeg(1) = sf(1)*omega(ip)
 
227
            sfomeg(2) = sf(2)*omega(im)
 
228
            pp3 = max(pp+p(3),rZero)
 
229
            chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
 
230
            if ( pp3.eq.rZero ) then
 
231
               chi(2) = dcmplx(-nh )
 
232
            else
 
233
               chi(2) = dcmplx( nh*p(1) , p(2) )/dsqrt(rTwo*pp*pp3)
 
234
            endif
 
235
 
 
236
            fi(1) = sfomeg(1)*chi(im)
 
237
            fi(2) = sfomeg(1)*chi(ip)
 
238
            fi(3) = sfomeg(2)*chi(im)
 
239
            fi(4) = sfomeg(2)*chi(ip)
 
240
 
 
241
         endif
 
242
 
 
243
      else
 
244
 
 
245
         if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
 
246
            sqp0p3 = 0d0
 
247
         else
 
248
            sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsf
 
249
         end if
 
250
         chi(1) = dcmplx( sqp0p3 )
 
251
         if ( sqp0p3.eq.rZero ) then
 
252
            chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
 
253
         else
 
254
            chi(2) = dcmplx( nh*p(1), p(2) )/sqp0p3
 
255
         endif
 
256
         if ( nh.eq.1 ) then
 
257
            fi(1) = dcmplx( rZero )
 
258
            fi(2) = dcmplx( rZero )
 
259
            fi(3) = chi(1)
 
260
            fi(4) = chi(2)
 
261
         else
 
262
            fi(1) = chi(2)
 
263
            fi(2) = chi(1)
 
264
            fi(3) = dcmplx( rZero )
 
265
            fi(4) = dcmplx( rZero )
 
266
         endif
 
267
      endif
 
268
c
 
269
      return
 
270
      end
 
271
 
 
272
 
 
273
      subroutine oxxxxx(p,fmass,nhel,nsf , fo)
 
274
c
 
275
c This subroutine computes a fermion wavefunction with the flowing-OUT
 
276
c fermion number.
 
277
c
 
278
c input:
 
279
c       real    p(0:3)         : four-momentum of fermion
 
280
c       real    fmass          : mass          of fermion
 
281
c       integer nhel = -1 or 1 : helicity      of fermion
 
282
c       integer nsf  = -1 or 1 : +1 for particle, -1 for anti-particle
 
283
c
 
284
c output:
 
285
c       complex fo(6)          : fermion wavefunction               <fo|
 
286
c
 
287
      implicit none
 
288
      double complex fo(6),chi(2)
 
289
      double precision p(0:3),sf(2),sfomeg(2),omega(2),fmass,
 
290
     &     pp,pp3,sqp0p3,sqm(0:1)
 
291
      integer nhel,nsf,nh,ip,im
 
292
 
 
293
      double precision rZero, rHalf, rTwo
 
294
      parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
 
295
 
 
296
c#ifdef HELAS_CHECK
 
297
c      double precision p2
 
298
c      double precision epsi
 
299
c      parameter( epsi = 2.0d-5 )
 
300
c      integer stdo
 
301
c      parameter( stdo = 6 )
 
302
c#endif
 
303
c
 
304
c#ifdef HELAS_CHECK
 
305
c      pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
 
306
c      if ( abs(p(0))+pp.eq.rZero ) then
 
307
c         write(stdo,*)
 
308
c     &        ' helas-error : p(0:3) in oxxxxx is zero momentum'
 
309
c      endif
 
310
c      if ( p(0).le.rZero ) then
 
311
c         write(stdo,*)
 
312
c     &        ' helas-error : p(0:3) in oxxxxx has non-positive energy'
 
313
c         write(stdo,*)
 
314
c     &        '         : p(0) = ',p(0)
 
315
c      endif
 
316
c      p2 = (p(0)-pp)*(p(0)+pp)
 
317
c      if ( abs(p2-fmass**2).gt.p(0)**2*epsi ) then
 
318
c         write(stdo,*)
 
319
c     &        ' helas-error : p(0:3) in oxxxxx has inappropriate mass'
 
320
c         write(stdo,*)
 
321
c     &        '             : p**2 = ',p2,' : fmass**2 = ',fmass**2
 
322
c      endif
 
323
c      if ( abs(nhel).ne.1 ) then
 
324
c         write(stdo,*) ' helas-error : nhel in oxxxxx is not -1,1'
 
325
c         write(stdo,*) '             : nhel = ',nhel
 
326
c      endif
 
327
c      if ( abs(nsf).ne.1 ) then
 
328
c         write(stdo,*) ' helas-error : nsf in oxxxxx is not -1,1'
 
329
c         write(stdo,*) '             : nsf = ',nsf
 
330
c      endif
 
331
c#endif
 
332
 
 
333
      fo(1) = dcmplx(p(0),p(3))*nsf
 
334
      fo(2) = dcmplx(p(1),p(2))*nsf
 
335
 
 
336
      nh = nhel*nsf
 
337
 
 
338
      if ( fmass.ne.rZero ) then
 
339
 
 
340
         pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
341
 
 
342
         if ( pp.eq.rZero ) then
 
343
 
 
344
            sqm(0) = dsqrt(abs(fmass)) ! possibility of negative fermion masses
 
345
            sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
 
346
            ip = -((1+nh)/2)
 
347
            im =  (1-nh)/2
 
348
 
 
349
            fo(3) = im     * sqm(im)
 
350
            fo(4) = ip*nsf * sqm(im)
 
351
            fo(5) = im*nsf * sqm(-ip)
 
352
            fo(6) = ip     * sqm(-ip)
 
353
 
 
354
         else
 
355
 
 
356
            pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
357
            sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
 
358
            sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
 
359
            omega(1) = dsqrt(p(0)+pp)
 
360
            omega(2) = fmass/omega(1)
 
361
            ip = (3+nh)/2
 
362
            im = (3-nh)/2
 
363
            sfomeg(1) = sf(1)*omega(ip)
 
364
            sfomeg(2) = sf(2)*omega(im)
 
365
            pp3 = max(pp+p(3),rZero)
 
366
            chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
 
367
            if ( pp3.eq.rZero ) then
 
368
               chi(2) = dcmplx(-nh )
 
369
            else
 
370
               chi(2) = dcmplx( nh*p(1) , -p(2) )/dsqrt(rTwo*pp*pp3)
 
371
            endif
 
372
 
 
373
            fo(3) = sfomeg(2)*chi(im)
 
374
            fo(4) = sfomeg(2)*chi(ip)
 
375
            fo(5) = sfomeg(1)*chi(im)
 
376
            fo(6) = sfomeg(1)*chi(ip)
 
377
 
 
378
         endif
 
379
 
 
380
      else
 
381
 
 
382
         if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
 
383
            sqp0p3 = 0d0
 
384
         else
 
385
            sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsf
 
386
         end if
 
387
         chi(1) = dcmplx( sqp0p3 )
 
388
         if ( sqp0p3.eq.rZero ) then
 
389
            chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
 
390
         else
 
391
            chi(2) = dcmplx( nh*p(1), -p(2) )/sqp0p3
 
392
         endif
 
393
         if ( nh.eq.1 ) then
 
394
            fo(3) = chi(1)
 
395
            fo(4) = chi(2)
 
396
            fo(5) = dcmplx( rZero )
 
397
            fo(6) = dcmplx( rZero )
 
398
         else
 
399
            fo(3) = dcmplx( rZero )
 
400
            fo(4) = dcmplx( rZero )
 
401
            fo(5) = chi(2)
 
402
            fo(6) = chi(1)
 
403
         endif
 
404
 
 
405
      endif
 
406
c
 
407
      return
 
408
      end
 
409
 
 
410
      subroutine oxxxso(p,fmass,nhel,nsf , fo)
 
411
c Identical to oxxxxx, except that fo returns only the spinor (without the momentum)
 
412
      implicit none
 
413
      double complex fo(4),chi(2)
 
414
      double precision p(0:3),sf(2),sfomeg(2),omega(2),fmass,
 
415
     &     pp,pp3,sqp0p3,sqm(0:1)
 
416
      integer nhel,nsf,nh,ip,im
 
417
 
 
418
      double precision rZero, rHalf, rTwo
 
419
      parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
 
420
 
 
421
c#ifdef HELAS_CHECK
 
422
c      double precision p2
 
423
c      double precision epsi
 
424
c      parameter( epsi = 2.0d-5 )
 
425
c      integer stdo
 
426
c      parameter( stdo = 6 )
 
427
c#endif
 
428
c
 
429
c#ifdef HELAS_CHECK
 
430
c      pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
 
431
c      if ( abs(p(0))+pp.eq.rZero ) then
 
432
c         write(stdo,*)
 
433
c     &        ' helas-error : p(0:3) in oxxxxx is zero momentum'
 
434
c      endif
 
435
c      if ( p(0).le.rZero ) then
 
436
c         write(stdo,*)
 
437
c     &        ' helas-error : p(0:3) in oxxxxx has non-positive energy'
 
438
c         write(stdo,*)
 
439
c     &        '         : p(0) = ',p(0)
 
440
c      endif
 
441
c      p2 = (p(0)-pp)*(p(0)+pp)
 
442
c      if ( abs(p2-fmass**2).gt.p(0)**2*epsi ) then
 
443
c         write(stdo,*)
 
444
c     &        ' helas-error : p(0:3) in oxxxxx has inappropriate mass'
 
445
c         write(stdo,*)
 
446
c     &        '             : p**2 = ',p2,' : fmass**2 = ',fmass**2
 
447
c      endif
 
448
c      if ( abs(nhel).ne.1 ) then
 
449
c         write(stdo,*) ' helas-error : nhel in oxxxxx is not -1,1'
 
450
c         write(stdo,*) '             : nhel = ',nhel
 
451
c      endif
 
452
c      if ( abs(nsf).ne.1 ) then
 
453
c         write(stdo,*) ' helas-error : nsf in oxxxxx is not -1,1'
 
454
c         write(stdo,*) '             : nsf = ',nsf
 
455
c      endif
 
456
c#endif
 
457
 
 
458
c$$$      fo(1) = dcmplx(p(0),p(3))*nsf
 
459
c$$$      fo(2) = dcmplx(p(1),p(2))*nsf
 
460
 
 
461
      nh = nhel*nsf
 
462
 
 
463
      if ( fmass.ne.rZero ) then
 
464
 
 
465
         pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
466
 
 
467
         if ( pp.eq.rZero ) then
 
468
 
 
469
            sqm(0) = dsqrt(abs(fmass)) ! possibility of negative fermion masses
 
470
            sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
 
471
            ip = -((1+nh)/2)
 
472
            im =  (1-nh)/2
 
473
 
 
474
            fo(1) = im     * sqm(im)
 
475
            fo(2) = ip*nsf * sqm(im)
 
476
            fo(3) = im*nsf * sqm(-ip)
 
477
            fo(4) = ip     * sqm(-ip)
 
478
 
 
479
         else
 
480
 
 
481
            pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
482
            sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
 
483
            sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
 
484
            omega(1) = dsqrt(p(0)+pp)
 
485
            omega(2) = fmass/omega(1)
 
486
            ip = (3+nh)/2
 
487
            im = (3-nh)/2
 
488
            sfomeg(1) = sf(1)*omega(ip)
 
489
            sfomeg(2) = sf(2)*omega(im)
 
490
            pp3 = max(pp+p(3),rZero)
 
491
            chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
 
492
            if ( pp3.eq.rZero ) then
 
493
               chi(2) = dcmplx(-nh )
 
494
            else
 
495
               chi(2) = dcmplx( nh*p(1) , -p(2) )/dsqrt(rTwo*pp*pp3)
 
496
            endif
 
497
 
 
498
            fo(1) = sfomeg(2)*chi(im)
 
499
            fo(2) = sfomeg(2)*chi(ip)
 
500
            fo(3) = sfomeg(1)*chi(im)
 
501
            fo(4) = sfomeg(1)*chi(ip)
 
502
 
 
503
         endif
 
504
 
 
505
      else
 
506
 
 
507
         if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
 
508
            sqp0p3 = 0d0
 
509
         else
 
510
            sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsf
 
511
         end if
 
512
         chi(1) = dcmplx( sqp0p3 )
 
513
         if ( sqp0p3.eq.rZero ) then
 
514
            chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
 
515
         else
 
516
            chi(2) = dcmplx( nh*p(1), -p(2) )/sqp0p3
 
517
         endif
 
518
         if ( nh.eq.1 ) then
 
519
            fo(1) = chi(1)
 
520
            fo(2) = chi(2)
 
521
            fo(3) = dcmplx( rZero )
 
522
            fo(4) = dcmplx( rZero )
 
523
         else
 
524
            fo(1) = dcmplx( rZero )
 
525
            fo(2) = dcmplx( rZero )
 
526
            fo(3) = chi(2)
 
527
            fo(4) = chi(1)
 
528
         endif
 
529
 
 
530
      endif
 
531
c
 
532
      return
 
533
      end
 
534
 
 
535
      subroutine pxxxxx(p,tmass,nhel,nst , tc)
 
536
 
 
537
c    CP3 2009.NOV
 
538
 
 
539
c This subroutine computes a PSEUDOR wavefunction.
 
540
c
 
541
c input:
 
542
c       real    p(0:3)         : four-momentum of tensor boson
 
543
c       real    tmass          : mass          of tensor boson
 
544
c       integer nhel           : helicity      of tensor boson
 
545
c                = -2,-1,0,1,2 : (0 is forbidden if tmass=0.0)
 
546
c       integer nst  = -1 or 1 : +1 for final, -1 for initial
 
547
c
 
548
c output:
 
549
c       complex tc(18)         : PSEUDOR  wavefunction    epsilon^mu^nu(t)
 
550
c
 
551
      implicit none
 
552
      double precision p(0:3), tmass
 
553
      integer nhel, nst
 
554
      double complex tc(18)
 
555
 
 
556
      double complex ft(6,4), ep(4), em(4), e0(4)
 
557
      double precision pt, pt2, pp, pzpt, emp, sqh, sqs
 
558
      integer i, j
 
559
 
 
560
      double precision rZero, rHalf, rOne, rTwo
 
561
      parameter( rZero = 0.0d0, rHalf = 0.5d0 )
 
562
      parameter( rOne = 1.0d0, rTwo = 2.0d0 )
 
563
 
 
564
 
 
565
      tc(3)=NHEL
 
566
      tc(1) = dcmplx(p(0),p(3))*nst
 
567
      tc(2) = dcmplx(p(1),p(2))*nst
 
568
 
 
569
      return
 
570
      end
 
571
 
 
572
      subroutine sxxxxx(p,nss , sc)
 
573
c
 
574
c This subroutine computes a complex SCALAR wavefunction.
 
575
c
 
576
c input:
 
577
c       real    p(0:3)         : four-momentum of scalar boson
 
578
c       integer nss  = -1 or 1 : +1 for final, -1 for initial
 
579
c
 
580
c output:
 
581
c       complex sc(3)          : scalar wavefunction                   s
 
582
c
 
583
      implicit none
 
584
      double complex sc(3)
 
585
      double precision p(0:3)
 
586
      integer nss
 
587
 
 
588
      double precision rOne
 
589
      parameter( rOne = 1.0d0 )
 
590
 
 
591
c#ifdef HELAS_CHECK
 
592
c      double precision p2
 
593
c      double precision epsi
 
594
c      parameter( epsi = 2.0d-5 )
 
595
c      double precision rZero
 
596
c      parameter( rZero = 0.0d0 )
 
597
c      integer stdo
 
598
c      parameter( stdo = 6 )
 
599
c#endif
 
600
c
 
601
c#ifdef HELAS_CHECK
 
602
c      if ( abs(p(0))+abs(p(1))+abs(p(2))+abs(p(3)).eq.rZero ) then
 
603
c         write(stdo,*)
 
604
c     &        ' helas-error : p(0:3) in sxxxxx is zero momentum'
 
605
c      endif
 
606
c      if ( p(0).le.rZero ) then
 
607
c         write(stdo,*)
 
608
c     &        ' helas-error : p(0:3) in sxxxxx has non-positive energy'
 
609
c         write(stdo,*)
 
610
c     &        '             : p(0) = ',p(0)
 
611
c      endif
 
612
c      p2 = p(0)**2-(p(1)**2+p(2)**2+p(3)**2)
 
613
c      if ( p2.lt.-p(0)**2*epsi ) then
 
614
c         write(stdo,*) ' helas-error : p(0:3) in sxxxxx is spacelike'
 
615
c         write(stdo,*) '             : p**2 = ',p2
 
616
c      endif
 
617
c      if ( abs(nss).ne.1 ) then
 
618
c         write(stdo,*) ' helas-error : nss in sxxxxx is not -1,1'
 
619
c         write(stdo,*) '             : nss = ',nss
 
620
c      endif
 
621
c#endif
 
622
 
 
623
      sc(3) = dcmplx( rOne )
 
624
      sc(1) = dcmplx(p(0),p(3))*nss
 
625
      sc(2) = dcmplx(p(1),p(2))*nss
 
626
c
 
627
      return
 
628
      end
 
629
 
 
630
      subroutine txxxxx(p,tmass,nhel,nst , tc)
 
631
c
 
632
c This subroutine computes a TENSOR wavefunction.
 
633
c
 
634
c input:
 
635
c       real    p(0:3)         : four-momentum of tensor boson
 
636
c       real    tmass          : mass          of tensor boson
 
637
c       integer nhel           : helicity      of tensor boson
 
638
c                = -2,-1,0,1,2 : (0 is forbidden if tmass=0.0)
 
639
c       integer nst  = -1 or 1 : +1 for final, -1 for initial
 
640
c
 
641
c output:
 
642
c       complex tc(18)         : tensor wavefunction    epsilon^mu^nu(t)
 
643
c
 
644
      implicit none
 
645
      double precision p(0:3), tmass
 
646
      integer nhel, nst
 
647
      double complex tc(18)
 
648
 
 
649
      double complex ft(6,4), ep(4), em(4), e0(4)
 
650
      double precision pt, pt2, pp, pzpt, emp, sqh, sqs
 
651
      integer i, j
 
652
 
 
653
      double precision rZero, rHalf, rOne, rTwo
 
654
      parameter( rZero = 0.0d0, rHalf = 0.5d0 )
 
655
      parameter( rOne = 1.0d0, rTwo = 2.0d0 )
 
656
 
 
657
      integer stdo
 
658
      parameter( stdo = 6 )
 
659
 
 
660
      sqh = sqrt(rHalf)
 
661
      sqs = sqrt(rHalf/3.d0)
 
662
 
 
663
      pt2 = p(1)**2 + p(2)**2
 
664
      pp = min(p(0),sqrt(pt2+p(3)**2))
 
665
      pt = min(pp,sqrt(pt2))
 
666
 
 
667
      ft(5,1) = dcmplx(p(0),p(3))*nst
 
668
      ft(6,1) = dcmplx(p(1),p(2))*nst
 
669
 
 
670
      if ( nhel.ge.0 ) then
 
671
c construct eps+
 
672
         if ( pp.eq.rZero ) then
 
673
            ep(1) = dcmplx( rZero )
 
674
            ep(2) = dcmplx( -sqh )
 
675
            ep(3) = dcmplx( rZero , nst*sqh )
 
676
            ep(4) = dcmplx( rZero )
 
677
         else
 
678
            ep(1) = dcmplx( rZero )
 
679
            ep(4) = dcmplx( pt/pp*sqh )
 
680
            if ( pt.ne.rZero ) then
 
681
               pzpt = p(3)/(pp*pt)*sqh
 
682
               ep(2) = dcmplx( -p(1)*pzpt , -nst*p(2)/pt*sqh )
 
683
               ep(3) = dcmplx( -p(2)*pzpt ,  nst*p(1)/pt*sqh )
 
684
            else
 
685
               ep(2) = dcmplx( -sqh )
 
686
               ep(3) = dcmplx( rZero , nst*sign(sqh,p(3)) )
 
687
            endif
 
688
         endif
 
689
      end if
 
690
 
 
691
      if ( nhel.le.0 ) then
 
692
c construct eps-
 
693
         if ( pp.eq.rZero ) then
 
694
            em(1) = dcmplx( rZero )
 
695
            em(2) = dcmplx( sqh )
 
696
            em(3) = dcmplx( rZero , nst*sqh )
 
697
            em(4) = dcmplx( rZero )
 
698
         else
 
699
            em(1) = dcmplx( rZero )
 
700
            em(4) = dcmplx( -pt/pp*sqh )
 
701
            if ( pt.ne.rZero ) then
 
702
               pzpt = -p(3)/(pp*pt)*sqh
 
703
               em(2) = dcmplx( -p(1)*pzpt , -nst*p(2)/pt*sqh )
 
704
               em(3) = dcmplx( -p(2)*pzpt ,  nst*p(1)/pt*sqh )
 
705
            else
 
706
               em(2) = dcmplx( sqh )
 
707
               em(3) = dcmplx( rZero , nst*sign(sqh,p(3)) )
 
708
            endif
 
709
         endif
 
710
      end if
 
711
 
 
712
      if ( abs(nhel).le.1 ) then
 
713
c construct eps0
 
714
         if ( pp.eq.rZero ) then
 
715
            e0(1) = dcmplx( rZero )
 
716
            e0(2) = dcmplx( rZero )
 
717
            e0(3) = dcmplx( rZero )
 
718
            e0(4) = dcmplx( rOne )
 
719
         else
 
720
            emp = p(0)/(tmass*pp)
 
721
            e0(1) = dcmplx( pp/tmass )
 
722
            e0(4) = dcmplx( p(3)*emp )
 
723
            if ( pt.ne.rZero ) then
 
724
               e0(2) = dcmplx( p(1)*emp )
 
725
               e0(3) = dcmplx( p(2)*emp )
 
726
            else
 
727
               e0(2) = dcmplx( rZero )
 
728
               e0(3) = dcmplx( rZero )
 
729
            endif
 
730
         end if
 
731
      end if
 
732
 
 
733
      if ( nhel.eq.2 ) then
 
734
         do j = 1,4
 
735
            do i = 1,4
 
736
               ft(i,j) = ep(i)*ep(j)
 
737
            end do
 
738
         end do
 
739
      else if ( nhel.eq.-2 ) then
 
740
         do j = 1,4
 
741
            do i = 1,4
 
742
               ft(i,j) = em(i)*em(j)
 
743
            end do
 
744
         end do
 
745
      else if (tmass.eq.0) then
 
746
         do j = 1,4
 
747
            do i = 1,4
 
748
               ft(i,j) = 0
 
749
            end do
 
750
         end do
 
751
      else if (tmass.ne.0) then
 
752
        if  ( nhel.eq.1 ) then
 
753
           do j = 1,4
 
754
              do i = 1,4
 
755
                 ft(i,j) = sqh*( ep(i)*e0(j) + e0(i)*ep(j) )
 
756
              end do
 
757
           end do
 
758
        else if ( nhel.eq.0 ) then
 
759
           do j = 1,4
 
760
              do i = 1,4
 
761
                 ft(i,j) = sqs*( ep(i)*em(j) + em(i)*ep(j)
 
762
     &                                + rTwo*e0(i)*e0(j) )
 
763
              end do
 
764
           end do
 
765
        else if ( nhel.eq.-1 ) then
 
766
           do j = 1,4
 
767
              do i = 1,4
 
768
                 ft(i,j) = sqh*( em(i)*e0(j) + e0(i)*em(j) )
 
769
              end do
 
770
           end do
 
771
        else
 
772
           write(stdo,*) 'invalid helicity in TXXXXX'
 
773
           stop
 
774
        end if
 
775
      end if
 
776
 
 
777
      tc(3) = ft(1,1)
 
778
      tc(4) = ft(1,2)
 
779
      tc(5) = ft(1,3)
 
780
      tc(6) = ft(1,4)
 
781
      tc(7) = ft(2,1)
 
782
      tc(8) = ft(2,2)
 
783
      tc(9) = ft(2,3)
 
784
      tc(10) = ft(2,4)
 
785
      tc(11) = ft(3,1)
 
786
      tc(12) = ft(3,2)
 
787
      tc(13) = ft(3,3)
 
788
      tc(14) = ft(3,4)
 
789
      tc(15) = ft(4,1)
 
790
      tc(16) = ft(4,2)
 
791
      tc(17) = ft(4,3)
 
792
      tc(18) = ft(4,4)
 
793
      tc(1) = ft(5,1)
 
794
      tc(2) = ft(6,1)
 
795
 
 
796
      return
 
797
      end
 
798
 
 
799
 
 
800
      subroutine vxxxxx(p,vmass,nhel,nsv , vc)
 
801
c
 
802
c This subroutine computes a VECTOR wavefunction.
 
803
c
 
804
c input:
 
805
c       real    p(0:3)         : four-momentum of vector boson
 
806
c       real    vmass          : mass          of vector boson
 
807
c       integer nhel = -1, 0, 1: helicity      of vector boson
 
808
c                                (0 is forbidden if vmass=0.0)
 
809
c       integer nsv  = -1 or 1 : +1 for final, -1 for initial
 
810
c
 
811
c output:
 
812
c       complex vc(6)          : vector wavefunction       epsilon^mu(v)
 
813
c
 
814
      implicit none
 
815
      double complex vc(6)
 
816
      double precision p(0:3),vmass,hel,hel0,pt,pt2,pp,pzpt,emp,sqh
 
817
      integer nhel,nsv,nsvahl
 
818
 
 
819
      double precision rZero, rHalf, rOne, rTwo
 
820
      parameter( rZero = 0.0d0, rHalf = 0.5d0 )
 
821
      parameter( rOne = 1.0d0, rTwo = 2.0d0 )
 
822
 
 
823
c#ifdef HELAS_CHECK
 
824
c      double precision p2
 
825
c      double precision epsi
 
826
c      parameter( epsi = 2.0d-5 )
 
827
c      integer stdo
 
828
c      parameter( stdo = 6 )
 
829
c#endif
 
830
c
 
831
c#ifdef HELAS_CHECK
 
832
c      pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
 
833
c      if ( abs(p(0))+pp.eq.rZero ) then
 
834
c         write(stdo,*)
 
835
c     &        ' helas-error : p(0:3) in vxxxxx is zero momentum'
 
836
c      endif
 
837
c      if ( p(0).le.rZero ) then
 
838
c         write(stdo,*)
 
839
c     &        ' helas-error : p(0:3) in vxxxxx has non-positive energy'
 
840
c         write(stdo,*)
 
841
c     &        '             : p(0) = ',p(0)
 
842
c      endif
 
843
c      p2 = (p(0)+pp)*(p(0)-pp)
 
844
c      if ( abs(p2-vmass**2).gt.p(0)**2*2.e-5 ) then
 
845
c         write(stdo,*)
 
846
c     &        ' helas-error : p(0:3) in vxxxxx has inappropriate mass'
 
847
c         write(stdo,*)
 
848
c     &        '             : p**2 = ',p2,' : vmass**2 = ',vmass**2
 
849
c      endif
 
850
c      if ( vmass.ne.rZero ) then
 
851
c         if ( abs(nhel).gt.1 ) then
 
852
c            write(stdo,*) ' helas-error : nhel in vxxxxx is not -1,0,1'
 
853
c            write(stdo,*) '             : nhel = ',nhel
 
854
c         endif
 
855
c      else
 
856
c         if ( abs(nhel).ne.1 ) then
 
857
c            write(stdo,*) ' helas-error : nhel in vxxxxx is not -1,1'
 
858
c            write(stdo,*) '             : nhel = ',nhel
 
859
c         endif
 
860
c      endif
 
861
c      if ( abs(nsv).ne.1 ) then
 
862
c         write(stdo,*) ' helas-error : nsv in vmxxxx is not -1,1'
 
863
c         write(stdo,*) '             : nsv = ',nsv
 
864
c      endif
 
865
c#endif
 
866
 
 
867
      sqh = dsqrt(rHalf)
 
868
      hel = dble(nhel)
 
869
      nsvahl = nsv*dabs(hel)
 
870
      pt2 = p(1)**2+p(2)**2
 
871
      pp = min(p(0),dsqrt(pt2+p(3)**2))
 
872
      pt = min(pp,dsqrt(pt2))
 
873
 
 
874
      vc(1) = dcmplx(p(0),p(3))*nsv
 
875
      vc(2) = dcmplx(p(1),p(2))*nsv
 
876
 
 
877
c#ifdef HELAS_CHECK
 
878
c nhel=4 option for scalar polarization
 
879
c      if( nhel.eq.4 ) then
 
880
c         if( vmass.eq.rZero ) then
 
881
c            vc(1) = rOne
 
882
c            vc(2) = p(1)/p(0)
 
883
c            vc(3) = p(2)/p(0)
 
884
c            vc(4) = p(3)/p(0)
 
885
c         else
 
886
c            vc(1) = p(0)/vmass
 
887
c            vc(2) = p(1)/vmass
 
888
c            vc(3) = p(2)/vmass
 
889
c            vc(4) = p(3)/vmass
 
890
c         endif
 
891
c         return
 
892
c      endif
 
893
c#endif
 
894
 
 
895
      if ( vmass.ne.rZero ) then
 
896
 
 
897
         hel0 = rOne-dabs(hel)
 
898
 
 
899
         if ( pp.eq.rZero ) then
 
900
 
 
901
            vc(3) = dcmplx( rZero )
 
902
            vc(4) = dcmplx(-hel*sqh )
 
903
            vc(5) = dcmplx( rZero , nsvahl*sqh )
 
904
            vc(6) = dcmplx( hel0 )
 
905
 
 
906
         else
 
907
 
 
908
            emp = p(0)/(vmass*pp)
 
909
            vc(3) = dcmplx( hel0*pp/vmass )
 
910
            vc(6) = dcmplx( hel0*p(3)*emp+hel*pt/pp*sqh )
 
911
            if ( pt.ne.rZero ) then
 
912
               pzpt = p(3)/(pp*pt)*sqh*hel
 
913
               vc(4) = dcmplx( hel0*p(1)*emp-p(1)*pzpt ,
 
914
     &                         -nsvahl*p(2)/pt*sqh       )
 
915
               vc(5) = dcmplx( hel0*p(2)*emp-p(2)*pzpt ,
 
916
     &                          nsvahl*p(1)/pt*sqh       )
 
917
            else
 
918
               vc(4) = dcmplx( -hel*sqh )
 
919
               vc(5) = dcmplx( rZero , nsvahl*sign(sqh,p(3)) )
 
920
            endif
 
921
 
 
922
         endif
 
923
 
 
924
      else
 
925
 
 
926
         pp = p(0)
 
927
         pt = sqrt(p(1)**2+p(2)**2)
 
928
         vc(3) = dcmplx( rZero )
 
929
         vc(6) = dcmplx( hel*pt/pp*sqh )
 
930
         if ( pt.ne.rZero ) then
 
931
            pzpt = p(3)/(pp*pt)*sqh*hel
 
932
            vc(4) = dcmplx( -p(1)*pzpt , -nsv*p(2)/pt*sqh )
 
933
            vc(5) = dcmplx( -p(2)*pzpt ,  nsv*p(1)/pt*sqh )
 
934
         else
 
935
            vc(4) = dcmplx( -hel*sqh )
 
936
            vc(5) = dcmplx( rZero , nsv*sign(sqh,p(3)) )
 
937
         endif
 
938
 
 
939
      endif
 
940
c
 
941
      return
 
942
      end
 
943
 
 
944
      subroutine boostx(p,q , pboost)
 
945
c
 
946
c This subroutine performs the Lorentz boost of a four-momentum.  The
 
947
c momentum p is assumed to be given in the rest frame of q.  pboost is
 
948
c the momentum p boosted to the frame in which q is given.  q must be a
 
949
c timelike momentum.
 
950
c
 
951
c input:
 
952
c       real    p(0:3)         : four-momentum p in the q rest  frame
 
953
c       real    q(0:3)         : four-momentum q in the boosted frame
 
954
c
 
955
c output:
 
956
c       real    pboost(0:3)    : four-momentum p in the boosted frame
 
957
c
 
958
      implicit none
 
959
      double precision p(0:3),q(0:3),pboost(0:3),pq,qq,m,lf
 
960
 
 
961
      double precision rZero
 
962
      parameter( rZero = 0.0d0 )
 
963
 
 
964
c#ifdef HELAS_CHECK
 
965
c      integer stdo
 
966
c      parameter( stdo = 6 )
 
967
c      double precision pp
 
968
c#endif
 
969
c
 
970
      qq = q(1)**2+q(2)**2+q(3)**2
 
971
 
 
972
c#ifdef HELAS_CHECK
 
973
c      if (abs(p(0))+abs(p(1))+abs(p(2))+abs(p(3)).eq.rZero) then
 
974
c         write(stdo,*)
 
975
c     &        ' helas-error : p(0:3) in boostx is zero momentum'
 
976
c      endif
 
977
c      if (abs(q(0))+qq.eq.rZero) then
 
978
c         write(stdo,*)
 
979
c     &        ' helas-error : q(0:3) in boostx is zero momentum'
 
980
c      endif
 
981
c      if (p(0).le.rZero) then
 
982
c         write(stdo,*)
 
983
c     &        ' helas-warn  : p(0:3) in boostx has not positive energy'
 
984
c         write(stdo,*)
 
985
c     &        '             : p(0) = ',p(0)
 
986
c      endif
 
987
c      if (q(0).le.rZero) then
 
988
c         write(stdo,*)
 
989
c     &        ' helas-error : q(0:3) in boostx has not positive energy'
 
990
c         write(stdo,*)
 
991
c     &        '             : q(0) = ',q(0)
 
992
c      endif
 
993
c      pp=p(0)**2-p(1)**2-p(2)**2-p(3)**2
 
994
c      if (pp.lt.rZero) then
 
995
c         write(stdo,*)
 
996
c     &        ' helas-warn  : p(0:3) in boostx is spacelike'
 
997
c         write(stdo,*)
 
998
c     &        '             : p**2 = ',pp
 
999
c      endif
 
1000
c      if (q(0)**2-qq.le.rZero) then
 
1001
c         write(stdo,*)
 
1002
c     &        ' helas-error : q(0:3) in boostx is not timelike'
 
1003
c         write(stdo,*)
 
1004
c     &        '             : q**2 = ',q(0)**2-qq
 
1005
c      endif
 
1006
c      if (qq.eq.rZero) then
 
1007
c         write(stdo,*)
 
1008
c     &   ' helas-warn  : q(0:3) in boostx has zero spacial components'
 
1009
c      endif
 
1010
c#endif
 
1011
 
 
1012
      if ( qq.ne.rZero ) then
 
1013
         pq = p(1)*q(1)+p(2)*q(2)+p(3)*q(3)
 
1014
         m = sqrt(max(q(0)**2-qq,1d-99))
 
1015
         lf = ((q(0)-m)*pq/qq+p(0))/m
 
1016
         pboost(0) = (p(0)*q(0)+pq)/m
 
1017
         pboost(1) =  p(1)+q(1)*lf
 
1018
         pboost(2) =  p(2)+q(2)*lf
 
1019
         pboost(3) =  p(3)+q(3)*lf
 
1020
      else
 
1021
         pboost(0) = p(0)
 
1022
         pboost(1) = p(1)
 
1023
         pboost(2) = p(2)
 
1024
         pboost(3) = p(3)
 
1025
      endif
 
1026
c
 
1027
      return
 
1028
      end
 
1029
 
 
1030
      subroutine momntx(energy,mass,costh,phi , p)
 
1031
c
 
1032
c This subroutine sets up a four-momentum from the four inputs.
 
1033
c
 
1034
c input:
 
1035
c       real    energy         : energy
 
1036
c       real    mass           : mass
 
1037
c       real    costh          : cos(theta)
 
1038
c       real    phi            : azimuthal angle
 
1039
c
 
1040
c output:
 
1041
c       real    p(0:3)         : four-momentum
 
1042
c
 
1043
      implicit none
 
1044
      double precision p(0:3),energy,mass,costh,phi,pp,sinth
 
1045
 
 
1046
      double precision rZero, rOne
 
1047
      parameter( rZero = 0.0d0, rOne = 1.0d0 )
 
1048
 
 
1049
c#ifdef HELAS_CHECK
 
1050
c      double precision rPi, rTwo
 
1051
c      parameter( rPi = 3.14159265358979323846d0, rTwo = 2.d0 )
 
1052
c      integer stdo
 
1053
c      parameter( stdo = 6 )
 
1054
c#endif
 
1055
c
 
1056
c#ifdef HELAS_CHECK
 
1057
c      if (energy.lt.mass) then
 
1058
c         write(stdo,*)
 
1059
c     &        ' helas-error : energy in momntx is less than mass'
 
1060
c         write(stdo,*)
 
1061
c     &        '             : energy = ',energy,' : mass = ',mass
 
1062
c      endif
 
1063
c      if (mass.lt.rZero) then
 
1064
c         write(stdo,*) ' helas-error : mass in momntx is negative'
 
1065
c         write(stdo,*) '             : mass = ',mass
 
1066
c      endif
 
1067
c      if (abs(costh).gt.rOne) then
 
1068
c         write(stdo,*) ' helas-error : costh in momntx is improper'
 
1069
c         write(stdo,*) '             : costh = ',costh
 
1070
c      endif
 
1071
c      if (phi.lt.rZero .or. phi.gt.rTwo*rPi) then
 
1072
c         write(stdo,*)
 
1073
c     &   ' helas-warn  : phi in momntx does not lie on 0.0 thru 2.0*pi'
 
1074
c         write(stdo,*)
 
1075
c     &   '             : phi = ',phi
 
1076
c      endif
 
1077
c#endif
 
1078
 
 
1079
      p(0) = energy
 
1080
 
 
1081
      if ( energy.eq.mass ) then
 
1082
 
 
1083
         p(1) = rZero
 
1084
         p(2) = rZero
 
1085
         p(3) = rZero
 
1086
 
 
1087
      else
 
1088
 
 
1089
         pp = sqrt((energy-mass)*(energy+mass))
 
1090
         sinth = sqrt((rOne-costh)*(rOne+costh))
 
1091
         p(3) = pp*costh
 
1092
         if ( phi.eq.rZero ) then
 
1093
            p(1) = pp*sinth
 
1094
            p(2) = rZero
 
1095
         else
 
1096
            p(1) = pp*sinth*cos(phi)
 
1097
            p(2) = pp*sinth*sin(phi)
 
1098
         endif
 
1099
 
 
1100
      endif
 
1101
c
 
1102
      return
 
1103
      end
 
1104
      subroutine rotxxx(p,q , prot)
 
1105
c
 
1106
c This subroutine performs the spacial rotation of a four-momentum.
 
1107
c the momentum p is assumed to be given in the frame where the spacial
 
1108
c component of q points the positive z-axis.  prot is the momentum p
 
1109
c rotated to the frame where q is given.
 
1110
c
 
1111
c input:
 
1112
c       real    p(0:3)         : four-momentum p in q(1)=q(2)=0 frame
 
1113
c       real    q(0:3)         : four-momentum q in the rotated frame
 
1114
c
 
1115
c output:
 
1116
c       real    prot(0:3)      : four-momentum p in the rotated frame
 
1117
c
 
1118
      implicit none
 
1119
      double precision p(0:3),q(0:3),prot(0:3),qt2,qt,psgn,qq,p1
 
1120
 
 
1121
      double precision rZero, rOne
 
1122
      parameter( rZero = 0.0d0, rOne = 1.0d0 )
 
1123
 
 
1124
c#ifdef HELAS_CHECK
 
1125
c      integer stdo
 
1126
c      parameter( stdo = 6 )
 
1127
c#endif
 
1128
c
 
1129
      prot(0) = p(0)
 
1130
 
 
1131
      qt2 = q(1)**2 + q(2)**2
 
1132
 
 
1133
c#ifdef HELAS_CHECK
 
1134
c      if (abs(p(0))+abs(p(1))+abs(p(2))+abs(p(3)).eq.rZero) then
 
1135
c         write(stdo,*)
 
1136
c     &        ' helas-error : p(0:3) in rotxxx is zero momentum'
 
1137
c      endif
 
1138
c      if (abs(q(0))+abs(q(3))+qt2.eq.rZero) then
 
1139
c         write(stdo,*)
 
1140
c     &        ' helas-error : q(0:3) in rotxxx is zero momentum'
 
1141
c      endif
 
1142
c      if (qt2+abs(q(3)).eq.rZero) then
 
1143
c         write(stdo,*)
 
1144
c     &   ' helas-warn  : q(0:3) in rotxxx has zero spacial momentum'
 
1145
c      endif
 
1146
c#endif
 
1147
 
 
1148
      if ( qt2.eq.rZero ) then
 
1149
          if ( q(3).eq.rZero ) then
 
1150
             prot(1) = p(1)
 
1151
             prot(2) = p(2)
 
1152
             prot(3) = p(3)
 
1153
          else
 
1154
             psgn = dsign(rOne,q(3))
 
1155
             prot(1) = p(1)*psgn
 
1156
             prot(2) = p(2)*psgn
 
1157
             prot(3) = p(3)*psgn
 
1158
          endif
 
1159
      else
 
1160
          qq = sqrt(qt2+q(3)**2)
 
1161
          qt = sqrt(qt2)
 
1162
          p1 = p(1)
 
1163
          prot(1) = q(1)*q(3)/qq/qt*p1 -q(2)/qt*p(2) +q(1)/qq*p(3)
 
1164
          prot(2) = q(2)*q(3)/qq/qt*p1 +q(1)/qt*p(2) +q(2)/qq*p(3)
 
1165
          prot(3) =          -qt/qq*p1               +q(3)/qq*p(3)
 
1166
      endif
 
1167
c
 
1168
      return
 
1169
      end
 
1170
 
 
1171
      subroutine mom2cx(esum,mass1,mass2,costh1,phi1 , p1,p2)
 
1172
c
 
1173
c This subroutine sets up two four-momenta in the two particle rest
 
1174
c frame.
 
1175
c
 
1176
c input:
 
1177
c       real    esum           : energy sum of particle 1 and 2
 
1178
c       real    mass1          : mass            of particle 1
 
1179
c       real    mass2          : mass            of particle 2
 
1180
c       real    costh1         : cos(theta)      of particle 1
 
1181
c       real    phi1           : azimuthal angle of particle 1
 
1182
c
 
1183
c output:
 
1184
c       real    p1(0:3)        : four-momentum of particle 1
 
1185
c       real    p2(0:3)        : four-momentum of particle 2
 
1186
c     
 
1187
      implicit none
 
1188
      double precision p1(0:3),p2(0:3),
 
1189
     &     esum,mass1,mass2,costh1,phi1,md2,ed,pp,sinth1
 
1190
 
 
1191
      double precision rZero, rHalf, rOne, rTwo
 
1192
      parameter( rZero = 0.0d0, rHalf = 0.5d0 )
 
1193
      parameter( rOne = 1.0d0, rTwo = 2.0d0 )
 
1194
 
 
1195
c#ifdef HELAS_CHECK
 
1196
c      double precision rPi
 
1197
c      parameter( rPi = 3.14159265358979323846d0 )
 
1198
c      integer stdo
 
1199
c      parameter( stdo = 6 )
 
1200
c#endif
 
1201
cc
 
1202
c#ifdef HELAS_CHECK
 
1203
c      if (esum.lt.mass1+mass2) then
 
1204
c         write(stdo,*)
 
1205
c     &        ' helas-error : esum in mom2cx is less than mass1+mass2'
 
1206
c         write(stdo,*)
 
1207
c     &        '             : esum = ',esum,
 
1208
c     &        ' : mass1+mass2 = ',mass1,mass2
 
1209
c      endif
 
1210
c      if (mass1.lt.rZero) then
 
1211
c         write(stdo,*) ' helas-error : mass1 in mom2cx is negative'
 
1212
c         write(stdo,*) '             : mass1 = ',mass1
 
1213
c      endif
 
1214
c      if (mass2.lt.rZero) then
 
1215
c         write(stdo,*) ' helas-error : mass2 in mom2cx is negative'
 
1216
c         write(stdo,*) '             : mass2 = ',mass2
 
1217
c      endif
 
1218
c      if (abs(costh1).gt.1.) then
 
1219
c         write(stdo,*) ' helas-error : costh1 in mom2cx is improper'
 
1220
c         write(stdo,*) '             : costh1 = ',costh1
 
1221
c      endif
 
1222
c      if (phi1.lt.rZero .or. phi1.gt.rTwo*rPi) then
 
1223
c         write(stdo,*)
 
1224
c     &   ' helas-warn  : phi1 in mom2cx does not lie on 0.0 thru 2.0*pi'
 
1225
c         write(stdo,*) 
 
1226
c     &   '             : phi1 = ',phi1
 
1227
c      endif
 
1228
c#endif
 
1229
 
 
1230
      md2 = (mass1-mass2)*(mass1+mass2)
 
1231
      ed = md2/esum
 
1232
      if ( mass1*mass2.eq.rZero ) then
 
1233
         pp = (esum-abs(ed))*rHalf
 
1234
      else
 
1235
         pp = sqrt(max((md2/esum)**2-rTwo*(mass1**2+mass2**2)+esum**2,1d-99))*rHalf
 
1236
      endif
 
1237
      sinth1 = sqrt((rOne-costh1)*(rOne+costh1))
 
1238
 
 
1239
      p1(0) = max((esum+ed)*rHalf,rZero)
 
1240
      p1(1) = pp*sinth1*cos(phi1)
 
1241
      p1(2) = pp*sinth1*sin(phi1)
 
1242
      p1(3) = pp*costh1
 
1243
 
 
1244
      p2(0) = max((esum-ed)*rHalf,rZero)
 
1245
      p2(1) = -p1(1)
 
1246
      p2(2) = -p1(2)
 
1247
      p2(3) = -p1(3)
 
1248
c
 
1249
      return
 
1250
      end
 
1251
      subroutine irxxxx(p,rmass,nhel,nsr , ri)
 
1252
c
 
1253
c This subroutine computes a Rarita-Schwinger wavefunction of spin-3/2
 
1254
c fermion with the flowing-IN fermion number.
 
1255
c
 
1256
c input:
 
1257
c       real    p(0:3)           : four-momentum of RS fermion
 
1258
c       real    rmass            : mass          of RS fermion
 
1259
c       integer nhel = -3,-1,1,3 : helicity      of RS fermion
 
1260
c                                  (1- and 1 is forbidden if rmass = 0)
 
1261
c       integer nsr  = -1 or 1   : +1 for particle, -1 for anti-particle
 
1262
c
 
1263
c output:
 
1264
c       complex ri(18)           : RS fermion wavefunction         |ri>v   
 
1265
c     
 
1266
c- by K.Mawatari - 2008/02/26
 
1267
c
 
1268
      implicit none
 
1269
      double precision p(0:3),rmass
 
1270
      integer nhel,nsr
 
1271
      double complex ri(18)
 
1272
 
 
1273
      double complex rc(6,4),ep(4),em(4),e0(4),fip(4),fim(4),chi(2)
 
1274
      double precision pp,pt2,pt,pzpt,emp, sf(2),sfomeg(2),omega(2),pp3,
 
1275
     &                 sqp0p3,sqm      
 
1276
      integer i,j,nsv,ip,im,nh
 
1277
 
 
1278
      double precision rZero, rHalf, rOne, rTwo, rThree, sqh,sq2,sq3
 
1279
      parameter( rZero = 0.0d0, rHalf = 0.5d0 )
 
1280
      parameter( rOne = 1.0d0, rTwo = 2.0d0, rThree = 3.0d0 )
 
1281
 
 
1282
c#ifdef HELAS_CHECK
 
1283
c      double precision p2
 
1284
c      double precision epsi
 
1285
c      parameter( epsi = 2.0d-5 )
 
1286
c      integer stdo
 
1287
c      parameter( stdo = 6 )
 
1288
c#endif
 
1289
c
 
1290
c#ifdef HELAS_CHECK
 
1291
c      pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
 
1292
c      if ( abs(p(0))+pp.eq.rZero ) then
 
1293
c         write(stdo,*)
 
1294
c     &        ' helas-error : p(0:3) in irxxxx is zero momentum'
 
1295
c      endif
 
1296
c      if ( p(0).le.rZero ) then
 
1297
c         write(stdo,*)
 
1298
c     &        ' helas-error : p(0:3) in irxxxx has non-positive energy'
 
1299
c         write(stdo,*)
 
1300
c     &        '             : p(0) = ',p(0)
 
1301
c      endif
 
1302
c      p2 = (p(0)-pp)*(p(0)+pp)
 
1303
c      if ( abs(p2-rmass**2).gt.p(0)**2*epsi ) then
 
1304
c         write(stdo,*)
 
1305
c     &        ' helas-error : p(0:3) in irxxxx has inappropriate mass'
 
1306
c         write(stdo,*)
 
1307
c     &        '             : p**2 = ',p2,' : rmass**2 = ',rmass**2
 
1308
c      endif
 
1309
c      if (abs(nhel).gt.3 .or. abs(nhel).eq.2 .or. abs(nhel).eq.0 ) then
 
1310
c         write(stdo,*) ' helas-error : nhel in irxxxx is not -3,-1,1,3'
 
1311
c         write(stdo,*) '             : nhel = ',nhel
 
1312
c      endif
 
1313
c      if (abs(nsr).ne.1) then
 
1314
c         write(stdo,*) ' helas-error : nsr in irxxxx is not -1,1'
 
1315
c         write(stdo,*) '             : nsr = ',nsr
 
1316
c      endif
 
1317
c#endif
 
1318
 
 
1319
      sqh = sqrt(rHalf)
 
1320
      sq2 = sqrt(rTwo)
 
1321
      sq3 = sqrt(rThree)
 
1322
 
 
1323
      pt2 = p(1)**2 + p(2)**2
 
1324
      pp = min(p(0),sqrt(pt2+p(3)**2))
 
1325
      pt = min(pp,sqrt(pt2))
 
1326
 
 
1327
      rc(5,1) = -1*dcmplx(p(0),p(3))*nsr
 
1328
      rc(6,1) = -1*dcmplx(p(1),p(2))*nsr
 
1329
 
 
1330
      nsv = -nsr ! nsv=+1 for final, -1 for initial
 
1331
 
 
1332
      if ( nhel.ge.1 ) then 
 
1333
c construct eps+
 
1334
         if ( pp.eq.rZero ) then
 
1335
            ep(1) = dcmplx( rZero )
 
1336
            ep(2) = dcmplx( -sqh )
 
1337
            ep(3) = dcmplx( rZero , nsv*sqh )
 
1338
            ep(4) = dcmplx( rZero )
 
1339
         else
 
1340
            ep(1) = dcmplx( rZero )
 
1341
            ep(4) = dcmplx( pt/pp*sqh )
 
1342
            if ( pt.ne.rZero ) then
 
1343
               pzpt = p(3)/(pp*pt)*sqh
 
1344
               ep(2) = dcmplx( -p(1)*pzpt , -nsv*p(2)/pt*sqh )
 
1345
               ep(3) = dcmplx( -p(2)*pzpt ,  nsv*p(1)/pt*sqh )
 
1346
            else
 
1347
               ep(2) = dcmplx( -sqh )
 
1348
               ep(3) = dcmplx( rZero , nsv*sign(sqh,p(3)) )
 
1349
            endif
 
1350
         endif
 
1351
      end if
 
1352
 
 
1353
      if ( nhel.le.-1 ) then 
 
1354
c construct eps-
 
1355
         if ( pp.eq.rZero ) then
 
1356
            em(1) = dcmplx( rZero )
 
1357
            em(2) = dcmplx( sqh )
 
1358
            em(3) = dcmplx( rZero , nsv*sqh )
 
1359
            em(4) = dcmplx( rZero )
 
1360
         else
 
1361
            em(1) = dcmplx( rZero )
 
1362
            em(4) = dcmplx( -pt/pp*sqh )
 
1363
            if ( pt.ne.rZero ) then
 
1364
               pzpt = -p(3)/(pp*pt)*sqh
 
1365
               em(2) = dcmplx( -p(1)*pzpt , -nsv*p(2)/pt*sqh )
 
1366
               em(3) = dcmplx( -p(2)*pzpt ,  nsv*p(1)/pt*sqh )
 
1367
            else
 
1368
               em(2) = dcmplx( sqh )
 
1369
               em(3) = dcmplx( rZero , nsv*sign(sqh,p(3)) )
 
1370
            endif
 
1371
         endif
 
1372
      end if
 
1373
 
 
1374
      if ( abs(nhel).le.1 ) then  
 
1375
c construct eps0
 
1376
         if ( pp.eq.rZero ) then
 
1377
            e0(1) = dcmplx( rZero )
 
1378
            e0(2) = dcmplx( rZero )
 
1379
            e0(3) = dcmplx( rZero )
 
1380
            e0(4) = dcmplx( rOne )
 
1381
         else
 
1382
            emp = p(0)/(rmass*pp)
 
1383
            e0(1) = dcmplx( pp/rmass )
 
1384
            e0(4) = dcmplx( p(3)*emp )
 
1385
            if ( pt.ne.rZero ) then
 
1386
               e0(2) = dcmplx( p(1)*emp )
 
1387
               e0(3) = dcmplx( p(2)*emp )
 
1388
            else
 
1389
               e0(2) = dcmplx( rZero )
 
1390
               e0(3) = dcmplx( rZero )
 
1391
            endif
 
1392
         end if
 
1393
      end if
 
1394
 
 
1395
      if ( nhel.ge.-1 ) then
 
1396
c constract spinor+ 
 
1397
       nh = nsr
 
1398
       if ( rmass.ne.rZero ) then
 
1399
         pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
1400
         if ( pp.eq.rZero ) then
 
1401
            sqm = dsqrt(rmass)
 
1402
            ip = (1+nh)/2
 
1403
            im = (1-nh)/2
 
1404
              fip(1) = ip     * sqm
 
1405
              fip(2) = im*nsr * sqm
 
1406
              fip(3) = ip*nsr * sqm
 
1407
              fip(4) = im     * sqm
 
1408
         else
 
1409
            sf(1) = dble(1+nsr+(1-nsr)*nh)*rHalf
 
1410
            sf(2) = dble(1+nsr-(1-nsr)*nh)*rHalf
 
1411
            omega(1) = dsqrt(p(0)+pp)
 
1412
            omega(2) = rmass/omega(1)
 
1413
            ip = (3+nh)/2
 
1414
            im = (3-nh)/2
 
1415
            sfomeg(1) = sf(1)*omega(ip)
 
1416
            sfomeg(2) = sf(2)*omega(im)
 
1417
            pp3 = max(pp+p(3),rZero)
 
1418
            chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
 
1419
            if ( pp3.eq.rZero ) then
 
1420
               chi(2) = dcmplx(-nh )
 
1421
            else
 
1422
               chi(2) = dcmplx( nh*p(1) , p(2) )/dsqrt(rTwo*pp*pp3)
 
1423
            endif
 
1424
              fip(1) = sfomeg(1)*chi(im)
 
1425
              fip(2) = sfomeg(1)*chi(ip)
 
1426
              fip(3) = sfomeg(2)*chi(im)
 
1427
              fip(4) = sfomeg(2)*chi(ip)
 
1428
         endif
 
1429
       else
 
1430
         sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsr
 
1431
         chi(1) = dcmplx( sqp0p3 )
 
1432
         if ( sqp0p3.eq.rZero ) then
 
1433
            chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
 
1434
         else
 
1435
            chi(2) = dcmplx( nh*p(1), p(2) )/sqp0p3
 
1436
         endif
 
1437
         if ( nh.eq.1 ) then
 
1438
            fip(1) = dcmplx( rZero )
 
1439
            fip(2) = dcmplx( rZero )
 
1440
            fip(3) = chi(1)
 
1441
            fip(4) = chi(2)
 
1442
         else
 
1443
            fip(1) = chi(2)
 
1444
            fip(2) = chi(1)
 
1445
            fip(3) = dcmplx( rZero )
 
1446
            fip(4) = dcmplx( rZero )
 
1447
         endif
 
1448
       endif
 
1449
      end if
 
1450
 
 
1451
      if ( nhel.le.1 ) then
 
1452
c constract spinor-
 
1453
       nh = -nsr
 
1454
       if ( rmass.ne.rZero ) then
 
1455
         pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
1456
         if ( pp.eq.rZero ) then
 
1457
            sqm = dsqrt(rmass)
 
1458
            ip = (1+nh)/2
 
1459
            im = (1-nh)/2
 
1460
              fim(1) = ip     * sqm
 
1461
              fim(2) = im*nsr * sqm
 
1462
              fim(3) = ip*nsr * sqm
 
1463
              fim(4) = im     * sqm
 
1464
         else
 
1465
            sf(1) = dble(1+nsr+(1-nsr)*nh)*rHalf
 
1466
            sf(2) = dble(1+nsr-(1-nsr)*nh)*rHalf
 
1467
            omega(1) = dsqrt(p(0)+pp)
 
1468
            omega(2) = rmass/omega(1)
 
1469
            ip = (3+nh)/2
 
1470
            im = (3-nh)/2
 
1471
            sfomeg(1) = sf(1)*omega(ip)
 
1472
            sfomeg(2) = sf(2)*omega(im)
 
1473
            pp3 = max(pp+p(3),rZero)
 
1474
            chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
 
1475
            if ( pp3.eq.rZero ) then
 
1476
               chi(2) = dcmplx(-nh )
 
1477
            else
 
1478
               chi(2) = dcmplx( nh*p(1) , p(2) )/dsqrt(rTwo*pp*pp3)
 
1479
            endif
 
1480
              fim(1) = sfomeg(1)*chi(im)
 
1481
              fim(2) = sfomeg(1)*chi(ip)
 
1482
              fim(3) = sfomeg(2)*chi(im)
 
1483
              fim(4) = sfomeg(2)*chi(ip)
 
1484
         endif
 
1485
       else
 
1486
         sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsr
 
1487
         chi(1) = dcmplx( sqp0p3 )
 
1488
         if ( sqp0p3.eq.rZero ) then
 
1489
            chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
 
1490
         else
 
1491
            chi(2) = dcmplx( nh*p(1), p(2) )/sqp0p3
 
1492
         endif
 
1493
         if ( nh.eq.1 ) then
 
1494
            fim(1) = dcmplx( rZero )
 
1495
            fim(2) = dcmplx( rZero )
 
1496
            fim(3) = chi(1)
 
1497
            fim(4) = chi(2)
 
1498
         else
 
1499
            fim(1) = chi(2)
 
1500
            fim(2) = chi(1)
 
1501
            fim(3) = dcmplx( rZero )
 
1502
            fim(4) = dcmplx( rZero )
 
1503
         endif
 
1504
       endif
 
1505
      end if
 
1506
 
 
1507
c spin-3/2 fermion wavefunction
 
1508
      if ( nhel.eq.3 ) then
 
1509
         do j = 1,4
 
1510
            do i = 1,4
 
1511
               rc(i,j) = ep(i)*fip(j)
 
1512
            end do
 
1513
         end do
 
1514
      else if ( nhel.eq.1 ) then
 
1515
         do j = 1,4
 
1516
            do i = 1,4
 
1517
              if     ( pt.eq.rZero .and. p(3).ge.0d0 ) then 
 
1518
               rc(i,j) =  sq2/sq3*e0(i)*fip(j) +rOne/sq3*ep(i)*fim(j)
 
1519
              elseif ( pt.eq.rZero .and. p(3).lt.0d0 ) then 
 
1520
               rc(i,j) =  sq2/sq3*e0(i)*fip(j) -rOne/sq3*ep(i)*fim(j)
 
1521
              else
 
1522
               rc(i,j) =  sq2/sq3*e0(i)*fip(j) 
 
1523
     &                  +rOne/sq3*ep(i)*fim(j) *dcmplx(P(1),nsr*P(2))/pt  
 
1524
              endif
 
1525
            end do
 
1526
         end do
 
1527
      else if ( nhel.eq.-1 ) then
 
1528
         do j = 1,4
 
1529
            do i = 1,4
 
1530
              if     ( pt.eq.rZero .and.p(3).ge.0d0 ) then 
 
1531
               rc(i,j) = rOne/sq3*em(i)*fip(j) +sq2/sq3*e0(i)*fim(j)
 
1532
              elseif ( pt.eq.rZero .and.p(3).lt.0d0 ) then 
 
1533
               rc(i,j) = rOne/sq3*em(i)*fip(j) -sq2/sq3*e0(i)*fim(j)
 
1534
              else
 
1535
               rc(i,j) = rOne/sq3*em(i)*fip(j) 
 
1536
     &                  + sq2/sq3*e0(i)*fim(j) *dcmplx(P(1),nsr*P(2))/pt  
 
1537
              endif
 
1538
            end do
 
1539
         end do
 
1540
      else
 
1541
         do j = 1,4
 
1542
            do i = 1,4
 
1543
              if     ( pt.eq.rZero .and. p(3).ge.0d0 ) then 
 
1544
               rc(i,j) =  em(i)*fim(j)
 
1545
              elseif ( pt.eq.rZero .and. p(3).lt.0d0 ) then 
 
1546
               rc(i,j) = -em(i)*fim(j)
 
1547
              else
 
1548
               rc(i,j) =  em(i)*fim(j) *dcmplx(P(1),nsr*P(2))/pt  
 
1549
              endif
 
1550
            end do
 
1551
         end do
 
1552
      end if
 
1553
 
 
1554
      ri(3) = rc(1,1)
 
1555
      ri(4) = rc(1,2)
 
1556
      ri(5) = rc(1,3)
 
1557
      ri(6) = rc(1,4)
 
1558
      ri(7) = rc(2,1)
 
1559
      ri(8) = rc(2,2)
 
1560
      ri(9) = rc(2,3)
 
1561
      ri(10) = rc(2,4)
 
1562
      ri(11) = rc(3,1)
 
1563
      ri(12) = rc(3,2)
 
1564
      ri(13) = rc(3,3)
 
1565
      ri(14) = rc(3,4)
 
1566
      ri(15) = rc(4,1)
 
1567
      ri(16) = rc(4,2)
 
1568
      ri(17) = rc(4,3)
 
1569
      ri(18) = rc(4,4)
 
1570
      ri(1) = rc(5,1)
 
1571
      ri(2) = rc(6,1)
 
1572
 
 
1573
      return
 
1574
      end
 
1575
      subroutine orxxxx(p,rmass,nhel,nsr , ro)
 
1576
c
 
1577
c This subroutine computes a Rarita-Schwinger wavefunction of spin-3/2
 
1578
c fermion with the flowing-IN fermion number.
 
1579
c
 
1580
c input:
 
1581
c       real    p(0:3)           : four-momentum of RS fermion
 
1582
c       real    rmass            : mass          of RS fermion
 
1583
c       integer nhel = -3,-1,1,3 : helicity      of RS fermion
 
1584
c                                  (1- and 1 is forbidden if rmass = 0)
 
1585
c       integer nsr  = -1 or 1   : +1 for particle, -1 for anti-particle
 
1586
c
 
1587
c output:
 
1588
c       complex ro(18)           : RS fermion wavefunction         |ro>v   
 
1589
c     
 
1590
c- by Y.Takaesu - 2011/01/11
 
1591
c
 
1592
      implicit none
 
1593
      double precision p(0:3),rmass
 
1594
      integer nhel,nsr
 
1595
      double complex ro(18),fipp(4),fimm(4)
 
1596
 
 
1597
      double complex rc(6,4),ep(4),em(4),e0(4),fop(4),fom(4),chi(2)
 
1598
      double precision pp,pt2,pt,pzpt,emp, sf(2),sfomeg(2),omega(2),pp3,
 
1599
     &                 sqp0p3,sqm(0:1)      
 
1600
      integer i,j,nsv,ip,im,nh
 
1601
 
 
1602
      double precision rZero, rHalf, rOne, rTwo, rThree, sqh,sq2,sq3
 
1603
      parameter( rZero = 0.0d0, rHalf = 0.5d0 )
 
1604
      parameter( rOne = 1.0d0, rTwo = 2.0d0, rThree = 3.0d0 )
 
1605
 
 
1606
c#ifdef HELAS_CHECK
 
1607
c      double precision p2
 
1608
c      double precision epsi
 
1609
c      parameter( epsi = 2.0d-5 )
 
1610
c      integer stdo
 
1611
c      parameter( stdo = 6 )
 
1612
c#endif
 
1613
c
 
1614
c#ifdef HELAS_CHECK
 
1615
c      pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
 
1616
c      if ( abs(p(0))+pp.eq.rZero ) then
 
1617
c         write(stdo,*)
 
1618
c     &        ' helas-error : p(0:3) in orxxxx is zero momentum'
 
1619
c      endif
 
1620
c      if ( p(0).le.rZero ) then
 
1621
c         write(stdo,*)
 
1622
c     &        ' helas-error : p(0:3) in orxxxx has non-positive energy'
 
1623
c         write(stdo,*)
 
1624
c     &        '             : p(0) = ',p(0)
 
1625
c      endif
 
1626
c      p2 = (p(0)-pp)*(p(0)+pp)
 
1627
c      if ( abs(p2-rmass**2).gt.p(0)**2*epsi ) then
 
1628
c         write(stdo,*)
 
1629
c     &        ' helas-error : p(0:3) in orxxxx has inappropriate mass'
 
1630
c         write(stdo,*)
 
1631
c     &        '             : p**2 = ',p2,' : rmass**2 = ',rmass**2
 
1632
c      endif
 
1633
c      if (abs(nhel).gt.3 .or. abs(nhel).eq.2 .or. abs(nhel).eq.0 ) then
 
1634
c         write(stdo,*) ' helas-error : nhel in orxxxx is not -3,-1,1,3'
 
1635
c         write(stdo,*) '             : nhel = ',nhel
 
1636
c      endif
 
1637
c      if (abs(nsr).ne.1) then
 
1638
c         write(stdo,*) ' helas-error : nsr in orxxxx is not -1,1'
 
1639
c         write(stdo,*) '             : nsr = ',nsr
 
1640
c      endif
 
1641
c#endif
 
1642
 
 
1643
      sqh = sqrt(rHalf)
 
1644
      sq2 = sqrt(rTwo)
 
1645
      sq3 = sqrt(rThree)
 
1646
 
 
1647
      pt2 = p(1)**2 + p(2)**2
 
1648
      pp = min(p(0),sqrt(pt2+p(3)**2))
 
1649
      pt = min(pp,sqrt(pt2))
 
1650
 
 
1651
      rc(5,1) = dcmplx(p(0),p(3))*nsr
 
1652
      rc(6,1) = dcmplx(p(1),p(2))*nsr
 
1653
 
 
1654
      nsv = nsr ! nsv=+1 for final, -1 for initial
 
1655
 
 
1656
      if ( nhel.ge.1 ) then 
 
1657
c construct eps+
 
1658
         if ( pp.eq.rZero ) then
 
1659
            ep(1) = dcmplx( rZero )
 
1660
            ep(2) = dcmplx( -sqh )
 
1661
            ep(3) = dcmplx( rZero , nsv*sqh )
 
1662
            ep(4) = dcmplx( rZero )
 
1663
         else
 
1664
            ep(1) = dcmplx( rZero )
 
1665
            ep(4) = dcmplx( pt/pp*sqh )
 
1666
            if ( pt.ne.rZero ) then
 
1667
               pzpt = p(3)/(pp*pt)*sqh
 
1668
               ep(2) = dcmplx( -p(1)*pzpt , -nsv*p(2)/pt*sqh )
 
1669
               ep(3) = dcmplx( -p(2)*pzpt ,  nsv*p(1)/pt*sqh )
 
1670
            else
 
1671
               ep(2) = dcmplx( -sqh )
 
1672
               ep(3) = dcmplx( rZero , nsv*sign(sqh,p(3)) )
 
1673
            endif
 
1674
         endif
 
1675
      end if
 
1676
 
 
1677
      if ( nhel.le.-1 ) then 
 
1678
c construct eps-
 
1679
         if ( pp.eq.rZero ) then
 
1680
            em(1) = dcmplx( rZero )
 
1681
            em(2) = dcmplx( sqh )
 
1682
            em(3) = dcmplx( rZero , nsv*sqh )
 
1683
            em(4) = dcmplx( rZero )
 
1684
         else
 
1685
            em(1) = dcmplx( rZero )
 
1686
            em(4) = dcmplx( -pt/pp*sqh )
 
1687
            if ( pt.ne.rZero ) then
 
1688
               pzpt = -p(3)/(pp*pt)*sqh
 
1689
               em(2) = dcmplx( -p(1)*pzpt , -nsv*p(2)/pt*sqh )
 
1690
               em(3) = dcmplx( -p(2)*pzpt ,  nsv*p(1)/pt*sqh )
 
1691
            else
 
1692
               em(2) = dcmplx( sqh )
 
1693
               em(3) = dcmplx( rZero , nsv*sign(sqh,p(3)) )
 
1694
            endif
 
1695
         endif
 
1696
      end if
 
1697
 
 
1698
      if ( abs(nhel).le.1 ) then  
 
1699
c construct eps0
 
1700
         if ( pp.eq.rZero ) then
 
1701
            e0(1) = dcmplx( rZero )
 
1702
            e0(2) = dcmplx( rZero )
 
1703
            e0(3) = dcmplx( rZero )
 
1704
            e0(4) = dcmplx( rOne )
 
1705
         else
 
1706
            emp = p(0)/(rmass*pp)
 
1707
            e0(1) = dcmplx( pp/rmass )
 
1708
            e0(4) = dcmplx( p(3)*emp )
 
1709
            if ( pt.ne.rZero ) then
 
1710
               e0(2) = dcmplx( p(1)*emp )
 
1711
               e0(3) = dcmplx( p(2)*emp )
 
1712
            else
 
1713
               e0(2) = dcmplx( rZero )
 
1714
               e0(3) = dcmplx( rZero )
 
1715
            endif
 
1716
         end if
 
1717
      end if
 
1718
 
 
1719
      if ( nhel.ge.-1 ) then
 
1720
c constract spinor+ 
 
1721
       nh = nsr
 
1722
 
 
1723
       if ( rmass.ne.rZero ) then
 
1724
 
 
1725
         pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
1726
 
 
1727
         if ( pp.eq.rZero ) then
 
1728
            
 
1729
            sqm(0) = dsqrt(abs(rmass)) ! possibility of negative fermion masses
 
1730
            sqm(1) = sign(sqm(0),rmass) ! possibility of negative fermion masses
 
1731
            ip = -((1+nh)/2)
 
1732
            im =  (1-nh)/2
 
1733
            
 
1734
            fop(1) = im     * sqm(im)
 
1735
            fop(2) = ip*nsr * sqm(im)
 
1736
            fop(3) = im*nsr * sqm(-ip)
 
1737
            fop(4) = ip     * sqm(-ip)
 
1738
            
 
1739
         else
 
1740
            
 
1741
            pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
1742
            sf(1) = dble(1+nsr+(1-nsr)*nh)*rHalf
 
1743
            sf(2) = dble(1+nsr-(1-nsr)*nh)*rHalf
 
1744
            omega(1) = dsqrt(p(0)+pp)
 
1745
            omega(2) = rmass/omega(1)
 
1746
            ip = (3+nh)/2
 
1747
            im = (3-nh)/2
 
1748
            sfomeg(1) = sf(1)*omega(ip)
 
1749
            sfomeg(2) = sf(2)*omega(im)
 
1750
            pp3 = max(pp+p(3),rZero)
 
1751
            chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
 
1752
            if ( pp3.eq.rZero ) then
 
1753
               chi(2) = dcmplx(-nh )
 
1754
            else
 
1755
               chi(2) = dcmplx( nh*p(1) , -p(2) )/dsqrt(rTwo*pp*pp3)
 
1756
            endif
 
1757
            
 
1758
            fop(1) = sfomeg(2)*chi(im)
 
1759
            fop(2) = sfomeg(2)*chi(ip)
 
1760
            fop(3) = sfomeg(1)*chi(im)
 
1761
            fop(4) = sfomeg(1)*chi(ip)
 
1762
 
 
1763
         endif
 
1764
         
 
1765
      else
 
1766
         
 
1767
         if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
 
1768
            sqp0p3 = 0d0
 
1769
         else
 
1770
            sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsr
 
1771
         end if
 
1772
         chi(1) = dcmplx( sqp0p3 )
 
1773
         if ( sqp0p3.eq.rZero ) then
 
1774
            chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
 
1775
         else
 
1776
            chi(2) = dcmplx( nh*p(1), -p(2) )/sqp0p3
 
1777
         endif
 
1778
         if ( nh.eq.1 ) then
 
1779
            fop(1) = chi(1)
 
1780
            fop(2) = chi(2)
 
1781
            fop(3) = dcmplx( rZero )
 
1782
            fop(4) = dcmplx( rZero )
 
1783
         else
 
1784
            fop(1) = dcmplx( rZero )
 
1785
            fop(2) = dcmplx( rZero )
 
1786
            fop(3) = chi(2)
 
1787
            fop(4) = chi(1)
 
1788
         endif
 
1789
       endif
 
1790
      endif
 
1791
 
 
1792
      if ( nhel.le.1 ) then
 
1793
c constract spinor+ 
 
1794
       nh = -nsr
 
1795
 
 
1796
      if ( rmass.ne.rZero ) then
 
1797
 
 
1798
         pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
1799
 
 
1800
         if ( pp.eq.rZero ) then
 
1801
            
 
1802
            sqm(0) = dsqrt(abs(rmass)) ! possibility of negative fermion masses
 
1803
            sqm(1) = sign(sqm(0),rmass) ! possibility of negative fermion masses
 
1804
            ip = -((1+nh)/2)
 
1805
            im =  (1-nh)/2
 
1806
            
 
1807
            fom(1) = im     * sqm(im)
 
1808
            fom(2) = ip*nsr * sqm(im)
 
1809
            fom(3) = im*nsr * sqm(-ip)
 
1810
            fom(4) = ip     * sqm(-ip)
 
1811
            
 
1812
         else
 
1813
            
 
1814
            pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
1815
            sf(1) = dble(1+nsr+(1-nsr)*nh)*rHalf
 
1816
            sf(2) = dble(1+nsr-(1-nsr)*nh)*rHalf
 
1817
            omega(1) = dsqrt(p(0)+pp)
 
1818
            omega(2) = rmass/omega(1)
 
1819
            ip = (3+nh)/2
 
1820
            im = (3-nh)/2
 
1821
            sfomeg(1) = sf(1)*omega(ip)
 
1822
            sfomeg(2) = sf(2)*omega(im)
 
1823
            pp3 = max(pp+p(3),rZero)
 
1824
            chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
 
1825
            if ( pp3.eq.rZero ) then
 
1826
               chi(2) = dcmplx(-nh )
 
1827
            else
 
1828
               chi(2) = dcmplx( nh*p(1) , -p(2) )/dsqrt(rTwo*pp*pp3)
 
1829
            endif
 
1830
            
 
1831
            fom(1) = sfomeg(2)*chi(im)
 
1832
            fom(2) = sfomeg(2)*chi(ip)
 
1833
            fom(3) = sfomeg(1)*chi(im)
 
1834
            fom(4) = sfomeg(1)*chi(ip)
 
1835
 
 
1836
         endif
 
1837
         
 
1838
      else
 
1839
         
 
1840
         if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
 
1841
            sqp0p3 = 0d0
 
1842
         else
 
1843
            sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsr
 
1844
         end if
 
1845
         chi(1) = dcmplx( sqp0p3 )
 
1846
         if ( sqp0p3.eq.rZero ) then
 
1847
            chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
 
1848
         else
 
1849
            chi(2) = dcmplx( nh*p(1), -p(2) )/sqp0p3
 
1850
         endif
 
1851
         if ( nh.eq.1 ) then
 
1852
            fom(1) = chi(1)
 
1853
            fom(2) = chi(2)
 
1854
            fom(3) = dcmplx( rZero )
 
1855
            fom(4) = dcmplx( rZero )
 
1856
         else
 
1857
            fom(1) = dcmplx( rZero )
 
1858
            fom(2) = dcmplx( rZero )
 
1859
            fom(3) = chi(2)
 
1860
            fom(4) = chi(1)
 
1861
         endif
 
1862
       endif 
 
1863
      endif
 
1864
      
 
1865
c spin-3/2 fermion wavefunction
 
1866
      if ( nhel.eq.3 ) then
 
1867
         do j = 1,4
 
1868
            do i = 1,4
 
1869
               rc(i,j) = ep(i)*fop(j)
 
1870
            end do
 
1871
         end do
 
1872
      else if ( nhel.eq.1 ) then
 
1873
         do j = 1,4
 
1874
            do i = 1,4
 
1875
              if     ( pt.eq.rZero .and. p(3).ge.0d0 ) then 
 
1876
               rc(i,j) =  sq2/sq3*e0(i)*fop(j)
 
1877
     &                    +rOne/sq3*ep(i)*fom(j)
 
1878
              elseif ( pt.eq.rZero .and. p(3).lt.0d0 ) then 
 
1879
               rc(i,j) =  sq2/sq3*e0(i)*fop(j)
 
1880
     &                    -rOne/sq3*ep(i)*fom(j)
 
1881
              else
 
1882
               rc(i,j) =  sq2/sq3*e0(i)*fop(j) 
 
1883
     &                  +rOne/sq3*ep(i)*fom(j)
 
1884
     &                   *dcmplx(P(1),-nsr*P(2))/pt  
 
1885
              endif
 
1886
            end do
 
1887
         end do
 
1888
      else if ( nhel.eq.-1 ) then
 
1889
         do j = 1,4
 
1890
            do i = 1,4
 
1891
              if     ( pt.eq.rZero .and.p(3).ge.0d0 ) then 
 
1892
               rc(i,j) = rOne/sq3*em(i)*fop(j)
 
1893
     &                   +sq2/sq3*e0(i)*fom(j)
 
1894
              elseif ( pt.eq.rZero .and.p(3).lt.0d0 ) then 
 
1895
               rc(i,j) = rOne/sq3*em(i)*fop(j)
 
1896
     &                   -sq2/sq3*e0(i)*fom(j)
 
1897
              else
 
1898
               rc(i,j) = rOne/sq3*em(i)*fop(j) 
 
1899
     &                  + sq2/sq3*e0(i)*fom(j)
 
1900
     &                   *dcmplx(P(1),-nsr*P(2))/pt  
 
1901
              endif
 
1902
            end do
 
1903
         end do
 
1904
      else
 
1905
         do j = 1,4
 
1906
            do i = 1,4
 
1907
              if     ( pt.eq.rZero .and. p(3).ge.0d0 ) then 
 
1908
               rc(i,j) =  em(i)*fom(j)
 
1909
              elseif ( pt.eq.rZero .and. p(3).lt.0d0 ) then 
 
1910
               rc(i,j) = -em(i)*fom(j)
 
1911
              else
 
1912
               rc(i,j) =  em(i)*fom(j)*dcmplx(P(1),-nsr*P(2))/pt  
 
1913
              endif
 
1914
            end do
 
1915
         end do
 
1916
      end if
 
1917
 
 
1918
      ro(3) = rc(1,1)
 
1919
      ro(4) = rc(1,2)
 
1920
      ro(5) = rc(1,3)
 
1921
      ro(6) = rc(1,4)
 
1922
      ro(7) = rc(2,1)
 
1923
      ro(8) = rc(2,2)
 
1924
      ro(9) = rc(2,3)
 
1925
      ro(10) = rc(2,4)
 
1926
      ro(11) = rc(3,1)
 
1927
      ro(12) = rc(3,2)
 
1928
      ro(13) = rc(3,3)
 
1929
      ro(14) = rc(3,4)
 
1930
      ro(15) = rc(4,1)
 
1931
      ro(16) = rc(4,2)
 
1932
      ro(17) = rc(4,3)
 
1933
      ro(18) = rc(4,4)
 
1934
      ro(1) = rc(5,1)
 
1935
      ro(2) = rc(6,1)
 
1936
 
 
1937
      return
 
1938
      end