BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKSKSK.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: EvtDToKSKSK.cc
11//
12// Description: D+ -> KS0 KS0 K+
13//
14// Modification history:
15//
16// Liaoyuan Dong Wed Nov 9 00:25:20 2022 Module created
17//------------------------------------------------------------------------
18#include "EvtGenBase/EvtPatches.hh"
19#include "EvtGenBase/EvtParticle.hh"
20#include "EvtGenBase/EvtGenKine.hh"
21#include "EvtGenBase/EvtPDL.hh"
22#include "EvtGenBase/EvtReport.hh"
23#include "EvtGenModels/EvtDToKSKSK.hh"
24#include "EvtGenBase/EvtComplex.hh"
25#include "EvtGenBase/EvtDecayTable.hh"
26#include <stdlib.h>
27#include <iostream>
28#include <cmath>
29using namespace std;
30
31EvtDToKSKSK::~EvtDToKSKSK() {}
32
33void EvtDToKSKSK::getName(std::string& model_name){
34 model_name="DToKSKSK";
35}
36
40
41void EvtDToKSKSK::init(){
42 checkNArg(0);
43 checkNDaug(3);
44 checkSpinParent(EvtSpinType::SCALAR);
45
46 mass[0] = 0.980;
47 mass[1] = 0.980;
48 width[0] = 0.0972;
49 width[1] = 0.0972;
50
51 rho[0] = 1.0;
52 rho[1] = 1.1141e+01;
53 phi[0] = 0.00;
54 phi[1] = -5.6220e+00;
55
56// spin[0] = 0;
57// spin[1] = 1;
58// spin[2] = 1;
59 modetype[0] = 12;
60 modetype[1] = 122;
61
62 //std::cout << "Initializing EvtDToKSKSK" << std::endl;
63 //for (int i=0; i<2; i++) {
64 // cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
65 //}
66
67 math_K = 0.493677;
68 GS1 = 0.636619783;
69 GS2 = 0.23459086;
70 GS3 = 0.1591549458;
71 GS4 = 0.078196955;
72
73 mD0 = 1.86966;
74 mK0 = 0.497611;
75 mKa = 0.493677;
76 mPi = 0.13957;
77 mK02 = 0.237616707;
78 mPi2 = 0.01947978;
79 mass_EtaP = 0.95778;
80 mass_Kaon = 0.49368;
81
82 math_pi = 3.1415926;
83 mass_Pion = 0.13957;
84 mass_Pion2 = 0.0194797849;
85 mass_2Pion = 0.27914;
86 math_2pi = 6.2831852;
87 rD2 = 25.0;
88 rRes2 = 9.0;
89 g2 = 0.23;
90
91 rho_omega = 0.00294;
92 phi_omega = -0.02;
93
94}
95
96void EvtDToKSKSK::initProbMax() {
97 setProbMax(878.9);
98}
99
100void EvtDToKSKSK::decay( EvtParticle *p){
101
102/* double maxprob = 0.0;
103 for(int ir=0;ir<=60000000;ir++){
104 p->initializePhaseSpace(getNDaug(),getDaugs());
105 EvtVector4R D1 = p->getDaug(0)->getP4();
106 EvtVector4R D2 = p->getDaug(1)->getP4();
107 EvtVector4R D3 = p->getDaug(2)->getP4();
108
109 double P1[4], P2[4], P3[4];
110 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
111 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
112 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
113
114 double value;
115 int g0[3]={1,1,1};
116 int spin[3]={0,1,0};
117 int nstates=2;
118 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
119 if (value<0) continue;
120 if(value>maxprob) {
121 maxprob=value;
122 cout << "ir= " << ir << endl;
123 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3] <<"};"<< endl;
124 cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3] <<"};"<< endl;
125 cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< endl;
126 cout << "MAX====> " << maxprob << endl;
127 }
128 }
129 printf("MAXprob = %.10f\n",maxprob);
130
131*/
132 p->initializePhaseSpace(getNDaug(),getDaugs());
133 EvtVector4R D1 = p->getDaug(0)->getP4();
134 EvtVector4R D2 = p->getDaug(1)->getP4();
135 EvtVector4R D3 = p->getDaug(2)->getP4();
136
137 double P1[4], P2[4], P3[4];
138 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
139 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
140 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
141
142 double value;
143 int g0[2]={1,1};
144 int spin[2]={0,0};
145 int nstates=2;
146 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
147 setProb(value);
148
149 return ;
150}
151
152void EvtDToKSKSK::Com_Multi(double a1[2], double a2[2], double res[2])
153{
154 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
155 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
156}
157void EvtDToKSKSK::Com_Divide(double a1[2], double a2[2], double res[2])
158{
159 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
160 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
161 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
162}
163//------------base---------------------------------
164double EvtDToKSKSK::SCADot(double a1[4], double a2[4])
165{
166 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
167 return _cal;
168}
169double EvtDToKSKSK::barrier(int l, double sa, double sb, double sc, double r, double mass)
170{
171 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
172 // if(q < 0) q = 1e-16;
173 double z;
174 z = q*r*r;
175 double sa0;
176 sa0 = mass*mass;
177 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
178 // if(q0 < 0) q0 = 1e-16;
179 double z0 = q0*r*r;
180 double F = 0.0;
181 if(l == 0) F = 1;
182 if(l == 1) F = sqrt((1+z0)/(1+z));
183 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
184 return F;
185}
186
187void EvtDToKSKSK::calt1(double daug1[4], double daug2[4], double t1[4])
188{
189 double p, pq, tmp;
190 double pa[4], qa[4];
191 for(int i=0; i<4; i++) {
192 pa[i] = daug1[i] + daug2[i];
193 qa[i] = daug1[i] - daug2[i];
194 }
195 p = SCADot(pa,pa);
196 pq = SCADot(pa,qa);
197 tmp = pq/p;
198 for(int i=0; i<4; i++) {
199 t1[i] = qa[i] - tmp*pa[i];
200 }
201}
202
203void EvtDToKSKSK::calt2(double daug1[4], double daug2[4], double t2[4][4])
204{
205 double p, r;
206 double pa[4], t1[4];
207 calt1(daug1,daug2,t1);
208 r = SCADot(t1,t1)/3.0;
209 for(int i=0; i<4; i++) {
210 pa[i] = daug1[i] + daug2[i];
211 }
212 p = SCADot(pa,pa);
213 for(int i=0; i<4; i++) {
214 for(int j=0; j<4; j++) {
215 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
216 }
217 }
218}
219
220 void EvtDToKSKSK::propagatorCBW(double mass, double width, double sx, double prop[2])
221{
222 double a[2], b[2];
223 a[0] = 1;
224 a[1] = 0;
225 b[0] = mass*mass-sx;
226 b[1] = -mass*width;
227 Com_Divide(a,b,prop);
228}
229
230
231double EvtDToKSKSK::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
232{
233 double widm = 0.;
234 double m = sqrt(sa);
235 double tmp = sb-sc;
236 double tmp1 = sa+tmp;
237 double q = fabs(0.25*tmp1*tmp1/sa-sb);
238 double tmp2 = mass2+tmp;
239 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
240 double z = q*r2;
241 double z0 = q0*r2;
242 double t = q/q0;
243
244 if(l == 0) {widm = sqrt(t)*mass/m;}
245 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
246 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
247 return widm;
248}
249double EvtDToKSKSK::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
250{
251 double widm = 0.;
252 double m = sqrt(sa);
253 double tmp = sb-sc;
254 double tmp1 = sa+tmp;
255 double q = fabs(0.25*tmp1*tmp1/sa-sb);
256 double tmp2 = mass2+tmp;
257 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
258 double z = q*r2;
259 double z0 = q0*r2;
260 double F = (1+z0)/(1+z);
261 double t = q/q0;
262 widm = t*sqrt(t)*mass/m*F;
263 return widm;
264
265}
266 void EvtDToKSKSK::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
267{
268 double a[2], b[2];
269 double mass2 = mass*mass;
270
271 a[0] = 1;
272 a[1] = 0;
273 b[0] = mass2-sa;
274 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
275 Com_Divide(a,b,prop);
276}
277
278void EvtDToKSKSK::propagatorFlatte(double mass, double width, double sa, double sb, double sc, double prop[2]){
279 double rhoab, rhoKK;
280 rhoab=1.0/sa*sqrt((sa-(0.547862+0.139570)*(0.547862+0.139570))*(sa-(0.547862-0.139570)*(0.547862-0.139570)));
281 rhoKK=1.0/sa*sqrt((sa-(0.497611+0.493677)*(0.497611+0.493677))*(sa-(0.497611-0.493677)*(0.497611-0.493677)));
282 double a[2], b[2];
283 a[0] = 1;
284 a[1] = 0;
285 b[0] = mass*mass - sa;
286 b[1] = - (0.324*0.324*rhoab+ 1.03*0.324*0.324*rhoKK);
287 Com_Divide(a,b,prop);
288}
289
290 void EvtDToKSKSK::propagatorGS(double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
291{
292 double a[2], b[2];
293 double mass2 = mass*mass;
294 double tmp = sb-sc;
295 double tmp1 = sa+tmp;
296 double q2 = fabs(0.25*tmp1*tmp1/sa-sb);
297 double tmp2 = mass2+tmp;
298 double q02 = fabs(0.25*tmp2*tmp2/mass2-sb);
299 double q = sqrt(q2);
300 double q0 = sqrt(q02);
301 double m = sqrt(sa);
302 double q03 = q0*q02;
303 double tmp3 = log(mass+2*q0)+0.0087501713;
304 double h = GS1*q/m*(log(m+2*q)+0.0087501713);
305 double h0 = GS1*q0/mass*tmp3;
306 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
307 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
308 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
309 a[0] = 1.0+d*width/mass;
310 a[1] = 0.0;
311 b[0] = mass2-sa+width*f;
312 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
313 Com_Divide(a,b,prop);
314}
315
316 void EvtDToKSKSK::PiPiSWASS(double sa, double sb, double sc, double prop[2]) {
317 double tmp = sb-sc;
318 double tmp2 = sa+tmp;
319 double qs = 0.25*tmp2*tmp2/sa-sb;
320 double q = sqrt(qs);
321 double a0 = -0.11/mass_Pion;
322 prop[0] = 1/(1+a0*a0*q*q);
323 prop[1] = a0*q/(1+a0*a0*q*q);
324}
325
326 void EvtDToKSKSK::KPiSLASS(double sa, double sb, double sc, double prop[2]) {
327 const double m1430 = 1.441;
328 const double sa0 = 2.076481;
329 const double w1430 = 0.193;
330 const double Lass1 = 0.25/sa0;
331 double tmp = sb-sc;
332 double tmp1 = sa0+tmp;
333 double q0 = fabs(Lass1*tmp1*tmp1-sb);
334 double tmp2 = sa+tmp;
335 double qs = 0.25*tmp2*tmp2/sa-sb;
336 double q = sqrt(qs);
337 double width = w1430*q*m1430/sqrt(sa*q0);
338 double temp_R = atan(m1430*width/(sa0-sa));
339 if(temp_R<0) temp_R += math_K;
340 double deltaR = -109.7*math_K/180.0 + temp_R;
341 double temp_F = atan(0.226*q/(2.0-3.8194*qs));
342 if(temp_F<0) temp_F += math_K;
343 double deltaF = 0.1*math_K/180.0 + temp_F;
344 double deltaS = deltaR + 2.0*deltaF;
345 double t1 = 0.96*sin(deltaF);
346 double t2 = sin(deltaR);
347 double CF[2], CS[2];
348 CF[0] = cos(deltaF);
349 CF[1] = sin(deltaF);
350 CS[0] = cos(deltaS);
351 CS[1] = sin(deltaS);
352 prop[0] = t1*CF[0] + t2*CS[0];
353 prop[1] = t1*CF[1] + t2*CS[1];
354}
355
356 void EvtDToKSKSK::Flatte_rhoab(double sa, double sb, double sc, double rho[2]){
357 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
358 if(q>0) {
359 rho[0]=2* sqrt(q/sa);
360 rho[1]=0;
361 }
362 else if(q<0){
363 rho[0]=0;
364 rho[1]=2*sqrt(-q/sa);
365 }
366}
367
368 void EvtDToKSKSK::propagatorKstr1430(double mass, double sx, double *sb, double *sc, double prop[2])
369{
370 double unit[2]={1.0};
371 double ci[2]={0,1};
372 double rho1[2];
373 Flatte_rhoab(sx,sb[0],sc[0],rho1);
374 double rho2[2];
375 Flatte_rhoab(sx,sb[1],sc[1],rho2);
376 double gKPi_Kstr1430=0.2990, gEtaPK_Kstr1430=0.0529;
377 double tmp1[2]={gKPi_Kstr1430,0};
378 double tmp11[2];
379 double tmp2[2]={gEtaPK_Kstr1430,0};
380 double tmp22[2];
381 Com_Multi(tmp1,rho1,tmp11);
382 Com_Multi(tmp2,rho2,tmp22);
383 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
384 double tmp31[2];
385 Com_Multi(tmp3, ci,tmp31);
386 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
387 Com_Divide( unit,tmp4, prop);
388}
389 void EvtDToKSKSK::rhoab(double sa, double sb, double sc, double res[2]) {
390 double tmp = sa+sb-sc;
391 double q = 0.25*tmp*tmp/sa-sb;
392 if(q>=0) {
393 res[0]=2.0*sqrt(q/sa);
394 res[1]=0.0;
395 } else {
396 res[0]=0.0;
397 res[1]=2.0*sqrt(-q/sa);
398 }
399}
400 void EvtDToKSKSK::rho4Pi(double sa, double res[2]) {
401 double temp = 1.0-0.3116765584/sa;
402 if(temp>=0) {
403 res[0]=sqrt(temp)/(1.0+exp(9.8-3.5*sa));
404 res[1]=0.0;
405 } else {
406 res[0]=0.0;
407 res[1]=sqrt(-temp)/(1.0+exp(9.8-3.5*sa));
408 }
409}
410
411 void EvtDToKSKSK::propagatorsigma500(double sa, double sb, double sc, double prop[2]) {
412 double f = 0.5843+1.6663*sa;
413 const double M = 0.9264;
414 const double mass2 = 0.85821696;
415 const double mpi2d2 = 0.00973989245;
416 double g1 = f*(sa-mpi2d2)/(mass2-mpi2d2)*exp((mass2-sa)/1.082);
417 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
418 rhoab(sa,sb,sc,rho1s);
419 rhoab(mass2,sb,sc,rho1M);
420 rho4Pi(sa,rho2s);
421 rho4Pi(mass2,rho2M);
422 Com_Divide(rho1s,rho1M,rho1);
423 Com_Divide(rho2s,rho2M,rho2);
424 double a[2], b[2];
425 a[0] = 1.0;
426 a[1] = 0.0;
427 b[0] = mass2-sa+M*(g1*rho1[1]+0.0024*rho2[1]);
428 b[1] = -M*(g1*rho1[0]+0.0024*rho2[0]);
429 Com_Divide(a,b,prop);
430}
431
432 void EvtDToKSKSK::getprop(double sa, double sb, double sc, double mass, double width, double prop[2]){
433 double prop1[2], prop2[2];
434 propagatorGS(mass,width,sa,sb,sc,9.0,prop1);
435 propagatorRBW(0.783,0.008,sa,sb,sc,3.0,1,prop2);
436 double coef_omega[2];
437 coef_omega[0] = rho_omega*cos(phi_omega),
438 coef_omega[1] = rho_omega*sin(phi_omega);
439 double one[2]; one[0] = 1; one[1] = 0;
440 double temp[2];
441 Com_Multi(coef_omega,prop2,temp);
442 temp[0] = one[0] + 0.783*0.783*temp[0];
443 temp[1] = one[1] + 0.783*0.783*temp[1];
444 Com_Multi(prop1,temp,prop);
445}
446
447double EvtDToKSKSK::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
448 double pR[4], pD[4];
449 double temp_PDF, v_re;
450 temp_PDF = 0.0;
451 v_re = 0.0;
452 double B[2], s1, s2, s3, sR, sD;
453 for(int i=0; i<4; i++){
454 pR[i] = P1[i] + P2[i];
455 pD[i] = pR[i] + P3[i];
456 }
457 s1 = SCADot(P1,P1);
458 s2 = SCADot(P2,P2);
459 s3 = SCADot(P3,P3);
460 sR = SCADot(pR,pR);
461 sD = SCADot(pD,pD);
462 int G[4][4];
463 for(int i=0; i!=4; i++){
464 for(int j=0; j!=4; j++){
465 if(i==j){
466 if(i==0) G[i][j] = 1;
467 else G[i][j] = -1;
468 }
469 else G[i][j] = 0;
470 }
471 }
472 if(Ang == 0){
473 B[0] = 1;
474 B[1] = 1;
475 temp_PDF = 1;
476 }
477 if(Ang == 1){
478 B[0] = barrier(1,sR,s1,s2,3.0,mass);
479 B[1] = barrier(1,sD,sR,s3,5.0,mD0);
480 double t1[4], T1[4];
481 calt1(P1,P2,t1);
482 calt1(pR,P3,T1);
483 temp_PDF = 0;
484 for(int i=0; i<4; i++){
485 temp_PDF += t1[i]*T1[i]*G[i][i];
486 }
487 }
488 if(Ang == 2){
489 B[0] = barrier(2,sR,s1,s2,3.0,mass);
490 B[1] = barrier(2,sD,sR,s3,5.0,mD0);
491 double t2[4][4], T2[4][4];
492 calt2(P1,P2,t2);
493 calt2(pR,P3,T2);
494 temp_PDF = 0;
495 for(int i=0; i<4; i++){
496 for(int j=0; j<4; j++){
497 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
498 }
499 }
500 }
501 v_re = temp_PDF*B[0]*B[1];
502 return v_re;
503}
504
505
506
507void EvtDToKSKSK::calEva(double* Ks01, double* Ks02, double* Kp, double *mass1, double *width1, double *amp, double *phase,int* g0,int* spin, int* modetype, int nstates, double & Result)
508{
509 double numEvents;
510 double P23[4], P13[4], P12[4];
511// double Ks01[4], Ks02[4], Kp[4], P23[4], P13[4], P12[4];
512 double cof[2], amp_PDF[2], PDF[2];
513 double s12, s23, s13;
514 for(int i=0; i<4; i++){
515 P12[i] = Ks01[i] + Ks02[i];
516 P13[i] = Ks01[i] + Kp[i];
517 P23[i] = Ks02[i]+ Kp[i];
518 }
519 s12 = SCADot(P12,P12);
520 s13 = SCADot(P13,P13);
521 s23 = SCADot(P23,P23);
522 double s1,s2,s3;
523 s1 = SCADot(Ks01,Ks01);
524 s2 = SCADot(Ks02,Ks02);
525 s3 = SCADot(Kp,Kp);
526 double pro[2], temp_PDF, amp_tmp[2];
527 double pro1[2], temp_PDF1,pro2[2], temp_PDF2;
528 double pro3[2], temp_PDF3,pro4[2], temp_PDF4;
529 double pro5[2], temp_PDF5,pro6[2], temp_PDF6;
530 double pro7[2], temp_PDF7,pro8[2], temp_PDF8;
531 double pro9[2], temp_PDF9,pro10[2], temp_PDF10;
532 double Amp_KPiS[2];
533 amp_PDF[0] = 0;
534 amp_PDF[1] = 0;
535 PDF[0] = 0;
536 PDF[1] = 0;
537 amp_tmp[0] = 0;
538 amp_tmp[1] = 0;
539 for(int i=0; i<3; i++) {
540 amp_tmp[0] = 0;
541 amp_tmp[1] = 0;
542 cof[0] = amp[i]*cos(phase[i]);
543 cof[1] = amp[i]*sin(phase[i]);
544 temp_PDF = 0;
545 if(modetype[i] == 12){
546 temp_PDF1 = DDalitz(Ks01, Kp, Ks02, spin[i], mass1[i]);
547 if(g0[i]==1) propagatorFlatte(mass1[i],width1[i],s13,s1,s3,pro1);
548
549 if(g0[i]==2) {
550 double skm2[2]={s1, mass_EtaP *mass_EtaP};
551 double spi2[2]={s2, mass_Kaon *mass_Kaon};
552 propagatorKstr1430(mass1[i],s12,skm2,spi2,pro1);
553 }
554 if(g0[i]==4) KPiSLASS(s12,s1,s2,pro1);
555 if(g0[i]==0){
556 pro1[0] = 1;
557 pro1[1] = 0;
558 }
559
560 temp_PDF2 = DDalitz(Ks02, Kp, Ks01, spin[i], mass1[i]);
561 if(g0[i]==1) propagatorFlatte(mass1[i],width1[i],s23,s2,s3,pro2);
562
563 if(g0[i]==2) {
564 double skm2[2]={s1, mass_EtaP *mass_EtaP};
565 double spi2[2]={s3, mass_Kaon *mass_Kaon};
566 propagatorKstr1430(mass1[i],s13,skm2,spi2,pro2);
567 }
568 if(g0[i]==4) KPiSLASS(s13,s1,s3,pro2);
569 if(g0[i]==0){
570 pro2[0] = 1;
571 pro2[1] = 0;
572 }
573 amp_tmp[0] = temp_PDF1*pro1[0] + temp_PDF2*pro2[0];
574 amp_tmp[1] = temp_PDF1*pro1[1] + temp_PDF2*pro2[1];
575 }
576
577
578 if(modetype[i] == 123){
579 temp_PDF3 = DDalitz(Ks01, Kp, Ks02, spin[i], mass1[i]);
580 if(g0[i]==1) propagatorGS(mass1[i],width1[i],s13,s1,s3,9.0,pro3);
581 if(g0[i]==2) {
582 double skm2[2]={s1, mass_EtaP *mass_EtaP};
583 double spi2[2]={s2, mass_Kaon *mass_Kaon};
584 propagatorKstr1430(mass1[i],s12,skm2,spi2,pro3);
585 }
586 if(g0[i]==4) KPiSLASS(s12,s1,s2,pro3);
587 if(g0[i]==0){
588 pro1[0] = 1;
589 pro1[1] = 0;
590 }
591
592 temp_PDF4 = DDalitz(Ks02, Kp, Ks01, spin[i], mass1[i]);
593 if(g0[i]==1) propagatorGS(mass1[i],width1[i],s23,s2,s3,9.0,pro4);
594
595 if(g0[i]==2) {
596 double skm2[2]={s1, mass_EtaP *mass_EtaP};
597 double spi2[2]={s3, mass_Kaon *mass_Kaon};
598 propagatorKstr1430(mass1[i],s13,skm2,spi2,pro2);
599 }
600 if(g0[i]==4) KPiSLASS(s13,s1,s3,pro2);
601 if(g0[i]==0){
602 pro2[0] = 1;
603 pro2[1] = 0;
604 }
605 amp_tmp[0] = temp_PDF3*pro3[0] + temp_PDF4*pro4[0];
606 amp_tmp[1] = temp_PDF3*pro3[1] + temp_PDF4*pro4[1];
607 }
608
609 if(modetype[i] == 122){
610 amp_tmp[0] =2;
611 amp_tmp[1] =0;
612 }
613// cout<<"pdf: "<<i<<", "<<amp_tmp[0]<<", "<<amp_tmp[1]<<endl;
614 Com_Multi(amp_tmp,cof,amp_PDF);
615 PDF[0] += amp_PDF[0];
616 PDF[1] += amp_PDF[1];
617
618// cout<<"PDF: "<<i<<", "<<amp_PDF[0]<<", "<<amp_PDF[1]<<endl;
619 }
620 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
621// cout<<"value: "<<value<<endl;
622 if(value <=0) value = 1e-20;
623
624 Result = value;
625}
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
mass
double mass
Definition CosmicGenerator.cxx:138
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
g1
TF1 * g1
Definition DataBase/tau_mode.c:66
a0
character *LEPTONflag integer iresonances real zeta5 real a0
Definition Ekhara/Ekhara-00-01-02/src/ekhara/common.h:4
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:252
unit
*******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
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
CS
TCrossPart * CS
Definition Mcgpj.cxx:51
t
TTree * t
Definition binning.cxx:23
EvtDToKSKSK::initProbMax
void initProbMax()
Definition EvtDToKSKSK.cc:96
EvtDToKSKSK::getName
void getName(std::string &name)
Definition EvtDToKSKSK.cc:33
EvtDToKSKSK::~EvtDToKSKSK
virtual ~EvtDToKSKSK()
Definition EvtDToKSKSK.cc:31
EvtDToKSKSK::clone
EvtDecayBase * clone()
Definition EvtDToKSKSK.cc:37
EvtDToKSKSK::decay
void decay(EvtParticle *p)
Definition EvtDToKSKSK.cc:100
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
s12
double double double double * s12
Definition qcdloop1.h:77
s23
double double double double double * s23
Definition qcdloop1.h:77
one
double precision pisqo6 one
Definition qlconstants.h:4
b
const double b
Definition slope.cxx:9