BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKpipi.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of the EvtGen package developed jointly
4// for the BaBar and CLEO collaborations. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/COPYRIGHT
8// Copyright (C) 1998 Caltech, UCSB
9//
10// Module: EvtDsToKpipi.cc
11//
12// Description: Routine to handle three-body decays of Ds+/Ds-
13//
14// Modification history:
15// Liaoyuan Dong Sun Jun 02 01:23:25 2024 Module created
16//------------------------------------------------------------------------
17#include "EvtGenBase/EvtPatches.hh"
18#include <fstream>
19#include <stdlib.h>
20#include "EvtGenBase/EvtParticle.hh"
21#include "EvtGenBase/EvtGenKine.hh"
22#include "EvtGenBase/EvtPDL.hh"
23#include "EvtGenBase/EvtReport.hh"
24#include "EvtGenBase/EvtResonance.hh"
25#include "EvtGenBase/EvtResonance2.hh"
26#include "EvtGenModels/EvtDsToKpipi.hh"
27#include <string>
28#include "EvtGenBase/EvtConst.hh"
29#include "EvtGenBase/EvtDecayTable.hh"
30using std::endl;
31
32EvtDsToKpipi::~EvtDsToKpipi() {}
33
34void EvtDsToKpipi::getName(std::string& model_name){
35 model_name="DsToKpipi";
36}
37
41
42void EvtDsToKpipi::init(){
43 // check that there are 0 arguments
44 checkNArg(0);
45 checkNDaug(3);
46 checkSpinParent(EvtSpinType::SCALAR);
47 checkSpinDaughter(0,EvtSpinType::SCALAR);
48 checkSpinDaughter(1,EvtSpinType::SCALAR);
49 checkSpinDaughter(2,EvtSpinType::SCALAR);
50
51 phi[0] = 0;
52 rho[0] = 1;
53 phi[1] = 0;
54 rho[1] = 0;
55 phi[2] = 0;
56 rho[2] = 0;
57 phi[3] = 0;
58 rho[3] = 0;
59 phi[4] = 0;
60 rho[4] = 0;
61 phi[5] = 0;
62 rho[5] = 0;
63 phi[6] = 0;
64 rho[6] = 0;
65 phi[7] = 0;
66 rho[7] = 0;
67
68 phi[1] = -3.47995752;
69 phi[2] = -1.249864467;
70 phi[3] = -0.3067504308;
71 phi[4] = 0.9229242379;
72 phi[5] = -3.278567926;
73 phi[6] = -0.6346408622;
74 phi[7] = 1.762377475;
75
76 rho[1] = 2.463898984;
77 rho[2] = 0.7361782665;
78 rho[3] = 1.90216812;
79 rho[4] = 2.140687169;
80 rho[5] = 0.914684056;
81 rho[6] = 1.116206881;
82 rho[7] = 1.483440555;
83
84 modetype[0]= 23;
85 modetype[1]= 23;
86 modetype[2]= 23;
87 modetype[3]= 23;
88 modetype[4]= 23;
89 modetype[5]= 13;
90 modetype[6]= 13;
91 modetype[7]= 13;
92
93 //cout << "DsToKpipi :" << endl;
94 //for (int i=0; i<8; i++) {
95 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
96 //}
97
98 width[0] = 0.1478;
99 width[1] = 0.400;
100 width[2] = 0.100;
101 width[3] = 0.265;
102 width[4] = 0.270;
103 width[5] = 0.0473;
104 width[6] = 0.232;
105 width[7] = 0.270;
106
107 mass[0] = 0.77526;
108 mass[1] = 1.465;
109 mass[2] = 0.965;
110 mass[3] = 1.35;
111 mass[4] = 1.425;
112 mass[5] = 0.89555;
113 mass[6] = 1.414;
114 mass[7] = 1.432787726;
115
116 mDsM = 1.9683;
117 mD = 1.86486;
118 mDs = 1.9683;
119// rRes = 9.0;
120 rD = 25.0;
121 metap = 0.95778;
122 mkstr = 0.89594;
123 mk0 = 0.497614;
124 mass_Kaon = 0.49368;
125 mass_Pion = 0.13957;
126 mass_Pi0 = 0.1349766;
127 math_pi = 3.1415926;
128 mKa2 = 0.24371994;
129 mPi2 = 0.01947978;
130 mPi = 0.13957;
131 mKa = 0.49368;
132 //mrho = 0.77549;
133 //Grho = 0.1491;
134 ma0 = 0.99;
135 Ga0 = 0.0756;
136 meta = 0.547862;
137
138 GS1 = 0.636619783;
139 GS2 = 0.01860182466;
140 GS3 = 0.1591549458; // 1/(2*math_2pi)
141 GS4 = 0.00620060822; // mass_Pion2/math_pi
142
143 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
144 for (int i=0; i<4; i++) {
145 for (int j=0; j<4; j++) {
146 G[i][j] = GG[i][j];
147 }
148 }
149}
150
151void EvtDsToKpipi::initProbMax() {
152 setProbMax(825.0);
153}
154
155void EvtDsToKpipi::decay( EvtParticle *p){
156/*
157 double maxprob = 0.0;
158 for(int ir=0;ir<=60000000;ir++){
159 p->initializePhaseSpace(getNDaug(),getDaugs());
160 EvtVector4R D1 = p->getDaug(0)->getP4();
161 EvtVector4R D2 = p->getDaug(1)->getP4();
162 EvtVector4R D3 = p->getDaug(2)->getP4();
163
164 double P1[4], P2[4], P3[4];
165 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
166 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
167 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
168
169 double value;
170 int g0[8]={1,1,4,2,500,1,1,2};
171 int spin[8]={1,1,0,0,0,1,1,0};
172 int nstates=8;
173 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
174 if (value<0) continue;
175 if(value>maxprob) {
176 maxprob=value;
177 cout << "ir= " << ir << endl;
178 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3] <<"};"<< endl;
179 cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3] <<"};"<< endl;
180 cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< endl;
181 cout << "MAX====> " << maxprob << endl;
182 }
183 }
184 printf("MAXprob = %.10f\n",maxprob);
185*/
186
187 p->initializePhaseSpace(getNDaug(),getDaugs());
188 EvtVector4R D1 = p->getDaug(0)->getP4();
189 EvtVector4R D2 = p->getDaug(1)->getP4();
190 EvtVector4R D3 = p->getDaug(2)->getP4();
191
192 double P1[4], P2[4], P3[4];
193 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
194 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
195 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
196
197 //P1[0] =0.863405 ;P1[1] =-0.239724 ;P1[2] =-0.665631 ; P1[3] =-0.0349068 ;
198 //P2[0] =0.209097 ;P2[1] =0.0992859 ;P2[2] =-0.119905 ; P2[3] =-0.00262863 ;
199 //P3[0] =0.941396 ;P3[1] =-0.270638 ;P3[2] =0.890782 ; P3[3] = -0.00293441 ;
200
201 double value;
202 int g0[8]={1,1,4,2,500,1,1,2};
203 int spin[8]={1,1,0,0,0,1,1,0};
204 int nstates=8;
205 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
206
207 setProb(value);
208
209 return ;
210}
211
212void EvtDsToKpipi::Com_Multi(double a1[2], double a2[2], double res[2])
213{
214 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
215 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
216}
217void EvtDsToKpipi::Com_Divide(double a1[2], double a2[2], double res[2])
218{
219 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
220 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
221 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
222}
223
224double EvtDsToKpipi::SCADot(double a1[4], double a2[4])
225{
226 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
227 return _cal;
228}
229double EvtDsToKpipi::barrier(int l, double sa, double sb, double sc, double r, double mass)
230{
231 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
232 if(q < 0) q = 1e-16;
233 double z;
234 z = q*r*r;
235 double sa0;
236 sa0 = mass*mass;
237 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
238 if(q0 < 0) q0 = 1e-16;
239 double z0 = q0*r*r;
240 double F = 0.0;
241 if(l == 0) F = 1;
242 if(l == 1) F = sqrt((1+z0)/(1+z));
243 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
244 return F;
245}
246void EvtDsToKpipi::calt1(double daug1[4], double daug2[4], double t1[4])
247{
248 double p, pq, tmp;
249 double pa[4], qa[4];
250 for(int i=0; i<4; i++) {
251 pa[i] = daug1[i] + daug2[i];
252 qa[i] = daug1[i] - daug2[i];
253 }
254 p = SCADot(pa,pa);
255 pq = SCADot(pa,qa);
256 tmp = pq/p;
257 for(int i=0; i<4; i++) {
258 t1[i] = qa[i] - tmp*pa[i];
259 }
260}
261void EvtDsToKpipi::calt2(double daug1[4], double daug2[4], double t2[4][4])
262{
263 double p, r;
264 double pa[4], t1[4];
265 calt1(daug1,daug2,t1);
266 r = SCADot(t1,t1)/3.0;
267 for(int i=0; i<4; i++) {
268 pa[i] = daug1[i] + daug2[i];
269 }
270 p = SCADot(pa,pa);
271 for(int i=0; i<4; i++) {
272 for(int j=0; j<4; j++) {
273 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
274 }
275 }
276}
277
278//-------------------prop--------------------------------------------
279
280void EvtDsToKpipi::propagatorCBW(double mass, double width, double sx, double prop[2])
281{
282 double a[2], b[2];
283 a[0] = 1;
284 a[1] = 0;
285 b[0] = mass*mass-sx;
286 b[1] = -mass*width;
287 Com_Divide(a,b,prop);
288}
289double EvtDsToKpipi::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
290{
291 double widm = 0.;
292 double m = sqrt(sa);
293 double tmp = sb-sc;
294 double tmp1 = sa+tmp;
295 double q = 0.25*tmp1*tmp1/sa-sb;
296 if(q<0) q = 1e-16;
297 double tmp2 = mass2+tmp;
298 double q0 = 0.25*tmp2*tmp2/mass2-sb;
299 if(q0<0) q0 = 1e-16;
300 double z = q*r2;
301 double z0 = q0*r2;
302 double t = q/q0;
303 if(l == 0) {widm = sqrt(t)*mass/m;}
304 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
305 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
306 return widm;
307}
308
309double EvtDsToKpipi::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
310{
311 double widm = 0.;
312 double m = sqrt(sa);
313 double tmp = sb-sc;
314 double tmp1 = sa+tmp;
315 double q = 0.25*tmp1*tmp1/sa-sb;
316 if(q<0) q = 1e-16;
317 double tmp2 = mass2+tmp;
318 double q0 = 0.25*tmp2*tmp2/mass2-sb;
319 if(q0<0) q0 = 1e-16;
320 double z = q*r2;
321 double z0 = q0*r2;
322 double F = (1+z0)/(1+z);
323 double t = q/q0;
324 widm = t*sqrt(t)*mass/m*F;
325 return widm;
326}
327void EvtDsToKpipi::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
328{
329 double a[2], b[2];
330 double mass2 = mass*mass;
331
332 a[0] = 1;
333 a[1] = 0;
334 b[0] = mass2-sa;
335 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
336 Com_Divide(a,b,prop);
337}
338void EvtDsToKpipi::propagatorFlatte(double mass, double width, double sa, double prop[2]){
339
340 double q2_Pi, q2_Ka;
341 double rhoPi[2], rhoKa[2];
342
343 q2_Pi = 0.25*sa-mPi*mPi;
344 q2_Ka = 0.25*sa-mKa*mKa;
345
346 if (q2_Pi > 0) {
347 rhoPi[0] = 2.0*sqrt(q2_Pi/sa);
348 rhoPi[1] = 0.0;
349 }
350 if (q2_Pi <= 0) {
351 rhoPi[0] = 0.0;
352 rhoPi[1] = 2.0*sqrt(-q2_Pi/sa);
353 }
354
355 if (q2_Ka > 0) {
356 rhoKa[0] = 2.0*sqrt(q2_Ka/sa);
357 rhoKa[1] = 0.0;
358 }
359 if (q2_Ka <= 0) {
360 rhoKa[0] = 0.0;
361 rhoKa[1] = 2.0*sqrt(-q2_Ka/sa);
362 }
363
364double a[2], b[2];
365 a[0] = 1;
366 a[1] = 0;
367 b[0] = mass*mass - sa + 0.165*rhoPi[1] + 0.69465*rhoKa[1];
368 b[1] = - (0.165*rhoPi[0] + 0.69465*rhoKa[0]);
369 Com_Divide(a,b,prop);
370}
371void EvtDsToKpipi::propagatorGS(double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
372{
373 double a[2], b[2];
374 double mass2 = mass*mass;
375 double tmp = sb-sc;
376 double tmp1 = sa+tmp;
377 double q2 = 0.25*tmp1*tmp1/sa-sb;
378 if(q2<0) q2 = 1e-16;
379
380 double tmp2 = mass2+tmp;
381 double q02 = 0.25*tmp2*tmp2/mass2-sb;
382 if(q02<0) q02 = 1e-16;
383
384 double q = sqrt(q2);
385 double q0 = sqrt(q02);
386 double m = sqrt(sa);
387 double q03 = q0*q02;
388 double tmp3 = log(mass+2*q0)+1.2760418309; // log(mass_2Pion) = 1.2760418309;
389
390 double h = GS1*q/m*(log(m+2*q)+1.2760418309);
391 double h0 = GS1*q0/mass*tmp3;
392 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
393 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
394 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
395
396 a[0] = 1.0+d*width/mass;
397 a[1] = 0.0;
398 b[0] = mass2-sa+width*f;
399 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
400 Com_Divide(a,b,prop);
401}
402
403void EvtDsToKpipi::Flatte_rhoab(double sa, double sb, double sc, double rho[2]){
404 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
405 if(q>0) {
406 rho[0]=2* sqrt(q/sa);
407 rho[1]=0;
408 }
409 else if(q<0){
410 rho[0]=0;
411 rho[1]=2*sqrt(-q/sa);
412 }
413}
414
415
416void EvtDsToKpipi::propagatorKstr1430(double mass, double sx, double *sb, double *sc, double prop[2]) //K*1430 Flatte
417{
418 double unit[2]={1.0};
419 double ci[2]={0,1};
420 double rho1[2];
421 Flatte_rhoab(sx,sb[0],sc[0],rho1);
422 double rho2[2];
423 Flatte_rhoab(sx,sb[1],sc[1],rho2);
424 double gKPi_Kstr1430=0.2990, gEtaPK_Kstr1430=0.0529;
425 double tmp1[2]={gKPi_Kstr1430,0};
426 double tmp11[2];
427 double tmp2[2]={gEtaPK_Kstr1430,0};
428 double tmp22[2];
429 Com_Multi(tmp1,rho1,tmp11);
430 Com_Multi(tmp2,rho2,tmp22);
431 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
432 double tmp31[2];
433 Com_Multi(tmp3, ci,tmp31);
434 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
435 Com_Divide( unit,tmp4, prop);
436}
437
438double EvtDsToKpipi::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
439 double pR[4], pD[4];
440 double temp_PDF, v_re;
441 temp_PDF = 0.0;
442 v_re = 0.0;
443 double B[2], s1, s2, s3, sR, sD;
444 for(int i=0; i<4; i++){
445 pR[i] = P1[i] + P2[i];
446 pD[i] = pR[i] + P3[i];
447 }
448 s1 = SCADot(P1,P1);
449 s2 = SCADot(P2,P2);
450 s3 = SCADot(P3,P3);
451 sR = SCADot(pR,pR);
452 sD = SCADot(pD,pD);
453 int G[4][4];
454 for(int i=0; i!=4; i++){
455 for(int j=0; j!=4; j++){
456 if(i==j){
457 if(i==0) G[i][j] = 1;
458 else G[i][j] = -1;
459 }
460 else G[i][j] = 0;
461 }
462 }
463 if(Ang == 0){
464 B[0] = 1;
465 B[1] = 1;
466 temp_PDF = 1;
467 }
468 if(Ang == 1){
469 B[0] = barrier(1,sR,s1,s2,3.0,mass);
470 B[1] = barrier(1,sD,sR,s3,5.0,mDsM);
471 //B[0] = Barrier(1,sR,s1,s2,9.0);
472 //B[1] = Barrier(1,sD,sR,s3,25.0);
473 double t1[4], T1[4];
474 calt1(P1,P2,t1);
475 calt1(pR,P3,T1);
476 temp_PDF = 0;
477 for(int i=0; i<4; i++){
478 temp_PDF += t1[i]*T1[i]*G[i][i];
479 }
480 }
481 if(Ang == 2){
482 B[0] = barrier(2,sR,s1,s2,3.0,mass);
483 B[1] = barrier(2,sD,sR,s3,5.0,mDsM);
484 double t2[4][4], T2[4][4];
485 calt2(P1,P2,t2);
486 calt2(pR,P3,T2);
487 temp_PDF = 0;
488 for(int i=0; i<4; i++){
489 for(int j=0; j<4; j++){
490 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
491 }
492 }
493 }
494 v_re = temp_PDF*B[0]*B[1];
495 return v_re;
496}
497
498
499void EvtDsToKpipi::rhoab(double sa, double sb, double sc, double res[2]) {
500 double tmp = sa+sb-sc;
501 double q = 0.25*tmp*tmp/sa-sb;
502 if(q>=0) {
503 res[0]=2.0*sqrt(q/sa);
504 res[1]=0.0;
505 } else {
506 res[0]=0.0;
507 res[1]=2.0*sqrt(-q/sa);
508 }
509}
510void EvtDsToKpipi::rho4Pi(double sa, double res[2]) {
511 double temp = 1.0-0.3116765584/sa; // 0.3116765584=0.13957*0.13957*16
512 if(temp>=0) {
513 res[0]=sqrt(temp)/(1.0+exp(9.8-3.5*sa));
514 res[1]=0.0;
515 } else {
516 res[0]=0.0;
517 res[1]=sqrt(-temp)/(1.0+exp(9.8-3.5*sa));//?????????????????????
518 }
519}
520
521void EvtDsToKpipi::propagatorsigma500(double sa, double sb, double sc, double prop[2]) {//for pipi_s wave
522 double f = 0.5843+1.6663*sa;
523 const double M = 0.9264;
524 const double mass2 = 0.85821696; // M*M
525 const double mpi2d2 = 0.00973989245;
526 double g1 = f*(sa-mpi2d2)/(mass2-mpi2d2)*exp((mass2-sa)/1.082);
527 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
528 rhoab(sa,sb,sc,rho1s);
529 rhoab(mass2,sb,sc,rho1M);
530 rho4Pi(sa,rho2s);
531 rho4Pi(mass2,rho2M);
532 Com_Divide(rho1s,rho1M,rho1);
533 Com_Divide(rho2s,rho2M,rho2);
534 double a[2], b[2];
535 a[0] = 1.0;
536 a[1] = 0.0;
537 b[0] = mass2-sa+M*(g1*rho1[1]+0.0024*rho2[1]);
538 b[1] = -M*(g1*rho1[0]+0.0024*rho2[0]);
539 Com_Divide(a,b,prop);
540}
541
542
543void EvtDsToKpipi::calEva(double* K, double* Pi1, double* Pi2, double *mass1, double *width1, double *amp, double *phase,int* g0,int* spin, int* modetype, int nstates, double & Result)
544{
545 //double Ks0[4], Pic[4], Pi0[4];
546 double rRes = 3.0;
547 double rRess = 9.0;
548 double P23[4], P13[4];
549 double cof[2], amp_PDF[2], PDF[2];
550// double scpi, sks02, sks01;
551// double s12, s13, s23;
552 double s13, s23;
553 for(int i=0; i<4; i++){
554 // P12[i] = Ks01[i] + Ks02[i];
555 P13[i] = K[i] + Pi2[i];
556 P23[i] = Pi1[i] + Pi2[i];
557 }
558 s13 = SCADot(P13,P13);
559 s23 = SCADot(P23,P23);
560 double s1,s2,s3;
561 s1 = SCADot(K,K);
562 s2 = SCADot(Pi1,Pi1);
563 s3 = SCADot(Pi2,Pi2);
564 double pro[2], temp_PDF, amp_tmp[2],temp_PDF1 ,temp_PDF2,pro1[2],pro2[2];
565// double mass1sq;
566 amp_PDF[0] = 0;
567 amp_PDF[1] = 0;
568 PDF[0] = 0;
569 PDF[1] = 0;
570 amp_tmp[0] = 0;
571 amp_tmp[1] = 0;
572 for(int i=0; i<nstates; i++) {
573 amp_tmp[0] = 0;
574 amp_tmp[1] = 0;
575 // mass1sq = mass1[i]*mass1[i];
576 cof[0] = amp[i]*cos(phase[i]);
577 cof[1] = amp[i]*sin(phase[i]);
578 temp_PDF = 0;
579
580 if(modetype[i] == 13){
581 temp_PDF = DDalitz(K, Pi2, Pi1, spin[i], mass1[i]);
582 if(g0[i]==1) propagatorRBW(mass1[i],width1[i],s13,mKa2,mPi2,rRess,spin[i],pro);
583 if(g0[i]==2) { // K*1430 Flatte
584 //double skm2[2]={s1, mass_EtaP *mass_EtaP};
585 double skm2[2]={s1, 0.95778*0.95778};
586 double spi2[2]={s3, mass_Kaon *mass_Kaon};
587 propagatorKstr1430(mass1[i],s13,skm2,spi2,pro);
588 }
589 if(g0[i]==0){
590 pro[0] = 1;
591 pro[1] = 0;
592 }
593 amp_tmp[0] = temp_PDF*pro[0];
594 amp_tmp[1] = temp_PDF*pro[1];
595 }
596
597 if(modetype[i] == 23){
598 temp_PDF = DDalitz(Pi1, Pi2, K, spin[i], mass1[i]);
599 if(g0[i]==4) propagatorFlatte(mass1[i],width1[i],s23,pro); // Only for f0(980)
600 if(g0[i]==3) propagatorCBW(mass1[i],width1[i],s23,pro);
601 if(g0[i]==2) propagatorRBW(mass1[i],width1[i],s23,mPi2,mPi2,rRess,spin[i],pro);
602 if(g0[i]==1) propagatorGS(mass1[i],width1[i],s23,mPi2,mPi2,9.0,pro);
603 if(g0[i]==500) propagatorsigma500(s23, s2, s3, pro);
604 if(g0[i]==0){
605 pro[0] = 1;
606 pro[1] = 0;
607 }
608 amp_tmp[0] = temp_PDF*pro[0];
609 amp_tmp[1] = temp_PDF*pro[1];
610 }
611 Com_Multi(amp_tmp,cof,amp_PDF);
612 PDF[0] += amp_PDF[0];
613 PDF[1] += amp_PDF[1];
614 //cout<<"PDF: "<<PDF[0]<<" "<<PDF[1]<<endl;
615 //cout<<"amp_PDF: "<<amp_PDF[0]<<" "<<amp_PDF[1]<<endl;
616 }
617
618 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
619 // cout<<"value = "<<value<<endl;
620 if(value <=0) value = 1e-20;
621 Result = value;
622 // cout<<"value = "<<value<<" Result = "<<Result<<endl;
623}
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
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
t
TTree * t
Definition binning.cxx:23
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:533
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
EvtDsToKpipi::decay
void decay(EvtParticle *p)
Definition EvtDsToKpipi.cc:155
EvtDsToKpipi::getName
void getName(std::string &name)
Definition EvtDsToKpipi.cc:34
EvtDsToKpipi::clone
EvtDecayBase * clone()
Definition EvtDsToKpipi.cc:38
EvtDsToKpipi::~EvtDsToKpipi
virtual ~EvtDsToKpipi()
Definition EvtDsToKpipi.cc:32
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
s23
double double double double double * s23
Definition qcdloop1.h:77
b
const double b
Definition slope.cxx:9