BOSS 7.0.8
BESIII Offline Software System
Loading...
Searching...
No Matches
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
16Initialize::Initialize()
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
47Initialize::~Initialize()
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
238Initialize::Initialize()
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 */
ICache::getMu2
static double getMu2()
Definition: cache.cpp:20
Kinem1
Definition: kinem.h:93
Kinem1::m1
double m1() const
Definition: kinem.h:100
Kinem2
Definition: kinem.h:105
Kinem2::m2
double m2() const
Definition: kinem.h:118
Kinem2::m1
double m1() const
Definition: kinem.h:117
Kinem2::p1
double p1() const
Definition: kinem.h:116
Kinem3
Definition: kinem.h:123
Kinem3::m1
double m1() const
Definition: kinem.h:140
Kinem3::p1
double p1() const
Definition: kinem.h:137
Kinem3::p3
double p3() const
Definition: kinem.h:139
Kinem3::p2
double p2() const
Definition: kinem.h:138
Kinem3::m2
double m2() const
Definition: kinem.h:141
Kinem3::m3
double m3() const
Definition: kinem.h:142
Kinem4
Definition: kinem.h:147
Kinem4::m1
double m1() const
Definition: kinem.h:170
Kinem4::m3
double m3() const
Definition: kinem.h:172
Kinem4::s12
double s12() const
Definition: kinem.h:174
Kinem4::p1
double p1() const
Definition: kinem.h:166
Kinem4::s23
double s23() const
Definition: kinem.h:175
Kinem4::m2
double m2() const
Definition: kinem.h:171
Kinem4::p3
double p3() const
Definition: kinem.h:168
Kinem4::p2
double p2() const
Definition: kinem.h:167
Kinem4::p4
double p4() const
Definition: kinem.h:169
Kinem4::m4
double m4() const
Definition: kinem.h:173
MinorBase::getmeps
static double getmeps()
Definition: minor.h:148
qlI1
ICache::Ival qlI1(const Kinem1 &k)
qlI2
ICache::Ival qlI2(const Kinem2 &k)
qlI4
ICache::Ival qlI4(const Kinem4 &k)
qlI3
ICache::Ival qlI3(const Kinem3 &k)
ICache::Ival::val
ncomplex val[3]
Definition: cache.h:57