Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pymael.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pymael.f
1 
2 C*********************************************************************
3 
4 C...PYMAEL
5 C...Auxiliary to PYSHOW and PYPTFS.
6 C...Matrix elements for gluon (or photon) emission from
7 C...a two-body state; to be used by the parton shower routine.
8 C...Here X_i = 2 E_i/E_cm, R_i = m_i/E_cm and
9 C...1/sigma_0 d(sigma)/d(x_1)d(x_2) =
10 C... = (alpha-strong/2 pi) * CF * PYMAEL,
11 C...i.e. normalization is such that one recovers the familiar
12 C...(X1**2+X2**2)/((1-X1)*(1-X2)) for the massless case.
13 C...Coupling structure:
14 C...NI = 6- 9 : eikonal soft-gluon expression (spin-independent)
15 C... = 11-14 : V -> q qbar (V = vector/axial vector colour singlet)
16 C... = 16-19 : q -> q V
17 C... = 21-24 : S -> q qbar (S = scalar/pseudoscalar colour singlet)
18 C... = 26-29 : q -> q S
19 C... = 31-34 : V -> ~q ~qbar (~q = squark)
20 C... = 36-39 : ~q -> ~q V
21 C... = 41-44 : S -> ~q ~qbar
22 C... = 46-49 : ~q -> ~q S
23 C... = 51-54 : chi -> q ~qbar (chi = neutralino/chargino)
24 C... = 56-59 : ~q -> q chi
25 C... = 61-64 : q -> ~q chi
26 C... = 66-69 : ~g -> q ~qbar
27 C... = 71-74 : ~q -> q ~g
28 C... = 76-79 : q -> ~q ~g
29 C... = 81-84 : (9/4)*(eikonal) for gg -> ~g ~g
30 C...Note that the order of the decay products is important.
31 C...In each set of four, the variants are ordered as:
32 C...ICOMBI = 1 : pure non-gamma5, i.e. vector/scalar/...
33 C... = 2 : pure gamma5, i.e. axial vector/pseudoscalar/....
34 C... = 3 : mixture alpha*(ICOMBI=1) + (1-alpha)*(ICOMBI=2)
35 C... = 4 : mixture (ICOMBI=1) +- (ICOMBI=2)
36 
37  FUNCTION pymael(NI,X1,X2,R1,R2,ALPHA)
38 
39 C...Double precision and integer declarations.
40  IMPLICIT DOUBLE PRECISION(a-h, o-z)
41  IMPLICIT INTEGER(i-n)
42 
43 C...Check input values. Return zero outside allowed phase space.
44  pymael=0d0
45  IF(x1.LE.2d0*r1.OR.x1.GE.1d0+r1**2-r2**2) RETURN
46  IF(x2.LE.2d0*r2.OR.x2.GE.1d0+r2**2-r1**2) RETURN
47  IF(x1+x2.LE.1d0+(r1+r2)**2) RETURN
48  IF((2d0-2d0*x1-2d0*x2+x1*x2+2d0*r1**2+2d0*r2**2)**2.GE.
49  &(x1**2-4d0*r1**2)*(x2**2-4d0*r2**2)) RETURN
50  alpcor=max(0d0,min(1d0,alpha))
51 
52 C...Initial values and flags.
53  iclass=ni/5
54  icombi=ni-5*iclass
55  isset1=0
56  isset2=0
57  isset4=0
58 
59 C... Phase space.
60  ps=sqrt((1d0-(r1+r2)**2)*(1d0-(r1-r2)**2))
61 
62 C...Eikonal expression; also acts as default.
63  IF(iclass.LE.1.OR.iclass.GE.17.OR.icombi.EQ.0) THEN
64  rlo=ps
65  IF(icombi.EQ.0.OR.icombi.EQ.1) THEN
66  anum=0d0
67  ELSEIF(icombi.EQ.2) THEN
68  anum=(2d0-x1-x2)**2
69  ELSEIF(icombi.EQ.3) THEN
70  anum=alpcor*(2d0-x1-x2)**2
71  ELSE
72  anum=0.5d0*(2d0-x1-x2)**2
73  ENDIF
74  rfo=ps*2d0*((x1+x2-1d0+anum-r1**2-r2**2)/
75  & ((1d0+r1**2-r2**2-x1)*(1d0+r2**2-r1**2-x2))-
76  & r1**2/(1d0+r2**2-r1**2-x2)**2-
77  & r2**2/(1d0+r1**2-r2**2-x1)**2)
78  icombi=0
79 
80 C...V -> q qbar (V = gamma*/Z0/W+-/...).
81  ELSEIF(iclass.EQ.2) THEN
82  IF(icombi.EQ.1.OR.icombi.EQ.3) THEN
83  rlo1=ps*(2-r1**2-r1**4+6*r1*r2-r2**2+2*r1**2*r2**2-r2**4)/2.d0
84  rfo1=-1.d0*(3+6*r1**2+r1**4-6*r1*r2+6*r1**3*r2-2*r2**2
85  & -6*r1**2*r2**2+6*r1*r2**3+r2**4-3*x1+6*r1*r2*x1
86  & +2*r2**2*x1+x1**2-2*r1**2*x1**2+3*r1**2*(2-x1-x2)
87  & +6*r1*r2*(2-x1-x2)-r2**2*(2-x1-x2)-2*x1*(2-x1-x2)
88  & -5*r1**2*x1*(2-x1-x2)+r2**2*x1*(2-x1-x2)+x1**2*(2-x1-x2)
89  & -3*(2-x1-x2)**2-3*r1**2*(2-x1-x2)**2+r2**2*(2-x1-x2)**2
90  & +2*x1*(2-x1-x2)**2+(2-x1-x2)**3-x2)/
91  & (-1+r1**2-r2**2+x2)**2
92  rfo1=rfo1-2*(-3+r1**2-6*r1*r2+6*r1**3*r2+3*r2**2-4*r1**2*r2**2
93  & +6*r1*r2**3+2*x1+3*r1**2*x1+r2**2*x1-x1**2-r1**2*x1**2
94  & -r2**2*x1**2+4*(2-x1-x2)+2*r1**2*(2-x1-x2)+3*r1*r2*(2-x1
95  & -x2)-r2**2*(2-x1-x2)-3*x1*(2-x1-x2)-2*r1**2*x1*(2-x1-x2)
96  & +x1**2*(2-x1-x2)-(2-x1-x2)**2-r1**2*(2-x1-x2)**2+r1*r2*(2
97  & -x1-x2)**2+x1*(2-x1-x2)**2)/
98  & (-1-r1**2+r2**2+x1)/(-1+r1**2-r2**2+x2)
99  rfo1=rfo1-1.d0*(-1+2*r1**2+r1**4+6*r1*r2+6*r1**3*r2-2*r2**2
100  & -6*r1**2*r2**2+6*r1*r2**3+r2**4-x1-2*r1**2*x1-6*r1*r2*x1
101  & +8*r2**2*x1+x1**2-2*r2**2*x1**2-r1**2*(2-x1-x2)+r2**2*(2
102  & -x1-x2)-r1**2*x1*(2-x1-x2)+r2**2*x1*(2-x1-x2)+x1**2*
103  & (2-x1-x2)+x2)/(-1-r1**2+r2**2+x1)**2
104  rfo1=rfo1/2.d0
105  isset1=1
106  ENDIF
107  IF(icombi.EQ.2.OR.icombi.EQ.3) THEN
108  rlo2=ps*(2-r1**2-r1**4-6*r1*r2-r2**2+2*r1**2*r2**2-r2**4)/2.d0
109  rfo2=-1*(3+6*r1**2+r1**4+6*r1*r2-6*r1**3*r2-2*r2**2
110  & -6*r1**2*r2**2-6*r1*r2**3+r2**4-3*x1-6*r1*r2*x1+2*r2**2*x1
111  & +x1**2-2*r1**2*x1**2+3*r1**2*(2-x1-x2)-6*r1*r2*(2-x1-x2)
112  & -r2**2*(2-x1-x2)-2*x1*(2-x1-x2)-5*r1**2*x1*(2-x1-x2)
113  & +r2**2*x1*(2-x1-x2)+x1**2*(2-x1-x2)-3*(2-x1-x2)**2
114  & -3*r1**2*(2-x1-x2)**2+r2**2*(2-x1-x2)**2+2*x1*(2-x1-x2)**2
115  & +(2-x1-x2)**3-x2)/(-1+r1**2-r2**2+x2)**2
116  rfo2=rfo2-2*(-3+r1**2+6*r1*r2-6*r1**3*r2+3*r2**2-4*r1**2*r2**2
117  & -6*r1*r2**3+2*x1+3*r1**2*x1+r2**2*x1-x1**2-r1**2*x1**2
118  & -r2**2*x1**2+4*(2-x1-x2)+2*r1**2*(2-x1-x2)-3*r1*r2*(2-x1
119  & -x2)-r2**2*(2-x1-x2)-3*x1*(2-x1-x2)-2*r1**2*x1*(2-x1-x2)
120  & +x1**2*(2-x1-x2)-(2-x1-x2)**2-r1**2*(2-x1-x2)**2-r1*r2*(2
121  & -x1-x2)**2+x1*(2-x1-x2)**2)/
122  & (-1-r1**2+r2**2+x1)/(-1+r1**2-r2**2+x2)
123  rfo2=rfo2-1*(-1+2*r1**2+r1**4-6*r1*r2-6*r1**3*r2-2*r2**2
124  & -6*r1**2*r2**2-6*r1*r2**3+r2**4-x1-2*r1**2*x1+6*r1*r2*x1
125  & +8*r2**2*x1+x1**2-2*r2**2*x1**2-r1**2*(2-x1-x2)+r2**2*(2-x1
126  & -x2)-r1**2*x1*(2-x1-x2)+r2**2*x1*(2-x1-x2)+x1**2*(2-x1-x2)
127  & +x2)/(-1-r1**2+r2**2+x1)**2
128  rfo2=rfo2/2.d0
129  isset2=1
130  ENDIF
131  IF(icombi.EQ.4) THEN
132  rlo4=ps*(2d0-r1**2-r1**4-r2**2+2d0*r1**2*r2**2-r2**4)/2d0
133  rfo4=(1-r1**4+6*r1**2*r2**2-r2**4+x1+3*r1**2*x1-9*r2**2*x1
134  & -3*x1**2-r1**2*x1**2+3*r2**2*x1**2+x1**3-x2-r1**2*x2
135  & +r2**2*x2-r1**2*x1*x2+r2**2*x1*x2+x1**2*x2)/
136  & (-1-r1**2+r2**2+x1)**2
137  rfo4=rfo4
138  & -2*(1+r1**2+r2**2-4*r1**2*r2**2+r1**2*x1+2*r2**2*x1-x1**2
139  & -r2**2*x1**2+2*r1**2*x2+r2**2*x2-3*x1*x2+x1**2*x2-x2**2
140  & -r1**2*x2**2+x1*x2**2)/
141  & (-1-r1**2+r2**2+x1)/(-1+r1**2-r2**2+x2)
142  rfo4=rfo4+(1-r1**4+6*r1**2*r2**2-r2**4-x1+r1**2*x1-r2**2*x1+x2
143  & -9*r1**2*x2+3*r2**2*x2+r1**2*x1*x2-r2**2*x1*x2-3*x2**2
144  & +3*r1**2*x2**2-r2**2*x2**2+x1*x2**2+x2**3)/
145  & (-1+r1**2-r2**2+x2)**2
146  rfo4=rfo4/2.d0
147  isset4=1
148  ENDIF
149 
150 C...q -> q V.
151  ELSEIF(iclass.EQ.3) THEN
152  IF(icombi.EQ.1.OR.icombi.EQ.3) THEN
153  rlo1=ps*(1d0-2d0*r1**2+r1**4+r2**2-6d0*r1*r2**2
154  & +r1**2*r2**2-2d0*r2**4)
155  rfo1=2*(-1+r1-2*r1**2+2*r1**3-r1**4+r1**5-r2**2+r1*r2**2
156  & -5*r1**2*r2**2+r1**3*r2**2-2*r1*r2**4+2*x1-2*r1*x1
157  & +2*r1**2*x1-2*r1**3*x1+2*r2**2*x1+5*r1*r2**2*x1
158  & +r1**2*r2**2*x1+2*r2**4*x1-x1**2+r1*x1**2-r2**2*x1**2+3*x2
159  & +4*r1**2*x2+r1**4*x2+2*r2**2*x2+2*r1**2*r2**2*x2-4*x1*x2
160  & -2*r1**2*x1*x2-r2**2*x1*x2+x1**2*x2-2*x2**2
161  & -2*r1**2*x2**2+x1*x2**2)/(1-r1**2+r2**2-x2)/(-2+x1+x2)
162  rfo1=rfo1+(2*r2**2+6*r1*r2**2-6*r1**2*r2**2+6*r1**3*r2**2
163  & +2*r2**4+6*r1*r2**4-r2**2*x1+r1**2*r2**2*x1-r2**4*x1+x2
164  & -r1**4*x2-3*r2**2*x2-6*r1*r2**2*x2+9*r1**2*r2**2*x2
165  & -2*r2**4*x2-x1*x2+r1**2*x1*x2-x2**2-3*r1**2*x2**2
166  & +2*r2**2*x2**2+x1*x2**2)/(-1+r1**2-r2**2+x2)**2
167  rfo1=rfo1+(-4-8*r1**2-4*r1**4+4*r2**2-4*r1**2*r2**2+8*r2**4
168  & +9*x1+10*r1**2*x1+r1**4*x1-3*r2**2*x1+6*r1*r2**2*x1
169  & +r1**2*r2**2*x1-2*r2**4*x1-6*x1**2-2*r1**2*x1**2+x1**3
170  & +7*x2+8*r1**2*x2+r1**4*x2-7*r2**2*x2+6*r1*r2**2*x2
171  & +r1**2*r2**2*x2-2*r2**4*x2-9*x1*x2-3*r1**2*x1*x2
172  & +2*r2**2*x1*x2+2*x1**2*x2-3*x2**2-r1**2*x2**2
173  & +2*r2**2*x2**2+x1*x2**2)/(-2+x1+x2)**2
174  isset1=1
175  ENDIF
176  IF(icombi.EQ.2.OR.icombi.EQ.3) THEN
177  rlo2=ps*(1d0-2d0*r1**2+r1**4+r2**2+6d0*r1*r2**2
178  & +r1**2*r2**2-2d0*r2**4)
179  rfo2=2*(1+r1+2*r1**2+2*r1**3+r1**4+r1**5+r2**2+r1*r2**2
180  & +5*r1**2*r2**2+r1**3*r2**2-2*r1*r2**4-2*x1-2*r1*x1
181  & -2*r1**2*x1-2*r1**3*x1-2*r2**2*x1+5*r1*r2**2*x1
182  & -r1**2*r2**2*x1-2*r2**4*x1+x1**2+r1*x1**2+r2**2*x1**2-3*x2
183  & -4*r1**2*x2-r1**4*x2-2*r2**2*x2-2*r1**2*r2**2*x2+4*x1*x2
184  & +2*r1**2*x1*x2+r2**2*x1*x2-x1**2*x2+2*x2**2+2*r1**2*x2**2
185  & -x1*x2**2)/(-1+r1**2-r2**2+x2)/(-2+x1+x2)
186  rfo2=rfo2+(2*r2**2-6*r1*r2**2-6*r1**2*r2**2-6*r1**3*r2**2
187  & +2*r2**4-6*r1*r2**4-r2**2*x1+r1**2*r2**2*x1-r2**4*x1+x2
188  & -r1**4*x2-3*r2**2*x2+6*r1*r2**2*x2+9*r1**2*r2**2*x2
189  & -2*r2**4*x2-x1*x2+r1**2*x1*x2-x2**2-3*r1**2*x2**2
190  & +2*r2**2*x2**2+x1*x2**2)/(-1+r1**2-r2**2+x2)**2
191  rfo2=rfo2+(-4-8*r1**2-4*r1**4+4*r2**2-4*r1**2*r2**2+8*r2**4+9*x1
192  & +10*r1**2*x1+r1**4*x1-3*r2**2*x1-6*r1*r2**2*x1
193  & +r1**2*r2**2*x1-2*r2**4*x1-6*x1**2-2*r1**2*x1**2+x1**3
194  & +7*x2+8*r1**2*x2+r1**4*x2-7*r2**2*x2-6*r1*r2**2*x2
195  & +r1**2*r2**2*x2-2*r2**4*x2-9*x1*x2-3*r1**2*x1*x2
196  & +2*r2**2*x1*x2+2*x1**2*x2-3*x2**2-r1**2*x2**2+2*r2**2*x2**2
197  & +x1*x2**2)/(-2+x1+x2)**2
198  isset2=1
199  ENDIF
200  IF(icombi.EQ.4) THEN
201  rlo4=ps*(1.d0-2.d0*r1**2+r1**4+r2**2+r1**2*r2**2-2.d0*r2**4)
202  rfo4=2*(1+2*r1**2+r1**4+r2**2+5*r1**2*r2**2-2*x1-2*r1**2*x1
203  & -2*r2**2*x1-r1**2*r2**2*x1-2*r2**4*x1+x1**2+r2**2*x1**2
204  & -3*x2-4*r1**2*x2-r1**4*x2-2*r2**2*x2-2*r1**2*r2**2*x2
205  & +4*x1*x2+2*r1**2*x1*x2+r2**2*x1*x2-x1**2*x2+2*x2**2
206  & +2*r1**2*x2**2-x1*x2**2)/(-1+r1**2-r2**2+x2)/(-2+x1+x2)
207  rfo4=rfo4+(2*r2**2-6*r1**2*r2**2+2*r2**4-r2**2*x1+r1**2*r2**2*x1
208  & -r2**4*x1+x2-r1**4*x2-3*r2**2*x2+9*r1**2*r2**2*x2
209  & -2*r2**4*x2-x1*x2+r1**2*x1*x2-x2**2-3*r1**2*x2**2
210  & +2*r2**2*x2**2+x1*x2**2)/(-1+r1**2-r2**2+x2)**2
211  rfo4=rfo4+(-4-8*r1**2-4*r1**4+4*r2**2-4*r1**2*r2**2+8*r2**4+9*x1
212  & +10*r1**2*x1+r1**4*x1-3*r2**2*x1+r1**2*r2**2*x1-2*r2**4*x1
213  & -6*x1**2-2*r1**2*x1**2+x1**3+7*x2+8*r1**2*x2+r1**4*x2
214  & -7*r2**2*x2+r1**2*r2**2*x2-2*r2**4*x2-9*x1*x2-3*r1**2*x1*x2
215  & +2*r2**2*x1*x2+2*x1**2*x2-3*x2**2-r1**2*x2**2+2*r2**2*x2**2
216  & +x1*x2**2)/(2-x1-x2)**2
217  isset4=1
218  ENDIF
219 
220 C...S -> q qbar (S = h0/H0/A0/H+-/...).
221  ELSEIF(iclass.EQ.4) THEN
222  IF(icombi.EQ.1.OR.icombi.EQ.3) THEN
223  rlo1=ps*(1d0-r1**2-r2**2-2d0*r1*r2)
224  rfo1=-(-1+r1**4-2*r1*r2-2*r1**3*r2-6*r1**2*r2**2-2*r1*r2**3
225  & +r2**4+x1-r1**2*x1+2*r1*r2*x1+3*r2**2*x1+x2+r1**2*x2
226  & -r2**2*x2-x1*x2)/(-1-r1**2+r2**2+x1)**2
227  & -2*(r1**2+r1**4-2*r1**3*r2+r2**2-6*r1**2*r2**2-2*r1*r2**3
228  & +r2**4-r1**2*x1+r1*r2*x1+2*r2**2*x1+2*r1**2*x2+r1*r2*x2
229  & -r2**2*x2-x1*x2)/(-1-r1**2+r2**2+x1)/(-1+r1**2-r2**2+x2)
230  & -(-1+r1**4-2*r1*r2-2*r1**3*r2-6*r1**2*r2**2-2*r1*r2**3
231  & +r2**4+x1-r1**2*x1+r2**2*x1+x2+3*r1**2*x2+2*r1*r2*x2
232  & -r2**2*x2-x1*x2)/(-1+r1**2-r2**2+x2)**2
233  isset1=1
234  ENDIF
235  IF(icombi.EQ.2.OR.icombi.EQ.3) THEN
236  rlo2=ps*(1d0-r1**2-r2**2+2d0*r1*r2)
237  rfo2=-(-1+r1**4+2*r1*r2+2*r1**3*r2-6*r1**2*r2**2+2*r1*r2**3
238  & +r2**4+x1-r1**2*x1-2*r1*r2*x1+3*r2**2*x1+x2+r1**2*x2
239  & -r2**2*x2-x1*x2)/(-1-r1**2+r2**2+x1)**2
240  & -(-1+r1**4+2*r1*r2+2*r1**3*r2-6*r1**2*r2**2+2*r1*r2**3
241  & +r2**4+x1-r1**2*x1+r2**2*x1+x2+3*r1**2*x2-2*r1*r2*x2
242  & -r2**2*x2-x1*x2)/(-1+r1**2-r2**2+x2)**2
243  & +2*(-r1**2-r1**4-2*r1**3*r2-r2**2+6*r1**2*r2**2
244  & -2*r1*r2**3-r2**4+r1**2*x1+r1*r2*x1-2*r2**2*x1
245  & -2*r1**2*x2+r1*r2*x2+r2**2*x2+x1*x2)/
246  & (-1-r1**2+r2**2+x1)/(-1+r1**2-r2**2+x2)
247  isset2=1
248  ENDIF
249  IF(icombi.EQ.4) THEN
250  rlo4=ps*(1d0-r1**2-r2**2)
251  rfo4=-(-1+r1**4-6*r1**2*r2**2+r2**4+x1-r1**2*x1+3*r2**2*x1+x2
252  & +r1**2*x2-r2**2*x2-x1*x2)/(-1-r1**2+r2**2+x1)**2
253  & -2*(r1**2+r1**4+r2**2-6*r1**2*r2**2+r2**4-r1**2*x1
254  & +2*r2**2*x1+2*r1**2*x2-r2**2*x2-x1*x2)/
255  & (-1-r1**2+r2**2+x1)/(-1+r1**2-r2**2+x2)
256  & -(-1+r1**4-6*r1**2*r2**2+r2**4+x1-r1**2*x1+r2**2*x1
257  & +x2+3*r1**2*x2-r2**2*x2-x1*x2)/(-1+r1**2-r2**2+x2)**2
258  isset4=1
259  ENDIF
260 
261 C...q -> q S.
262  ELSEIF(iclass.EQ.5) THEN
263  IF(icombi.EQ.1.OR.icombi.EQ.3) THEN
264  rlo1=ps*(1d0+r1**2-r2**2+2d0*r1)
265  rfo1=(4-4*r1**2+4*r2**2-3*x1-2*r1*x1+r1**2*x1-r2**2*x1-5*x2
266  & -2*r1*x2+r1**2*x2-r2**2*x2+x1*x2+x2**2)/(-2+x1+x2)**2
267  & +2*(3-r1-5*r1**2-r1**3+3*r2**2+r1*r2**2-2*x1-r1*x1
268  & +r1**2*x1-4*x2+2*r1**2*x2-r2**2*x2+x1*x2+x2**2)/
269  & (1-r1**2+r2**2-x2)/(-2+x1+x2)
270  & +(2-2*r1-6*r1**2-2*r1**3+2*r2**2-2*r1*r2**2-x1+r1**2*x1
271  & -r2**2*x1-3*x2+2*r1*x2+3*r1**2*x2-r2**2*x2+x1*x2+x2**2)/
272  & (-1+r1**2-r2**2+x2)**2
273  isset1=1
274  ENDIF
275  IF(icombi.EQ.2.OR.icombi.EQ.3) THEN
276  rlo2=ps*(1d0+r1**2-r2**2-2d0*r1)
277  rfo2=(4-4*r1**2+4*r2**2-3*x1+2*r1*x1+r1**2*x1-r2**2*x1-5*x2
278  & +2*r1*x2+r1**2*x2-r2**2*x2+x1*x2+x2**2)/(-2+x1+x2)**2
279  & +2*(3+r1-5*r1**2+r1**3+3*r2**2-r1*r2**2-2*x1+r1*x1
280  & +r1**2*x1-4*x2+2*r1**2*x2-r2**2*x2+x1*x2+x2**2)/
281  & (1-r1**2+r2**2-x2)/(-2+x1+x2)
282  & +(2+2*r1-6*r1**2+2*r1**3+2*r2**2+2*r1*r2**2-x1+r1**2*x1
283  & -r2**2*x1-3*x2-2*r1*x2+3*r1**2*x2-r2**2*x2+x1*x2+x2**2)/
284  & (-1+r1**2-r2**2+x2)**2
285  isset2=1
286  ENDIF
287  IF(icombi.EQ.4) THEN
288  rlo4=ps*(1d0+r1**2-r2**2)
289  rfo4=(4-4*r1**2+4*r2**2-3*x1+r1**2*x1-r2**2*x1-5*x2+r1**2*x2
290  & -r2**2*x2+x1*x2+x2**2)/(-2+x1+x2)**2
291  & +2*(3-5*r1**2+3*r2**2-2*x1+r1**2*x1-4*x2+2*r1**2*x2
292  & -r2**2*x2+x1*x2+x2**2)/(1-r1**2+r2**2-x2)/(-2+x1+x2)
293  & +(2-6*r1**2+2*r2**2-x1+r1**2*x1-r2**2*x1-3*x2+3*r1**2*x2
294  & -r2**2*x2+x1*x2+x2**2)/(-1+r1**2-r2**2+x2)**2
295  isset4=1
296  ENDIF
297 
298 C...V -> ~q ~qbar (~q = squark).
299  ELSEIF(iclass.EQ.6) THEN
300  rlo1=ps*(1d0-2d0*r1**2+r1**4-2d0*r2**2-2d0*r1**2*r2**2+r2**4)
301  rfo1=2d0*3d0+(1+r1**2+r2**2-x1)*(4*r1**2-x1**2)/
302  & (-1-r1**2+r2**2+x1)**2
303  & -2d0*(-1-3*r1**2-r2**2+x1+x1**2/2+x2-x1*x2/2)/
304  & (-1-r1**2+r2**2+x1)
305  & +(1+r1**2+r2**2-x2)*(4*r2**2-x2**2)
306  & /(-1+r1**2-r2**2+x2)**2
307  & -2d0*(-1-r1**2-3*r2**2+x1+x2-x1*x2/2+x2**2/2)/
308  & (-1+r1**2-r2**2+x2)
309  & -(-4*r1**2-4*r1**4-4*r2**2-8*r1**2*r2**2-4*r2**4+2*x1
310  & +6*r1**2*x1+6*r2**2*x1-2*x1**2+2*x2+6*r1**2*x2+6*r2**2*x2
311  & -4*x1*x2-2*r1**2*x1*x2-2*r2**2*x1*x2+x1**2*x2-2*x2**2
312  & +x1*x2**2)/(-1-r1**2+r2**2+x1)/(-1+r1**2-r2**2+x2)
313  isset1=1
314 
315 C...~q -> ~q V.
316  ELSEIF(iclass.EQ.7) THEN
317  rlo1=ps*(1d0-2d0*r1**2+r1**4-2d0*r2**2-2d0*r1**2*r2**2+r2**4)
318  rfo1=16*r2**2+8*(4*r2**2+2*r2**2*x1+x2+r1**2*x2+r2**2*x2-x1*x2
319  & -2*x2**2)/(3*(-1+r1**2-r2**2+x2))+8*(1+r1**2+r2**2-x2)*
320  & (4*r2**2-x2**2)/(3*(-1+r1**2-r2**2+x2)**2)+8*(x1+x2)*
321  & (-1-2*r1**2-r1**4-2*r2**2+2*r1**2*r2**2-r2**4+2*x1
322  & +2*r1**2*x1+2*r2**2*x1-x1**2+2*x2+2*r1**2*x2+2*r2**2*x2
323  & -2*x1*x2-x2**2)/(3*(-2+x1+x2)**2)+8*(-1-r1**2+r2**2-x1)*
324  & (2*r2**2*x1+x2+r1**2*x2+r2**2*x2-x1*x2-x2**2)/
325  & (3*(-1+r1**2-r2**2+x2)*(-2+x1+x2))+8*(1+2*r1**2+r1**4
326  & +2*r2**2-2*r1**2*r2**2+r2**4-2*x1-2*r1**2*x1-4*r2**2*x1
327  & +x1**2-3*x2-3*r1**2*x2-3*r2**2*x2+3*x1*x2+2*x2**2)/
328  & (3*(-2+x1+x2))
329  rfo1=3d0*rfo1/8d0
330  isset1=1
331 
332 C...S -> ~q ~qbar.
333  ELSEIF(iclass.EQ.8) THEN
334  rlo1=ps
335  rfo1=(-1-2*r1**2-r1**4-2*r2**2+2*r1**2*r2**2-r2**4+2*x1
336  & +2*r1**2*x1+2*r2**2*x1-x1**2-r2**2*x1**2+2*x2+2*r1**2*x2
337  & +2*r2**2*x2-3*x1*x2-r1**2*x1*x2-r2**2*x1*x2+x1**2*x2-x2**2
338  & -r1**2*x2**2+x1*x2**2)/
339  & (1+r1**2-r2**2-x1)**2/(-1+r1**2-r2**2+x2)**2
340  rfo1=2d0*rfo1
341  isset1=1
342 
343 C...~q -> ~q S.
344  ELSEIF(iclass.EQ.9) THEN
345  rlo1=ps
346  rfo1=(-1-r1**2-r2**2+x2)/(-1+r1**2-r2**2+x2)**2
347  & +(1+r1**2-r2**2+x1)/(-1+r1**2-r2**2+x2)/(-2+x1+x2)
348  & -(x1+x2)/(-2+x1+x2)**2
349  isset1=1
350 
351 C...chi -> q ~qbar (chi = neutralino/chargino).
352  ELSEIF(iclass.EQ.10) THEN
353  IF(icombi.EQ.1.OR.icombi.EQ.3) THEN
354  rlo1=ps*(1d0+r1**2-r2**2+2d0*r1)
355  rfo1=(2*r1+x1)*(-1-r1**2-r2**2+x1)/(-1-r1**2+r2**2+x1)**2
356  & +2*(-1-r1**2-2*r1**3-r2**2-2*r1*r2**2+3*x1/2+r1*x1
357  & -r1**2*x1/2-r2**2*x1/2+x2+r1*x2+r1**2*x2-x1*x2/2)/
358  & (-1-r1**2+r2**2+x1)/(-1+r1**2-r2**2+x2)
359  & +(2-2*r1-6*r1**2-2*r1**3+2*r2**2-2*r1*r2**2-x1+r1**2*x1
360  & -r2**2*x1-3*x2+2*r1*x2+3*r1**2*x2-r2**2*x2+x1*x2+x2**2)/
361  & (-1+r1**2-r2**2+x2)**2
362  isset1=1
363  ENDIF
364  IF(icombi.EQ.2.OR.icombi.EQ.3) THEN
365  rlo2=ps*(1d0-2d0*r1+r1**2-r2**2)
366  rfo2=(2*r1-x1)*(1+r1**2+r2**2-x1)/(-1-r1**2+r2**2+x1)**2
367  & +2*(-1-r1**2+2*r1**3-r2**2+2*r1*r2**2+3*x1/2-r1*x1
368  & -r1**2*x1/2-r2**2*x1/2+x2-r1*x2+r1**2*x2-x1*x2/2)/
369  & (-1-r1**2+r2**2+x1)/(-1+r1**2-r2**2+x2)
370  & +(2+2*r1-6*r1**2+2*r1**3+2*r2**2+2*r1*r2**2-x1+r1**2*x1
371  & -r2**2*x1-3*x2-2*r1*x2+3*r1**2*x2-r2**2*x2+x1*x2+x2**2)/
372  & (-1+r1**2-r2**2+x2)**2
373  isset2=1
374  ENDIF
375  IF(icombi.EQ.4) THEN
376  rlo4=ps*(1+r1**2-r2**2)
377  rfo4=x1*(-1-r1**2-r2**2+x1)/(-1-r1**2+r2**2+x1)**2
378  & +2d0*(-1-r1**2-r2**2+3*x1/2-r1**2*x1/2-r2**2*x1/2
379  & +x2+r1**2*x2-x1*x2/2)/
380  & (-1-r1**2+r2**2+x1)/(-1+r1**2-r2**2+x2)
381  & +(2-6*r1**2+2*r2**2-x1+r1**2*x1-r2**2*x1-3*x2+3*r1**2*x2
382  & -r2**2*x2+x1*x2+x2**2)/(-1+r1**2-r2**2+x2)**2
383  isset4=1
384  ENDIF
385 
386 C...~q -> q chi.
387  ELSEIF(iclass.EQ.11) THEN
388  IF(icombi.EQ.1.OR.icombi.EQ.3) THEN
389  rlo1=ps*(1d0-(r1+r2)**2)
390  rfo1=(1+r1**2+2*r1*r2+r2**2-x1-x2)*(x1+x2)/(-2+x1+x2)**2
391  & -(-1+r1**4-2*r1*r2-2*r1**3*r2-6*r1**2*r2**2-2*r1*r2**3
392  & +r2**4+x1-r1**2*x1+r2**2*x1+x2+3*r1**2*x2+2*r1*r2*x2
393  & -r2**2*x2-x1*x2)/(-1+r1**2-r2**2+x2)**2
394  & +(-1-2*r1**2-r1**4-2*r1*r2-2*r1**3*r2+2*r1*r2**3+r2**4
395  & +x1+r1**2*x1-2*r1*r2*x1-3*r2**2*x1+2*r1**2*x2-2*r2**2*x2
396  & +x1*x2+x2**2)/(-1+r1**2-r2**2+x2)/(-2+x1+x2)
397  isset1=1
398  ENDIF
399  IF(icombi.EQ.2.OR.icombi.EQ.3) THEN
400  rlo2=ps*(1d0-(r1-r2)**2)
401  rfo2=(1+r1**2-2*r1*r2+r2**2-x1-x2)*(x1+x2)/
402  & (-2+x1+x2)**2
403  & -(-1+r1**4+2*r1*r2+2*r1**3*r2-6*r1**2*r2**2+2*r1*r2**3
404  & +r2**4+x1-r1**2*x1+r2**2*x1+x2+3*r1**2*x2-2*r1*r2*x2
405  & -r2**2*x2-x1*x2)/(-1+r1**2-r2**2+x2)**2
406  & +(-1-2*r1**2-r1**4+2*r1*r2+2*r1**3*r2-2*r1*r2**3+r2**4
407  & +x1+r1**2*x1+2*r1*r2*x1-3*r2**2*x1+2*r1**2*x2-2*r2**2*x2
408  & +x1*x2+x2**2)/(-1+r1**2-r2**2+x2)/(-2+x1+x2)
409  isset2=1
410  ENDIF
411  IF(icombi.EQ.4) THEN
412  rlo4=ps*(1d0-r1**2-r2**2)
413  rfo4=(1+r1**2+r2**2-x1-x2)*(x1+x2)/(-2+x1+x2)**2
414  & -(-1+r1**4-6*r1**2*r2**2+r2**4+x1-r1**2*x1+r2**2*x1+x2
415  & +3*r1**2*x2-r2**2*x2-x1*x2)/
416  & (-1+r1**2-r2**2+x2)**2
417  & -(-1-2*r1**2-r1**4+r2**4+x1+r1**2*x1-3*r2**2*x1
418  & +2*r1**2*x2-2*r2**2*x2+x1*x2+x2**2)/
419  & (2-x1-x2)/(-1+r1**2-r2**2+x2)
420  isset4=1
421  ENDIF
422 
423 C...q -> ~q chi.
424  ELSEIF(iclass.EQ.12) THEN
425  IF(icombi.EQ.1.OR.icombi.EQ.3) THEN
426  rlo1=ps*(1d0-r1**2+r2**2+2d0*r2)
427  rfo1=(2*r2+x2)*(-1-r1**2-r2**2+x2)/(-1+r1**2-r2**2+x2)**2
428  & +(4+4*r1**2-4*r2**2-5*x1-r1**2*x1-2*r2*x1+r2**2*x1+x1**2
429  & -3*x2-r1**2*x2-2*r2*x2+r2**2*x2+x1*x2)/
430  & (-2+x1+x2)**2-2*(-1-r1**2+r2+r1**2*r2-r2**2-r2**3+x1
431  & +r2*x1+r2**2*x1+2*x2+r1**2*x2-x1*x2/2-x2**2/2)/
432  & (2-x1-x2)/(-1+r1**2-r2**2+x2)
433  isset1=1
434  END IF
435  IF(icombi.EQ.2.OR.icombi.EQ.3) THEN
436  rlo2=ps*(1d0-r1**2+r2**2-2d0*r2)
437  rfo2=(2*r2-x2)*(1+r1**2+r2**2-x2)/(-1+r1**2-r2**2+x2)**2
438  & +(4+4*r1**2-4*r2**2-5*x1-r1**2*x1+2*r2*x1+r2**2*x1+x1**2
439  & -3*x2-r1**2*x2+2*r2*x2+r2**2*x2+x1*x2)/
440  & (-2+x1+x2)**2-2*(-1-r1**2-r2-r1**2*r2-r2**2+r2**3+x1
441  & -r2*x1+r2**2*x1+2*x2+r1**2*x2-x1*x2/2-x2**2/2)/
442  & (2-x1-x2)/(-1+r1**2-r2**2+x2)
443  isset2=1
444  END IF
445  IF(icombi.EQ.4) THEN
446  rlo4=ps*(1d0-r1**2+r2**2)
447  rfo4=x2*(-1-r1**2-r2**2+x2)/(-1+r1**2-r2**2+x2)**2
448  & +(4+4*r1**2-4*r2**2-5*x1-r1**2*x1+r2**2*x1+x1**2
449  & -3*x2-r1**2*x2+r2**2*x2+x1*x2)/
450  & (-2+x1+x2)**2-2*(-1-r1**2-r2**2+x1+r2**2*x1+2*x2
451  & +r1**2*x2-x1*x2/2-x2**2/2)/
452  & (2-x1-x2)/(-1+r1**2-r2**2+x2)
453  isset4=1
454  END IF
455 
456 C...~g -> q ~qbar.
457  ELSEIF(iclass.EQ.13) THEN
458  IF(icombi.EQ.1.OR.icombi.EQ.3) THEN
459  rlo1=ps*(1d0+r1**2-r2**2+2d0*r1)
460  rfo1=4*(2*r1+x1)*(-1-r1**2-r2**2+x1)/(3*(-1-r1**2+r2**2+x1)**2)
461  & -(-1-r1**2-2*r1**3-r2**2-2*r1*r2**2+3*x1/2+r1*x1-r1**2*x1/2
462  & -r2**2*x1/2+x2+r1*x2+r1**2*x2-x1*x2/2)/(3*(-1-r1**2+r2**2
463  & +x1)*(-1+r1**2-r2**2+x2))-3*(-1+r1-r1**2-r1**3-r2**2
464  & +r1*r2**2+2*x1+r2**2*x1-x1**2/2+x2+r1*x2+r1**2*x2-x1*x2/2)/
465  & ((-1-r1**2+r2**2+x1)*(2-x1-x2))+3*(4-4*r1**2+4*r2**2-3*x1
466  & -2*r1*x1+r1**2*x1-r2**2*x1-5*x2-2*r1*x2+r1**2*x2-r2**2*x2
467  & +x1*x2+x2**2)/(-2+x1+x2)**2+3*(3-r1-5*r1**2-r1**3+3*r2**2
468  & +r1*r2**2-2*x1-r1*x1+r1**2*x1-4*x2+2*r1**2*x2-r2**2*x2
469  & +x1*x2+x2**2)/((1-r1**2+r2**2-x2)*(-2+x1+x2))+4*(2-2*r1
470  & -6*r1**2-2*r1**3+2*r2**2-2*r1*r2**2-x1+r1**2*x1-r2**2*x1
471  & -3*x2+2*r1*x2+3*r1**2*x2-r2**2*x2+x1*x2+x2**2)/
472  & (3*(-1+r1**2-r2**2+x2)**2)
473  rfo1=3d0*rfo1/4d0
474  isset1=1
475  ENDIF
476  IF(icombi.EQ.2.OR.icombi.EQ.3) THEN
477  rlo2=ps*(1d0+r1**2-r2**2-2d0*r1)
478  rfo2=4*(2*r1-x1)*(1+r1**2+r2**2-x1)/(3*(-1-r1**2+r2**2+x1)**2)
479  & -3*(-1-r1-r1**2+r1**3-r2**2-r1*r2**2+2*x1+r2**2*x1-x1**2/2
480  & +x2-r1*x2+r1**2*x2-x1*x2/2)/((-1-r1**2+r2**2+x1)*(2-x1-x2))
481  & +(2+2*r1**2-4*r1**3+2*r2**2-4*r1*r2**2-3*x1+2*r1*x1
482  & +r1**2*x1+r2**2*x1-2*x2+2*r1*x2-2*r1**2*x2+x1*x2)/
483  & (6*(-1-r1**2+r2**2+x1)*(-1+r1**2-r2**2+x2))+3*(4-4*r1**2
484  & +4*r2**2-3*x1+2*r1*x1+r1**2*x1-r2**2*x1-5*x2+2*r1*x2
485  & +r1**2*x2-r2**2*x2+x1*x2+x2**2)/(-2+x1+x2)**2+3*(3+r1
486  & -5*r1**2+r1**3+3*r2**2-r1*r2**2-2*x1+r1*x1+r1**2*x1-4*x2
487  & +2*r1**2*x2-r2**2*x2+x1*x2+x2**2)/
488  & ((1-r1**2+r2**2-x2)*(-2+x1+x2))+4*(2+2*r1-6*r1**2+2*r1**3
489  & +2*r2**2+2*r1*r2**2-x1+r1**2*x1-r2**2*x1-3*x2-2*r1*x2
490  & +3*r1**2*x2-r2**2*x2+x1*x2+x2**2)/
491  & (3*(-1+r1**2-r2**2+x2)**2)
492  rfo2=3d0*rfo2/4d0
493  isset2=1
494  ENDIF
495  IF(icombi.EQ.4) THEN
496  rlo4=ps*(1d0+r1**2-r2**2)
497  rfo4=8*x1*(-1-r1**2-r2**2+x1)/(3*(-1-r1**2+r2**2+x1)**2)-6*(-1
498  & -r1**2-r2**2+2*x1+r2**2*x1-x1**2/2+x2+r1**2*x2-x1*x2/2)/
499  & ((-1-r1**2+r2**2+x1)*(2-x1-x2))+(2+2*r1**2+2*r2**2-3*x1
500  & +r1**2*x1+r2**2*x1-2*x2-2*r1**2*x2+x1*x2)/(3*(-1-r1**2
501  & +r2**2+x1)*(-1+r1**2-r2**2+x2))+6*(4-4*r1**2+4*r2**2-3*x1
502  & +r1**2*x1-r2**2*x1-5*x2+r1**2*x2-r2**2*x2+x1*x2+x2**2)/
503  & (-2+x1+x2)**2+6*(3-5*r1**2+3*r2**2-2*x1+r1**2*x1-4*x2
504  & +2*r1**2*x2-r2**2*x2+x1*x2+x2**2)/
505  & ((1-r1**2+r2**2-x2)*(-2+x1+x2))+8*(2-6*r1**2+2*r2**2-x1
506  & +r1**2*x1-r2**2*x1-3*x2+3*r1**2*x2-r2**2*x2+x1*x2+x2**2)/
507  & (3*(-1+r1**2-r2**2+x2)**2)
508  rfo4=3d0*rfo4/8d0
509  isset4=1
510  ENDIF
511 
512 C...~q -> q ~g.
513  ELSEIF(iclass.EQ.14) THEN
514  IF(icombi.EQ.1.OR.icombi.EQ.3) THEN
515  rlo1=ps*(1-r1**2-r2**2-2d0*r1*r2)
516  rfo1=64*(1+r1**2+2*r1*r2+r2**2-x1-x2)*(x1+x2)/(9*(-2+x1+x2)**2)
517  & -16*(-1+r1**4-2*r1*r2-2*r1**3*r2-6*r1**2*r2**2-2*r1*r2**3
518  & +r2**4+x1-r1**2*x1+2*r1*r2*x1+3*r2**2*x1+x2+r1**2*x2
519  & -r2**2*x2-x1*x2)/(-1-r1**2+r2**2+x1)**2-16*(r1**2+r1**4
520  & -2*r1**3*r2+r2**2-6*r1**2*r2**2-2*r1*r2**3+r2**4
521  & -r1**2*x1+r1*r2*x1+2*r2**2*x1+2*r1**2*x2+r1*r2*x2-r2**2*x2
522  & -x1*x2)/((-1-r1**2+r2**2+x1)*(-1+r1**2-r2**2+x2))
523  & -64*(-1+r1**4-2*r1*r2-2*r1**3*r2-6*r1**2*r2**2-2*r1*r2**3
524  & +r2**4+x1-r1**2*x1+r2**2*x1+x2+3*r1**2*x2+2*r1*r2*x2
525  & -r2**2*x2-x1*x2)/(9*(-1+r1**2-r2**2+x2)**2)
526  & +8*(-1+r1**4-2*r1*r2+2*r1**3*r2-2*r2**2-2*r1*r2**3-r2**4
527  & -2*r1**2*x1+2*r2**2*x1+x1**2+x2-3*r1**2*x2-2*r1*r2*x2
528  & +r2**2*x2+x1*x2)/((-1-r1**2+r2**2+x1)*(-2+x1+x2))
529  rfo1=rfo1
530  & +8*(-1-2*r1**2-r1**4-2*r1*r2-2*r1**3*r2+2*r1*r2**3+r2**4
531  & +x1+r1**2*x1-2*r1*r2*x1-3*r2**2*x1+2*r1**2*x2-2*r2**2*x2
532  & +x1*x2+x2**2)/(9*(2-x1-x2)*(-1+r1**2-r2**2+x2))
533  rfo1=9d0*rfo1/64d0
534  isset1=1
535  ENDIF
536  IF(icombi.EQ.2.OR.icombi.EQ.3) THEN
537  rlo2=ps*(1-r1**2-r2**2+2d0*r1*r2)
538  rfo2=64*(1+r1**2-2*r1*r2+r2**2-x1-x2)*(x1+x2)/(9*(-2+x1+x2)**2)
539  & -16*(-1+r1**4+2*r1*r2+2*r1**3*r2-6*r1**2*r2**2+2*r1*r2**3
540  & +r2**4+x1-r1**2*x1-2*r1*r2*x1+3*r2**2*x1+x2+r1**2*x2
541  & -r2**2*x2-x1*x2)/(-1-r1**2+r2**2+x1)**2-64*(-1+r1**4
542  & +2*r1*r2+2*r1**3*r2-6*r1**2*r2**2+2*r1*r2**3+r2**4+x1
543  & -r1**2*x1+r2**2*x1+x2+3*r1**2*x2-2*r1*r2*x2-r2**2*x2
544  & -x1*x2)/(9*(-1+r1**2-r2**2+x2)**2)+16*(-r1**2-r1**4
545  & -2*r1**3*r2-r2**2+6*r1**2*r2**2-2*r1*r2**3-r2**4+r1**2*x1
546  & +r1*r2*x1-2*r2**2*x1-2*r1**2*x2+r1*r2*x2+r2**2*x2+x1*x2)/
547  & ((-1-r1**2+r2**2+x1)*(-1+r1**2-r2**2+x2))
548  rfo2=rfo2
549  & +8*(-1+r1**4+2*r1*r2-2*r1**3*r2-2*r2**2+2*r1*r2**3-r2**4
550  & -2*r1**2*x1+2*r2**2*x1+x1**2+x2-3*r1**2*x2+2*r1*r2*x2
551  & +r2**2*x2+x1*x2)/((-1-r1**2+r2**2+x1)*(-2+x1+x2))
552  & +8*(-1-2*r1**2-r1**4+2*r1*r2+2*r1**3*r2-2*r1*r2**3
553  & +r2**4+x1+r1**2*x1+2*r1*r2*x1-3*r2**2*x1+2*r1**2*x2
554  & -2*r2**2*x2+x1*x2+x2**2)/(9*(2-x1-x2)*(-1+r1**2-r2**2+x2))
555  rfo2=9d0*rfo2/64d0
556  isset2=1
557  ENDIF
558  IF(icombi.EQ.4) THEN
559  rlo4=ps*(1-r1**2-r2**2)
560  rfo4=128*(1+r1**2+r2**2-x1-x2)*(x1+x2)/(9*(-2+x1+x2)**2)-32*(-1
561  & +r1**4-6*r1**2*r2**2+r2**4+x1-r1**2*x1+3*r2**2*x1+x2
562  & +r1**2*x2-r2**2*x2-x1*x2)/(-1-r1**2+r2**2+x1)**2
563  & -32*(r1**2+r1**4+r2**2-6*r1**2*r2**2+r2**4-r1**2*x1
564  & +2*r2**2*x1+2*r1**2*x2-r2**2*x2-x1*x2)/
565  & ((-1-r1**2+r2**2+x1)*(-1+r1**2-r2**2+x2))-128*(-1+r1**4
566  & -6*r1**2*r2**2+r2**4+x1-r1**2*x1+r2**2*x1+x2+3*r1**2*x2
567  & -r2**2*x2-x1*x2)/(9*(-1+r1**2-r2**2+x2)**2)
568  & +16*(-1+r1**4-2*r2**2-r2**4-2*r1**2*x1+2*r2**2*x1+x1**2
569  & +x2-3*r1**2*x2+r2**2*x2+x1*x2)/
570  & ((-1-r1**2+r2**2+x1)*(-2+x1+ x2))
571  rfo4=rfo4+16*(-1-2*r1**2-r1**4+r2**4+x1+r1**2*x1-3*r2**2*x1
572  & +2*r1**2*x2-2*r2**2*x2+x1*x2+x2**2)/
573  & (9*(1-r1**2+r2**2-x2)*(-2+x1+x2))
574  rfo4=9d0*rfo4/128d0
575  isset4=1
576  ENDIF
577 
578 C...q -> ~q ~g.
579  ELSEIF(iclass.EQ.15) THEN
580  IF(icombi.EQ.1.OR.icombi.EQ.3) THEN
581  rlo1=ps*(1d0-r1**2+r2**2+2d0*r2)
582  rfo1=32*(2*r2+x2)*(-1-r1**2-r2**2+x2)/(9*(-1+r1**2-r2**2+x2)**2)
583  & +8*(-1-r1**2-2*r1**2*r2-r2**2-2*r2**3+x1+r2*x1+r2**2*x1
584  & +3*x2/2-r1**2*x2/2+r2*x2-r2**2*x2/2-x1*x2/2)/
585  & ((-1-r1**2+r2**2+x1)*(-1+r1**2-r2**2+x2))+8*(2+2*r1**2-2*r2
586  & -2*r1**2*r2-6*r2**2-2*r2**3-3*x1-r1**2*x1+2*r2*x1
587  & +3*r2**2*x1+x1**2-x2-r1**2*x2+r2**2*x2+x1*x2)/
588  & (-1-r1**2+r2**2+x1)**2+32*(4+4*r1**2-4*r2**2-5*x1
589  & -r1**2*x1-2*r2*x1+r2**2*x1+x1**2-3*x2-r1**2*x2-2*r2*x2
590  & +r2**2*x2+x1*x2)/(9*(-2+x1+x2)**2)
591  rfo1=rfo1+8*(3+3*r1**2-r2+r1**2*r2-5*r2**2-r2**3-4*x1-r1**2*x1
592  & +2*r2**2*x1+x1**2-2*x2-r2*x2+r2**2*x2+x1*x2)/
593  & ((-1-r1**2+r2**2+x1)*(2-x1-x2))+8*(-1-r1**2+r2+r1**2*r2
594  & -r2**2-r2**3+x1+r2*x1+r2**2*x1+2*x2+r1**2*x2-x1*x2/2
595  & -x2**2/2)/(9*(2-x1-x2)*(-1+r1**2-r2**2+x2))
596  rfo1=9d0*rfo1/32d0
597  isset1=1
598  END IF
599  IF(icombi.EQ.2.OR.icombi.EQ.3) THEN
600  rlo2=ps*(1d0-r1**2+r2**2-2d0*r2)
601  rfo2=32*(2*r2-x2)*(1+r1**2+r2**2-x2)/(9*(-1+r1**2-r2**2+x2)**2)
602  & +8*(-1-r1**2+2*r1**2*r2-r2**2+2*r2**3+x1-r2*x1+r2**2*x1
603  & +3*x2/2-r1**2*x2/2-r2*x2-r2**2*x2/2-x1*x2/2)/
604  & ((-1-r1**2+r2**2+x1)*(-1+r1**2-r2**2+x2))+8*(2+2*r1**2+2*r2
605  & +2*r1**2*r2-6*r2**2+2*r2**3-3*x1-r1**2*x1-2*r2*x1
606  & +3*r2**2*x1+x1**2-x2-r1**2*x2+r2**2*x2+x1*x2)/
607  & (-1-r1**2+r2**2+x1)**2+8*(3+3*r1**2+r2-r1**2*r2-5*r2**2
608  & +r2**3-4*x1-r1**2*x1+2*r2**2*x1+x1**2-2*x2+r2*x2+r2**2*x2
609  & +x1*x2)/((-1-r1**2+r2**2+x1)*(2-x1-x2))
610  rfo2=rfo2+32*(4+4*r1**2-4*r2**2-5*x1-r1**2*x1+2*r2*x1+r2**2*x1
611  & +x1**2-3*x2-r1**2*x2+2*r2*x2+r2**2*x2+x1*x2)/
612  & (9*(-2+x1+x2)**2)+8*(-1-r1**2-r2-r1**2*r2-r2**2+r2**3+x1
613  & -r2*x1+r2**2*x1+2*x2+r1**2*x2-x1*x2/2-x2**2/2)/
614  & (9*(2-x1-x2)*(-1+r1**2-r2**2+x2))
615  rfo2=9d0*rfo2/32d0
616  isset2=1
617  END IF
618  IF(icombi.EQ.4) THEN
619  rlo4=ps*(1d0-r1**2+r2**2)
620  rfo4=64*x2*(-1-r1**2-r2**2+x2)/(9*(-1+r1**2-r2**2+x2)**2)
621  & +16*(-1-r1**2-r2**2+x1+r2**2*x1+3*x2/2-r1**2*x2/2
622  & -r2**2*x2/2-x1*x2/2)/
623  & ((-1-r1**2+r2**2+x1)*(-1+r1**2-r2**2+x2))+16*(3+3*r1**2
624  & -5*r2**2-4*x1-r1**2*x1+2*r2**2*x1+x1**2-2*x2+r2**2*x2
625  & +x1*x2)/((-1-r1**2+r2**2+x1)*(2-x1-x2))
626  & +64*(4+4*r1**2-4*r2**2-5*x1-r1**2*x1+r2**2*x1+x1**2-3*x2
627  & -r1**2*x2+r2**2*x2+x1*x2)/(9*(-2+x1+x2)**2)
628  rfo4=rfo4+16*(2+2*r1**2-6*r2**2-3*x1-r1**2*x1+3*r2**2*x1+x1**2
629  & -x2-r1**2*x2+r2**2*x2+x1*x2)/(-1-r1**2+r2**2+x1)**2
630  & +16*(-1-r1**2-r2**2+x1+r2**2*x1+2*x2+r1**2*x2-x1*x2/2
631  & -x2**2/2)/(9*(2-x1-x2)*(-1+r1**2-r2**2+x2))
632  rfo4=9d0*rfo4/64d0
633  isset4=1
634  END IF
635 
636 C...g -> ~g ~g. Use (9/4)*eikonal. May be changed in the future.
637  ELSEIF(iclass.EQ.16) THEN
638  rlo=ps
639  IF(icombi.EQ.0.OR.icombi.EQ.1) THEN
640  anum=0d0
641  ELSEIF(icombi.EQ.2) THEN
642  anum=(2d0-x1-x2)**2
643  ELSEIF(icombi.EQ.3) THEN
644  anum=alpcor*(2d0-x1-x2)**2
645  ELSE
646  anum=0.5d0*(2d0-x1-x2)**2
647  ENDIF
648  rfo=ps*2d0*((x1+x2-1d0+anum-r1**2-r2**2)/
649  & ((1d0+r1**2-r2**2-x1)*(1d0+r2**2-r1**2-x2))-
650  & r1**2/(1d0+r2**2-r1**2-x2)**2-
651  & r2**2/(1d0+r1**2-r2**2-x1)**2)
652  rfo=9d0*rfo/4d0
653  icombi=0
654  ENDIF
655 
656 C...Find relevant LO and FO expression.
657  IF(icombi.EQ.0) THEN
658  ELSEIF(icombi.EQ.1.AND.isset1.EQ.1) THEN
659  rlo=rlo1
660  rfo=rfo1
661  ELSEIF(icombi.EQ.2.AND.isset2.EQ.1) THEN
662  rlo=rlo2
663  rfo=rfo2
664  ELSEIF(icombi.EQ.3.AND.isset1.EQ.1.AND.isset2.EQ.1) THEN
665  rlo=alpcor*rlo1+(1d0-alpcor)*rlo2
666  rfo=alpcor*rfo1+(1d0-alpcor)*rfo2
667  ELSEIF(isset4.EQ.1) THEN
668  rlo=rlo4
669  rfo=rfo4
670  ELSEIF(icombi.EQ.4.AND.isset1.EQ.1.AND.isset2.EQ.1) THEN
671  rlo=0.5d0*(rlo1+rlo2)
672  rfo=0.5d0*(rfo1+rfo2)
673  ELSEIF(isset1.EQ.1) THEN
674  rlo=rlo1
675  rfo=rfo1
676  ELSE
677  CALL pyerrm(16,'(PYMAEL:) not implemented ME code')
678  rlo=1d0
679  rfo=0d0
680  ENDIF
681 
682 C...Output.
683  pymael=rfo/rlo
684 
685  RETURN
686  END