StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
WeakShowerMEs.cc
1 // WeakShowerMEs.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Function definitions (not found in the header) for the
7 // WeakShowerMEs class.
8 
9 #include "Pythia8/WeakShowerMEs.h"
10 
11 namespace Pythia8 {
12 
13 //==========================================================================
14 
15 // The WeakShowerMEs class.
16 
17 //--------------------------------------------------------------------------
18 
19 // Calculate the 2 to 2 ME uG -> uG, with an overall factor of
20 // g_s^4 / 9 was removed.
21 
22 double WeakShowerMEs::getTchanneluGuGME(double sHat, double tHat,
23  double uHat) {
24 
25  double sH2 = sHat* sHat;
26  double sH3 = sH2 * sHat;
27  double tH2 = tHat * tHat;
28  double uH2 = uHat * uHat;
29  return (18.*sH3*uHat - 4*tH2*uH2 + 9*sHat*uH2*(tHat + 2*uHat) +
30  sH2*(-4.*tH2 + 9*tHat*uHat + 18*uH2))/(sHat*tH2*uHat);
31 
32 }
33 
34 //--------------------------------------------------------------------------
35 
36 // Calculate the 2 to 2 ME ud -> ud, with an overall factor of
37 // g_s^4 / 9 was removed.
38 
39 double WeakShowerMEs::getTchannelududME(double sHat, double tHat,
40  double uHat) {
41 
42  double sH2 = sHat * sHat;
43  double tH2 = tHat * tHat;
44  double uH2 = uHat * uHat;
45  return 4.*(sH2+uH2)/tH2;
46 
47 }
48 
49 //--------------------------------------------------------------------------
50 
51 // Calculate the 2 to 2 ME uu -> uu, with an overall factor of
52 // g_s^4 / 9 was removed.
53 
54 double WeakShowerMEs::getTchanneluuuuME(double sHat, double tHat,
55  double uHat) {
56 
57  double sH2 = sHat * sHat;
58  double tH2 = tHat * tHat;
59  double uH2 = uHat * uHat;
60  return 4./2.*((sH2+uH2)/tH2 + (sH2+tH2)/uH2 - 2.*sH2/(3.*tHat*uHat));
61 
62 }
63 
64 //--------------------------------------------------------------------------
65 
66 // Calculate the 2 to 3 ME uG -> uGZ, with an overall factor of
67 // \frac{(g_s^4 * 4\pi \alpha_{em} * (lU^2 + rU^2))}{
68 // 9 * \cos^2(\theta_W) \sin^2(\theta_W)}$ was removed.
69 // p1 = incoming quark, p2 = incoming gluon, p3 = outgoing gluon,
70 // p4 = outgoing Z, p5 = outgoing quark.
71 
72 double WeakShowerMEs::getTchanneluGuGZME(Vec4 p1, Vec4 p2, Vec4 p3,
73  Vec4 p4, Vec4 p5) {
74 
75  double p12 = p1*p2;
76  double p13 = p1*p3;
77  double p14 = p1*p4;
78  double p23 = p2*p3;
79  double p24 = p2*p4;
80  double p44 = p4*p4;
81  double p12P = (p1 + p2).m2Calc();
82  double p13M = (p1 - p3).m2Calc();
83  double p14M = (p1 - p4).m2Calc();
84  double p23M = (p2 - p3).m2Calc();
85  double p25M = (p2 - p5).m2Calc();
86  double p35P = (p3 + p5).m2Calc();
87  double p45P = (p4 + p5).m2Calc();
88  double p44S = p44*p44;
89 
90  double dia1 = (-4*p12*(p12* (p44 - 2*p24) -p23 * (p44 - 2*p24) +2*p13*p24))
91  / pow2(p12P * p45P);
92 
93  double dia2 = -(((p12 - p13 - p23) * (p12 * (p44 - 2*p14)
94  - p13 * (p44 - 2*p14) + 2*p14*p23)) / (p12P * p13M * pow2(p45P)));
95 
96  double dia3 = (8*(p44 - 2*p12) * p12 * (p12 - p23 - p24))
97  / (pow2(p12P) * p35P * p45P);
98 
99  double dia4 = -(8*pow3(p12) + 2*pow2(p12) * (p44 - 6*p13 - 6*p14 - 6*p23
100  - 4*p24) + 4*pow2(p13) * (p44 - 2*p14 - p24) - 2*p14*p23
101  * (-p44 + 2*p14 + 2*p23 + 4*p24) + 2*p13 * (-4*pow2(p14)
102  + 2*p14 * (p44 - 3*p23 - 3*p24) + 2*p23 * (p44 - p24) + p44*p24)
103  + p12* (p44S + 4*pow2(p13) + 4*pow2(p14) + 4*pow2(p23) - 4*p14
104  * (p44 - 4*p23 - 3*p24) - 4*p13* (p44 - 4*p14 - 4*p23 - 3*p24)
105  - 2*p44*p24 + 4*p23*p24)) / (4. * p12P * p13M * p25M * p45P);
106 
107  double dia5 = (-9*(2*pow3(p12) + (p13 + p23) * (-(p14* p23)
108  + p13 * (p44 - 2*p14 - p24)) + pow2(p12)* (2*p44 - 4*p13 - 3*p14
109  - 4*p23 - 3*p24) + p12 * (2*pow2(p13) + p23 * (-p44 + 4*p14
110  + 2*p23 + 3*p24) + p13 * (-2*p44 + 5*p14 + 4*p23 + 4*p24))))
111  / (p12P * p23M* pow2(p45P));
112 
113  double dia6 = (-4*(2*pow3(p12) + pow2(p12)* (7*p44 - 4*p13 - 6*p14
114  - 2*p23 - 4*p24) + p12 * (2*pow2(p13) + 4*pow2(p14) - 2*p44*p23
115  - 3*p44*p24 + 2*p23* p24 + 2*pow2(p24) + p14 * (-8*p44 + 6*p23 + 2*p24)
116  + p13* (-7*p44 + 10*p14 + 6*p23 + 6*p24)) + 2*(pow2(p13) * (p44
117  - 2*p14 - p24) + p14* (p14* (p44 - 2*p23) + p23 * (p44 - p24) + p44*p24)
118  + p13 * (-2*pow2(p14) + p23 * (p44 - 2*p24) + p14* (2*p44 - 3*p23
119  - 2*p24) + (p44 - p24)* p24)))) / (p12P * p14M * p35P * p45P);
120 
121  double dia7 = ((p12 - p13 - p14) * (p44 - 2*p23) * (p12 - p13 - p14
122  - p23 - p24)) / (p12P * p14M * p25M * p45P);
123 
124  double dia8 = (9*(2*pow2((p1*p2))* (p44 - 10*p14 - 16*p23 + 4*p24)
125  + p12* (p44S + 20*pow2(p14) - 6*p44*p23 + 16*pow2(p23)
126  + p13* (-18*p44 + 28*p14 + 32*p23 - 8*p24) + 2*p44*p24
127  + 16*p23* p24 - 8*pow2(p24) + p14* (-8*p44 + 52*p23 + 4*p24))
128  + 2*(4*pow2(p13)* (p44 - p14 - 2*p23) - p14*p23* (-3*p44 + 10*p14
129  + 8*p23 + 8*p24) + p13* (-4*pow2(p14) - 8*pow2(p23) + 4*p23* (p44 - p24)
130  + 2*p14* (2*p44 - 8*p23 + p24) + p24* (p44 + 4*p24)))))
131  / (8. * p12P * p23M * p14M * p45P);
132 
133  double dia9 = (-4*p13* (2*pow2(p12) + p12* (p44 - 4*p13 - 2*p14
134  - 4*p23 - 2*p24) + 2*(p13 + p23)* (p13 + p14 + p23 + p24)))
135  / (pow2(p13M) * pow2(p45P));
136 
137  double dia10 = (4*pow3(p12) - 4*pow2(p12) * (3*p13 + 2*p14 + 2*p23 + p24)
138  + p12* (p44S + 8*pow2(p13) + 4*pow2(p14) + 4*pow2(p23)
139  - 4*p14* (p44 - 3*p23 - 3*p24) - 2*p44*p24 + 4*p23*p24
140  + p13* (-6*p44 + 20*p14 + 4*p23 + 12*p24)) + 2*(p14*p23* (p44
141  - 2*p14 - 2*p23 - 4*p24) + 2*pow2(p13)* (p44 - 2*p14 - 2*p24)
142  + p13* (-4*pow2(p14) + 2*p14* (p44 - 2*p23 - 3*p24)
143  + 2*p23* (p44 - p24) + p44*p24))) / (4. * p12P * p13M * p35P* p45P);
144 
145  double dia11 = (4*p13* (p44 + 2*p13)* (p44 + 2*p12 - 2*p14 - 2*p24))
146  / (pow2(p13M)* p25M * p45P);
147 
148  double dia12 = (-9*(-2*pow3(p12) + 2*pow3(p13) - 2*p14* pow2(p23)
149  + pow2(p12)* (p44 + 6*p13 - 2*p14 + 4*p23 + 2*p24) + p13*p23
150  * (-p44 - 2*p14 + 2*p23 + 6*p24) + pow2(p13)* (p44 + 4*p23 + 6*p24)
151  - p12* (6*pow2(p13) + p23* (p44 - 4*p14 + 2*p23 + 2*p24)
152  + p13* (-2*p14 + 8*(p23 + p24))))) / (2.* p13M * p23M * pow2(p45P));
153 
154  double dia13 = -((p12 - p13 - p14)* (p44 - 2*p23)* (p44 - 2*p24))
155  / (2. * p13M * p14M * p35P * p45P);
156 
157  double dia14 = (2*(-8*pow2(p13)* (p44 - p14) - 2*pow2(p12)
158  * (p44 - 2*p14 + 2*p23 - 2*p24) + 2*p14* (-p44S - 2*pow2(p23)
159  + p23 * (p44 + 2*p14 - 2*p24) + 2*p44*p24) + p12* (p44S
160  - 4*pow2(p14) + 4*pow2(p23) + 2*p13* (3*p44 - 6*p14 - 2*p23 - 2*p24)
161  - 4*pow2(p24)) + p13* (-2*p44S + 8*pow2(p14) + 2*(p44 + 2*p23)* p24
162  + 4*pow2(p24) + p14* (-8*p44 + 8*p23 + 4*p24))))
163  / (p13M * p14M * p25M * p45P);
164 
165  double dia15 = (-9*(8*pow3(p12) + 2*pow2(p12)* (p44 - 8*p13 - 2*p14
166  - 4*p23 - 8*p24) + p12* (3*p44S + 8*pow2(p13) - 4*pow2(p14)
167  - 6*p44*p23 - 10*p44*p24 + 24*p23* p24 + 8*pow2(p24) + 4*p14
168  * (p23 + 3*p24) + 2*p13 * (3*p44 - 2*p14 + 8*p23 + 12*p24))
169  + 2*(p14*p23 * (p44 + 2*p14 - 8*p24) + 4*pow2(p13)* (p14
170  + p23 - p24) + p13* (4*pow2(p14) - 4*pow2(p23)
171  + 4*p23* (p44 - 4*p24) - 10*p14* p24 + (3*p44 - 4*p24)* p24))))
172  / (8. * p13M * p23M * p14M * p45P);
173 
174  double dia16 = (4*p12*p23* (p44 + 2*p12 - 2*p14 - 2*p24))
175  / (pow2(p12P) * pow2(p35P));
176 
177  double dia17 = ((p12 - p13 - p14)* (p44 - 2*p24)
178  * (p12 - p13 - p14 - p23 - p24)) / (p12P * p13M * p25M * p35P);
179 
180  double dia18 = (9*(-20*pow3(p12) + 4*pow2(p12)* (4*p44 + p13 + 2*p14
181  + p23) + 2*p14*p23* (p44 + 2*p14 - 4*p23 - 2*p24) + 8*pow2(p13)
182  * (p44 - p14 - p24) + p13* (-8*pow2(p14) + 8*p44*p23
183  + 4*p14* (2*p44 - 2*p23 - 5*p24) + 6*p44*p24 - 4*pow2(p24))
184  + p12* (p44S - 4*pow2(p14) - 14*p44*p23 + 8*pow2(p23)
185  - 8*p14* (p44 - p23 - 2*p24) - 12*p44*p24 + 12*p23* p24
186  + 4*pow2(p24) + p13 * (-22*p44 + 20*p14 - 8*p23 + 24*p24))))
187  / (8. * p12P * p23M * p35P * p45P);
188 
189  double dia19 = (-8*p13* (p44 + 2*p12 - 2*p14 - 2*p24)*(p12 - p14 - p24))
190  / (p12P * p14M * pow2(p35P));
191 
192  double dia20 = (-8*pow3(p12) - p12* (4*pow2(p13) + 4*pow2(p14)
193  - 4*p14* (p44 - 3*p23 - 4*p24) - 4*p13* (p44 - 4*p14 - 3*p23 - 4*p24)
194  + (p44 - 2*p24)* (p44 - 2*p23 - 2*p24)) - 2*pow2(p12)
195  * (p44 - 6*p13 - 6*p14 - 4*p23 - 6*p24) + 2*p14*p23
196  * (-p44 + 2*p14 + 2*p24) + pow2(p13) * (-4*p44 + 8*p14 + 4*p24)
197  + p13* (8*pow2(p14) - 4*p14* (p44 - 3*p23 - 3*p24) - 2*(p44 - 2*p24)
198  * (2*p23 + p24))) / (4. * p12P * p14M * p25M* p35P);
199 
200  double dia21 = (-9*(14*pow3(p12) + pow2(p12)* (7*p44 - 10*p13 - 30*p14
201  - 10*p23 - 22*p24) + 2*p12 * (2*pow2(p13) + 12*pow2(p14) - p44*p23
202  - 3*p44*p24 + 2*p23* p24 + 4*pow2(p24) + p13* (-4*p44 + 10*p14
203  + 6*p23 + 5*p24) + p14* (-5*p44 + 7*p23 + 16*p24)) + 2*(2*pow2(p13)
204  * (p44 - 2*p14 - p24) + p13* (-4*pow2(p14) + 2*p23* (p44 - 2*p24)
205  + 2*p14* (p44 - 3*p23 - 2*p24) + p44*p24) + p14* (-4*pow2(p14)
206  + 2*p14* (p44 - p23 - 3*p24) + p23* (p44 - 2*p24)
207  + 2*(p44 - p24)* p24)))) / (4. * p12P * p23M * p14M* p35P);
208 
209  double dia22 = (-8*p13*p23* (-p12 + p23 + p24))
210  / (pow2(p13M)* pow2(p25M));
211 
212  double dia23 = (-9*(2*pow3(p12) - pow2(p12)* (p44 + 10*p13 + 6*p14)
213  + p12* (p44S + 8*pow2(p13) + 4*pow2(p14) + 3*p44*p23 - 2*pow2(p23)
214  - p44*p24 - 4*p23 * p24 - 2*pow2(p24) + 2*p13* (6*p14 + 7*p23 + 2*p24)
215  + p14* (-2*p44 + 4*p23 + 6*p24)) - 2*(4*pow3(p13) + 2*pow2(p13)
216  * (p44 + 2*p14 + p23 - p24) + p14*p23* (2*p14 - p23 + p24)
217  + p13 * (p44S + 4*pow2(p14) + (-2*p44 + p23)* p24 - pow2(p24)
218  + p14* (-2*p44 + 4*p23 + 2*p24))))) / (4. * p13M * p23M * p25M * p45P);
219 
220  double dia24 = (-2*pow2(p12)* (p44 + p13 + p23 - p24) + p12* (2*pow2(p13)
221  - p13* (p44 - 6*p14) + (p44 + 2*p14 + 2*p23 - 2*p24)* (p23 + p24))
222  + 2*(pow2(p13)* (p44 - 2*p14 - p24) - p14*p23 * (p23 + p24)
223  + p13* (-2*pow2(p14) + p14 * (p44 - p23 - 2*p24) + p24* (p23 + p24))))
224  / (2. * p13M * p14M * p25M* p35P);
225 
226  double dia25 = (8*p12* (p44 + 2*p12 - 2*p23 - 2*p24)*(p12 - p23 - p24))
227  / (p13M * p14M * pow2(p25M));
228 
229  double dia26 = (9*(4*pow3(p12) + 2*pow2(p13)* (p44 - 2*p14 - 3*p24)
230  - 2*pow2(p12)* (p13 + 6*p23 + 2*p24) + p12* (p44S + 6*pow2(p13)
231  - 2*p44*p14 + 4*pow2(p14) - 2*p44*p23 + 8*pow2(p23) + 3*p13
232  * (p44 + 2*p14 - 2*p23 - 2*p24) - 2*p44*p24 + 8*p23*p24)
233  + p14* (-p44S + 2*p14* (p44 - 2*p23 - 2*p24) + 4*p23* p24
234  + 4*pow2(p24)) - p13 * (p44S + 4*pow2(p14) - 4*pow2(p23)
235  + 2*p23 * (p44 - 6*p24) + 2*p44*p24 - 8*pow2(p24)
236  + p14* (-4*p44 + 2*p23 + 8*p24)))) / (4. * p13M * p23M * p14M * p25M);
237 
238  double dia27 = (-9*(3*pow3(p12) - pow3(p13) + p14 * pow2(p23)
239  + pow2(p12)* (p44 - 7*p13 - 3*p14 - 6*p23 - 2*p24) - p13*p23
240  * (2*p44 + p23 - 2*p24) - pow2(p13)* (p14 + 2*p23 - 2*p24)
241  + p12 * (5*pow2(p13) + p13* (p44 + 4*p14 + 8*p23)
242  + p23 * (p44 + 2*p14 + 3*p23 + 2*p24)))) / (pow2(p23M) * pow2(p45P));
243 
244  double dia28 = (9*(2*pow3(p12) - pow2(p12)* (p44 + 4*p13 + 10*p14
245  + 8*p23 - 2*p24) + 2*pow2(p13)* (p44 - 2*p14 - p24) + p14* (p44S
246  + 2*p44*p23 + 2*pow2(p23) - 2*p14* (p44 + 6*p23) - 4*p44*p24)
247  + p12* (-p44S + 2*pow2(p13) + 8*pow2(p14) + p13 * (-5*p44 + 14*p14)
248  - 3*p44*p23 - 2*pow2(p23) + 2*p14* (p44 + 10*p23) + 6*p44*p24
249  + 2*p23* p24 - 4*pow2(p24)) + p13* (-4*pow2(p14) + pow2(p44 - 2*p24)
250  + 2*p23* (p44 + p24) + p14* (-6*p23 + 4*p24))))
251  / (4. * p23M * p14M * p35P * p45P);
252 
253  double dia29 = (-9*(6*pow3(p12) - 2*pow2(p13)* (p44 - 2*p14 + 3*p24)
254  - pow2(p12)* (p44 + 12*p13 + 2*p14 + 8*p23 + 10*p24) + p14* (p44S
255  - 2*p14* (p44 - 6*p23) - 6*p44*p23 - 2*pow2(p23) - 4*p44*p24)
256  + p13* (p44S + 4*pow2(p14) - 8*pow2(p23) - 2*p44* p24 - 4*pow2(p24)
257  - 2*p14* (2*p44 - 5*p23 + 4*p24) - 2*p23* (2*p44 + 5*p24))
258  + p12* (6* pow2(p13) - 4*pow2(p14) + 3*p44*p23 + 2*pow2(p23)
259  + 6*p23* p24 + 4*pow2(p24) + p14 * (6*p44 - 4*p23 + 4*p24)
260  + p13* (7*p44 - 2*p14 + 16*p23 + 16*p24))))
261  / (4. * p23M * p14M * p25M * p45P);
262 
263  double dia30 = (-9*(4*pow3(p12) + 2*pow2(p12)* (p44 - 4*p13 + 2*p14
264  - 6*p24) + p12* (p44S + 4*pow2(p13) - 8*pow2(p14) - 6*p44*p23
265  + 12*pow2(p23) - 6*p44 *p24 + 4*p23* p24 + 8*pow2(p24) + 2*p14
266  * (p44 - 4*p23 + 2*p24) + 8*p13* (p44 - p14 + 2*p23 + 2*p24))
267  - 2*(-2*p14* (p44 + 2*p14 - 3*p23)* p23 + pow2(p13)* (p44 - 2*p14
268  + 2*p24) + p13 * (-2*pow2(p14) + 8*pow2(p23) - 3*p23* (p44 - 2*p24)
269  + p24* (-p44 + 4*p24) + p14* (p44 + 6*p24)))))
270  / (2. * pow2(p23M) * p14M * p45P);
271 
272  double dia31 = (-2*(p13* (p44 - 2*p14) + p14* (p44 + 2*p12 - 2*p14
273  - 2*p23 - 2*p24))* (p44 + 2*p12 - 2*p14 - 2*p24))
274  / (pow2(p14M) * pow2(p35P));
275 
276  double dia32 = -((p44 - 2*p14)* (p12* (p44 - 2*p14) - p13* (p44 - 2*p14)
277  + 2*p14*p23)) / (2. * pow2(p14M) * p25M * p35P);
278 
279  double dia33 = (-9*(4*pow3(p14) + 2*pow2(p12)* (p44 + 2*p14)
280  - 4*pow2(p14)* (p44 + p23 - 3*p24) - 2*p44*p13* p24
281  + p14* (p44S + 2*p44*p23 + (-6*p44 + 4*p13)* p24 + 8*pow2(p24))
282  + p12 * (2*p13* (p44 - 2*p14) + p44*(3*p44 + 2*p23 - 2*p24)
283  - 2*p14* (3*p44 + 2*p23 + 6*p24)))) / (4. * p23M * pow2(p14M) * p35P);
284 
285  double dia34 = (-4*(p12 - p23 - p24)* (p44*p12-2*p14*p24))
286  / (pow2(p14M)* pow2(p25M));
287 
288  double dia35 = (-9*(2*pow2(p12)* (p44 + 2*p14) - 2*p13 * (p44 - 2*p14)
289  * (p44 - 2*p14 + 2*p23 + p24) + 2*p14* (4*pow2(p23) + 8*p23* p24
290  + p24* (-p44 + 2*p14 + 4*p24)) + p12* (2*p13* (p44 - 2*p14)
291  + p44*(p44 - 2*p23 - 2*p24) - 2*p14* (p44 + 6*p23 + 6*p24))))
292  / (4. * p23M * pow2(p14M) * p25M);
293 
294  double dia36 = (-9*(4*pow2(p12)* (p44 + 2*p14) - p13 * (p44 - 2*p14)
295  * (p44 - 2*p14 + 8*p23 + 4*p24) + p14* (-p44S - 4*pow2(p14)
296  + 16*pow2(p23) - 4*p23 * (p44 - 4*p24) - 8*p44*p24 + 16*pow2(p24)
297  + 4*p14* (p44 + 2*p23 + 4*p24)) + p12* (4*p13 * (p44 - 2*p14)
298  - 4*pow2(p14) + p44*(3*p44 + 4*p23 - 4*p24) - 4*p14* (p44 + 6*p23
299  + 6*p24)))) / (4. * pow2(p23M) * pow2(p14M));
300 
301  return dia1 + dia2 + dia3 + dia4 + dia5 + dia6 + dia7 + dia8 + dia9
302  + dia10 + dia11 + dia12 + dia13 + dia14 + dia15 + dia16 + dia17 + dia18
303  + dia19 + dia20 + dia21 + dia22 + dia23 + dia24 + dia25 + dia26 + dia27
304  + dia28 + dia29 + dia30 + dia31 + dia32 + dia33 + dia34 + dia35 + dia36;
305 
306 }
307 
308 //--------------------------------------------------------------------------
309 
310 // Calculate the 2 to 3 ME ud -> udZ, with the coupling between Z and d
311 // set to zero. An overall factor of
312 // \frac{(g_s^4 * 4\pi \alpha_{em} * (lU^2 + rU^2))}{
313 // 9 * \cos^2(\theta_W) \sin^2(\theta_W)}$ was removed.
314 // p1 = incoming u, p2 = incoming d, p3 = outgoing Z,
315 // p4 = outgoing d, p5 = outgoing u.
316 
317 double WeakShowerMEs::getTchannelududZME(Vec4 p1,Vec4 p2,Vec4 p3,
318  Vec4 p4,Vec4 p5) {
319 
320  double p12 = p1*p2;
321  double p13 = p1*p3;
322  double p14 = p1*p4;
323  double p23 = p2*p3;
324  double p24 = p2*p4;
325  double p33 = p3*p3;
326  double p13M = (p1 - p3).m2Calc();
327  double p24M = (p2 - p4).m2Calc();
328  double p35P = (p3 + p5).m2Calc();
329  double p33S = p33*p33;
330 
331  double dia1 = (-4*(2*pow3(p12) + pow2(p12) * (p33 - 2*p13 - 4*p14
332  - 2*p23 - 4*p24) + p14* (2*p14* p23 - (p33 - 2*p23)* p24)
333  + p12 * (2*pow2(p14) + 2*p24* (p13 + p23 + p24)
334  + p14* (p33 + 2*p13 + 4*p24)))) / pow2(p24M* p35P);
335 
336  double dia2 = (-2*(4*pow3(p12) + 4*pow2(p12)
337  * (p33 - 2*p14 - 3*p23) + p12 * (p33S - 4*pow2(p13) + 4*pow2(p14)
338  - 6*p33*p23 + 8*pow2(p23) + 4*p13 * (p23 - p24) - 4*p33*p24 + 4*p23
339  * p24 + 4*pow2(p24) + 4*p14* (p33 + 4*p23 + 4*p24)) + 2*(-2*pow2(p14)
340  * p23 + p13* (p33 + 2*p13 - 2*p24) * p24 + p14* (-4*pow2(p23) + p23
341  * (p33 - 6*p13 - 6*p24) + 2*(p33 - p13 - 2*p24)* p24))))
342  / (p13M * pow2(p24M) * p35P);
343 
344  double dia3 = (-2*(2*pow2(p12)* (p33 + 2*p13)
345  - 2*p33*p14* (p23 + p24) + 4*pow2(p13)* (2*p23 + p24) + p12
346  * (-4*pow2(p13) + p33*(p33 + 2*p14 - 2*p23) - 4*p13* (p14 + 3*p23
347  + 2*p24)) + p13* (8*pow2(p23) + 2*p24 * (-p33 + 2*p14 + 2*p24) + p23
348  * (-4*p33 + 4*p14 + 8*p24)))) / pow2(p13M * p24M);
349 
350  return dia1 + dia2 + dia3;
351 
352 }
353 
354 //==========================================================================
355 
356 } // end namespace Pythia8