BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
eemmg-lib/src/integral.cpp
Go to the documentation of this file.
1/*
2 * integral.h - scalar integrals wrappers
3 *
4 * this file is part of PJFry library
5 * Copyright 2011 Valery Yundin
6 */
7
8#include "integral.h"
9#include "minor.h"
10#include "cache.h"
11#include <stdio.h>
12
13static Initialize pjinit=Initialize();
14
15#ifdef USE_QCDLOOP
17{
18 printf("PJFRY init\n");
19 const double dbl_min=std::numeric_limits<double>::min();
20
21 F77_FUNC(qlinit,QLINIT)();
22
23 if (qlprec.xalogm < dbl_min) {
24 qlprec.xalogm=dbl_min;
25 qlprec.xalog2=sqrt(dbl_min);
26 #ifndef NDEBUG
27 printf("Set xalogm to normalized value %e\n", qlprec.xalogm);
28 printf("Set xalog2 to normalized value %e\n", qlprec.xalog2);
29 #endif
30 }
31 if (qlprec.xclogm < dbl_min) {
32 qlprec.xclogm=dbl_min;
33 qlprec.xclog2=sqrt(dbl_min);
34 #ifndef NDEBUG
35 printf("Set xclogm to normalized value %e\n", qlprec.xclogm);
36 printf("Set xclog2 to normalized value %e\n", qlprec.xclog2);
37 #endif
38 }
39 if (qlflag.lwarn) {
40 qlflag.lwarn=0;
41 #ifndef NDEBUG
42 printf("Disable FF warnings %d\n",qlflag.lwarn);
43 #endif
44 }
45}
46
48{
49 F77_FUNC(ffexi,FFEXI)();
50}
51
52ICache::Ival qlI1(const Kinem1& k)
53{
54 int ep;
55 ICache::Ival ivalue;
56
57 double m1=k.m1();
58 double mu2=ICache::getMu2();
59
60#ifdef USE_F2C
61 std::complex<double> rslt;
62 ep=0;
63 F77_FUNC(qli1,QLI1)(&rslt, &m1, &mu2, &ep);
64 ivalue.val[0]=rslt;
65 ep=-1;
66 F77_FUNC(qli1,QLI1)(&rslt, &m1, &mu2, &ep);
67 ivalue.val[1]=rslt;
68 ep=-2;
69 F77_FUNC(qli1,QLI1)(&rslt, &m1, &mu2, &ep);
70 ivalue.val[2]=rslt;
71#else
72 ep=0;
73 ivalue.val[0]=F77_FUNC(qli1,QLI1)(&m1, &mu2, &ep);
74 ep=-1;
75 ivalue.val[1]=F77_FUNC(qli1,QLI1)(&m1, &mu2, &ep);
76 ep=-2;
77 ivalue.val[2]=F77_FUNC(qli1,QLI1)(&m1, &mu2, &ep);
78#endif
79
80 return ivalue;
81}
82
83ICache::Ival qlI2(const Kinem2& k)
84{
85 int ep;
86 ICache::Ival ivalue;
87
88 double p1=k.p1();
89 double m1=k.m1();
90 double m2=k.m2();
91 double mu2=ICache::getMu2();
92
93#ifdef USE_F2C
94 std::complex<double> rslt;
95 ep=0;
96 F77_FUNC(qli2,QLI2)(&rslt, &p1, &m1, &m2, &mu2, &ep);
97 ivalue.val[0]=rslt;
98 ep=-1;
99 F77_FUNC(qli2,QLI2)(&rslt, &p1, &m1, &m2, &mu2, &ep);
100 ivalue.val[1]=rslt;
101 ep=-2;
102 F77_FUNC(qli2,QLI2)(&rslt, &p1, &m1, &m2, &mu2, &ep);
103 ivalue.val[2]=rslt;
104#else
105 ep=0;
106 ivalue.val[0]=F77_FUNC(qli2,QLI2)(&p1, &m1, &m2, &mu2, &ep);
107 ep=-1;
108 ivalue.val[1]=F77_FUNC(qli2,QLI2)(&p1, &m1, &m2, &mu2, &ep);
109 ep=-2;
110 ivalue.val[2]=F77_FUNC(qli2,QLI2)(&p1, &m1, &m2, &mu2, &ep);
111#endif
112
113 return ivalue;
114}
115
116ICache::Ival qlI3(const Kinem3& k)
117{
118 int ep;
119 ICache::Ival ivalue;
120
121 double p1=k.p1();
122 double p2=k.p2();
123 double p3=k.p3();
124 double m1=k.m1();
125 double m2=k.m2();
126 double m3=k.m3();
127 double mu2=ICache::getMu2();
128
129#ifdef USE_F2C
130 std::complex<double> rslt;
131 ep=0;
132 F77_FUNC(qli3,QLI3)(&rslt, &p1, &p2, &p3, &m1, &m2, &m3, &mu2, &ep);
133 ivalue.val[0]=rslt;
134 ep=-1;
135 F77_FUNC(qli3,QLI3)(&rslt, &p1, &p2, &p3, &m1, &m2, &m3, &mu2, &ep);
136 ivalue.val[1]=rslt;
137 ep=-2;
138 F77_FUNC(qli3,QLI3)(&rslt, &p1, &p2, &p3, &m1, &m2, &m3, &mu2, &ep);
139 ivalue.val[2]=rslt;
140#else
141 ep=0;
142 ivalue.val[0]=F77_FUNC(qli3,QLI3)(&p1, &p2, &p3, &m1, &m2, &m3, &mu2, &ep);
143 ep=-1;
144 ivalue.val[1]=F77_FUNC(qli3,QLI3)(&p1, &p2, &p3, &m1, &m2, &m3, &mu2, &ep);
145 ep=-2;
146 ivalue.val[2]=F77_FUNC(qli3,QLI3)(&p1, &p2, &p3, &m1, &m2, &m3, &mu2, &ep);
147#endif
148
149 return ivalue;
150}
151
152ICache::Ival qlI4(const Kinem4& k)
153{
154 int ep;
155 ICache::Ival ivalue;
156
157 double p1=k.p1();
158 double p2=k.p2();
159 double p3=k.p3();
160 double p4=k.p4();
161 double s12=k.s12();
162 double s23=k.s23();
163 double m1=k.m1();
164 double m2=k.m2();
165 double m3=k.m3();
166 double m4=k.m4();
167 double mu2=ICache::getMu2();
168
169 if (s12==0. || s23==0.) {
170 if (p1!=0. && p3!=0.) {
171# ifdef USE_F2C
172 std::complex<double> rslt;
173 ep=0;
174 F77_FUNC(qli4,QLI4)(&rslt, &s12, &p2, &s23, &p4, &p1, &p3, &m1, &m3, &m2, &m4, &mu2, &ep);
175 ivalue.val[0]=rslt;
176 ep=-1;
177 F77_FUNC(qli4,QLI4)(&rslt, &s12, &p2, &s23, &p4, &p1, &p3, &m1, &m3, &m2, &m4, &mu2, &ep);
178 ivalue.val[1]=rslt;
179 ep=-2;
180 F77_FUNC(qli4,QLI4)(&rslt, &s12, &p2, &s23, &p4, &p1, &p3, &m1, &m3, &m2, &m4, &mu2, &ep);
181 ivalue.val[2]=rslt;
182# else
183 ep=0;
184 ivalue.val[0]=F77_FUNC(qli4,QLI4)(&s12, &p2, &s23, &p4, &p1, &p3, &m1, &m3, &m2, &m4, &mu2, &ep);
185 ep=-1;
186 ivalue.val[1]=F77_FUNC(qli4,QLI4)(&s12, &p2, &s23, &p4, &p1, &p3, &m1, &m3, &m2, &m4, &mu2, &ep);
187 ep=-2;
188 ivalue.val[2]=F77_FUNC(qli4,QLI4)(&s12, &p2, &s23, &p4, &p1, &p3, &m1, &m3, &m2, &m4, &mu2, &ep);
189# endif
190 } else if (p2!=0. && p4!=0.) {
191# ifdef USE_F2C
192 std::complex<double> rslt;
193 ep=0;
194 F77_FUNC(qli4,QLI4)(&rslt, &s23, &p3, &s12, &p1, &p2, &p4, &m2, &m4, &m3, &m1, &mu2, &ep);
195 ivalue.val[0]=rslt;
196 ep=-1;
197 F77_FUNC(qli4,QLI4)(&rslt, &s23, &p3, &s12, &p1, &p2, &p4, &m2, &m4, &m3, &m1, &mu2, &ep);
198 ivalue.val[1]=rslt;
199 ep=-2;
200 F77_FUNC(qli4,QLI4)(&rslt, &s23, &p3, &s12, &p1, &p2, &p4, &m2, &m4, &m3, &m1, &mu2, &ep);
201 ivalue.val[2]=rslt;
202# else
203 ep=0;
204 ivalue.val[0]=F77_FUNC(qli4,QLI4)(&s23, &p3, &s12, &p1, &p2, &p4, &m2, &m4, &m3, &m1, &mu2, &ep);
205 ep=-1;
206 ivalue.val[1]=F77_FUNC(qli4,QLI4)(&s23, &p3, &s12, &p1, &p2, &p4, &m2, &m4, &m3, &m1, &mu2, &ep);
207 ep=-2;
208 ivalue.val[2]=F77_FUNC(qli4,QLI4)(&s23, &p3, &s12, &p1, &p2, &p4, &m2, &m4, &m3, &m1, &mu2, &ep);
209# endif
210 } else { assert(0); }
211 } else {
212# ifdef USE_F2C
213 std::complex<double> rslt;
214 ep=0;
215 F77_FUNC(qli4,QLI4)(&rslt, &p1, &p2, &p3, &p4, &s12, &s23, &m1, &m2, &m3, &m4, &mu2, &ep);
216 ivalue.val[0]=rslt;
217 ep=-1;
218 F77_FUNC(qli4,QLI4)(&rslt, &p1, &p2, &p3, &p4, &s12, &s23, &m1, &m2, &m3, &m4, &mu2, &ep);
219 ivalue.val[1]=rslt;
220 ep=-2;
221 F77_FUNC(qli4,QLI4)(&rslt, &p1, &p2, &p3, &p4, &s12, &s23, &m1, &m2, &m3, &m4, &mu2, &ep);
222 ivalue.val[2]=rslt;
223# else
224 ep=0;
225 ivalue.val[0]=F77_FUNC(qli4,QLI4)(&p1, &p2, &p3, &p4, &s12, &s23, &m1, &m2, &m3, &m4, &mu2, &ep);
226 ep=-1;
227 ivalue.val[1]=F77_FUNC(qli4,QLI4)(&p1, &p2, &p3, &p4, &s12, &s23, &m1, &m2, &m3, &m4, &mu2, &ep);
228 ep=-2;
229 ivalue.val[2]=F77_FUNC(qli4,QLI4)(&p1, &p2, &p3, &p4, &s12, &s23, &m1, &m2, &m3, &m4, &mu2, &ep);
230# endif
231 }
232 return ivalue;
233}
234
235#endif /* USE_QCDLOOP */
236
237#ifdef USE_ONELOOP
239{
240 printf("PJFRY init\n");
241
242 double mu=sqrt(ICache::getMu2());
243 F77_FUNC_(avh_olo_mu_set,AVH_OLO_MU_SET)(&mu);
244
245 double thrs=Minor5::getmeps();
246 F77_FUNC_(avh_olo_onshell,AVH_OLO_ONSHELL)(&thrs);
247}
248
249ICache::Ival qlI1(const Kinem1& k)
250{
251 ICache::Ival ivalue;
252
253 double m1=k.m1();
254 std::complex<double> rslt[3];
255 F77_FUNC_(avh_olo_a0m,AVH_OLO_A0M)(rslt, &m1);
256
257 ivalue.val[0]=rslt[0];
258 ivalue.val[1]=rslt[1];
259 ivalue.val[2]=rslt[2];
260 return ivalue;
261}
262
263ICache::Ival qlI2(const Kinem2& k)
264{
265 ICache::Ival ivalue;
266
267 double p1=k.p1();
268 double m1=k.m1();
269 double m2=k.m2();
270 std::complex<double> rslt[3];
271 F77_FUNC_(avh_olo_b0m,AVH_OLO_B0M)(rslt, &p1, &m1, &m2);
272
273 ivalue.val[0]=rslt[0];
274 ivalue.val[1]=rslt[1];
275 ivalue.val[2]=rslt[2];
276 return ivalue;
277}
278
279ICache::Ival qlI3(const Kinem3& k)
280{
281 ICache::Ival ivalue;
282
283 double p1=k.p1();
284 double p2=k.p2();
285 double p3=k.p3();
286 double m1=k.m1();
287 double m2=k.m2();
288 double m3=k.m3();
289 std::complex<double> rslt[3];
290 F77_FUNC_(avh_olo_c0m,AVH_OLO_C0M)(rslt, &p1, &p2, &p3, &m1, &m2, &m3);
291
292 ivalue.val[0]=rslt[0];
293 ivalue.val[1]=rslt[1];
294 ivalue.val[2]=rslt[2];
295 return ivalue;
296}
297
298ICache::Ival qlI4(const Kinem4& k)
299{
300 ICache::Ival ivalue;
301
302 double p1=k.p1();
303 double p2=k.p2();
304 double p3=k.p3();
305 double p4=k.p4();
306 double s12=k.s12();
307 double s23=k.s23();
308 double m1=k.m1();
309 double m2=k.m2();
310 double m3=k.m3();
311 double m4=k.m4();
312 std::complex<double> rslt[3];
313 if (s12==0. || s23==0.) {
314 if (p1!=0. && p3!=0.) {
315 F77_FUNC_(avh_olo_d0m,AVH_OLO_D0M)(rslt, &s12, &p2, &s23, &p4, &p1, &p3, &m1, &m3, &m2, &m4);
316 } else if (p2!=0. && p4!=0.) {
317 F77_FUNC_(avh_olo_d0m,AVH_OLO_D0M)(rslt, &s23, &p3, &s12, &p1, &p2, &p4, &m2, &m4, &m3, &m1);
318 } else { assert(0); }
319 } else {
320 F77_FUNC_(avh_olo_d0m,AVH_OLO_D0M)(rslt, &p1, &p2, &p3, &p4, &s12, &s23, &m1, &m2, &m3, &m4);
321 }
322
323 ivalue.val[0]=rslt[0];
324 ivalue.val[1]=rslt[1];
325 ivalue.val[2]=rslt[2];
326 return ivalue;
327}
328
329#endif /* USE_QCDLOOP1 */
double p1[4]
double p2[4]
static double getMu2()
double m1() const
double m2() const
double m1() const
double p1() const
double m1() const
double p1() const
double p3() const
double p2() const
double m2() const
double m3() const
double m1() const
double m3() const
double s12() const
double p1() const
double s23() const
double m2() const
double p3() const
double p2() const
double p4() const
double m4() const
static double getmeps()
ICache::Ival qlI1(const Kinem1 &k)
ICache::Ival qlI2(const Kinem2 &k)
ICache::Ival qlI4(const Kinem4 &k)
ICache::Ival qlI3(const Kinem3 &k)
double double double double double * m3
Definition qcdloop1.h:76
double double double * p4
Definition qcdloop1.h:77
double double double double * s12
Definition qcdloop1.h:77
double double double double double * s23
Definition qcdloop1.h:77
double double double double double double double double double * m4
Definition qcdloop1.h:77
double int * ep
Definition qcdloop1.h:74
double * mu2
Definition qcdloop1.h:74
double * m1
Definition qcdloop1.h:75
double double * p3
Definition qcdloop1.h:76
double double * m2
Definition qcdloop1.h:75