a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
PVfunctions.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2012 HEPfit Collaboration
3 *
4 *
5 * For the licensing terms see doc/COPYING.
6 */
7
8#include <iostream>
9#include <stdexcept>
10#include <cmath>
11#include <algorithm>
12#include <gsl/gsl_complex.h>
13#include <gsl/gsl_sf.h>
14#include "PVfunctions.h"
15
16
17PVfunctions::PVfunctions(const bool bExtraMinusSign)
18{
19 if (bExtraMinusSign) ExtraMinusSign = -1.0;
20 else ExtraMinusSign = 1.0;
21}
22
23double PVfunctions::A0(const double mu2, const double m2) const
24{
25#ifdef USE_LOOPTOOLS
26 return ( ExtraMinusSign * myLT.PV_A0(mu2, m2) );
27#else
28 if ( mu2<=0.0 || m2<0.0 )
29 throw std::runtime_error("PVfunctions::A0(): Invalid argument!");
30
31 const double Tolerance = 1.0e-10;
32 const bool m2zero = (m2 <= mu2*Tolerance);
33
34 if ( m2zero )
35 return ( 0.0 );
36 else
37 return ( ExtraMinusSign * m2*(-log(m2/mu2)+1.0) );
38#endif
39}
40
41gslpp::complex PVfunctions::B0(const double mu2, const double p2,
42 const double m02, const double m12) const
43{
44#ifdef USE_LOOPTOOLS
45 return myLT.PV_B0(mu2, p2, m02, m12);
46#else
47 if ( mu2<=0.0 || p2<0.0 || m02<0.0 || m12<0.0 )
48 throw std::runtime_error("PVfunctions::B0(): Invalid argument!");
49
50 const double Tolerance = 1.0e-8;
51 const double maxM2 = std::max(p2, std::max(m02, m12));
52 const bool p2zero = (p2 <= maxM2*Tolerance);
53 const bool m02zero = (m02 <= maxM2*Tolerance);
54 const bool m12zero = (m12 <= maxM2*Tolerance);
55 const bool m02_eq_m12 = (fabs(m02 - m12) <= (m02 + m12)*Tolerance);
56
57 gslpp::complex B0(0.0, 0.0, false);
58 if ( p2zero ) {
59 if ( !m02zero && !m12zero ) {
60 if ( m02_eq_m12 )
61 B0 = - log(m02/mu2);
62 else
63 B0 = - m02/(m02-m12)*log(m02/mu2) + m12/(m02-m12)*log(m12/mu2) + 1.0;
64 } else if ( m02zero && !m12zero )
65 B0 = - log(m12/mu2) + 1.0;
66 else if ( !m02zero && m12zero )
67 B0 = - log(m02/mu2) + 1.0;
68 else if ( m02zero && m12zero )
69 throw std::runtime_error("PVfunctions::B0(): IR divergent! (vanishes in DR)");
70 else
71 throw std::runtime_error("PVfunctions::B0(): Undefined!");
72 } else {
73 if ( !m02zero && !m12zero ) {
74 double m0 = sqrt(m02), m1 = sqrt(m12);
75 double Lambda = sqrt( fabs((p2-m02-m12)*(p2-m02-m12) - 4.0*m02*m12) );
76 double R;
77 if ( p2 > (m0-m1)*(m0-m1) && p2 < (m0+m1)*(m0+m1) ) {
78 if ( p2-m02-m12 > 0.0 ) {
79 if ( p2-m02-m12 > Lambda*Tolerance )
80 R = - Lambda/p2*(atan(Lambda/(p2-m02-m12)) - M_PI);
81 else
82 R = - Lambda/p2*(M_PI/2.0 - M_PI);
83 } else {
84 if ( - (p2-m02-m12) > Lambda*Tolerance )
85 R = - Lambda/p2*atan(Lambda/(p2-m02-m12));
86 else
87 R = - Lambda/p2*( -M_PI/2.0 );
88 }
89 } else
90 R = Lambda/p2*log( fabs((p2-m02-m12+Lambda)/2.0/m0/m1) );
91 B0 = - log(m0*m1/mu2) + (m02-m12)/2.0/p2*log(m12/m02) - R + 2.0;
92 if ( p2 > (m0+m1)*(m0+m1) )
93 B0 += M_PI*Lambda/p2*gslpp::complex::i();// imaginary part
94 } else if ( (m02zero && !m12zero) || (!m02zero && m12zero) ) {
95 double M2;
96 if (!m02zero) M2 = m02;
97 else M2 = m12;
98 B0 = - log(M2/mu2) + 2.0;
99 if ( p2 < M2 )
100 B0 += - (1.0 - M2/p2)*log(1.0 - p2/M2);
101 else if ( p2 > M2 ) {
102 B0 += - (1.0 - M2/p2)*log(p2/M2 - 1.0);
103 B0 += (1.0 - M2/p2)*M_PI*gslpp::complex::i();// imaginary part
104 } else
105 B0 += 0.0;
106 } else if ( m02zero && m12zero ) {
107 if ( p2 < 0.0 )
108 B0 = - log(-p2/mu2) + 2.0;
109 else {
110 B0 = - log(p2/mu2) + 2.0;
111 B0 += M_PI*gslpp::complex::i();// imaginary part
112 }
113 } else
114 throw std::runtime_error("PVfunctions::B0(): Undefined!");
115 }
116 return B0;
117#endif
118}
119
120gslpp::complex PVfunctions::B1(const double mu2, const double p2,
121 const double m02, const double m12) const
122{
123#ifdef USE_LOOPTOOLS
124 return myLT.PV_B1(mu2, p2, m02, m12);
125#else
126 if ( mu2<=0.0 || p2<0.0 || m02<0.0 || m12<0.0 )
127 throw std::runtime_error("PVfunctions::B1(): Invalid argument!");
128
129 const double Tolerance = 1.0e-8;
130 const double maxM2 = std::max(p2, std::max(m02, m12));
131 const bool p2zero = (p2 <= maxM2*Tolerance);
132 const bool m02zero = (m02 <= maxM2*Tolerance);
133 const bool m12zero = (m12 <= maxM2*Tolerance);
134 const bool m02_eq_m12 = (fabs(m02 - m12) <= (m02 + m12)*Tolerance);
135
136 gslpp::complex B1(0.0, 0.0, false);
137 double DeltaM2 = m02 - m12;
138 if ( p2zero ) {
139 if ( !m02zero && !m12zero ) {
140 if ( m02_eq_m12 )
141 B1.real() = 1.0/2.0*log(m12/mu2);
142 else {
143 double F0 = - log(m12/m02);
144 double F1 = - 1.0 + m02/DeltaM2*F0;
145 double F2 = - 1.0/2.0 + m02/DeltaM2*F1;
146 B1.real() = 1.0/2.0*( log(m12/mu2) + F2 );
147 }
148 } else if ( m02zero && !m12zero )
149 B1.real() = 1.0/2.0*log(m12/mu2) - 1.0/4.0;
150 else if ( !m02zero && m12zero )
151 B1.real() = 1.0/2.0*log(m02/mu2) - 1.0/4.0;
152 else
153 B1 = 0.0;
154 } else
155 B1 = -1.0/2.0/p2*(- ExtraMinusSign*A0(mu2,m02)
156 + ExtraMinusSign*A0(mu2,m12)
157 + (DeltaM2 + p2)*B0(mu2,p2,m02,m12));
158 return B1;
159#endif
160}
161
162gslpp::complex PVfunctions::B11(const double mu2, const double p2,
163 const double m02, const double m12) const
164{
165#ifdef USE_LOOPTOOLS
166 return myLT.PV_B11(mu2, p2, m02, m12);
167#else
168 if ( mu2<=0.0 || p2<0.0 || m02<0.0 || m12<0.0 )
169 throw std::runtime_error("PVfunctions::B11(): Invalid argument!");
170
171 const double Tolerance = 1.0e-8;
172 const double maxM2 = std::max(p2, std::max(m02, m12));
173 const bool p2zero = (p2 <= maxM2*Tolerance);
174 const bool m02zero = (m02 <= maxM2*Tolerance);
175 const bool m12zero = (m12 <= maxM2*Tolerance);
176 const bool m02_eq_m12 = (fabs(m02 - m12) <= (m02 + m12)*Tolerance);
177
178 gslpp::complex B11(0.0, 0.0, false);
179 double DeltaM2 = m02 - m12;
180 if ( p2zero ) {
181 if ( !m02zero && !m12zero ) {
182 if ( m02_eq_m12 )
183 B11.real() = - 1.0/3.0*log(m12/mu2);
184 else {
185 double F0 = - log(m12/m02);
186 double F1 = - 1.0 + m02/DeltaM2*F0;
187 double F2 = - 1.0/2.0 + m02/DeltaM2*F1;
188 double F3 = - 1.0/3.0 + m02/DeltaM2*F2;
189 B11.real() = - 1.0/3.0*( log(m12/mu2) + F3 );
190 }
191 } else if ( m02zero && !m12zero )
192 B11.real() = - 1.0/3.0*log(m12/mu2) + 1.0/9.0;
193 else if ( !m02zero && m12zero )
194 B11.real() = - 1.0/3.0*log(m02/mu2) + 1.0/9.0;
195 else
196 throw std::runtime_error("PVfunctions::B11(): Undefined!");
197 } else {
198 double Lambdabar2 = (p2-m02-m12)*(p2-m02-m12) - 4.0*m02*m12;
199 B11 = - (3.0*(m02 + m12) - p2)/18.0/p2
200 + (DeltaM2 + p2)/3.0/p2/p2*(-ExtraMinusSign)*A0(mu2,m02)
201 - (DeltaM2 + 2.0*p2)/3.0/p2/p2*(-ExtraMinusSign)*A0(mu2,m12)
202 + (Lambdabar2 + 3.0*p2*m02)/3.0/p2/p2*B0(mu2,p2,m02,m12);
203 }
204 return B11;
205#endif
206}
207
208gslpp::complex PVfunctions::B00(const double mu2, const double p2,
209 const double m02, const double m12) const
210{
211#ifdef USE_LOOPTOOLS
212 return myLT.PV_B00(mu2, p2, m02, m12);
213#else
214 if ( mu2<=0.0 || p2<0.0 || m02<0.0 || m12<0.0 )
215 throw std::runtime_error("PVfunctions::B00(): Invalid argument!");
216
217 const double Tolerance = 1.0e-8;
218 const double maxM2 = std::max(p2, std::max(m02, m12));
219 const bool p2zero = (p2 <= maxM2*Tolerance);
220 const bool m02zero = (m02 <= maxM2*Tolerance);
221 const bool m12zero = (m12 <= maxM2*Tolerance);
222 const bool m02_eq_m12 = (fabs(m02 - m12) <= (m02 + m12)*Tolerance);
223
224 gslpp::complex B00(0.0, 0.0, false);
225 if ( p2zero ) {
226 if ( !m02zero && !m12zero ) {
227 if( m02_eq_m12 )
228 B00 = m02/2.0*(- log(m02/mu2) + 1.0);
229 else
230 B00 = 1.0/4.0*(m02 + m12)*(- log(sqrt(m02)*sqrt(m12)/mu2) + 3.0/2.0)
231 - (m02*m02 + m12*m12)/8.0/(m02 - m12)*log(m02/m12);
232 } else if ( (!m02zero && m12zero) || (m02zero && !m12zero) ) {
233 double M2;
234 if ( !m02zero ) M2 = m02;
235 else M2 = m12;
236 B00 = M2/4.0*(- log(M2/mu2) + 3.0/2.0);
237 } else
238 B00 = 0.0;
239 } else {
240 if ( !m02zero && !m12zero ) {
241 if ( m02_eq_m12 )
242 B00 = (6.0*m02 - p2)/18.0 + ExtraMinusSign*A0(mu2,m02)/6.0
243 - (p2 - 4.0*m02)/12.0*B0(mu2,p2,m02,m12);
244 else {
245 double DeltaM2 = m02 - m12;
246 double Lambdabar2 = (p2-m02-m12)*(p2-m02-m12) - 4.0*m02*m12;
247 B00 = (3.0*(m02 + m12) - p2)/18.0
248 - (DeltaM2 + p2)/12.0/p2*(-ExtraMinusSign)*A0(mu2,m02)
249 + (DeltaM2 - p2)/12.0/p2*(-ExtraMinusSign)*A0(mu2,m12)
250 - Lambdabar2*B0(mu2,p2,m02,m12)/12.0/p2;
251 }
252 } else if ( (!m02zero && m12zero) || (m02zero && !m12zero) ) {
253 double M2;
254 if ( !m02zero ) M2 = m02;
255 else M2 = m12;
256 B00 = (3.0*M2 - p2)/18.0 - (M2 + p2)/12.0/p2*(-ExtraMinusSign)*A0(mu2,M2)
257 - (M2 - p2)*(M2 - p2)/12.0/p2*B0(mu2,p2,M2,0.0);
258 } else
259 B00 = - p2/18.0 - p2/12.0*B0(mu2,p2,0.0,0.0);
260 }
261 return B00;
262#endif
263}
264
265gslpp::complex PVfunctions::Bf(const double mu2, const double p2,
266 const double m02, const double m12) const
267{
268 if ( mu2<=0.0 || p2<0.0 || m02<0.0 || m12<0.0 )
269 throw std::runtime_error("PVfunctions::Bf(): Invalid argument!");
270
271 gslpp::complex Bf(0.0, 0.0, false);
272 Bf = 2.0*(B11(mu2,p2,m02,m12) + B1(mu2,p2,m02,m12));
273 return Bf;
274}
275
276gslpp::complex PVfunctions::B0p(const double muIR2, const double p2,
277 const double m02, const double m12) const
278{
279#ifdef USE_LOOPTOOLS
280 return myLT.PV_B0p(muIR2, p2, m02, m12);
281#else
282 if ( muIR2<=0.0 || p2<0.0 || m02<0.0 || m12<0.0 )
283 throw std::runtime_error("PVfunctions::B0p(): Invalid argument!");
284
285 const double Tolerance = 1.0e-8;
286 const double maxM2 = std::max(p2, std::max(m02, m12));
287 const bool p2zero = (p2 <= maxM2*Tolerance);
288 const bool m02zero = (m02 <= maxM2*Tolerance);
289 const bool m12zero = (m12 <= maxM2*Tolerance);
290 const bool m02_eq_m12 = (fabs(m02 - m12) <= (m02 + m12)*Tolerance);
291
292 gslpp::complex B0p(0.0, 0.0, false);
293 if ( p2zero ) {
294 if ( !m02zero && !m12zero ) {
295 if ( m02_eq_m12 )
296 B0p = 1.0/6.0/m02;
297 else {
298 double DeltaM2 = m02 - m12;
299 B0p = (m02 + m12)/2.0/pow(DeltaM2,2.0)
300 + m02*m12/pow(DeltaM2,3.0)*log(m12/m02);
301 }
302 } else if ( !m02zero && m12zero )
303 B0p = 1.0/2.0/m02;
304 else if ( m02zero && !m12zero )
305 B0p = 1.0/2.0/m12;
306 else
307 throw std::runtime_error("PVfunctions::B0p(): Undefined!");
308 } else {
309 if ( !m02zero && !m12zero ) {
310 double m0 = sqrt(m02), m1 = sqrt(m12);
311 double Lambda = sqrt( fabs((p2-m02-m12)*(p2-m02-m12) - 4.0*m02*m12) );
312 double Rprime;
313 if ( p2 > (m0-m1)*(m0-m1) && p2 < (m0+m1)*(m0+m1) ) {
314 if ( p2-m02-m12 > 0.0 ) {
315 if ( p2-m02-m12 > Lambda*Tolerance )
316 Rprime = ((p2 - m02 - m12)/Lambda + Lambda/p2)
317 *(atan(Lambda/(p2-m02-m12)) - M_PI);
318 else
319 Rprime = ((p2 - m02 - m12)/Lambda + Lambda/p2)
320 *(M_PI/2.0 - M_PI);
321 } else {
322 if ( - (p2-m02-m12) > Lambda*Tolerance )
323 Rprime = ((p2 - m02 - m12)/Lambda + Lambda/p2)
324 *atan(Lambda/(p2-m02-m12));
325 else
326 Rprime = ((p2 - m02 - m12)/Lambda + Lambda/p2)
327 *( -M_PI/2.0 );
328 }
329 } else
330 Rprime = ((p2 - m02 - m12)/Lambda - Lambda/p2)
331 *log( fabs((p2-m02-m12+Lambda)/2.0/m0/m1) );
332 B0p = - (m02 - m12)/2.0/p2/p2*log(m12/m02) - (Rprime + 1.0)/p2;
333 if ( p2 > (m0+m1)*(m0+m1) )
334 B0p += M_PI/p2*((p2 - m02 - m12)/Lambda - Lambda/p2)
335 *gslpp::complex::i();// imaginary part
336 } else if ( (m02zero && !m12zero) || (!m02zero && m12zero) ) {
337 double M2;
338 if ( !m02zero ) M2 = m02;
339 else M2 = m12;
340 if ( p2 < M2 )
341 B0p = - M2/p2/p2*log(1.0 - p2/M2) - 1.0/p2;
342 else if ( p2 > M2 ) {
343 B0p = - M2/p2/p2*log(p2/M2 - 1.0) - 1.0/p2;
344 B0p += M2/p2/p2*M_PI*gslpp::complex::i();// imaginary part
345 } else /* p2=M2 */
346 B0p = 1.0/2.0/M2*(log(M2/muIR2) - 2.0);
347 } else if ( m02zero && m12zero )
348 B0p = - 1.0/p2;
349 else
350 throw std::runtime_error("PVfunctions::B0p(): Undefined!");
351 }
352 return B0p;
353#endif
354}
355
356gslpp::complex PVfunctions::B1p(const double mu2, const double p2,
357 const double m02, const double m12) const
358{
359#ifdef USE_LOOPTOOLS
360 return myLT.PV_B1p(mu2, p2, m02, m12);
361#else
362 if ( mu2<=0.0 || p2<0.0 || m02<0.0 || m12<0.0 )
363 throw std::runtime_error("PVfunctions::B1p(): Invalid argument!");
364
365 const double Tolerance = 1.0e-8;
366 const double maxM2 = std::max(p2, std::max(m02, m12));
367 const bool p2zero = (p2 <= maxM2*Tolerance);
368
369 gslpp::complex B1p(0.0, 0.0, false);
370 if ( p2zero )
371 throw std::runtime_error("PVfunctions::B1p(): Undefined!");
372 else {
373 double DeltaM2 = m02 - m12;
374 B1p = - ( 2.0*B1(mu2,p2,m02,m12) + B0(mu2,p2,m02,m12)
375 + (DeltaM2 + p2)*B0p(mu2,p2,m02,m12) )/2.0/p2;
376 }
377 return B1p;
378#endif
379}
380
381gslpp::complex PVfunctions::B11p(const double mu2, const double p2,
382 const double m02, const double m12) const
383{
384#ifdef USE_LOOPTOOLS
385 return myLT.PV_B11p(mu2, p2, m02, m12);
386#else
387 if ( mu2<=0.0 || p2<0.0 || m02<0.0 || m12<0.0 )
388 throw std::runtime_error("PVfunctions::B11p(): Invalid argument!");
389
390 const double Tolerance = 1.0e-8;
391 const double maxM2 = std::max(p2, std::max(m02, m12));
392 const bool p2zero = (p2 <= maxM2*Tolerance);
393
394 double p4 = p2*p2, p6=p2*p2*p2;
395 double DeltaM2 = m02 - m12;
396 gslpp::complex B11p(0.0, 0.0, false);
397 if ( p2zero )
398 throw std::runtime_error("PVfunctions::B11p(): Undefined!");
399 else {
400 double Lambdabar2 = (p2-m02-m12)*(p2-m02-m12) - 4.0*m02*m12;
401 B11p = (m02 + m12)/6.0/p4 - (2.0*DeltaM2 + p2)/3.0/p6*(-ExtraMinusSign)*A0(mu2,m02)
402 + 2.0*(DeltaM2 + p2)/3.0/p6*(-ExtraMinusSign)*A0(mu2,m12)
403 - (2.0*DeltaM2*DeltaM2 + p2*m02 - 2.0*p2*m12)/3.0/p6*B0(mu2,p2,m02,m12)
404 + (Lambdabar2 + 3.0*p2*m02)/3.0/p4*B0p(mu2,p2,m02,m12);
405 }
406 return B11p;
407#endif
408}
409
410gslpp::complex PVfunctions::B00p(const double mu2, const double p2,
411 const double m02, const double m12) const
412{
413#ifdef USE_LOOPTOOLS
414 return myLT.PV_B00p(mu2, p2, m02, m12);
415#else
416 if ( mu2<=0.0 || p2<0.0 || m02<0.0 || m12<0.0 )
417 throw std::runtime_error("PVfunctions::B00p(): Invalid argument!");
418
419 const double Tolerance = 1.0e-8;
420 const double maxM2 = std::max(p2, std::max(m02, m12));
421 const bool p2zero = (p2 <= maxM2*Tolerance);
422 const bool m02zero = (m02 <= maxM2*Tolerance);
423 const bool m12zero = (m12 <= maxM2*Tolerance);
424 const bool m02_eq_m12 = (fabs(m02 - m12) <= (m02 + m12)*Tolerance);
425
426 gslpp::complex B00p(0.0, 0.0, false);
427 double DeltaM2 = m02 - m12;
428 if ( p2zero ) {
429 if ( m02_eq_m12 )
430 B00p = - 1.0/18.0 - 1.0/12.0*B0(mu2,0.0,m02,m12)
431 + (m02 + m12)/6.0*B0p(mu2,0.0,m02,m12);
432 else if ( m02zero || m12zero )
433 B00p = - 1.0/18.0 - 1.0/12.0*B0(mu2,0.0,m02,m12)
434 + (m02 + m12)/6.0*B0p(mu2,0.0,m02,m12) - 1.0/72.0;
435 else
436 B00p = - 1.0/18.0 - 1.0/12.0*B0(mu2,0.0,m02,m12)
437 + (m02 + m12)/6.0*B0p(mu2,0.0,m02,m12)
438 - 1.0/24.0
439 *( (m02*m02 + 10.0*m02*m12 + m12*m12)/3.0/DeltaM2/DeltaM2
440 + 2.0*m02*m12*(m02 + m12)/DeltaM2/DeltaM2/DeltaM2*log(m12/m02) );
441 } else {
442 double Lambdabar2 = (p2-m02-m12)*(p2-m02-m12) - 4.0*m02*m12;
443 if ( m02_eq_m12 )
444 B00p = - 1.0/18.0 - B0(mu2,p2,m02,m12)/12.0
445 - Lambdabar2/12.0/p2*B0p(mu2,p2,m02,m12);
446 else
447 B00p = - 1.0/18.0
448 + DeltaM2/12.0/p2/p2*(- ExtraMinusSign*A0(mu2,m02)
449 + ExtraMinusSign*A0(mu2,m12)
450 + DeltaM2*B0(mu2,p2,m02,m12))
451 - B0(mu2,p2,m02,m12)/12.0
452 - Lambdabar2/12.0/p2*B0p(mu2,p2,m02,m12);
453 }
454 return B00p;
455#endif
456}
457
458gslpp::complex PVfunctions::Bfp(const double mu2, const double p2,
459 const double m02, const double m12) const
460{
461 if ( mu2<=0.0 || p2<0.0 || m02<0.0 || m12<0.0 )
462 throw std::runtime_error("PVfunctions::Bfp(): Invalid argument!");
463
464 gslpp::complex Bfp(0.0, 0.0, false);
465 Bfp = 2.0*(B11p(mu2,p2,m02,m12) + B1p(mu2,p2,m02,m12));
466 return Bfp;
467}
468
469gslpp::complex PVfunctions::C0(const double p1, const double p2, const double p1p22,
470 const double m02, const double m12, const double m22) const
471{
472#ifdef USE_LOOPTOOLS
473 return ( ExtraMinusSign * myLT.PV_C0(p1,p2,p1p22,m02,m12,m22) );
474#else
475 throw std::runtime_error("PVfunctions::C0(p1,p2,p1p22,m02,m12,m22): Only available with LoopTools.");
476#endif
477} //AG:added
478
479gslpp::complex PVfunctions::C0(const double p2,
480 const double m02, const double m12, const double m22) const
481{
482#ifdef USE_LOOPTOOLS
483 return ( ExtraMinusSign * myLT.PV_C0(p2, m02, m12, m22) );
484#else
485 if ( p2<0.0 || m02<0.0 || m12<0.0 || m22<0.0 )
486 throw std::runtime_error("PVfunctions::C0(): Invalid argument!");
487
488 const double Tolerance = 1.0e-8;
489 const double maxM2 = std::max(p2, std::max(m02, std::max(m12, m22)));
490 const bool p2zero = (p2 <= maxM2*Tolerance);
491 const bool m02zero = (m02 <= maxM2*Tolerance);
492 const bool m12zero = (m12 <= maxM2*Tolerance);
493 const bool m22zero = (m22 <= maxM2*Tolerance);
494 bool diff01 = (fabs(m02 - m12) > (m02 + m12)*Tolerance);
495 bool diff12 = (fabs(m12 - m22) > (m12 + m22)*Tolerance);
496 bool diff20 = (fabs(m22 - m02) > (m22 + m02)*Tolerance);
497
498 gslpp::complex C0(0.0, 0.0, false);
499 if ( p2zero ) {
500 if ( !m02zero && !m12zero && !m22zero ) {
501 if ( diff01 && diff12 && diff20 )
502 return ( - ( m12/(m02 - m12)*log(m12/m02)
503 - m22/(m02 - m22)*log(m22/m02) )/(m12 - m22) );
504 else if ( !diff01 && diff12 && diff20 )
505 return ( - (- m02 + m22 - m22*log(m22/m02))/(m02 - m22)/(m02 - m22) );
506 else if ( diff01 && !diff12 && diff20 )
507 return ( - ( m02 - m12 + m02*log(m12/m02))/(m02 - m12)/(m02 - m12) );
508 else if ( diff01 && diff12 && !diff20 )
509 return ( - (- m02 + m12 - m12*log(m12/m02))/(m02 - m12)/(m02 - m12) );
510 else
511 return ( 1.0/2.0/m02 );
512 }
513 } else {
514 if ( !diff20 && diff01 ) {
515 double epsilon = 1.0e-12;
516 gsl_complex tmp = gsl_complex_rect(1.0 - 4.0*m02/p2, epsilon);
517 tmp = gsl_complex_sqrt(tmp);
518 gslpp::complex tmp_complex(GSL_REAL(tmp), GSL_IMAG(tmp), false);
519 gslpp::complex x0 = 1.0 - (m02 - m12)/p2;
520 gslpp::complex x1 = (1.0 + tmp_complex)/2.0;
521 gslpp::complex x2 = (1.0 - tmp_complex)/2.0;
522 gslpp::complex x3 = m02/(m02 - m12);
523
524 if ( x0==x1 || x0==x2 || x0==x3)
525 throw std::runtime_error("PVfunctions::C0(): Undefined-2!");
526
527 gslpp::complex arg[6];
528 arg[0] = (x0 - 1.0)/(x0 - x1);
529 arg[1] = x0/(x0 - x1);
530 arg[2] = (x0 - 1.0)/(x0 - x2);
531 arg[3] = x0/(x0 - x2);
532 arg[4] = (x0 - 1.0)/(x0 - x3);
533 arg[5] = x0/(x0 - x3);
534
535 gslpp::complex Li2[6];
536 for (int i=0; i<6; i++) {
537 gsl_sf_result re, im;
538 gsl_sf_complex_dilog_xy_e(arg[i].real(), arg[i].imag(), &re, &im);
539 Li2[i].real() = re.val;
540 Li2[i].imag() = im.val;
541
542 /* Check the sizes of errors */
543 //std::cout << "re.val=" << re.val << " re.err=" << re.err << std::endl;
544 //std::cout << "im.val=" << im.val << " im.err=" << im.err << std::endl;
545 }
546 C0 = - 1.0/p2*( Li2[0] - Li2[1] + Li2[2] - Li2[3] - Li2[4] + Li2[5]);
547 } else if ( !m02zero && !m22zero && diff20 && m12zero ) {
548 double epsilon = 1.0e-12;
549 double tmp_real = pow((m02+m22-p2),2.0) - 4.0*m02*m22;
550 gsl_complex tmp = gsl_complex_rect(tmp_real, epsilon);
551 tmp = gsl_complex_sqrt(tmp);
552 gslpp::complex tmp_complex(GSL_REAL(tmp), GSL_IMAG(tmp), false);
553 gslpp::complex x1 = (p2 - m02 + m22 + tmp_complex)/2.0/p2;
554 gslpp::complex x2 = (p2 - m02 + m22 - tmp_complex)/2.0/p2;
555
556 if ( x1==0.0 || x1==1.0 || x2==0.0 || x2==1.0 )
557 throw std::runtime_error("PVfunctions::C0(): Undefined-3!");
558
559 gslpp::complex arg1 = (x1 - 1.0)/x1;
560 gslpp::complex arg2 = x2/(x2 - 1.0);
561 gsl_complex arg1_tmp = gsl_complex_rect(arg1.real(), arg1.imag());
562 gsl_complex arg2_tmp = gsl_complex_rect(arg2.real(), arg2.imag());
563 C0.real() = - 1.0/p2*( GSL_REAL(gsl_complex_log(arg1_tmp))
564 *GSL_REAL(gsl_complex_log(arg2_tmp))
565 - GSL_IMAG(gsl_complex_log(arg1_tmp))
566 *GSL_IMAG(gsl_complex_log(arg2_tmp)) );
567 C0.imag() = - 1.0/p2*( GSL_REAL(gsl_complex_log(arg1_tmp))
568 *GSL_IMAG(gsl_complex_log(arg2_tmp))
569 + GSL_IMAG(gsl_complex_log(arg1_tmp))
570 *GSL_REAL(gsl_complex_log(arg2_tmp)) );
571 } else if ( m02zero && !m12zero && !m22zero ) {
572 gslpp::complex arg[2];
573 arg[0] = 1.0 - m22/m12;
574 arg[1] = 1.0 - (-p2+m22)/m12;
575 gslpp::complex Li2[2];
576 for (int i=0; i<2; i++) {
577 gsl_sf_result re, im;
578 gsl_sf_complex_dilog_xy_e(arg[i].real(), arg[i].imag(), &re, &im);
579 Li2[i].real() = re.val;
580 Li2[i].imag() = im.val;
581 }
582 C0 = 1./(-p2)*(Li2[0]-Li2[1]);
583 } else if ( !m02zero && !m12zero && !m22zero && diff01 && diff12 ) {
584 double x0 = 1.0 - (m02-m12)/p2;
585 double x1 = -(-p2+m02-m22-sqrt(fabs((m02+m22-p2)*(m02+m22-p2) - 4.*m02*m22)))/p2/2.0;
586 double x2 = -(-p2+m02-m22+sqrt(fabs((m02+m22-p2)*(m02+m22-p2) - 4.*m02*m22)))/p2/2.0;
587 double x3 = m22/(m22-m12);
588
589 gslpp::complex arg[6];
590 arg[0] = (x0-1.0)/(x0-x1);
591 arg[1] = x0/(x0-x1);
592 arg[2] = (x0-1.0)/(x0-x2);
593 arg[3] = x0/(x0-x2);
594 arg[4] = (x0-1.0)/(x0-x3);
595 arg[5] = x0/(x0-x3);
596
597 gslpp::complex Li2[6];
598 for (int i=0; i<2; i++) {
599 gsl_sf_result re, im;
600 gsl_sf_complex_dilog_xy_e(arg[i].real(), arg[i].imag(), &re, &im);
601 Li2[i].real() = re.val;
602 Li2[i].imag() = im.val;
603 }
604
605 C0 = -1.0/p2*(Li2[0] - Li2[1] + Li2[2] - Li2[3] - Li2[4] + Li2[5]);
606 } else
607 throw std::runtime_error("PVfunctions::C0(): Undefined-4!");
608 }
609 return ( - ExtraMinusSign * C0 );
610#endif
611}
612
613double PVfunctions::C11(const double m12, const double m22, const double m32) const
614{
615 if ( m12<=0.0 || m22<=0.0 || m32<=0.0 )
616 throw std::runtime_error("PVfunctions::C11(): Argument is not positive!");
617
618 const double Tolerance = 2.5e-3;
619 double C11;
620
621 if ( 2.0*fabs(m12-m22) > (m12+m22)*Tolerance && 2.0*fabs(m12-m32) > (m12+m32)*Tolerance && 2.0*fabs(m22-m32) > (m22+m32)*Tolerance ) {
622 C11=(-m12*(m12-m22)*(m12-m32)*(m22-m32)
623 +m12*m12*m22*(2.0*m12-m22)*log(m12/m22)
624 +m12*m12*m32*(-2.0*m12+m32)*log(m12/m32)
625 +m22*m32*(-2.0*m12+m22)*(-2.0*m12+m32)*log(m22/m32))
626 /(2.0*(m12-m22)*(m12-m22)*(m12-m32)*(m12-m32)*(m22-m32));
627 }
628 else if ( 2.0*fabs(m12-m22) > (m12+m22)*Tolerance && 2.0*fabs(m12-m32) > (m12+m32)*Tolerance && 2.0*fabs(m22-m32) <= (m22+m32)*Tolerance ) {
629 C11=-((-12.0*m12*m12+8.0*m12*(m22+m32)-(m22+m32)*(m22+m32)+8.0*m12*m12*log((2.0*m12)/(m22+m32)))
630 /pow(-2.0*m12+m22+m32,3));
631 }
632 else if ( 2.0*fabs(m12-m22) > (m12+m22)*Tolerance && 2.0*fabs(m12-m32) <= (m12+m32)*Tolerance && 2.0*fabs(m22-m32) > (m22+m32)*Tolerance ) {
633 C11=(3.0*m12*m12-8.0*m12*m22+4.0*m22*m22+6.0*m12*m32-8.0*m22*m32+3.0*m32*m32+8.0*m22*(m12-m22+m32)*log((2.0*m22)/(m12+m32)))
634 /(2.0*pow(m12-2.0*m22+m32,3));
635 }
636 else if ( 2.0*fabs(m12-m22) <= (m12+m22)*Tolerance && 2.0*fabs(m12-m32) <= (m12+m32)*Tolerance && 2.0*fabs(m22-m32) <= (m22+m32)*Tolerance ) {
637 C11=1/(m12+m22+m32);
638 }
639 else {
640 C11=(3.0*m12*m12+6.0*m12*m22+3.0*m22*m22-8.0*m12*m32-8.0*m22*m32+4.0*m32*m32-8.0*m32*(m12+m22-m32)*log((m12+m22)/(2.0*m32)))
641 /(2.0*pow(m12+m22-2.0*m32,3));
642 }
643
644 return C11;
645}
646
647double PVfunctions::C12(const double m12, const double m22, const double m32) const
648{
649 if ( m12<=0.0 || m22<=0.0 || m32<=0.0 )
650 throw std::runtime_error("PVfunctions::C12(): Argument is not positive!");
651
652 const double Tolerance = 2.5e-3;
653 double C12;
654
655 if ( 2.0*fabs(m12-m22) > (m12+m22)*Tolerance && 2.0*fabs(m12-m32) > (m12+m32)*Tolerance && 2.0*fabs(m22-m32) > (m22+m32)*Tolerance ) {
656 C12=((m12-m22)*(m12-m32)*(m22-m32)*m32
657 +m12*m12*(m22*m22*log(m12/m22) +m32*(-2.0*m22+m32)*log(m12/m32))
658 +m22*m22*(2.0*m12-m32)*m32*log(m22/m32))
659 /(2.0*(m12-m22)*(m12-m32)*(m12-m32)*(m22-m32)*(m22-m32));
660 }
661 else if ( 2.0*fabs(m12-m22) > (m12+m22)*Tolerance && 2.0*fabs(m12-m32) > (m12+m32)*Tolerance && 2.0*fabs(m22-m32) <= (m22+m32)*Tolerance ) {
662 C12=-(-12.0*m12*m12+8.0*m12*(m22+m32)-(m22+m32)*(m22+m32)+8.0*m12*m12*log((2.0*m12)/(m22+m32)))
663 /(2.0*pow(-2.0*m12+m22+m32,3));
664 }
665 else if ( 2.0*fabs(m12-m22) > (m12+m22)*Tolerance && 2.0*fabs(m12-m32) <= (m12+m32)*Tolerance && 2.0*fabs(m22-m32) > (m22+m32)*Tolerance ) {
666 C12=(m12*m12-8.0*m12*m22+12.0*m22*m22+2.0*m12*m32-8.0*m22*m32+m32*m32-8.0*m22*m22*log((2.0*m22)/(m12+m32)))
667 /(2.0*pow(m12-2.0*m22+m32,3));
668 }
669 else if ( 2.0*fabs(m12-m22) <= (m12+m22)*Tolerance && 2.0*fabs(m12-m32) <= (m12+m32)*Tolerance && 2.0*fabs(m22-m32) <= (m22+m32)*Tolerance ) {
670 C12=1.0/(2.0*(m12+m22+m32));
671 }
672 else {
673 C12=(m12*m12+2.0*m12*m22+m22*m22-4.0*m32*m32-4.0*(m12+m22)*m32*log((m12+m22)/(2.0*m32)))
674 /pow(m12+m22-2.0*m32,3);
675 }
676
677 return C12;
678}
679
680gslpp::complex PVfunctions::D0(const double s, const double t, const double m02,
681 const double m12, const double m22, const double m32) const
682{
683 if ( m02<0.0 || m12<0.0 || m22<0.0 || m32<0.0 )
684 throw std::runtime_error("PVfunctions::D0(): Invalid argument!");
685
686 const double Tolerance = 1.0e-8;
687 const double maxM2 = std::max(s, std::max(t, std::max(m02, std::max(m12, std::max(m22, m32)))));
688 const bool szero = (s <= maxM2*Tolerance);
689 const bool tzero = (t <= maxM2*Tolerance);
690 const bool m02zero = (m02 <= maxM2*Tolerance);
691 const bool m12zero = (m12 <= maxM2*Tolerance);
692 const bool m22zero = (m22 <= maxM2*Tolerance);
693 const bool m32zero = (m32 <= maxM2*Tolerance);
694 bool diff01 = (fabs(m02 - m12) > (m02 + m12)*Tolerance);
695 bool diff02 = (fabs(m02 - m22) > (m02 + m22)*Tolerance);
696 bool diff03 = (fabs(m02 - m32) > (m02 + m32)*Tolerance);
697 bool diff23 = (fabs(m22 - m32) > (m22 + m32)*Tolerance);
698
699 if ( szero && tzero ) {
700 if ( diff01 )
701 return ( ( ExtraMinusSign * C0(0.0, m02, m22, m32)
702 - ExtraMinusSign * C0(0.0, m12, m22, m32) ) / (m02 - m12) );
703 else {
704 if ( !m02zero && !m12zero && !m22zero && !m32zero ) {
705 if ( diff02 && diff03 && diff23 )
706 return ( - 1.0/(m02 - m22)/(m02 - m32)
707 + m22/(m02 - m22)/(m02 - m22)/(m32 - m22)*log(m22/m02)
708 + m32/(m02 - m32)/(m02 - m32)/(m22 - m32)*log(m32/m02) );
709 else if ( !diff02 && diff03 && diff23 )
710 return ( (m02*m02 - m32*m32 + 2.0*m02*m32*log(m32/m02))
711 /2.0/m02/(m02 - m32)/(m02 - m32)/(m02 - m32) );
712 else if ( diff02 && !diff03 && diff23 )
713 return ( (m02*m02 - m22*m22 + 2.0*m02*m22*log(m22/m02))
714 /2.0/m02/(m02 - m22)/(m02 - m22)/(m02 - m22) );
715 else if ( diff02 && diff03 && !diff23 )
716 return ( ( - 2.0*m02 + 2.0*m22 - (m02 + m22)*log(m22/m02))
717 /(m02 - m22)/(m02 - m22)/(m02 - m22) );
718 else
719 return ( 1.0/6.0/m02/m02 );
720 } else
721 throw std::runtime_error("PVfunctions::D0(): Undefined!");
722 }
723 }
724
725#ifdef USE_LOOPTOOLS
726 return myLT.PV_D0(s, t, m02, m12, m22, m32);
727#else
728 gslpp::complex D0(0.0, 0.0, false);
729
730 if ( s>0.0 && t<0.0 && !m02zero && m12zero && !diff02 && !m32zero
731 && m02!=m32 && t-m32+m02!=0.0 && t-m32!=0.0 ) {
732 //D0(s,t; m02, 0.0, m02, m32)
733
734 /*
735 double x1, x2;
736 if ( s >= 4.0*m02 ) {
737 x1 = (1.0 - sqrt(1.0 - 4.0*m02/s))/2.0;
738 x2 = (1.0 + sqrt(1.0 - 4.0*m02/s))/2.0;
739 } else {
740 throw std::runtime_error("PVfunctions::D0(): Undefined!");
741 }
742 double x3 = m32/(m32 - m02);
743 double x4 = (t - m32)/(t - m32 + m02);
744 double d4 = 1.0 - 4.0*m02*t*(t - m32 + m02)/(s*(t - m32)*(t - m32));
745 double x1tilde, x2tilde;
746 if ( d4 >= 0.0 ) {
747 x1tilde = x4/2.0*(1.0 - sqrt(d4));
748 x2tilde = x4/2.0*(1.0 + sqrt(d4));
749 } else {
750 throw std::runtime_error("PVfunctions::D0(): Undefined!");
751 }
752 */
753
754 /* Write codes! */
755
756 throw std::runtime_error("PVfunctions::D0(): Undefined!");
757 } else if ( s>0.0 && t<0.0 && !m02zero && m12zero && !diff02 && m32zero ) {
758 //D0(s,t; m02, 0.0, m02, 0.0)
759
760 throw std::runtime_error("PVfunctions::D0(): Undefined!");
761 } else
762 throw std::runtime_error("PVfunctions::D0(): Undefined!");
763
764 return D0;
765#endif
766}
767
768gslpp::complex PVfunctions::D00(const double s, const double t, const double m02,
769 const double m12, const double m22, const double m32) const
770{
771 if ( m02<0.0 || m12<0.0 || m22<0.0 || m32<0.0 )
772 throw std::runtime_error("PVfunctions::D00(): Invalid argument!");
773
774 const double Tolerance = 1.0e-8;
775 const double maxM2 = std::max(s, std::max(t, std::max(m02, std::max(m12, std::max(m22, m32)))));
776 const bool szero = (s <= maxM2*Tolerance);
777 const bool tzero = (t <= maxM2*Tolerance);
778 const bool m02zero = (m02 <= maxM2*Tolerance);
779 const bool m12zero = (m12 <= maxM2*Tolerance);
780 const bool m22zero = (m22 <= maxM2*Tolerance);
781 const bool m32zero = (m32 <= maxM2*Tolerance);
782 bool diff01 = (fabs(m02 - m12) > (m02 + m12)*Tolerance);
783 bool diff02 = (fabs(m02 - m22) > (m02 + m22)*Tolerance);
784 bool diff03 = (fabs(m02 - m32) > (m02 + m32)*Tolerance);
785 bool diff23 = (fabs(m22 - m32) > (m22 + m32)*Tolerance);
786
787 if ( szero && tzero ) {
788 if ( diff01 )
789 return ( 0.25/(m02 - m12)*(m02*ExtraMinusSign*C0(0.0, m02, m22, m32)
790 - m12*ExtraMinusSign*C0(0.0, m12, m22, m32)) );
791 else {
792 if ( !m02zero && !m12zero && !m22zero && !m32zero ) {
793 if ( diff02 && diff03 && diff23 )
794 return ( m02/4.0/(m02 - m22)/(m32 - m02)
795 + m22*m22/4.0/(m02 - m22)/(m02 - m22)/(m32 - m22)*log(m22/m02)
796 + m32*m32/4.0/(m02 - m32)/(m02 - m32)/(m22 - m32)*log(m32/m02) );
797 else if ( !diff02 && diff03 && diff23 )
798 return ( ( - m02*m02 + 4.0*m02*m32 - 3.0*m32*m32
799 + 2.0*m32*m32*log(m32/m02) )
800 /8.0/(m02 - m32)/(m02 - m32)/(m02 - m32) );
801 else if ( diff02 && !diff03 && diff23 )
802 return ( ( - m02*m02 + 4.0*m02*m22 - 3.0*m22*m22
803 + 2.0*m22*m22*log(m22/m02) )
804 /8.0/(m02 - m22)/(m02 - m22)/(m02 - m22) );
805 else if ( diff02 && diff03 && !diff23 )
806 return ( ( - m02*m02 + m22*m22 - 2.0*m02*m22*log(m22/m02) )
807 /4.0/(m02 - m22)/(m02 - m22)/(m02 - m22) );
808 else
809 return ( - 1.0/12.0/m02 );
810 } else
811 throw std::runtime_error("PVfunctions::D00(): Undefined!");
812 }
813 }
814
815#ifdef USE_LOOPTOOLS
816 return myLT.PV_D00(s, t, m02, m12, m22, m32);
817#else
818 gslpp::complex D00(0.0, 0.0, false);
819
820 throw std::runtime_error("PVfunctions::D00(): Undefined!");
821
822 return D00;
823#endif
824}
825
F0
Test Observable.
Definition: NPSMEFT6dtopquark.h:785
PVfunctions::C12
double C12(const double m12, const double m22, const double m32) const
.
Definition: PVfunctions.cpp:647
PVfunctions::B1
gslpp::complex B1(const double mu2, const double p2, const double m02, const double m12) const
.
Definition: PVfunctions.cpp:120
PVfunctions::D00
gslpp::complex D00(const double s, const double t, const double m02, const double m12, const double m22, const double m32) const
.
Definition: PVfunctions.cpp:768
PVfunctions::myLT
LoopToolsWrapper myLT
An object of type LoopToolsWrapper.
Definition: PVfunctions.h:389
PVfunctions::B00
gslpp::complex B00(const double mu2, const double p2, const double m02, const double m12) const
.
Definition: PVfunctions.cpp:208
PVfunctions::B0
gslpp::complex B0(const double mu2, const double p2, const double m02, const double m12) const
.
Definition: PVfunctions.cpp:41
PVfunctions::B00p
gslpp::complex B00p(const double mu2, const double p2, const double m02, const double m12) const
.
Definition: PVfunctions.cpp:410
PVfunctions::C11
double C11(const double m12, const double m22, const double m32) const
.
Definition: PVfunctions.cpp:613
PVfunctions::PVfunctions
PVfunctions(const bool bExtraMinusSign)
Constructor.
Definition: PVfunctions.cpp:17
PVfunctions::B11
gslpp::complex B11(const double mu2, const double p2, const double m02, const double m12) const
.
Definition: PVfunctions.cpp:162
PVfunctions::B0p
gslpp::complex B0p(const double muIR2, const double p2, const double m02, const double m12) const
.
Definition: PVfunctions.cpp:276
PVfunctions::Bfp
gslpp::complex Bfp(const double mu2, const double p2, const double m02, const double m12) const
.
Definition: PVfunctions.cpp:458
PVfunctions::B11p
gslpp::complex B11p(const double mu2, const double p2, const double m02, const double m12) const
.
Definition: PVfunctions.cpp:381
PVfunctions::Bf
gslpp::complex Bf(const double mu2, const double p2, const double m02, const double m12) const
.
Definition: PVfunctions.cpp:265
PVfunctions::A0
double A0(const double mu2, const double m2) const
.
Definition: PVfunctions.cpp:23
PVfunctions::C0
gslpp::complex C0(const double p2, const double m02, const double m12, const double m22) const
.
Definition: PVfunctions.cpp:479
PVfunctions::ExtraMinusSign
double ExtraMinusSign
An overall factor for the one-point and three-point functions, initialized in PVfunctions().
Definition: PVfunctions.h:386
PVfunctions::B1p
gslpp::complex B1p(const double mu2, const double p2, const double m02, const double m12) const
.
Definition: PVfunctions.cpp:356
PVfunctions::D0
gslpp::complex D0(const double s, const double t, const double m02, const double m12, const double m22, const double m32) const
.
Definition: PVfunctions.cpp:680
s
Test Observable.
t
Test Observable.
Li2
cd Li2(cd x)
Definition: hpl.h:1011