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