BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKKpi.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: EvtDToKKpi.cc
11// the necessary file: Evt.hh
12//
13// Description: D+ -> K- K+ pi+,
14//
15//
16// Modification history:
17// Liaoyuan Dong Fri Dec 02 23:35:56 CST 2022 Module created
18//------------------------------------------------------------------------
19
20#include "EvtGenBase/EvtPatches.hh"
21#include <fstream>
22#include <stdlib.h>
23#include "EvtGenBase/EvtParticle.hh"
24#include "EvtGenBase/EvtGenKine.hh"
25#include "EvtGenBase/EvtPDL.hh"
26#include "EvtGenBase/EvtReport.hh"
27#include "EvtGenBase/EvtResonance.hh"
28#include "EvtGenBase/EvtResonance2.hh"
29#include "EvtGenModels/EvtDToKKpi.hh"
30#include <string>
31#include "EvtGenBase/EvtConst.hh"
32#include "EvtGenBase/EvtDecayTable.hh"
33using std::endl;
34
35 EvtDToKKpi::~EvtDToKKpi() {}
36
37void EvtDToKKpi::getName(std::string& model_name){
38 model_name="DToKKpi";
39}
40
44
45void EvtDToKKpi::init(){
46 // check that there are 0 arguments
47 checkNArg(0);
48 checkNDaug(3);
49 checkSpinParent(EvtSpinType::SCALAR);
50 checkSpinDaughter(0,EvtSpinType::SCALAR);
51 checkSpinDaughter(1,EvtSpinType::SCALAR);
52 checkSpinDaughter(2,EvtSpinType::SCALAR);
53
54 mass[0] = 0.89581;//Ksta0
55 mass[1] = 1.019461;//phi1020
56 mass[2] = 0.965;//f0(980)
57 mass[3] = 1.474;//a0(1450)
58 mass[4] = 1.68;//phi1680
59 mass[5] = 1.4324;//K21430
60 mass[6] = 0.89581;//KPi Swave
61 mass[7] = 1.3183;//a21320
62 width[0] = 0.0474;//Ksta0
63 width[1] = 0.004266;//phi1020
64 width[2] = 0.0400;//f0(980)
65 width[3] = 0.265;//a01450
66 width[4] = 0.15;//phi1680
67 width[5] = 0.109;//K21430
68 width[6] = 0.0474;//KPi Swave
69 width[7] = 0.107;//a01320
70
71 rho[0] = 1;//K*892
72 phi[0] = 0;//
73
74 rho[1] = 1.121155232058;// phi1020
75 phi[1] = 3.322105212534;//
76
77 rho[2] = 1.214264427529;// f0980
78 phi[2] = -5.321965906572;
79
80 rho[3] = 1.086952630367;// a01450
81 phi[3] = -4.213398459878;
82
83 rho[4] = 1.121071578335;// phi1680
84 phi[4] =-1.055595448689;//
85
86 rho[5] =-1.473195648049;// K21430
87 phi[5] =5.841826108667;
88
89 rho[6] = 5.396864431453;// KPISW
90 phi[6] =-0.360227442474;
91
92 rho[7] = 0.436480711396;// a21320
93 phi[7] = 0.778423040673;
94
95 spin[0] = 1; //K*892
96 spin[1] = 1; //phi1020
97 spin[2] = 0; //f0980
98 spin[3] = 0; //a01450
99 spin[4] = 1; //phi1680
100 spin[5] = 2; //K21430
101 spin[6] = 0; //KPISW
102 spin[7] = 2; //a21320
103
104 modetype[0] = 13;// K*892
105 modetype[1] = 12;// phi1020
106 modetype[2] = 12;// f0980
107 modetype[3] = 12;// a01450
108 modetype[4] = 12;// phi1680
109 modetype[5] = 13;// K21430 error
110 modetype[6] = 13;// KPiSW
111 modetype[7] = 12;// a21320 error
112
113/*
114 std::cout << "EvtDToKKpi (May 01, 2023) ==> Initialization" << std::endl;
115 for (int i=0; i<8; i++) {
116 cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
117 }
118*/
119
120 mDp = 1.86966;
121 mK0 = 0.497611;
122 mKa = 0.49368;
123 mPi = 0.13957;
124 mK02 = 0.237616707;
125 mPi2 = 0.01947978;
126 mass_EtaP = 0.95778;
127 mass_Kaon = 0.49368;
128
129 math_pi = 3.1415926;
130 mass_Pion = 0.13957;
131 mass_Pion2 = 0.0194797849;
132 mass_2Pion = 0.27914;
133 math_2pi = 6.2831852;
134 rD2 = 25.0; // 5*5
135 rRes2 = 9.0; // 3*3
136 //g1 = 0.5468;
137 g2 = 0.23;
138 GS1 = 0.636619783;
139 GS2 = 0.01860182466;
140 GS3 = 0.1591549458;
141 GS4 = 0.00620060822;
142
143
144 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
145 for (int i=0; i<4; i++) {
146 for (int j=0; j<4; j++) {
147 G[i][j] = GG[i][j];
148 }
149 }
150}
151
152void EvtDToKKpi::initProbMax() {
153 setProbMax(11700.0);//MAXprob = 11690.0512849074
154}
155
156
157void EvtDToKKpi::decay( EvtParticle *p){
158//==================================================get max============================================
159/*
160 double maxprob = 0.0;
161 for(int ir=0;ir<=60000000;ir++){
162 p->initializePhaseSpace(getNDaug(),getDaugs());
163 EvtVector4R D1 = p->getDaug(0)->getP4();
164 EvtVector4R D2 = p->getDaug(1)->getP4();
165 EvtVector4R D3 = p->getDaug(2)->getP4();
166
167 double P1[4], P2[4], P3[4];
168 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
169 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
170 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
171
172 double value;
173 int nstates=8;
174 int g0[8]={1,1,4,1,1,1,3,1};
175 double r0[8]={3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0};
176 double r1[8]={5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0};
177
178 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value,r0,r1);
179
180 if (value<0) continue;
181 if(value>maxprob) {
182 maxprob=value;
183 cout << "ir= " << ir << endl;
184 cout << "MAX====> " << maxprob << endl;
185 }
186 }
187 printf("MAXprob = %.10f\n",maxprob);
188*/
189//==================================================get max============================================
190
191
192 p->initializePhaseSpace(getNDaug(),getDaugs());
193 EvtVector4R D1 = p->getDaug(0)->getP4();
194 EvtVector4R D2 = p->getDaug(1)->getP4();
195 EvtVector4R D3 = p->getDaug(2)->getP4();
196
197 double P1[4], P2[4], P3[4];
198 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
199 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
200 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
201
202 //Check
203 //P1[0] = 0.694086; P1[1] = 0.343301; P1[2] = 0.263795; P1[3] = -0.224934;
204 //P2[0] = 0.719665; P2[1] = -0.165392; P2[2] = -0.495843; P2[3] = -0.0314011;
205 //P3[0] = 0.474568; P3[1] = 0.0149532 ; P3[2] = 0.405429; P3[3] = 0.202826;
206
207 double value;
208 int nstates=8;
209 int g0[8]={1,1,4,1,1,1,3,1};
210 double r0[8]={3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0};
211 double r1[8]={5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0};
212
213 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value,r0,r1);
214 setProb(value);
215 return ;
216
217}
218
219
220void EvtDToKKpi::Com_Multi(double a1[2], double a2[2], double res[2])
221{
222 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
223 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
224}
225void EvtDToKKpi::Com_Divide(double a1[2], double a2[2], double res[2])
226{
227 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
228 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
229 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
230}
231//------------base---------------------------------
232double EvtDToKKpi::SCADot(double a1[4], double a2[4])
233{
234 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
235 return _cal;
236}
237double EvtDToKKpi::barrier(int l, double sa, double sb, double sc, double r, double mass)
238{
239 double q = fabs((sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb);
240 double z;
241 z = q*r*r;
242 double sa0;
243 sa0 = mass*mass;
244 double q0 = fabs((sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb);
245 double z0 = q0*r*r;
246 double F = 0.0;
247 if(l == 0) F = 1;
248 if(l == 1) F = sqrt((1+z0)/(1+z));
249 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
250 return F;
251}
252//------------------spin-------------------------------------------
253void EvtDToKKpi::calt1(double daug1[4], double daug2[4], double t1[4])
254{
255 double p, pq, tmp;
256 double pa[4], qa[4];
257 for(int i=0; i<4; i++) {
258 pa[i] = daug1[i] + daug2[i];
259 qa[i] = daug1[i] - daug2[i];
260 }
261 p = SCADot(pa,pa);
262 pq = SCADot(pa,qa);
263 tmp = pq/p;
264 for(int i=0; i<4; i++) {
265 t1[i] = qa[i] - tmp*pa[i];
266 }
267}
268void EvtDToKKpi::calt2(double daug1[4], double daug2[4], double t2[4][4])
269{
270 double p, r;
271 double pa[4], t1[4];
272 calt1(daug1,daug2,t1);
273 r = SCADot(t1,t1)/3.0;
274 for(int i=0; i<4; i++) {
275 pa[i] = daug1[i] + daug2[i];
276 }
277 p = SCADot(pa,pa);
278 for(int i=0; i<4; i++) {
279 for(int j=0; j<4; j++) {
280 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
281 }
282 }
283}
284//-------------------prop--------------------------------------------
285void EvtDToKKpi::propagator(double mass2, double mass, double width, double sx, double prop[2])
286{
287 double a[2], b[2];
288 a[0] = 1;
289 a[1] = 0;
290 b[0] = mass2-sx;
291 b[1] = -mass*width;
292 Com_Divide(a,b,prop);
293}
294double EvtDToKKpi::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
295{
296 double widm = 0.;
297 double m = sqrt(sa);
298 double tmp = sb-sc;
299 double tmp1 = sa+tmp;
300 double q = fabs(0.25*tmp1*tmp1/sa-sb);
301 double tmp2 = mass2+tmp;
302 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
303 double z = q*r2;
304 double z0 = q0*r2;
305 double t = q/q0;
306 if(l == 0) {widm = sqrt(t)*mass/m;}
307 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
308 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
309 return widm;
310}
311double EvtDToKKpi::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
312{
313 double widm = 0.;
314 double m = sqrt(sa);
315 double tmp = sb-sc;
316 double tmp1 = sa+tmp;
317 double q = fabs(0.25*tmp1*tmp1/sa-sb);
318 double tmp2 = mass2+tmp;
319 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
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 EvtDToKKpi::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
328{
329
330 double a[2], b[2];
331 double mass2 = mass*mass;
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}
338
339void EvtDToKKpi::propagatorFlatte(double mass, double width, double sa, double prop[2]){
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
364 double 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}
371
372void EvtDToKKpi::rhoab(double sa, double sb, double sc, double res[2]) {
373 double tmp = sa+sb-sc;
374 double q = 0.25*tmp*tmp/sa-sb;
375 if(q>=0) {
376 res[0]=2.0*sqrt(q/sa);
377 res[1]=0.0;
378 } else {
379 res[0]=0.0;
380 res[1]=2.0*sqrt(-q/sa);
381 }
382}
383void EvtDToKKpi::rho4Pi(double sa, double res[2]) {
384 double temp = 1.0-0.3116765584/sa; // 0.3116765584=0.13957*0.13957*16
385 if(temp>=0) {
386 res[0]=sqrt(temp)/(1.0+exp(9.8-3.5*sa));
387 res[1]=0.0;
388 } else {
389 res[0]=0.0;
390 res[1]=sqrt(-temp)/(1.0+exp(9.8-3.5*sa));//?????????????????????
391 }
392}
393
394
395void EvtDToKKpi::propagatorsigma500(double sa, double sb, double sc, double prop[2]) {
396 double f = 0.5843+1.6663*sa;
397 const double M = 0.9264;//M(f0500)
398 const double mass2 = 0.85821696; // M*M
399 const double mpi2d2 = 0.00973989245;//mass_Pion2/2 = 0.0194797849/2
400 double g1 = f*(sa-mpi2d2)/(mass2-mpi2d2)*exp((mass2-sa)/1.082);
401 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
402 rhoab(sa,sb,sc,rho1s);
403 rhoab(mass2,sb,sc,rho1M);
404 rho4Pi(sa,rho2s);
405 rho4Pi(mass2,rho2M);
406 Com_Divide(rho1s,rho1M,rho1);
407 Com_Divide(rho2s,rho2M,rho2);
408 double a[2], b[2];
409 a[0] = 1.0;
410 a[1] = 0.0;
411 b[0] = mass2-sa+M*(g1*rho1[1]+0.0024*rho2[1]);
412 b[1] = -M*(g1*rho1[0]+0.0024*rho2[0]);
413 Com_Divide(a,b,prop);
414}
415void EvtDToKKpi::Flatte_rhoab(double sa, double sb, double rho[2])
416{
417 double q = 1.0-(4*sb/sa);
418
419 if(q>0) {
420 rho[0]=sqrt(q);
421 rho[1]=0;
422 }
423 else if(q<0){
424 rho[0]=0;
425 rho[1]=sqrt(-q);
426 }
427}
428
429void EvtDToKKpi::propagator980(double mass, double sx, double *sb, double prop[2])
430{
431
432 double gpipi1 = 2.0/3.0*0.165;
433 double gpipi2 = 1.0/3.0*0.165;
434 double gKK1 = 0.5*0.69465;
435 double gKK2 = 0.5*0.69465;
436
437 double unit[2]={1.0};
438 double ci[2]={0,1};
439 double rho1[2];
440 Flatte_rhoab(sx,sb[0],rho1);
441 double rho2[2];
442 Flatte_rhoab(sx,sb[1],rho2);
443 double rho3[2];
444 Flatte_rhoab(sx,sb[2],rho3);
445 double rho4[2];
446 Flatte_rhoab(sx,sb[3],rho4);
447
448 double tmp1[2]={gpipi1,0};
449 double tmp11[2];
450 double tmp2[2]={gpipi2,0};
451 double tmp22[2];
452 double tmp3[2]={gKK1,0};
453 double tmp33[2];
454 double tmp4[2]={gKK2,0};
455 double tmp44[2];
456
457 Com_Multi(tmp1,rho1,tmp11);
458 Com_Multi(tmp2,rho2,tmp22);
459 Com_Multi(tmp3,rho3,tmp33);
460 Com_Multi(tmp4,rho4,tmp44);
461
462 double tmp5[2]={tmp11[0]+tmp22[0]+tmp33[0]+tmp44[0],tmp11[1]+tmp22[1]+tmp33[1]+tmp44[1]};
463 double tmp51[2];
464 Com_Multi(tmp5, ci,tmp51);
465 double tmp6[2]={mass*mass-sx-tmp51[0], -1.0*tmp51[1]};
466
467 Com_Divide( unit,tmp6, prop);
468
469}
470
471void EvtDToKKpi::KPiSLASS(double sa, double sb, double sc, double prop[2]) {
472 const double m1430 = 1.441;
473 const double sa0 = 2.076481; // m1430*m1430;
474 const double w1430 = 0.193;
475 const double Lass1 = 0.25/sa0;
476 double tmp = sb-sc;
477 double tmp1 = sa0+tmp;
478 double q0 = fabs(Lass1*tmp1*tmp1-sb);
479 double tmp2 = sa+tmp;
480 double qs = 0.25*tmp2*tmp2/sa-sb;
481 double q = sqrt(qs);
482 double width = w1430*q*m1430/sqrt(sa*q0);
483 double temp_R = atan(m1430*width/(sa0-sa));
484 if(temp_R<0) temp_R += math_pi;
485 double deltaR = -109.7*math_pi/180.0 + temp_R;
486 double temp_F = atan(0.226*q/(2.0-3.8194*qs)); // 2.0*0.113 = 0.226; -33.8*0.113 = -3.8194
487 if(temp_F<0) temp_F += math_pi;
488 double deltaF = 0.1*math_pi/180.0 + temp_F;
489 double deltaS = deltaR + 2.0*deltaF;
490 double t1 = 0.96*sin(deltaF);
491 double t2 = sin(deltaR);
492 double CF[2], CS[2];
493 CF[0] = cos(deltaF);
494 CF[1] = sin(deltaF);
495 CS[0] = cos(deltaS);
496 CS[1] = sin(deltaS);
497 prop[0] = t1*CF[0] + t2*CS[0];
498 prop[1] = t1*CF[1] + t2*CS[1];
499}
500//------------GS---used by rho----------------------------
501void EvtDToKKpi::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
502{
503 double a[2], b[2];
504 double tmp = sb-sc;
505 double tmp1 = sa+tmp;
506 double q2 = fabs(0.25*tmp1*tmp1/sa-sb);
507
508 double tmp2 = mass2+tmp;
509 double q02 = fabs(0.25*tmp2*tmp2/mass2-sb);
510
511 double q = sqrt(q2);
512 double q0 = sqrt(q02);
513 double m = sqrt(sa);
514 double q03 = q0*q02;
515 double tmp3 = log(mass+2*q0)+1.2926305904; // log(mpi0+mpip) = -1.2926305904;
516
517 double h = GS1*q/m*(log(m+2*q)+1.2926305904);
518 double h0 = GS1*q0/mass*tmp3;
519 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
520 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
521 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
522
523 a[0] = 1.0+d*width/mass;
524 a[1] = 0.0;
525 b[0] = mass2-sa+width*f;
526 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
527 Com_Divide(a,b,prop);
528}
529
530double EvtDToKKpi::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
531 double pR[4], pD[4];
532 double temp_PDF, v_re;
533 temp_PDF = 0.0;
534 v_re = 0.0;
535 double B[2], s1, s2, s3, sR, sD;
536 for(int i=0; i<4; i++){
537 pR[i] = P1[i] + P2[i];
538 pD[i] = pR[i] + P3[i];
539 }
540 s1 = SCADot(P1,P1);
541 s2 = SCADot(P2,P2);
542 s3 = SCADot(P3,P3);
543 sR = SCADot(pR,pR);
544 sD = SCADot(pD,pD);
545 int G[4][4];
546 for(int i=0; i!=4; i++){
547 for(int j=0; j!=4; j++){
548 if(i==j){
549 if(i==0) G[i][j] = 1;
550 else G[i][j] = -1;
551 }
552 else G[i][j] = 0;
553 }
554 }
555 if(Ang == 0){
556 B[0] = 1;
557 B[1] = 1;
558 temp_PDF = 1;
559 }
560 if(Ang == 1){
561 B[0] = barrier(1,sR,s1,s2,3.0,mass);
562 B[1] = barrier(1,sD,sR,s3,5.0,mDp);
563 double t1[4], T1[4];
564 calt1(P1,P2,t1);
565 calt1(pR,P3,T1);
566 temp_PDF = 0;
567 for(int i=0; i<4; i++){
568 temp_PDF += t1[i]*T1[i]*G[i][i];
569 }
570 }
571 if(Ang == 2){
572 B[0] = barrier(2,sR,s1,s2,3.0,mass);
573 B[1] = barrier(2,sD,sR,s3,5.0,mDp);
574 double t2[4][4], T2[4][4];
575 calt2(P1,P2,t2);
576 calt2(pR,P3,T2);
577 temp_PDF = 0;
578 for(int i=0; i<4; i++){
579 for(int j=0; j<4; j++){
580 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
581 }
582 }
583 }
584 v_re = temp_PDF*B[0]*B[1];
585 return v_re;
586}
587
588void EvtDToKKpi::calEva(double* K01, double* K02, double* Pi, double *mass1, double *width1, double *amp, double *phase,int* g0,int* spin, int* modetype, int nstates, double & Result , double*r0 ,double*r1)
589{
590
591 double P12[4], P23[4], P13[4];
592 double cof[2], amp_PDF[2], PDF[2];
593 double Amp_KPiS[2];
594 double s1, s2, s3;
595 double s12, s13, s23;
596 for(int i=0; i<4; i++){
597 P12[i] = K01[i] + K02[i];
598 P13[i] = K01[i] + Pi[i];
599 P23[i] = K02[i] + Pi[i];
600 }
601 s1 = SCADot(K01,K01);
602 s2 = SCADot(K02,K02);
603 s3 = SCADot(Pi,Pi);
604
605 s12 = SCADot(P12,P12);
606 s13 = SCADot(P13,P13);
607 s23 = SCADot(P23,P23);
608
609 double pro[2], temp_PDF, amp_tmp[2];
610 double mass1sq;
611 amp_PDF[0] = 0;
612 amp_PDF[1] = 0;
613 PDF[0] = 0;
614 PDF[1] = 0;
615 amp_tmp[0] = 0;
616 amp_tmp[1] = 0;
617 for(int i=0; i<nstates; i++) {
618 amp_tmp[0] = 0;
619 amp_tmp[1] = 0;
620 mass1sq = mass1[i]*mass1[i];
621 cof[0] = amp[i]*cos(phase[i]);
622 cof[1] = amp[i]*sin(phase[i]);
623 temp_PDF = 0;
624
625 if(modetype[i] == 12){
626 temp_PDF = DDalitz(K01, K02, Pi, spin[i], mass1[i]);
627 if(g0[i]==1) propagatorRBW(mass1[i],width1[i],s12,s1,s2,rRes2,spin[i],pro);
628 if(g0[i]==4) propagatorFlatte(mass1[i],width1[i],s12,pro); // Only for f0(980)
629 amp_tmp[0] = temp_PDF*pro[0];
630 amp_tmp[1] = temp_PDF*pro[1];
631 }
632
633 if(modetype[i] == 13){
634 temp_PDF = DDalitz(K01, Pi, K02, spin[i], mass1[i]);
635 if(g0[i]==1) propagatorRBW(mass1[i],width1[i],s13,s1,s3,rRes2,spin[i],pro);
636 if(g0[i]==3) KPiSLASS(s13,s1,s3,pro);
637
638 amp_tmp[0] = temp_PDF*pro[0];
639 amp_tmp[1] = temp_PDF*pro[1];
640 }
641
642 Com_Multi(amp_tmp,cof,amp_PDF);
643 PDF[0] += amp_PDF[0];
644 PDF[1] += amp_PDF[1];
645 }
646 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
647 if(value <=0) value = 1e-20;
648 Result = value;
649 // printf("%s %8.15f\n","AAABBBvalue = ",value);
650 }
651
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
CS
TCrossPart * CS
Definition Mcgpj.cxx:51
t
TTree * t
Definition binning.cxx:23
EvtDToKKpi::getName
void getName(std::string &name)
Definition EvtDToKKpi.cc:37
EvtDToKKpi::initProbMax
void initProbMax()
Definition EvtDToKKpi.cc:152
EvtDToKKpi::clone
EvtDecayBase * clone()
Definition EvtDToKKpi.cc:41
EvtDToKKpi::decay
void decay(EvtParticle *p)
Definition EvtDToKKpi.cc:157
EvtDToKKpi::init
void init()
Definition EvtDToKKpi.cc:45
EvtDToKKpi::~EvtDToKKpi
virtual ~EvtDToKKpi()
Definition EvtDToKKpi.cc:35
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
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
b
const double b
Definition slope.cxx:9