BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKppipi.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: EvtDToKppipi.cc
11//
12// Description: D+ -> K+ pi+ pi-
13// See https://indico.ihep.ac.cn/event/18455/contributions/146055/attachments/74328/91168/DCS_D%2BtoK%2Bpi%2Bpi-.pdf
14//
15// Modification history:
16// Liaoyuan Dong Feb 17 23:19:33 2024 Module updated
17// Oct 10 00:03:38 2023 Module created
18//------------------------------------------------------------------------
20#include <fstream>
21#include <stdlib.h>
24#include "EvtGenBase/EvtPDL.hh"
29#include <string>
30#include <vector>
31#include <math.h>
32#include "TMath.h"
35using std::endl;
36
38
39void EvtDToKppipi::getName(std::string& model_name){
40 model_name="DToKppipi";
41}
42
46
48 // check that there are 0 arguments
49 checkNArg(0);
50 checkNDaug(3);
55
56 pi180inv = 1.0/EvtConst::radToDegrees;
57 phi[0] = 0.0;
58 phi[1] = 3.365002324407079470347526;
59 phi[2] = 0.6621583377125883629332748;
60 phi[3] = 0.4339644541377705166951273;
61 phi[4] = 2.496788137814539787484591;
62 phi[5] = 5.337638327720493514050304;
63
64 rho[0] = 1.0;
65 rho[1] = 2.176692936463668459623477;
66 rho[2] = 1.815490493394110060876301;
67 rho[3] = 0.3832709830473302048403639;
68 rho[4] = 1.12572937386981664076302;
69 rho[5] = -1.190394848723933307610423;
70
71 modetype[0]= 13; /////K*(892)
72 modetype[1]= 13; /////K*(1410)
73 modetype[2]= 13; /////K*2(1430)
74 modetype[3]= 23; /////PiPi
75 modetype[4]= 23; /////rho(770)
76 modetype[5]= 23; /////f2(1270)
77
78 //std::cout << "Initializing EvtDToKppipi" << std::endl;
79 //for (int i=0; i<9; i++) {
80 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
81 //}
82
83 width[0] = 0.0473000000000000017652546;
84 width[1] = 0.2320000000000000117683641;
85 width[2] = 0.1089999999999999996669331;
86 width[3] = 0.0200000000000000004163336;
87 width[4] = 0.1474000000000000032418512;
88 width[5] = 0.1867000000000000048405724;
89
90 mass[0] = 0.8955499999999999571898002;
91 mass[1] = 1.4139999999999999236166559;
92 mass[2] = 1.4323999999999998955502178;
93 mass[3] = 0.9000000000000000222044605;
94 mass[4] = 0.7752599999999999491606673;
95 mass[5] = 1.2755000000000000781597009;
96
97 mDp = 1.86966;
98 mK0 = 0.497611;
99 mKa = 0.49368;
100 mPi = 0.13957;
101 meta = 0.54775;
102 mK02 = 0.237616707;
103 mPi2 = 0.01947978;
104 mass_EtaP = 0.95778;
105 mass_Kaon = 0.49368;
106
107 math_pi = 3.1415926;
108 mass_Pion = 0.13957;
109 mass_Pion2 = 0.0194797849;
110 mass_2Pion = 0.27914;
111 math_2pi = 6.2831852;
112 rD2 = 25.0; // 5*5
113 rRes2 = 9.0; // 3*3
114 g1 = 0.5468;
115 g2 = 0.23;
116 GS1 = 0.636619783;
117 GS2 = 0.01860182466;
118 GS3 = 0.1591549458;
119 GS4 = 0.00620060822;
120
121 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
122 for (int i=0; i<4; i++) {
123 for (int j=0; j<4; j++) {
124 G[i][j] = GG[i][j];
125 }
126 }
127}
128
130 setProbMax(717.25);
131}
132
134/*
135 double maxprob = 0.0;
136 for(int ir=0;ir<=60000000;ir++){
137 p->initializePhaseSpace(getNDaug(),getDaugs());
138 EvtVector4R D1 = p->getDaug(0)->getP4();
139 EvtVector4R D2 = p->getDaug(1)->getP4();
140 EvtVector4R D3 = p->getDaug(2)->getP4();
141
142 double P1[4], P2[4], P3[4];
143 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
144 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
145 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
146
147 double value;
148 int nstates=6;
149 int g0[6]={1,1,1,6,1,2};
150 int spin[6]={1,1,2,0,1,2};
151 double r0[6]={3.0,3.0,3.0,3.0,3.0,3.0};
152 double r1[6]={5.0,5.0,5.0,5.0,5.0,5.0};
153
154 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value,r0,r1);
155
156 if (value<0) continue;
157 if(value>maxprob) {
158 maxprob=value;
159 cout << "ir= " << ir << endl;
160// cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3] <<"};"<< endl;
161// cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3] <<"};"<< endl;
162// cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< endl;
163 cout << "MAX====> " << maxprob << endl;
164 }
165 }
166 printf("MAXprob = %.10f\n",maxprob);
167*/
169 EvtVector4R D1 = p->getDaug(0)->getP4();
170 EvtVector4R D2 = p->getDaug(1)->getP4();
171 EvtVector4R D3 = p->getDaug(2)->getP4();
172
173 double P1[4], P2[4], P3[4];
174 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
175 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
176 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
177
178 double value;
179 int nstates=6;
180 int g0[6]={1,1,1,6,1,2};
181 int spin[6]={1,1,2,0,1,2};
182 double r0[6]={3.0,3.0,3.0,3.0,3.0,3.0};
183 double r1[6]={5.0,5.0,5.0,5.0,5.0,5.0};
184
185 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value,r0,r1);
186 setProb(value);
187 return ;
188
189}
190
191void EvtDToKppipi::Com_Multi(double a1[2], double a2[2], double res[2])
192{
193 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
194 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
195}
196void EvtDToKppipi::Com_Divide(double a1[2], double a2[2], double res[2])
197{
198 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
199 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
200 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
201}
202double EvtDToKppipi::CalRho4pi(double_t s)
203{
204 if(s>=1.){
205 return sqrt((s-16.*mass_Pion*mass_Pion)/s);
206 }
207 else{
208 double_t s0 = 1.2274+0.00370909/(s*s) - (0.111203)/(s) - 6.39017*s +16.8358*s*s - 21.8845*s*s*s + 11.3153*s*s*s*s;
209 double_t gam = s0*sqrt(1.0-(16.0*mass_Pion*mass_Pion));
210
211 return gam;
212 }
213}
214
215//------------base---------------------------------
216double EvtDToKppipi::SCADot(double a1[4], double a2[4])
217{
218 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
219 return _cal;
220}
221double EvtDToKppipi::barrier(int l, double sa, double sb, double sc, double r, double mass)
222{
223 double q = fabs((sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb);
224 double z;
225 z = q*r*r;
226 double sa0;
227 sa0 = mass*mass;
228 double q0 = fabs((sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb);
229 double z0 = q0*r*r;
230 double F = 0.0;
231 if(l == 0) F = 1;
232 if(l == 1) F = sqrt((1+z0)/(1+z));
233 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
234 return F;
235}
236//------------------spin-------------------------------------------
237void EvtDToKppipi::calt1(double daug1[4], double daug2[4], double t1[4])
238{
239 double p, pq, tmp;
240 double pa[4], qa[4];
241 for(int i=0; i<4; i++) {
242 pa[i] = daug1[i] + daug2[i];
243 qa[i] = daug1[i] - daug2[i];
244 }
245 p = SCADot(pa,pa);
246 pq = SCADot(pa,qa);
247 tmp = pq/p;
248 for(int i=0; i<4; i++) {
249 t1[i] = qa[i] - tmp*pa[i];
250 }
251}
252void EvtDToKppipi::calt2(double daug1[4], double daug2[4], double t2[4][4])
253{
254 double p, r;
255 double pa[4], t1[4];
256 calt1(daug1,daug2,t1);
257 r = SCADot(t1,t1)/3.0;
258 for(int i=0; i<4; i++) {
259 pa[i] = daug1[i] + daug2[i];
260 }
261 p = SCADot(pa,pa);
262 for(int i=0; i<4; i++) {
263 for(int j=0; j<4; j++) {
264 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
265 }
266 }
267}
268//-------------------prop--------------------------------------------
269void EvtDToKppipi::propagator(double mass2, double mass, double width, double sx, double prop[2])
270{
271 double a[2], b[2];
272 a[0] = 1;
273 a[1] = 0;
274 b[0] = mass2-sx;
275 b[1] = -mass*width;
276 Com_Divide(a,b,prop);
277}
278double EvtDToKppipi::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
279{
280 double widm = 0.;
281 double m = sqrt(sa);
282 double tmp = sb-sc;
283 double tmp1 = sa+tmp;
284 double q = fabs(0.25*tmp1*tmp1/sa-sb);
285 double tmp2 = mass2+tmp;
286 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
287 double z = q*r2;
288 double z0 = q0*r2;
289 double t = q/q0;
290 if(l == 0) {widm = sqrt(t)*mass/m;}
291 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
292 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
293 return widm;
294}
295double EvtDToKppipi::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
296{
297 double widm = 0.;
298 double m = sqrt(sa);
299 double tmp = sb-sc;
300 double tmp1 = sa+tmp;
301 double q = fabs(0.25*tmp1*tmp1/sa-sb);
302 double tmp2 = mass2+tmp;
303 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
304 double z = q*r2;
305 double z0 = q0*r2;
306 double F = (1+z0)/(1+z);
307 double t = q/q0;
308 widm = t*sqrt(t)*mass/m*F;
309 return widm;
310}
311void EvtDToKppipi::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
312{
313 double a[2], b[2];
314 double mass2 = mass*mass;
315 a[0] = 1;
316 a[1] = 0;
317 b[0] = mass2-sa;
318 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
319 Com_Divide(a,b,prop);
320}
321
322void EvtDToKppipi::propagatorFlatte(double mass, double width, double sa, double prop[2]){
323 double q2_Pi, q2_Ka;
324 double rhoPi[2], rhoKa[2];
325
326 q2_Pi = 0.25*sa-mPi*mPi;
327 q2_Ka = 0.25*sa-mKa*mKa;
328
329 if (q2_Pi > 0) {
330 rhoPi[0] = 2.0*sqrt(q2_Pi/sa);
331 rhoPi[1] = 0.0;
332 }
333 if (q2_Pi <= 0) {
334 rhoPi[0] = 0.0;
335 rhoPi[1] = 2.0*sqrt(-q2_Pi/sa);
336 }
337
338 if (q2_Ka > 0) {
339 rhoKa[0] = 2.0*sqrt(q2_Ka/sa);
340 rhoKa[1] = 0.0;
341 }
342 if (q2_Ka <= 0) {
343 rhoKa[0] = 0.0;
344 rhoKa[1] = 2.0*sqrt(-q2_Ka/sa);
345 }
346
347 double a[2], b[2];
348 a[0] = 1;
349 a[1] = 0;
350 b[0] = mass*mass - sa + 0.165*rhoPi[1] + 0.69465*rhoKa[1];
351 b[1] = - (0.165*rhoPi[0] + 0.69465*rhoKa[0]);
352 Com_Divide(a,b,prop);
353}
354
355void EvtDToKppipi::rhoab(double sa, double sb, double sc, double res[2]) {
356 double tmp = sa+sb-sc;
357 double q = 0.25*tmp*tmp/sa-sb;
358 if(q>=0) {
359 res[0]=2.0*sqrt(q/sa);
360 res[1]=0.0;
361 } else {
362 res[0]=0.0;
363 res[1]=2.0*sqrt(-q/sa);
364 }
365}
366void EvtDToKppipi::rho4Pi(double sa, double res[2]) {
367 double temp = 1.0-0.3116765584/sa; // 0.3116765584=0.13957*0.13957*16
368 if(temp>=0) {
369 res[0]=sqrt(temp)/(1.0+exp(9.8-3.5*sa));
370 res[1]=0.0;
371 } else {
372 res[0]=0.0;
373 res[1]=sqrt(-temp)/(1.0+exp(9.8-3.5*sa));
374 }
375}
376
377void EvtDToKppipi::propagatorsigma500(double sa, double sb, double sc, double prop[2]) {
378 double f = 0.5843+1.6663*sa;
379 const double M = 0.9264;//M(f0500)
380 const double mass2 = 0.85821696; // M*M
381 const double mpi2d2 = 0.00973989245;//mass_Pion2/2 = 0.0194797849/2
382 double g1 = f*(sa-mpi2d2)/(mass2-mpi2d2)*exp((mass2-sa)/1.082);
383 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
384 rhoab(sa,sb,sc,rho1s);
385 rhoab(mass2,sb,sc,rho1M);
386 rho4Pi(sa,rho2s);
387 rho4Pi(mass2,rho2M);
388 Com_Divide(rho1s,rho1M,rho1);
389 Com_Divide(rho2s,rho2M,rho2);
390 double a[2], b[2];
391 a[0] = 1.0;
392 a[1] = 0.0;
393 b[0] = mass2-sa+M*(g1*rho1[1]+0.0024*rho2[1]);
394 b[1] = -M*(g1*rho1[0]+0.0024*rho2[0]);
395 Com_Divide(a,b,prop);
396}
397void EvtDToKppipi::Flatte_rhoab(double sa, double sb, double rho[2])
398{
399 double q = 1.0-(4*sb/sa);
400
401 if(q>0) {
402 rho[0]=sqrt(q);
403 rho[1]=0;
404 }
405 else if(q<0){
406 rho[0]=0;
407 rho[1]=sqrt(-q);
408 }
409}
410
411void EvtDToKppipi::propagator980(double mass, double sx, double *sb, double prop[2])
412{
413 double gpipi1 = 2.0/3.0*0.165;
414 double gpipi2 = 1.0/3.0*0.165;
415 double gKK1 = 0.5*0.69465;
416 double gKK2 = 0.5*0.69465;
417
418 double unit[2]={1.0};
419 double ci[2]={0,1};
420 double rho1[2];
421 Flatte_rhoab(sx,sb[0],rho1);
422 double rho2[2];
423 Flatte_rhoab(sx,sb[1],rho2);
424 double rho3[2];
425 Flatte_rhoab(sx,sb[2],rho3);
426 double rho4[2];
427 Flatte_rhoab(sx,sb[3],rho4);
428
429 double tmp1[2]={gpipi1,0};
430 double tmp11[2];
431 double tmp2[2]={gpipi2,0};
432 double tmp22[2];
433 double tmp3[2]={gKK1,0};
434 double tmp33[2];
435 double tmp4[2]={gKK2,0};
436 double tmp44[2];
437
438 Com_Multi(tmp1,rho1,tmp11);
439 Com_Multi(tmp2,rho2,tmp22);
440 Com_Multi(tmp3,rho3,tmp33);
441 Com_Multi(tmp4,rho4,tmp44);
442
443 double tmp5[2]={tmp11[0]+tmp22[0]+tmp33[0]+tmp44[0],tmp11[1]+tmp22[1]+tmp33[1]+tmp44[1]};
444 double tmp51[2];
445 Com_Multi(tmp5, ci,tmp51);
446 double tmp6[2]={mass*mass-sx-tmp51[0], -1.0*tmp51[1]};
447
448 Com_Divide( unit,tmp6, prop);
449
450}
451
452void EvtDToKppipi::KPiSLASS(double sa, double sb, double sc, double prop[2])
453{
454 const double m1430 = 1.441;
455 const double sa0 = 2.076481; // m1430*m1430;
456 const double w1430 = 0.193;
457 const double Lass1 = 0.25/sa0;
458 double tmp = sb-sc;
459 double tmp1 = sa0+tmp;
460 double q0 = fabs(Lass1*tmp1*tmp1-sb);
461 double tmp2 = sa+tmp;
462 double qs = 0.25*tmp2*tmp2/sa-sb;
463 double q = sqrt(qs);
464 double width = w1430*q*m1430/sqrt(sa*q0);
465 double temp_R = atan(m1430*width/(sa0-sa));
466 if(temp_R<0) temp_R += math_pi;
467 double deltaR = -109.7*math_pi/180.0 + temp_R;
468 double temp_F = atan(0.226*q/(2.0-3.8194*qs)); // 2.0*0.113 = 0.226; -33.8*0.113 = -3.8194
469 if(temp_F<0) temp_F += math_pi;
470 double deltaF = 0.1*math_pi/180.0 + temp_F;
471 double deltaS = deltaR + 2.0*deltaF;
472 double t1 = 0.96*sin(deltaF);
473 double t2 = sin(deltaR);
474 double CF[2], CS[2];
475 CF[0] = cos(deltaF);
476 CF[1] = sin(deltaF);
477 CS[0] = cos(deltaS);
478 CS[1] = sin(deltaS);
479 prop[0] = t1*CF[0] + t2*CS[0];
480 prop[1] = t1*CF[1] + t2*CS[1];
481}
482//------------GS---used by rho----------------------------
483void EvtDToKppipi::propagatorGS(double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
484{
485 double a[2], b[2];
486 double mass2 = mass*mass;
487 double tmp = sb-sc;
488 double tmp1 = sa+tmp;
489 double q2 = fabs(0.25*tmp1*tmp1/sa-sb);
490
491 double tmp2 = mass2+tmp;
492 double q02 = fabs(0.25*tmp2*tmp2/mass2-sb);
493
494 double q = sqrt(q2);
495 double q0 = sqrt(q02);
496 double m = sqrt(sa);
497 double q03 = q0*q02;
498 double tmp3 = log(mass+2*q0)+1.2760418309; // log(mpi0+mpip) = -1.2926305904;
499
500 double h = GS1*q/m*(log(m+2*q)+1.2760418309);
501 double h0 = GS1*q0/mass*tmp3;
502 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
503 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
504 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
505
506 a[0] = 1.0+d*width/mass;
507 a[1] = 0.0;
508 b[0] = mass2-sa+width*f;
509 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
510 Com_Divide(a,b,prop);
511}
512
513double EvtDToKppipi::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
514 double pR[4], pD[4];
515 double temp_PDF, v_re;
516 temp_PDF = 0.0;
517 v_re = 0.0;
518 double B[2], s1, s2, s3, sR, sD;
519 for(int i=0; i<4; i++){
520 pR[i] = P1[i] + P2[i];
521 pD[i] = pR[i] + P3[i];
522 }
523 s1 = SCADot(P1,P1);
524 s2 = SCADot(P2,P2);
525 s3 = SCADot(P3,P3);
526 sR = SCADot(pR,pR);
527 sD = SCADot(pD,pD);
528 int G[4][4];
529 for(int i=0; i!=4; i++){
530 for(int j=0; j!=4; j++){
531 if(i==j){
532 if(i==0) G[i][j] = 1;
533 else G[i][j] = -1;
534 }
535 else G[i][j] = 0;
536 }
537 }
538 if(Ang == 0){
539 B[0] = 1;
540 B[1] = 1;
541 temp_PDF = 1;
542 }
543 if(Ang == 1){
544 B[0] = barrier(1,sR,s1,s2,3.0,mass);
545 B[1] = barrier(1,sD,sR,s3,5.0,mDp);
546 double t1[4], T1[4];
547 calt1(P1,P2,t1);
548 calt1(pR,P3,T1);
549 temp_PDF = 0;
550 for(int i=0; i<4; i++){
551 temp_PDF += t1[i]*T1[i]*G[i][i];
552 }
553 }
554 if(Ang == 2){
555 B[0] = barrier(2,sR,s1,s2,3.0,mass);
556 B[1] = barrier(2,sD,sR,s3,5.0,mDp);
557 double t2[4][4], T2[4][4];
558 calt2(P1,P2,t2);
559 calt2(pR,P3,T2);
560 temp_PDF = 0;
561 for(int i=0; i<4; i++){
562 for(int j=0; j<4; j++){
563 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
564 }
565 }
566 }
567 v_re = temp_PDF*B[0]*B[1];
568 return v_re;
569}
570
571void EvtDToKppipi::rhoMTX(int i, int j, double s, double Rho[2]){
572
573 double rhoijx;
574 double rhoijy;
575 double mpi = 0.13957;
576 if(i==j && i==0 ){
577 double m2 = 0.13957*0.13957;
578 if((1-(4*m2)/s)>0){
579 rhoijx = sqrt(1.0f - (4*m2)/s);
580 rhoijy = 0;
581 }else{
582 rhoijy = sqrt((4*m2)/s - 1.0f);
583 rhoijx = 0;
584 }
585 }
586 if(i==j && i==1 ){
587 double m2 = 0.49368*0.49368;
588 if((1-(4*m2)/s)>0){
589 rhoijx = sqrt(1.0f - (4*m2)/s);
590 rhoijy = 0;
591 }else{
592 rhoijy = sqrt((4*m2)/s - 1.0f);
593 rhoijx = 0;
594 }
595 }
596 if(i==j && i==2){
597 rhoijx = CalRho4pi(s);
598 rhoijy = 0;
599 }
600 if(i==j && i==3){
601 double m2 = 0.547862*0.547862;
602 if((1-(4*m2)/s)>0){
603 rhoijx = sqrt(1.0f - (4*m2)/s);
604 rhoijy = 0;
605 }else{
606 rhoijy = sqrt((4*m2)/s - 1.0f);
607 rhoijx = 0;
608 }
609 }
610 if(i==j && i==4){
611 double m_1 = 0.547862;
612 double m_2 = 0.95778;
613 double mp2 = (m_1+m_2)*(m_1+m_2);
614 double mm2 = (m_1-m_2)*(m_1-m_2);
615 if((1-mp2/s)>0){
616 rhoijx = sqrt(1.0f - mp2/s);
617 rhoijy = 0;
618 }else{
619 rhoijy = sqrt(mp2/s - 1.0f);
620 rhoijx = 0;
621 }
622 }
623
624 if(i!=j){
625 rhoijx = 0;
626 rhoijy = 0;
627 }
628 Rho[0] = rhoijx;
629 Rho[1] = rhoijy;
630
631}
632
633
634/* Kij = (sum_a gigj/(ma^2-s) + fij(1-s0)/(s-s0))((s-0.5sampi2)(1-sa0)/(s-sa0)) */
635void EvtDToKppipi::KMTX(int i, int j, double s,double KM[2]){
636
637 double Kijx;
638 double Kijy;
639 double mpi = 0.13957;
640 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206};
641
642 double g1[5] = { 0.22889,-0.55377, 0.00000,-0.39899,-0.34639};
643 double g2[5] = { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503};
644 double g3[5] = { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681};
645 double g4[5] = { 0.33650, 0.40907, 0.85679, 0.19906,-0.00984};
646 double g5[5] = { 0.18171,-0.17558,-0.79658,-0.00355, 0.22358};
647
648 double f1[5] = { 0.23399, 0.15044,-0.20545, 0.32825, 0.35412};
649
650 double upimag[5] = { 0,0,0,0,0};
651
652 for(int k=0; k<5; k++){
653 upimag[k] = 0;
654 }
655 double ss0 = -3.92637;
656 double sA = 1.0;//v1
657 double sA0 = -0.15;
658
659 if(i==0||j==0){
660 Kijx = (g1[i]*g1[j]/(m[0]*m[0]-s) + g2[i]*g2[j]/(m[1]*m[1]-s) + g3[i]*g3[j]/(m[2]*m[2]-s) + g4[i]*g4[j]/(m[3]*m[3]-s) + g5[i]*g5[j]/(m[4]*m[4]-s)+f1[j]*(1-ss0)/(s-ss0))*(1-sA0)/(s-sA0)*(s-sA*mpi*mpi*0.5);
661 Kijy = (g1[i]*g1[j]*upimag[0] + g2[i]*g2[j]*upimag[1] + g3[i]*g3[j]*upimag[2] + g4[i]*g4[j]*upimag[3] + g5[i]*g5[j]*upimag[4])*(1-sA0)/(s-sA0)*(s-sA*mpi*mpi*0.5);
662 }
663
664 else{
665 Kijx = (g1[i]*g1[j]/(m[0]*m[0]-s) + g2[i]*g2[j]/(m[1]*m[1]-s) + g3[i]*g3[j]/(m[2]*m[2]-s) + g4[i]*g4[j]/(m[3]*m[3]-s) + g5[i]*g5[j]/(m[4]*m[4]-s))*(1-sA0)/(s-sA0)*(s-sA*mpi*mpi*0.5);
666 Kijy = (g1[i]*g1[j]*upimag[0] + g2[i]*g2[j]*upimag[1] + g3[i]*g3[j]*upimag[2] + g4[i]*g4[j]*upimag[3] + g5[i]*g5[j]*upimag[4])*(1-sA0)/(s-sA0)*(s-sA*mpi*mpi*0.5);
667 }
668
669 KM[0] = Kijx;
670 KM[1] = Kijy;
671}
672
673
674void EvtDToKppipi::IMTX(int i, int j, double IMTX[2]){
675
676 double Iijx;
677 double Iijy;
678 if(i==j){
679 Iijx = 1;
680 Iijy = 0;
681 }else{
682 Iijx = 0;
683 Iijy = 0;
684 }
685 IMTX[0] = Iijx;
686 IMTX[1] = Iijy;
687
688}
689
690/* F = I - i*K*rho */
691void EvtDToKppipi::FMTX(double Kijx, double Kijy, double rhojjx, double rhojjy, int i, int j, double FM[2]){
692
693 double Fijx;
694 double Fijy;
695
696 double tmpx = Kijx*rhojjx - Kijy*rhojjy;
697 double tmpy = Kijy*rhojjx + Kijx*rhojjy;
698
699 double imtx[2];
700 IMTX(i,j,imtx);
701 Fijx = imtx[0]+tmpy;
702 Fijy = -tmpx;
703
704 FM[0] = Fijx;
705 FM[1] = Fijy;
706
707}
708
709void EvtDToKppipi::PVTR(int ID, double s, double PV[2]){
710
711 double VPix;
712 double VPiy;
713 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206};
714
715 double g[5][5] = {{ 0.22889,-0.55377, 0.00000,-0.39899,-0.34639},
716 { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503},
717 { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681},
718 { 0.33650, 0.40907, 0.85679, 0.19906,-0.00984},
719 { 0.18171,-0.17558,-0.79658,-0.00355, 0.22358}};
720
721 double betax[5], betay[5], fprodx[5], fprody[5];
722
723 betax[0] = 8.5*cos( 68.5*3.1415926/180.0); betay[0] = 8.5*sin( 68.5*3.1415926/180.0);
724 betax[1] = 12.2*cos( 24.0*3.1415926/180.0); betay[1] = 12.2*sin( 24.0*3.1415926/180.0);
725 betax[2] = 29.2*cos( -0.1*3.1415926/180.0); betay[2] = 29.2*sin( -0.1*3.1415926/180.0);
726 betax[3] = 10.8*cos(-51.9*3.1415926/180.0); betay[3] = 10.8*sin(-51.9*3.1415926/180.0);
727 betax[4] = 0.0; betay[4] = 0.0;
728
729 fprodx[0] = 8.0*cos(-126.0*3.1415926/180.0); fprody[0] = 8.0*sin(-126.0*3.1415926/180.0);
730 fprodx[1] = 26.3*cos(-152.3*3.1415926/180.0); fprody[1] = 26.3*sin(-152.3*3.1415926/180.0);
731 fprodx[2] = 33.0*cos( -93.2*3.1415926/180.0); fprody[2] = 33.0*sin( -93.2*3.1415926/180.0);
732 fprodx[3] = 26.2*cos(-121.4*3.1415926/180.0); fprody[3] = 26.2*sin(-121.4*3.1415926/180.0);
733 fprodx[4] = 0.0; fprody[4] = 0.0;
734
735 double V0x = 0.0, V0y = 0.0, V1x = 0.0, V1y = 0.0;
736 double s0_prod = -0.07;
737
738 for(int k=0;k<5;k++) {
739 V0x += betax[k]*g[k][ID]/(m[k]*m[k]-s);
740 V0y += betay[k]*g[k][ID]/(m[k]*m[k]-s);
741 }
742 V1x += (1.-s0_prod)/(s-s0_prod)*fprodx[ID];
743 V1y += (1.-s0_prod)/(s-s0_prod)*fprody[ID];
744
745 VPix = V0x+V1x;
746 VPiy = V0y+V1y;
747
748 PV[0] = VPix;
749 PV[1] = VPiy;
750
751}
752
753/* inverse for Matrix (I - i*rho*K) using LUP */
754void EvtDToKppipi::FINVMTX(double s, double *FINVx, double *FINVy){
755
756 int P[5] = { 0,1,2,3,4};
757
758 double Fx[5][5];
759 double Fy[5][5];
760
761 double Ux[5][5];
762 double Uy[5][5];
763 double Lx[5][5];
764 double Ly[5][5];
765
766 double UIx[5][5];
767 double UIy[5][5];
768 double LIx[5][5];
769 double LIy[5][5];
770
771 double rho[2];
772 double KM[2];
773 for(int k=0; k<5; k++){
774 rhoMTX(k,k,s,rho);
775 double rhokkx = rho[0];
776 double rhokky = rho[1];
777 Ux[k][k] = rhokkx;
778 Uy[k][k] = rhokky;
779 for(int l=k; l<5; l++){
780 KMTX(k,l,s,KM);
781 double Kklx = KM[0];
782 double Kkly = KM[1];
783 Lx[k][l] = Kklx;
784 Ly[k][l] = Kkly;
785 Lx[l][k] = Lx[k][l];
786 Ly[l][k] = Ly[k][l];
787 }
788 }
789
790 double AA[2];
791 for(int k=0; k<5; k++){
792 for(int l=0; l<5; l++){
793 FMTX(Lx[k][l],Ly[k][l],Ux[l][l],Uy[l][l],k,l,AA);
794 double Fklx = AA[0];
795 double Fkly = AA[1];
796 Fx[k][l] = Fklx;
797 Fy[k][l] = Fkly;
798 }
799 }
800
801 for(int k=0; k<5; k++){
802 double tmprM = (Fx[k][k]*Fx[k][k]+Fy[k][k]*Fy[k][k]);
803 int tmpID = 0;
804 for(int l=k; l<5; l++){
805 double tmprF = (Fx[l][k]*Fx[l][k]+Fy[l][k]*Fy[l][k]);
806 if(tmprM<=tmprF){
807 tmprM = tmprF;
808 tmpID = l;
809 }
810 }
811
812 int tmpP = P[k];
813 P[k] = P[tmpID];
814 P[tmpID] = tmpP;
815
816 for(int l=0; l<5; l++){
817
818 double tmpFx = Fx[k][l];
819 double tmpFy = Fy[k][l];
820
821 Fx[k][l] = Fx[tmpID][l];
822 Fy[k][l] = Fy[tmpID][l];
823
824 Fx[tmpID][l] = tmpFx;
825 Fy[tmpID][l] = tmpFy;
826
827 }
828 for(int l=k+1; l<5; l++){
829 double rFkk = Fx[k][k]*Fx[k][k] + Fy[k][k]*Fy[k][k];
830 double Fxlk = Fx[l][k];
831 double Fylk = Fy[l][k];
832 double Fxkk = Fx[k][k];
833 double Fykk = Fy[k][k];
834 Fx[l][k] = (Fxlk*Fxkk + Fylk*Fykk)/rFkk;
835 Fy[l][k] = (Fylk*Fxkk - Fxlk*Fykk)/rFkk;
836 for(int m=k+1; m<5; m++){
837 Fx[l][m] = Fx[l][m] - (Fx[l][k]*Fx[k][m] - Fy[l][k]*Fy[k][m]);
838 Fy[l][m] = Fy[l][m] - (Fx[l][k]*Fy[k][m] + Fy[l][k]*Fx[k][m]);
839 }
840 }
841 }
842
843 for(int k=0; k<5; k++){
844 for(int l=0; l<5 ;l++){
845 if(k==l){
846 Lx[k][k] = 1;
847 Ly[k][k] = 0;
848 Ux[k][k] = Fx[k][k];
849 Uy[k][k] = Fy[k][k];
850 }
851 if(k>l){
852 Lx[k][l] = Fx[k][l];
853 Ly[k][l] = Fy[k][l];
854 Ux[k][l] = 0;
855 Uy[k][l] = 0;
856 }
857 if(k<l){
858 Ux[k][l] = Fx[k][l];
859 Uy[k][l] = Fy[k][l];
860 Lx[k][l] = 0;
861 Ly[k][l] = 0;
862 }
863 }
864 }
865
866 for(int k=0; k<5; k++){
867
868 LIx[k][k] = 1;
869 LIy[k][k] = 0;
870
871 double rUkk = Ux[k][k]*Ux[k][k] + Uy[k][k]*Uy[k][k];
872 UIx[k][k] = Ux[k][k]/rUkk;
873 UIy[k][k] = -1.0f * Uy[k][k]/rUkk ;
874
875 for(int l=(k+1); l<5; l++){
876 LIx[k][l] = 0;
877 LIy[k][l] = 0;
878 UIx[l][k] = 0;
879 UIy[l][k] = 0;
880 }
881 for(int l=(k-1); l>=0; l--){ // U-1
882 double sx = 0;
883 double c_sx = 0;
884 double sy = 0;
885 double c_sy = 0;
886 for(int m=l+1; m<=k; m++){
887 sx = sx - c_sx;
888 double sx_tmp = sx + Ux[l][m]*UIx[m][k] - Uy[l][m]*UIy[m][k];
889 c_sx = (sx_tmp - sx) - (Ux[l][m]*UIx[m][k] - Uy[l][m]*UIy[m][k]);
890 sx = sx_tmp;
891
892 sy = sy - c_sy;
893 double sy_tmp = sy + Ux[l][m]*UIy[m][k] + Uy[l][m]*UIx[m][k];
894 c_sy = (sy_tmp - sy) - (Ux[l][m]*UIy[m][k] + Uy[l][m]*UIx[m][k]);
895 sy = sy_tmp;
896 }
897 UIx[l][k] = -1.0f * (UIx[l][l]*sx - UIy[l][l]*sy);
898 UIy[l][k] = -1.0f * (UIy[l][l]*sx + UIx[l][l]*sy);
899 }
900
901 for(int l=k+1; l<5; l++){ // L-1
902 double sx = 0;
903 double c_sx = 0;
904 double sy = 0;
905 double c_sy = 0;
906 for(int m=k; m<l; m++){
907 sx = sx - c_sx;
908 double sx_tmp = sx + Lx[l][m]*LIx[m][k] - Ly[l][m]*LIy[m][k];
909 c_sx = (sx_tmp - sx) - (Lx[l][m]*LIx[m][k] - Ly[l][m]*LIy[m][k]);
910 sx = sx_tmp;
911
912 sy = sy - c_sy;
913 double sy_tmp = sy + Lx[l][m]*LIy[m][k] + Ly[l][m]*LIx[m][k];
914 c_sy = (sy_tmp - sy) - (Lx[l][m]*LIy[m][k] + Ly[l][m]*LIx[m][k]);
915 sy = sy_tmp;
916 }
917 LIx[l][k] = -1.0f * sx;
918 LIy[l][k] = -1.0f * sy;
919 }
920 }
921 for(int m=0; m<5; m++){
922 double resX = 0;
923 double c_resX = 0;
924 double resY = 0;
925 double c_resY = 0;
926 for(int k=0; k<5; k++){
927 for(int l=0; l<5; l++){
928 double Plm = 0;
929 if(P[l] == m) Plm = 1;
930
931 resX = resX - c_resX;
932 double resX_tmp = resX + (UIx[0][k]*LIx[k][l] - UIy[0][k]*LIy[k][l])*Plm;
933 c_resX = (resX_tmp - resX) - ((UIx[0][k]*LIx[k][l] - UIy[0][k]*LIy[k][l])*Plm);
934 resX = resX_tmp;
935
936 resY = resY - c_resY;
937 double resY_tmp = resY + (UIx[0][k]*LIy[k][l] + UIy[0][k]*LIx[k][l])*Plm;
938 c_resY = (resY_tmp - resY) - ((UIx[0][k]*LIy[k][l] + UIy[0][k]*LIx[k][l])*Plm);
939 resY = resY_tmp;
940 }
941 }
942 FINVx[m] = resX;
943 FINVy[m] = resY;
944 }
945}
946
947void EvtDToKppipi::Fvector(double sa, double s0, double Fv[2]){
948
949 double outputx = 0;
950 double outputy = 0;
951
952 double FINVx[5] = {0,0,0,0,0};
953 double FINVy[5] = {0,0,0,0,0};
954
955 FINVMTX(sa,FINVx,FINVy);
956
957 double resx = 0;
958 double c_resx = 0;
959 double resy = 0;
960 double c_resy = 0;
961 double pv[2];
962 for(int j=0; j<5; j++){
963 PVTR(j,sa,pv);
964 double Plx = pv[0];
965 double Ply = pv[1];
966 resx = resx - c_resx;
967 double resx_tmp = resx + (FINVx[j]*Plx - FINVy[j]*Ply);
968 c_resx = (resx_tmp - resx) - (FINVx[j]*Plx - FINVy[j]*Ply);
969 resx = resx_tmp;
970
971 resy = resy - c_resy;
972 double resy_tmp = resy + (FINVx[j]*Ply + FINVy[j]*Plx);
973 c_resy = (resy_tmp - resy) - (FINVx[j]*Ply + FINVy[j]*Plx);
974 resy = resy_tmp;
975 }
976 outputx = resx;
977 outputy = resy;
978 Fv[0] = outputx;
979 Fv[1] = outputy;
980
981}
982void EvtDToKppipi::calEva(double* Kp, double* Pip, double* Pim, double *mass1, double *width1, double *amp, double *phase,int* g0,int* spin, int* modetype, int nstates, double & Result , double*r0 ,double*r1)
983{
984
985 double P12[4], P23[4], P13[4];
986 double cof[2], amp_PDF[2], PDF[2];
987 double Amp_KPiS[2];
988 double Fv[2];
989 double s1, s2, s3;
990 double Realpipis, Imagpipis;
991
992 double s12, s13, s23;
993 for(int i=0; i<4; i++){
994 P12[i] = Kp[i] + Pip[i];
995 P13[i] = Kp[i] + Pim[i];
996 P23[i] = Pip[i] + Pim[i];
997 }
998 s1 = SCADot(Kp,Kp);
999 s2 = SCADot(Pip,Pip);
1000 s3 = SCADot(Pim,Pim);
1001
1002 s12 = SCADot(P12,P12);
1003 s13 = SCADot(P13,P13);
1004 s23 = SCADot(P23,P23);
1005 //cout<<"Spipi = "<<s23<<"\n"<<endl;
1006 double pro[2], temp_PDF, pro1[2], temp_PDF1, pro2[2], temp_PDF2, amp_tmp[2];
1007 double mass1sq;
1008 amp_PDF[0] = 0;
1009 amp_PDF[1] = 0;
1010 PDF[0] = 0;
1011 PDF[1] = 0;
1012 amp_tmp[0] = 0;
1013 amp_tmp[1] = 0;
1014 for(int i=0; i<nstates; i++) {
1015 amp_tmp[0] = 0;
1016 amp_tmp[1] = 0;
1017 mass1sq = mass1[i]*mass1[i];
1018 cof[0] = amp[i]*cos(phase[i]);
1019 cof[1] = amp[i]*sin(phase[i]);
1020 temp_PDF = 0;
1021
1022
1023 if(modetype[i] == 12){
1024 temp_PDF = DDalitz(Kp, Pip, Pim, spin[i], mass1[i]);
1025 if(g0[i]==0){
1026 pro[0] = 1;
1027 pro[1] = 0;
1028 }
1029
1030 amp_tmp[0] = temp_PDF*pro[0];
1031 amp_tmp[1] = temp_PDF*pro[1];
1032 }
1033
1034 if(modetype[i] == 13){
1035 temp_PDF = DDalitz(Kp, Pim, Pip, spin[i], mass1[i]);
1036 if(g0[i]==1) propagatorRBW(mass1[i],width1[i],s13,s1,s3,rRes2,spin[i],pro);
1037 if(g0[i]==4) KPiSLASS(s13,s1,s3,pro);
1038 if(g0[i]==0){
1039 pro[0] = 1;
1040 pro[1] = 0;
1041 }
1042
1043 amp_tmp[0] = temp_PDF*pro[0];
1044 amp_tmp[1] = temp_PDF*pro[1];
1045
1046 }
1047
1048 if(modetype[i] == 23){
1049 temp_PDF = DDalitz(Pip, Pim, Kp, spin[i], mass1[i]);
1050 if(g0[i]==6) {
1051 Fvector(s23,-0.07,Fv);
1052 Realpipis = Fv[0];
1053 Imagpipis = Fv[1];
1054 /*Amp_PiPiS = K_matrix(s23);
1055 Realpipis = real(Amp_PiPiS);
1056 Imagpipis = imag(Amp_PiPiS);*/
1057 //cout<<"Realpipis = "<<Realpipis<<"\n"<<endl;
1058 //cout<<"Imagpipis = "<<Imagpipis<<"\n"<<endl;
1059 //cout<<"spipi = "<<s23<<"\n"<<endl;
1060 amp_tmp[0] = temp_PDF*Realpipis;
1061 amp_tmp[1] = temp_PDF*Imagpipis;
1062 }
1063 if(g0[i]==4) propagatorFlatte(mass1[i],width1[i],s23,pro); // Only for f0(980)
1064 if(g0[i]==1) propagatorGS(mass1[i],width1[i],s23,s2,s3,rRes2,pro);
1065 if(g0[i]==2) propagatorRBW(mass1[i],width1[i],s23,s2,s3,rRes2,spin[i],pro);
1066 if(g0[i]==0){
1067 pro[0] = 1;
1068 pro[1] = 0;
1069 }
1070 if(g0[i]==500) propagatorsigma500(s23, s2, s3, pro);
1071 if(g0[i]!=6) amp_tmp[0] = temp_PDF*pro[0];
1072 if(g0[i]!=6) amp_tmp[1] = temp_PDF*pro[1];
1073 }
1074
1075 if(modetype[i] == 132){
1076 KPiSLASS(s13,s1,s3,Amp_KPiS);
1077 amp_tmp[0] = Amp_KPiS[0];
1078 amp_tmp[1] = Amp_KPiS[1];
1079 }
1080
1081 Com_Multi(amp_tmp,cof,amp_PDF);
1082 PDF[0] += amp_PDF[0];
1083 PDF[1] += amp_PDF[1];
1084
1085
1086 }
1087 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
1088 if(value <=0) value = 1e-20;
1089 Result = value;
1090 }
double sin(const BesAngle a)
Definition BesAngle.h:210
double cos(const BesAngle a)
Definition BesAngle.h:213
double P(RecMdcKalTrack *trk)
double mass
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
TFile * f1
TF1 * g1
int ID[no]
EvtComplex exp(const EvtComplex &c)
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in !Latex Output unit
Definition FoamA.h:90
const double mpi
Definition Gam4pikp.cxx:47
XmlRpcServer s
****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
TCrossPart * CS
Definition Mcgpj.cxx:51
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA ** Rho(770) and Omega(782) are taken from CMD-2 F_pi fit *(hep-ex/9904027)
TTree * t
Definition binning.cxx:23
static const double radToDegrees
Definition EvtConst.hh:30
void decay(EvtParticle *p)
virtual ~EvtDToKppipi()
void getName(std::string &name)
EvtDecayBase * clone()
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
void checkSpinParent(EvtSpinType::spintype sp)
void setProbMax(double prbmx)
void checkNDaug(int d1, int d2=-1)
EvtId * getDaugs()
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
void setProb(double prob)
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
double get(int i) const
double double double double * s12
Definition qcdloop1.h:77
double double double double double * s23
Definition qcdloop1.h:77
double double * m2
Definition qcdloop1.h:75
const double b
Definition slope.cxx:9