BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKpipipi.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of models developed at BES collaboration
4// based on the EvtGen framework. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/BesCopyright
8// Copyright (A) 2006 Ping Rong-Gang @IHEP
9//
10// Module: EvtD0ToKpipipi.cc
11// the necessary file: EvtD0ToKpipipi.hh
12//
13// Description: D0 -> K- pi+ pi+ pi-,
14// see PHYSICAL REVIEW D 95, 072010 (2017)
15//
16// Modification history:
17//
18// Liaoyuan Dong Jan.15, 2020 Module created
19//
20//------------------------------------------------------------------------
21#include "EvtGenBase/EvtPatches.hh"
22#include "EvtGenBase/EvtParticle.hh"
23#include "EvtGenBase/EvtGenKine.hh"
24#include "EvtGenBase/EvtPDL.hh"
25#include "EvtGenBase/EvtReport.hh"
26#include "EvtGenModels/EvtD0ToKpipipi.hh"
27#include "EvtGenBase/EvtComplex.hh"
28#include "EvtGenBase/EvtDecayTable.hh"
29#include <stdlib.h>
30
31EvtD0ToKpipipi::~EvtD0ToKpipipi() {}
32
33void EvtD0ToKpipipi::getName(std::string& model_name){
34 model_name="D0ToKpipipi";
35}
36
37EvtDecayBase* EvtD0ToKpipipi::clone(){
38 return new EvtD0ToKpipipi;
39}
40
41void EvtD0ToKpipipi::init(){
42 checkNArg(0);
43 checkNDaug(4);
44 checkSpinParent(EvtSpinType::SCALAR);
45/*
46 checkSpinDaughter(0,EvtSpinType::SCALAR);
47 checkSpinDaughter(1,EvtSpinType::SCALAR);
48 checkSpinDaughter(3,EvtSpinType::SCALAR);
49 checkSpinDaughter(4,EvtSpinType::SCALAR);
50*/
51 //std::cout << "Initializing EvtD0ToKpipipi" << std::endl;
52
53 width[0] = 0.09;
54 width[1] = 0.044183653178315;
55 width[2] = 0.541879469380012;
56 width[3] = 0.148423336450619;
57 mass[0] = 1.272;
58 mass[1] = 0.894781734682169;
59 mass[2] = 1.3622013558915;
60 mass[3] = 0.779143408171384;
61
62 phi[0] = 2.34794687054858;
63 rho[0] = 0.0759345115620669;
64 phi[1] = -2.24641399153466;
65 rho[1] = 0.0383327604903577;
66 phi[2] = 2.48955684856045;
67 rho[2] = 0.0931445480476023;
68
69 phi[3] = 0;
70 rho[3] = 1;
71
72 phi[4] = -2.10558220063012;
73 rho[4] = 0.347041869435286;
74 phi[5] = 1.47445088061872;
75 phi[6] = 3.00243265559304;
76 rho[5] = 0.00965088341753795;
77 rho[6] = 0.120536507325731;
78 phi[7] = -2.45477499325158;
79 rho[7] = 0.101419048440676;
80 phi[8] = -1.35809992343491;
81 rho[8] = 4.28149643321317;
82 phi[9] = -2.45149221243198;
83 rho[9] = 0.339492272598394;
84 phi[10] = -0.17419389225461;
85 rho[10] = -0.143619437541254;
86
87 phi[11] = -2.08744386934208;
88 rho[11] = 0.296286583716349;
89 phi[12] = 0.;
90 rho[12] = 0.;
91 phi[13] = -0.432190571560873;
92 rho[13] = 0.657344690733276;
93 phi[14] = -1.39790294886865;
94 rho[14] = 1.71208007006123;
95 phi[15] = 1.58945300476228;
96 rho[15] = 3.58248347683687;
97 phi[16] = 2.58249107256307;
98 rho[16] = -1.10728829503506;
99 phi[17] = -0.163623135170955;
100 rho[17] = 1.70863070178363;
101 phi[18] = -0.134699023080211;
102 rho[18] = 0.567531283682344;
103 phi[19] = -2.12670610368279;
104 rho[19] = 0.276571752504914;
105 phi[20] = -1.3352622107357;
106 rho[20] = 0.416634203151278;
107
108 phi[21] = -2.91571684221842;
109 rho[21] = 0.423062298489176;
110 phi[22] = 2.4544220004327;
111 rho[22] = 1.4017194038459;
112 phi[23] = -2.23388390670423;
113 rho[23] = 4.11110400629068;
114
115 //for (int i=0; i<24; i++) {
116 // cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
117 //}
118
119 mD = 1.86486;
120 rRes = 3.0;
121 rD = 5.0;
122 metap = 0.95778;
123 mkstr = 0.89594;
124 mk0 = 0.497614;
125 mass_Kaon = 0.49368;
126 mass_Pion = 0.13957;
127 mass_Pi0 = 0.1349766;
128 math_pi = 3.1415926;
129
130 pi = 3.1415926;
131 mpi = 0.13957;
132 g1 = 0.5468;
133 g2 = 0.23;
134
135 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
136 int EE[4][4][4][4] =
137 { { {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
138 {{0,0,0,0}, {0,0,0,0}, {0,0,0,1}, {0,0,-1,0}},
139 {{0,0,0,0}, {0,0,0,-1}, {0,0,0,0}, {0,1,0,0} },
140 {{0,0,0,0}, {0,0,1,0}, {0,-1,0,0}, {0,0,0,0} } },
141 { {{0,0,0,0}, {0,0,0,0}, {0,0,0,-1}, {0,0,1,0} },
142 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
143 {{0,0,0,1}, {0,0,0,0}, {0,0,0,0}, {-1,0,0,0}},
144 {{0,0,-1,0}, {0,0,0,0}, {1,0,0,0}, {0,0,0,0} } },
145 { {{0,0,0,0}, {0,0,0,1}, {0,0,0,0}, {0,-1,0,0}},
146 {{0,0,0,-1}, {0,0,0,0}, {0,0,0,0}, {1,0,0,0} },
147 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
148 {{0,1,0,0}, {-1,0,0,0}, {0,0,0,0}, {0,0,0,0} } },
149 { {{0,0,0,0}, {0,0,-1,0}, {0,1,0,0}, {0,0,0,0} },
150 {{0,0,1,0}, {0,0,0,0}, {-1,0,0,0}, {0,0,0,0} },
151 {{0,-1,0,0}, {1,0,0,0}, {0,0,0,0}, {0,0,0,0} },
152 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} } } };
153 for (int i=0; i<4; i++) {
154 for (int j=0; j<4; j++) {
155 G[i][j] = GG[i][j];
156 for (int k=0; k<4; k++) {
157 for (int l=0; l<4; l++) {
158 E[i][j][k][l] = EE[i][j][k][l];
159 }
160 }
161 }
162 }
163}
164
165void EvtD0ToKpipipi::initProbMax() {
166 setProbMax(720.0);
167}
168
169void EvtD0ToKpipipi::decay( EvtParticle *p){
170/*
171 double maxprob = 0.0;
172 for(int ir=0;ir<=60000000;ir++){
173 p->initializePhaseSpace(getNDaug(),getDaugs());
174 EvtVector4R Km0 = p->getDaug(0)->getP4();
175 EvtVector4R pi1 = p->getDaug(1)->getP4();
176 EvtVector4R pi2 = p->getDaug(2)->getP4();
177 EvtVector4R pi3 = p->getDaug(3)->getP4();
178 double Km[4],Pip1[4],Pip2[4],Pim[4];
179 Km[0] = Km0.get(0); Pip1[0] = pi1.get(0); Pip2[0] = pi2.get(0); Pim[0] = pi3.get(0);
180 Km[1] = Km0.get(1); Pip1[1] = pi1.get(1); Pip2[1] = pi2.get(1); Pim[1] = pi3.get(1);
181 Km[2] = Km0.get(2); Pip1[2] = pi1.get(2); Pip2[2] = pi2.get(2); Pim[2] = pi3.get(2);
182 Km[3] = Km0.get(3); Pip1[3] = pi1.get(3); Pip2[3] = pi2.get(3); Pim[3] = pi3.get(3);
183 double Prob = calPDF(Km, Pip1, Pip2, Pim);
184 if(Prob>maxprob) {
185 maxprob=Prob;
186 std::cout << "Max PDF = " << ir << " prob= " << Prob << std::endl;
187 }
188 }
189 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
190*/
191 p->initializePhaseSpace(getNDaug(),getDaugs());
192 EvtVector4R Km0 = p->getDaug(0)->getP4();
193 EvtVector4R pi1 = p->getDaug(1)->getP4();
194 EvtVector4R pi2 = p->getDaug(2)->getP4();
195 EvtVector4R pi3 = p->getDaug(3)->getP4();
196
197 double Km[4],Pip1[4],Pip2[4],Pim[4];
198 Km[0] = Km0.get(0); Pip1[0] = pi1.get(0); Pip2[0] = pi2.get(0); Pim[0] = pi3.get(0);
199 Km[1] = Km0.get(1); Pip1[1] = pi1.get(1); Pip2[1] = pi2.get(1); Pim[1] = pi3.get(1);
200 Km[2] = Km0.get(2); Pip1[2] = pi1.get(2); Pip2[2] = pi2.get(2); Pim[2] = pi3.get(2);
201 Km[3] = Km0.get(3); Pip1[3] = pi1.get(3); Pip2[3] = pi2.get(3); Pim[3] = pi3.get(3);
202 double prob = calPDF(Km, Pip1, Pip2, Pim);
203 setProb(prob);
204 return;
205}
206
207double EvtD0ToKpipipi::calPDF(double Km[], double Pip1[], double Pip2[], double Pim[]) {
208 Km[0] = sqrt(mass_Kaon*mass_Kaon + Km[1]*Km[1] + Km[2]*Km[2] + Km[3]*Km[3]);
209 Pip1[0] = sqrt(mass_Pion*mass_Pion + Pip1[1]*Pip1[1] + Pip1[2]*Pip1[2] + Pip1[3]*Pip1[3]);
210 Pip2[0] = sqrt(mass_Pion*mass_Pion + Pip2[1]*Pip2[1] + Pip2[2]*Pip2[2] + Pip2[3]*Pip2[3]);
211 Pim[0] = sqrt(mass_Pion*mass_Pion + Pim[1]*Pim[1] + Pim[2]*Pim[2] + Pim[3]*Pim[3]);
212
213 EvtComplex PDF[24];
214 int g[3];
215 g[0] = 1; g[1] = 1;
216
217 //-----------D->K*rho--------
218 g[2] = 0;
219 PDF[0] = D2VV(Km,Pip1,Pip2,Pim,g) + D2VV(Km,Pip2,Pip1,Pim,g);
220 g[2] = 1;
221 PDF[1] = D2VV(Km,Pip1,Pip2,Pim,g) + D2VV(Km,Pip2,Pip1,Pim,g);
222 g[2] = 2;
223 PDF[2] = D2VV(Km,Pip1,Pip2,Pim,g) + D2VV(Km,Pip2,Pip1,Pim,g);
224 //-----------D->K*rho finish--
225 //----------D->a1K------------
226 g[2] = 0;
227 PDF[3] = D2AP_A2VP(Km,Pip2,Pip1,Pim,g,2) + D2AP_A2VP(Km,Pip1,Pip2,Pim,g,2);
228 g[2] = 2;
229 PDF[4] = D2AP_A2VP(Km,Pip2,Pip1,Pim,g,2) + D2AP_A2VP(Km,Pip1,Pip2,Pim,g,2);
230 //----------D->a1K finish-----
231 //--D->K1Pi--K1->K*Pi---------
232 g[2] = 0;
233 PDF[5] = D2AP_A2VP(Pip2,Pim,Km,Pip1,g,1) + D2AP_A2VP(Pip1,Pim,Km,Pip2,g,1);
234 g[2] = 2;
235 PDF[6] = D2AP_A2VP(Pip2,Pim,Km,Pip1,g,1) + D2AP_A2VP(Pip1,Pim,Km,Pip2,g,1);
236 //--D->K1Pi--K1->K*Pi-finish--
237 //--D->K1Pi--K1->rhoK---------
238 g[2] = 0;
239 PDF[7] = D2AP_A2VP(Pip2,Km,Pip1,Pim,g,3) + D2AP_A2VP(Pip1,Km,Pip2,Pim,g,3);
240 //--D->K1Pi--K1->rhoK-finish--
241 //--D->K1Pi--K1->K*0(1430)Pi--
242 PDF[8] = D2AP_A2SP(Pip2,Pim,Km,Pip1,1) + D2AP_A2SP(Pip1,Pim,Km,Pip2,1);
243 //--------res finish----------
244 //----------------non-res------------------
245 //--------KPiSrho------------------
246 PDF[9] = D2VS(Pip1,Pim,Km,Pip2,1,2) + D2VS(Pip2,Pim,Km,Pip1,1,2);
247 //--------K*PiPiS-----------------
248 PDF[10] = D2VS(Km,Pip1,Pip2,Pim,1,1) + D2VS(Km,Pip2,Pip1,Pim,1,1);
249 //--------K*PiP-------------------
250 PDF[11] = D2PP_P2VP(Pip2,Pim,Km,Pip1,1) + D2PP_P2VP(Pip1,Pim,Km,Pip2,1);
251 //--------rhoKA--------------------
252 g[0] = 1; g[1] = 0;
253 g[2] = 0;
254 PDF[12] = 0;
255 g[2] = 2;
256 PDF[13] = D2AP_A2VP(Pip2,Km,Pip1,Pim,g,3) + D2AP_A2VP(Pip1,Km,Pip2,Pim,g,3);
257 //-------PHSP-----------------------
258 PDF[14] = PHSP(Km,Pip1) + PHSP(Km,Pip2);
259 //------KPiVPiPiV-----------------------
260 g[0] = 0; g[1] = 0; g[2] = 0;
261 PDF[15] = D2VV(Km,Pip1,Pip2,Pim,g) + D2VV(Km,Pip2,Pip1,Pim,g);
262 //------KPiVPiPiS----------------------
263 PDF[16] = D2VS(Km,Pip1,Pip2,Pim,0,1) + D2VS(Km,Pip2,Pip1,Pim,0,1);
264 //------KPiSPiPiV----------------------
265 PDF[17] = D2VS(Pip1,Pim,Km,Pip2,0,2) + D2VS(Pip2,Pim,Km,Pip1,0,2);
266 //-------------------------------------
267 PDF[18] = D2PP_P2VP(Pip2,Km,Pip1,Pim,2) + D2PP_P2VP(Pip1,Km,Pip2,Pim,2);
268 //-------------------------------------
269 PDF[19] = D2VP_V2VP(Pip2,Pim,Km,Pip1,1) + D2VP_V2VP(Pip1,Pim,Km,Pip2,1);
270 //-------------------------------------
271 PDF[20] = D2VP_V2VP(Pip2,Km,Pip1,Pim,2) + D2VP_V2VP(Pip1,Km,Pip2,Pim,2);
272 //-------------------------------------
273 PDF[21] = D2TS(Km,Pip1,Pip2,Pim,1) + D2TS(Km,Pip2,Pip1,Pim,1);
274 PDF[22] = D2TS(Pip1,Pim,Km,Pip2,2) + D2TS(Pip2,Pim,Km,Pip1,2);
275 PDF[23] = D2AP_A2SP(Km,Pip2,Pip1,Pim,2) + D2AP_A2SP(Km,Pip1,Pip2,Pim,2);
276//------------------------------------------
277 EvtComplex cof;
278 EvtComplex pdf, module;
279 pdf = EvtComplex(0,0);
280 for(int i=0; i!=24; i++){
281 cof = EvtComplex(rho[i]*cos(phi[i]),rho[i]*sin(phi[i]));
282 //cout<<i << " " << (cof*PDF[i])<<" : "<<cof<<" "<<PDF[i]<<endl;
283 pdf = pdf + cof*PDF[i];
284 }
285 module = conj(pdf)*pdf;
286 double value;
287 value = real(module);
288 return (value <= 0) ? 1e-20 : value;
289}
290
291EvtComplex EvtD0ToKpipipi::KPiSFormfactor(double sa, double sb, double sc, double r)
292{
293 double m1430 = 1.463;
294 double sa0 = m1430*m1430;
295 double w1430 = 0.233;
296 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
297 if(q0<0) q0 = 1e-16;
298 double qs = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
299 double q = sqrt(qs);
300 double width = w1430*q*m1430/sqrt(sa*q0);
301 double temp_R = atan(m1430*width/(sa0-sa));
302 if(temp_R<0) temp_R += math_pi;
303 double deltaR = -5.31 + temp_R;
304 double temp_F = atan(2*1.07*q/(2+(-1.8)*1.07*qs));
305 if(temp_F<0) temp_F += math_pi;
306 double deltaF = 2.33 + temp_F;
307 EvtComplex cR(cos(deltaR), sin(deltaR));
308 EvtComplex cF(cos(deltaF), sin(deltaF));
309 EvtComplex amp = 0.8*sin(deltaF)*cF + sin(deltaR)*cR*cF*cF;
310 return amp;
311}
312
313EvtComplex EvtD0ToKpipipi::D2VV(double P1[], double P2[], double P3[], double P4[], int g[])
314{
315 double t1V1[4], t1V2[4], t1D[4], t2D[4][4];
316 double temp_PDF = 0;
317 EvtComplex amp_PDF(0,0);
318 EvtComplex pro[2];
319 //---------use g[0],g[1] to define res or non-res, g[2] represent S, P or D wave
320 double sa[3], sb[3], sc[3], B[3];
321 double pV1[4], pV2[4], pD[4];
322 for(int i=0; i!=4; i++){
323 pV1[i] = P1[i] + P2[i];
324 pV2[i] = P3[i] + P4[i];
325 pD[i] = pV1[i] + pV2[i];
326 }
327 sa[0] = dot(pV1,pV1);
328 sb[0] = dot(P1,P1);
329 sc[0] = dot(P2,P2);
330 sa[1] = dot(pV2,pV2);
331 sb[1] = dot(P3,P3);
332 sc[1] = dot(P4,P4);
333 sa[2] = dot(pD,pD);
334 sb[2] = sa[0];
335 sc[2] = sa[1];
336 if(g[1] == 1){
337 pro[1] = propagatorGS(mass[3],width[3],sa[1],sb[1],sc[1],rRes,1);
338 }
339 if(g[0] == 1){
340 pro[0] = propagatorRBW(mass[1],width[1],sa[0],sb[0],sc[0],rRes,1);
341 }
342 if(g[0] == 0) pro[0] = 1;
343 if(g[1] == 0) pro[1] = 1;
344 B[0] = barrier(1,sa[0],sb[0],sc[0],rRes);
345 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
346 calt1(P1,P2,t1V1);
347 calt1(P3,P4,t1V2);
348 if(g[2] == 0){
349 for(int i=0; i!=4; i++){
350 temp_PDF += (G[i][i])*t1V1[i]*t1V2[i];
351 }
352 B[2] = 1;
353 }
354 if(g[2] == 1){
355 calt1(pV1,pV2,t1D);
356 for(int i=0; i!=4; i++){
357 for(int j=0; j!=4; j++){
358 for(int k=0; k!=4; k++){
359 for(int l=0; l!=4; l++){
360 temp_PDF += E[i][j][k][l]*pD[i]*t1D[j]*t1V1[k]*t1V2[l]*(G[i][i])*(G[j][j])*(G[l][l])*(G[k][k]);
361 }
362 }
363 }
364 }
365 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
366 }
367 if(g[2] == 2){
368 calt2(pV1,pV2,t2D);
369 for(int i=0; i!=4; i++){
370 for(int j=0; j!=4; j++){
371 temp_PDF += t2D[i][j]*t1V1[i]*t1V2[j]*(G[i][i])*(G[j][j]);
372 }
373 }
374 B[2] = barrier(2,sa[2],sb[2],sc[2],rD);
375 }
376 amp_PDF = temp_PDF*B[0]*B[1]*B[2]*pro[0]*pro[1];
377 return amp_PDF;
378}
379
380EvtComplex EvtD0ToKpipipi::D2AP_A2VP(double P1[], double P2[], double P3[], double P4[], int g[], int flag)
381{
382 //flag = 1, V = K*, A = K1(1270); flag = 2, V = rho, A = a1(1260);
383 //flag = 3, V = rho, A = K1(1270);
384 double temp_PDF = 0;
385 EvtComplex amp_PDF(0,0);
386 EvtComplex pro[2];
387 double t1V[4], t1D[4], t2A[4][4];
388 double sa[3], sb[3], sc[3], B[3];
389 double pV[4],pA[4],pD[4];
390 for(int i=0; i!=4; i++){
391 pV[i] = P3[i] + P4[i];
392 pA[i] = pV[i] + P2[i];
393 pD[i] = pA[i] + P1[i];
394 }
395 sa[0] = dot(pV,pV);
396 sb[0] = dot(P3,P3);
397 sc[0] = dot(P4,P4);
398 sa[1] = dot(pA,pA);
399 sb[1] = sa[0];
400 sc[1] = dot(P2,P2);
401 sa[2] = dot(pD,pD);
402 sb[2] = sa[1];
403 sc[2] = dot(P1,P1);
404 if(g[0] == 1){
405 if(flag == 1) pro[0] = propagatorRBW(mass[1],width[1],sa[0],sb[0],sc[0],rRes,1);
406 if(flag == 2 || flag == 3) pro[0] = propagatorGS(mass[3],width[3],sa[0],sb[0],sc[0],rRes,1);
407 }
408 if(g[1] == 1){
409 if(flag == 1) pro[1] = propogator(mass[0],width[0],sa[1]);
410 if(flag == 2) pro[1] = propagatorRBW(mass[2],width[2],sa[1],sb[1],sc[1],rRes,g[2]);
411 if(flag == 3) pro[1] = propogator(mass[0],width[0],sa[1]);
412 }
413 if(g[0] == 0) pro[0] = 1;
414 if(g[1] == 0) pro[1] = 1;
415 B[0] = barrier(1,sa[0],sb[0],sc[0],rRes);
416 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
417 calt1(P3,P4,t1V);
418 calt1(pA,P1,t1D);
419 if(g[2] == 0){
420 for(int i=0; i!=4; i++){
421 for(int j=0; j!=4; j++){
422 temp_PDF += t1D[i]*(G[i][j] - pA[i]*pA[j]/sa[1])*t1V[j]*(G[i][i])*(G[j][j]);
423 }
424 }
425 B[1] = 1;
426 }
427 if(g[2] == 2){
428 calt2(pV,P2,t2A);
429 for(int i=0; i!=4; i++){
430 for(int j=0; j!=4; j++){
431 temp_PDF += t1D[i]*t2A[i][j]*t1V[j]*(G[i][i])*(G[j][j]);
432 }
433 }
434 B[1] = barrier(2,sa[1],sb[1],sc[1],rRes);
435 }
436 amp_PDF = temp_PDF*B[0]*B[1]*B[2]*pro[0]*pro[1];
437 return amp_PDF;
438}
439EvtComplex EvtD0ToKpipipi::D2AP_A2SP(double P1[], double P2[], double P3[], double P4[], int flag)
440{
441 //flag = 1, S = K*; flag = 2, S = rho
442 double temp_PDF = 0;
443 EvtComplex amp_PDF(0,0);
444 EvtComplex pro;
445 double sa[3], sb[3], sc[3], B[3];
446 double t1D[4], t1A[4];
447 double pS[4], pA[4], pD[4];
448 for(int i=0; i!=4; i++){
449 pS[i] = P3[i] + P4[i];
450 pA[i] = pS[i] + P2[i];
451 pD[i] = pA[i] + P1[i];
452 }
453 sa[0] = dot(pS,pS);
454 sb[0] = dot(P3,P3);
455 sc[0] = dot(P4,P4);
456 sa[1] = dot(pA,pA);
457 sb[1] = sa[0];
458 sc[1] = dot(P2,P2);
459 sa[2] = dot(pD,pD);
460 sb[2] = sa[1];
461 sc[2] = dot(P1,P1);
462 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
463 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
464 calt1(pA,P1,t1D);
465 calt1(pS,P2,t1A);
466 for(int i=0; i!=4; i++){
467 temp_PDF += t1D[i]*t1A[i]*(G[i][i]);
468 }
469 amp_PDF = temp_PDF*B[1]*B[2];
470 if(flag == 1) amp_PDF *= KPiSFormfactor(sa[0],sb[0],sc[0],rRes);
471 return amp_PDF;
472}
473EvtComplex EvtD0ToKpipipi::D2PP_P2VP(double P1[], double P2[], double P3[], double P4[], int flag)
474{
475 //flag = 1, V = K*; flag = 2, V = rho
476 double temp_PDF = 0;
477 EvtComplex amp_PDF(0,0);
478 EvtComplex pro;
479 double sa[3], sb[3], sc[3], B[3];
480 double t1V[4];
481 double pV[4], pP[4], pD[4];
482 for(int i=0; i!=4; i++){
483 pV[i] = P3[i] + P4[i];
484 pP[i] = pV[i] + P2[i];
485 pD[i] = pP[i] + P1[i];
486 }
487 sa[0] = dot(pV,pV);
488 sb[0] = dot(P3,P3);
489 sc[0] = dot(P4,P4);
490 sa[1] = dot(pP,pP);
491 sb[1] = sa[0];
492 sc[1] = dot(P2,P2);
493 sa[2] = dot(pD,pD);
494 sb[2] = sa[1];
495 sc[2] = dot(P1,P1);
496 B[0] = barrier(1,sa[0],sb[0],sc[0],rRes);
497 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
498 if(flag == 1) pro = propagatorRBW(mass[1],width[1],sa[0],sb[0],sc[0],rRes,1);
499 if(flag == 2) pro = propagatorGS(mass[3],width[3],sa[0],sb[0],sc[0],rRes,1);
500 calt1(P3,P4,t1V);
501 for(int i=0; i!=4; i++){
502 temp_PDF += P2[i]*t1V[i]*(G[i][i]);
503 }
504 amp_PDF = temp_PDF*B[0]*B[1]*pro;
505 return amp_PDF;
506}
507EvtComplex EvtD0ToKpipipi::D2VP_V2VP(double P1[], double P2[], double P3[], double P4[], int flag)
508{
509 //flag = 1, (K*Pi)V; flag = 2, (rhoK)V
510 double temp_PDF = 0;
511 EvtComplex amp_PDF(0,0);
512 EvtComplex pro;
513 double sa[3], sb[3], sc[3], B[3];
514 double pV1[4], pV2[4], qV1[4], qV2[4], pD[4];
515 for(int i=0; i!=4; i++){
516 pV2[i] = P3[i] + P4[i];
517 qV2[i] = P3[i] - P4[i];
518 pV1[i] = pV2[i] + P2[i];
519 qV1[i] = pV2[i] - P2[i];
520 pD[i] = pV1[i] + P1[i];
521 }
522 for(int i=0; i!=4; i++){
523 for(int j=0; j!=4; j++){
524 for(int k=0; k!=4; k++){
525 for(int l=0; l!=4; l++){
526 temp_PDF += E[i][j][k][l]*pV1[i]*qV1[j]*P1[k]*qV2[l]*(G[i][i])*(G[j][j])*(G[k][k])*(G[l][l]);
527 }
528 }
529 }
530 }
531 sa[0] = dot(pV2,pV2);
532 sb[0] = dot(P3,P3);
533 sc[0] = dot(P4,P4);
534 sa[1] = dot(pV1,pV1);
535 sb[1] = sa[0];
536 sc[1] = dot(P2,P2);
537 sa[2] = dot(pD,pD);
538 sb[2] = sa[1];
539 sc[2] = dot(P1,P1);
540 if(flag == 1) pro = propagatorRBW(mass[1],width[1],sa[0],sb[0],sc[0],rRes,1);
541 if(flag == 2) pro = propagatorGS(mass[3],width[3],sa[0],sb[0],sc[0],rRes,1);
542 B[0] = barrier(1,sa[0],sb[0],sc[0],rRes);
543 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
544 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
545 amp_PDF = temp_PDF*B[0]*B[1]*B[2]*pro;
546 return amp_PDF;
547}
548EvtComplex EvtD0ToKpipipi::D2VS(double P1[], double P2[], double P3[], double P4[], int g, int flag)
549{
550 //flag = 1, V = K*; flag = 2, V = rho
551 double temp_PDF = 0;
552 EvtComplex amp_PDF(0,0);
553 EvtComplex pro;
554 double sa[3], sb[3], sc[3], B[3];
555 double t1D[4], t1V[4];
556 double pS[4], pV[4], pD[4];
557 for(int i=0; i!=4; i++){
558 pS[i] = P3[i] + P4[i];
559 pV[i] = P1[i] + P2[i];
560 pD[i] = pS[i] + pV[i];
561 }
562 sa[0] = dot(pS,pS);
563 sb[0] = dot(P3,P3);
564 sc[0] = dot(P4,P4);
565 sa[1] = dot(pV,pV);
566 sb[1] = dot(P1,P1);
567 sc[1] = dot(P2,P2);
568 sa[2] = dot(pD,pD);
569 sb[2] = sa[0];
570 sc[2] = sa[1];
571 if(g == 1){
572 if(flag == 2) pro = propagatorGS(mass[3],width[3],sa[1],sb[1],sc[1],rRes,1);
573 if(flag == 1) pro = propagatorRBW(mass[1],width[1],sa[1],sb[1],sc[1],rRes,1);
574 }
575 if(g == 0) pro = 1;
576 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
577 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
578 calt1(P1,P2,t1V);
579 calt1(pS,pV,t1D);
580 for(int i=0; i!=4; i++){
581 temp_PDF += G[i][i]*t1D[i]*t1V[i];
582 }
583 amp_PDF = temp_PDF*B[1]*B[2]*pro;
584 if(flag == 2) amp_PDF *= KPiSFormfactor(sa[0],sb[0],sc[0],rRes);
585 return amp_PDF;
586}
587EvtComplex EvtD0ToKpipipi::D2TS(double P1[], double P2[], double P3[], double P4[], int flag)
588{
589 //flag = 1, T = K*; flag = 2, T = rho
590 double temp_PDF = 0;
591 EvtComplex amp_PDF(0,0);
592 double sa[3], sb[3], sc[3], B[3];
593 double t2D[4][4], t2T[4][4];
594 double pS[4], pT[4], pD[4];
595 for(int i=0; i!=4; i++){
596 pS[i] = P3[i] + P4[i];
597 pT[i] = P1[i] + P2[i];
598 pD[i] = pT[i] + pS[i];
599 }
600 sa[0] = dot(pT,pT);
601 sb[0] = dot(P1,P1);
602 sc[0] = dot(P2,P2);
603 sa[1] = dot(pS,pS);
604 sb[1] = dot(P3,P3);
605 sc[1] = dot(P4,P4);
606 sa[2] = dot(pD,pD);
607 sb[2] = sa[0];
608 sc[2] = sa[1];
609 B[0] = barrier(2,sa[0],sb[0],sc[0],rRes);
610 B[2] = barrier(2,sa[2],sb[2],sc[2],rD);
611 calt2(P1,P2,t2T);
612 calt2(pT,pS,t2D);
613 for(int i=0; i!=4; i++){
614 for(int j=0; j!=4; j++){
615 temp_PDF += t2D[i][j]*t2T[j][i]*(G[i][i])*(G[j][j]);
616 }
617 }
618 amp_PDF = temp_PDF*B[0]*B[2];
619 if(flag == 2) amp_PDF *= KPiSFormfactor(sa[1],sb[1],sc[1],rRes);
620 return amp_PDF;
621}
622EvtComplex EvtD0ToKpipipi::PHSP(double Km[], double Pip[])
623{
624 EvtComplex amp_PDF(0,0);
625 double sa,sb,sc;
626 double KPi[4];
627 for(int i=0; i!=4; i++){
628 KPi[i] = Km[i] + Pip[i];
629 }
630 sa = dot(KPi,KPi);
631 sb = dot(Km,Km);
632 sc = dot(Pip,Pip);
633 amp_PDF = KPiSFormfactor(sa,sb,sc,rRes);
634 return amp_PDF;
635}
636double EvtD0ToKpipipi::calDalEva(double P1[], double P2[], double P3[])
637{
638 //Ks Pi0 Pi0 exchange Pi0 in this case
639 //CLEOc model
640 EvtComplex PDF[7], cof, pdf, module;
641 //Ks Pi0 Pi0
642 double P[3][4];
643 for(int i=0; i<4; i++){
644 P[0][i] = P1[i];
645 P[1][i] = P2[i];
646 P[2][i] = P3[i];
647 }
648 PDF[0] = Spin_factor(P[1],P[2],P[0],0) + Spin_factor(P[2],P[1],P[0],0);
649 PDF[1] = Spin_factor(P[1],P[2],P[0],10) + Spin_factor(P[2],P[1],P[0],10);
650 PDF[2] = Spin_factor(P[1],P[2],P[0],100) + Spin_factor(P[2],P[1],P[0],100);
651 //-----------------
652 PDF[3] = Spin_factor(P[0],P[1],P[2],1) + Spin_factor(P[0],P[2],P[1],1);
653 PDF[4] = Spin_factor(P[0],P[1],P[2],11) + Spin_factor(P[0],P[2],P[1],11);
654 //-----------------
655 PDF[5] = Spin_factor(P[1],P[2],P[0],2) + Spin_factor(P[2],P[1],P[0],2);
656 PDF[6] = Spin_factor(P[0],P[1],P[2],12) + Spin_factor(P[0],P[2],P[1],12);
657
658 rho[0] = 4.2551e-01; phi[0] = 3.7370e+00;
659 rho[1] = -1.0086e+00; phi[1] = 6.9153e+00;
660 rho[2] = 6.5290e-01; phi[2] = 2.8041e+00;
661 rho[3] = 1.0; phi[3] = 0.0;
662 rho[4] = -2.5079e-01; phi[4] = 6.3740e-01;
663 rho[5] = -4.1640e-01; phi[5] = 7.7850e+00;
664 rho[6] = 2.1281e-01; phi[6] = 5.3474e+00;
665
666 pdf = EvtComplex(0,0);
667 for(int i=0; i<7; i++){
668 cof = EvtComplex(rho[i]*cos(phi[i]),rho[i]*sin(phi[i]));
669 pdf += cof*PDF[i];
670 }
671 module = conj(pdf)*pdf;
672 double value;
673 value = real(module);
674 return (value <= 0) ? 1e-20 : value;
675}
676EvtComplex EvtD0ToKpipipi::Spin_factor(double P1[], double P2[], double P3[], int spin)
677{
678 // D-> R P3, R->P1 P2
679 double R[4], s[3], sp2, sD, B[2];
680 for(int i=0; i<4; i++){
681 R[i] = P1[i] + P2[i];
682 }
683 s[0] = dot(R,R);
684 s[1] = dot(P1, P1);
685 s[2] = dot(P2, P2);
686 sp2 = dot(P3,P3);
687 sD = mD*mD;
688 EvtComplex amp(0,0), prop;
689 prop = getProp(s, spin);
690 double tmp(0.);
691 if(spin%10 == 0){
692 amp = prop;
693 }
694 else if(spin%10 == 1){
695 tmp = 0;
696 double T1[4], t1[4];
697 calt1(R,P3,T1);
698 calt1(P1,P2,t1);
699 B[0] = barrier(1,s[0],s[1],s[2],3.0);
700 B[1] = barrier(1,sD,s[0],sp2,3.0);
701 for(int i=0; i<4; i++){
702 tmp += T1[i]*t1[i]*G[i][i];
703 }
704 amp = tmp*prop*B[0]*B[1];
705 }
706 else if(spin%10 ==2){
707 tmp = 0;
708 double T2[4][4], t2[4][4];
709 calt2(R,P3,T2);
710 calt2(P1,P2,t2);
711 B[0] = barrier(2,s[0],s[1],s[2],3.0);
712 B[1] = barrier(2,sD,s[0],sp2,3.0);
713 for(int i=0; i<4; i++){
714 for(int j=0; j<4; j++){
715 tmp += T2[i][j]*t2[j][i]*G[j][j]*G[i][i];
716 }
717 }
718 amp = tmp*prop*B[0]*B[1];
719 }
720 else{
721 cout<<"Only S, P, D wave allowed"<<endl;
722 }
723 return amp;
724}
725EvtComplex EvtD0ToKpipipi::getProp(double s[], int flag)
726{
727 //Ks Pi0 Pi0
728 //0 980
729 //10 f0_1370
730 //20 f0_1500
731 //2 f2_1270
732 //1 K*892
733 //12 K*2_1430
734 //11 K*1680
735 //100 f0(500)
736 EvtComplex prop(1,0);
737 double snpi, scpi, snk, sck;
738 snpi = mass_Pi0*mass_Pi0;
739 scpi = mass_Pion*mass_Pion;
740 sck = mass_Kaon*mass_Kaon;
741 snk = mk0*mk0;
742 if(flag == 0){
743 double mass, gpipi, gkk;
744 mass = 0.965;
745 gpipi = 0.406;
746 gkk = 2*gpipi;
747 EvtComplex ci(0,1);
748 EvtComplex rhokk, rhopipi;
749 rhokk = 0.5*(rhofactor(s[0],sck)+rhofactor(s[0],snk));
750 rhopipi = 1.0/3.0 * (2.0 * rhofactor(s[0],scpi) + rhofactor(s[0],snpi));
751 prop = 1.0/(mass*mass - s[0] - ci*(gpipi*gpipi*rhopipi+gkk*gkk*rhokk));
752 }
753 if(flag == 10){
754 prop = propagatorRBW(1.35,0.265,s[0],s[1],s[2],3.0,0);
755 }
756 if(flag == 20){
757 prop = propagatorRBW(1.505,0.109,s[0],s[1],s[2],3.0,0);
758 }
759 if(flag == 1){
760 prop = propagatorRBW(0.896,0.0503,s[0],s[1],s[2],3.0,1);
761 }
762 if(flag == 11){
763 prop = propagatorRBW(1.717,0.322,s[0],s[1],s[2],3.0,1);
764 }
765 if(flag == 2){
766 prop = propagatorRBW(1.2751,0.185,s[0],s[1],s[2],3.0,2);
767 }
768 if(flag == 12){
769 prop = propagatorRBW(1.4324,0.109,s[0],s[1],s[2],3.0,2);
770 }
771 if(flag == 100){
772 EvtComplex pole(0.470,-0.220);
773 prop = 1.0/(s[0] - pole*pole);
774 }
775 return prop;
776}
777EvtComplex EvtD0ToKpipipi::rhofactor(double sx, double sdau)
778{
779 EvtComplex one(1,0);
780 EvtComplex ci(0,1);
781 EvtComplex res;
782 double tmp;
783 tmp = 1 - 4*sdau/sx;
784 if(tmp>0) res = one*sqrt(tmp);
785 if(tmp<0) res = ci*sqrt(fabs(tmp));
786 return res;
787}
788EvtComplex EvtD0ToKpipipi::propogator(double mass, double width, double sx)const
789{
790 EvtComplex ci(0,1);
791 EvtComplex prop=1.0/(mass*mass-sx-ci*mass*width);
792 return prop;
793}
794double EvtD0ToKpipipi::wid(double mass, double sa, double sb, double sc, double r, int l)const
795{
796 double widm(0.), q(0.), q0(0.);
797 double sa0 = mass*mass;
798 double m = sqrt(sa);
799 q = Qabcs(sa,sb,sc);
800 q0 = Qabcs(sa0,sb,sc);
801 double z,z0;
802 z = q*r*r;
803 z0 = q0*r*r;
804 double t = q/q0;
805 double F(0.);
806 if(l == 0) F = 1;
807 if(l == 1) F = sqrt((1+z0)/(1+z));
808 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
809 widm = pow(t,l+0.5)*mass/m*F*F;
810 return widm;
811}
812double EvtD0ToKpipipi::h(double m, double q)const
813{
814 double h(0.);
815 h = 2/pi*q/m*log((m+2*q)/(2*mpi));
816 return h;
817}
818double EvtD0ToKpipipi::dh(double mass, double q0)const
819{
820 double dh = h(mass,q0)*(1.0/(8*q0*q0)-1.0/(2*mass*mass))+1.0/(2*pi*mass*mass);
821 return dh;
822}
823double EvtD0ToKpipipi::f(double mass, double sx, double q0, double q)const
824{
825 double m = sqrt(sx);
826 double f = mass*mass/(pow(q0,3))*(q*q*(h(m,q)-h(mass,q0))+(mass*mass-sx)*q0*q0*dh(mass,q0));
827 return f;
828}
829double EvtD0ToKpipipi::d(double mass, double q0)const
830{
831 double d = 3.0/pi*mpi*mpi/(q0*q0)*log((mass+2*q0)/(2*mpi))+mass/(2*pi*q0) - (mpi*mpi*mass)/(pi*pow(q0,3));
832 return d;
833}
834EvtComplex EvtD0ToKpipipi::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r, int l)const
835{
836 EvtComplex ci(0,1);
837 EvtComplex prop=1.0/(mass*mass-sa-ci*mass*width*wid(mass,sa,sb,sc,r,l));
838 return prop;
839}
840EvtComplex EvtD0ToKpipipi::propagatorGS(double mass, double width, double sa, double sb, double sc, double r, int l)const
841{
842 EvtComplex ci(0,1);
843 double q = Qabcs(sa,sb,sc);
844 double sa0 = mass*mass;
845 double q0 = Qabcs(sa0,sb,sc);
846 q = sqrt(q);
847 q0 = sqrt(q0);
848 EvtComplex prop = (1+d(mass,q0)*width/mass)/(mass*mass-sa+width*f(mass,sa,q0,q)-ci*mass*width*wid(mass,sa,sb,sc,r,l));
849 return prop;
850}
851double EvtD0ToKpipipi::Flatte_rhoab(double sa, double sb, double sc)const
852{
853 double q = Qabcs(sa,sb,sc);
854 double rho = sqrt(q/sa);
855 return rho;
856}
857EvtComplex EvtD0ToKpipipi::propagatorFlatte(double mass, double width, double sx, double *sb, double *sc)const
858{
859 EvtComplex ci(0,1);
860 double rho1 = Flatte_rhoab(sx,sb[0],sc[0]);
861 double rho2 = Flatte_rhoab(sx,sb[1],sc[1]);
862 EvtComplex prop = 1.0/(mass*mass-sx-ci*(g1*g1*rho1+g2*g2*rho2));
863 return prop;
864}
865double EvtD0ToKpipipi::dot(double *a1, double *a2)const
866{
867 double dot = 0;
868 for(int i=0; i!=4; i++){
869 dot += a1[i]*a2[i]*G[i][i];
870 }
871 return dot;
872}
873double EvtD0ToKpipipi::Qabcs(double sa, double sb, double sc)const
874{
875 double Qabcs = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
876 if(Qabcs < 0) Qabcs = 1e-16;
877 return Qabcs;
878}
879double EvtD0ToKpipipi::barrier(double l, double sa, double sb, double sc, double r)const
880{
881 double q = Qabcs(sa,sb,sc);
882 q = sqrt(q);
883 double z = q*r;
884 z = z*z;
885 double F = 1;
886 if(l > 2) F = 0;
887 if(l == 0)
888 F = 1;
889 if(l == 1){
890 F = sqrt((2*z)/(1+z));
891 }
892 if(l == 2){
893 F = sqrt((13*z*z)/(9+3*z+z*z));
894 }
895 return F;
896}
897void EvtD0ToKpipipi::calt1(double daug1[], double daug2[], double t1[]) const
898{
899 double p, pq;
900 double pa[4], qa[4];
901 for(int i=0; i!=4; i++){
902 pa[i] = daug1[i] + daug2[i];
903 qa[i] = daug1[i] - daug2[i];
904 }
905 p = dot(pa,pa);
906 pq = dot(pa,qa);
907 for(int i=0; i!=4; i++){
908 t1[i] = qa[i] - pq/p*pa[i];
909 }
910}
911void EvtD0ToKpipipi::calt2(double daug1[], double daug2[], double t2[][4]) const
912{
913 double p,r;
914 double pa[4], t1[4];
915 calt1(daug1,daug2,t1);
916 r = dot(t1,t1);
917 for(int i=0; i!=4; i++){
918 pa[i] = daug1[i] + daug2[i];
919 }
920 p = dot(pa,pa);
921 for(int i=0; i!=4; i++){
922 for(int j=0; j!=4; j++){
923 t2[i][j] = t1[i]*t1[j] - 1.0/3*r*(G[i][j]-pa[i]*pa[j]/p);
924 }
925 }
926}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
P
double P(RecMdcKalTrack *trk)
Definition: CalibEventSelect.cxx:87
mass
double mass
Definition: CosmicGenerator.cxx:138
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition: Evt3Rank3C.cc:175
s
XmlRpcServer s
Definition: HelloServer.cpp:11
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition: KKsem.h:33
t
TTree * t
Definition: binning.cxx:23
EvtD0ToKpipipi::clone
EvtDecayBase * clone()
Definition: EvtD0ToKpipipi.cc:37
EvtD0ToKpipipi::getName
void getName(std::string &name)
Definition: EvtD0ToKpipipi.cc:33
EvtD0ToKpipipi::decay
void decay(EvtParticle *p)
Definition: EvtD0ToKpipipi.cc:169
EvtD0ToKpipipi::~EvtD0ToKpipipi
virtual ~EvtD0ToKpipipi()
Definition: EvtD0ToKpipipi.cc:31
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition: EvtDecayBase.cc:520
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition: EvtDecayBase.cc:297
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition: EvtDecayBase.cc:504
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition: EvtDecayBase.hh:65
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition: EvtDecayBase.cc:482
EvtDecayProb::setProb
void setProb(double prob)
Definition: EvtDecayProb.hh:34
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:120
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:84
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cc:1070
EvtVector4R::get
double get(int i) const
Definition: EvtVector4R.hh:179
rb::R
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition: TUtil.h:27
one
double precision pisqo6 one
Definition: qlconstants.h:4