BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKSKpi0.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: EvtDsToKSKpi0.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/EvtDsToKSKpi0.hh"
27#include <string>
28#include "EvtGenBase/EvtConst.hh"
29#include "EvtGenBase/EvtDecayTable.hh"
30using std::endl;
31
32EvtDsToKSKpi0::~EvtDsToKSKpi0() {}
33
34void EvtDsToKSKpi0::getName(std::string& model_name){
35 model_name="DsToKSKpi0";
36}
37
41
42void EvtDsToKSKpi0::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] = 0;
53 phi[1] = 0;
54 rho[1] = 1;
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[0] = 5.2793e+00;
64 phi[2] = -6.4411e+00;
65 phi[3] = 2.2286e-01;
66 phi[5] = -2.5771e+00;
67
68 rho[0] = 8.5648e-01;
69 rho[2] = 6.7257e-01;
70 rho[3] = 1.8923e+00;
71 rho[5] = 2.3549e+00;
72
73 modetype[0]= 12;
74 modetype[1]= 13;
75 modetype[2]= 23;
76 modetype[3]= 13;
77 modetype[4]= 23;
78 modetype[5]= 12;
79
80 //cout << "DsToKSKpi0 :" << endl;
81 //for (int i=0; i<6; i++) {
82 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
83 //}
84
85 width[0] = 0.0732;
86 width[1] = 0.0473;
87 width[2] = 0.0503;
88 width[3] = 0.232;
89 width[4] = 0.1767;
90 width[5] = 0.0993;
91
92 mass[0] = 1.0007;
93 mass[1] = 0.89555;
94 mass[2] = 0.89176;
95 mass[3] = 1.414;
96 mass[4] = 1.600;
97 mass[5] = 1.8154;
98
99 mD = 1.86486;
100 mDs = 1.9683;
101 rRes = 9.0;
102 rD = 5.0;
103 metap = 0.95778;
104 mkstr = 0.89594;
105 mk0 = 0.497614;
106 mass_Kaon = 0.49368;
107 mass_Pion = 0.13957;
108 mass_Pi0 = 0.1349766;
109 math_pi = 3.1415926;
110
111 GS1 = 0.636619783;
112 GS2 = 0.01860182466;
113 GS3 = 0.1591549458; // 1/(2*math_2pi)
114 GS4 = 0.00620060822; // mass_Pion2/math_pi
115
116 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
117 for (int i=0; i<4; i++) {
118 for (int j=0; j<4; j++) {
119 G[i][j] = GG[i][j];
120 }
121 }
122}
123
124void EvtDsToKSKpi0::initProbMax() {
125 setProbMax(1871.0);
126}
127
128void EvtDsToKSKpi0::decay( EvtParticle *p){
129/*
130 double maxprob = 0.0;
131 for(int ir=0;ir<=60000000;ir++){
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 //value = calDalEva(P1, P2, P3);
144 int g0[6]={1,1,1,1,1,1};
145 int spin[6]={0,1,1,1,1,0};
146
147 int nstates=6;
148 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
149 if (value<0) continue;
150 if(value>maxprob) {
151 maxprob=value;
152 cout << "ir= " << ir << endl;
153 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3] <<"};"<<" "<< sqrt(P1[0]*P1[0]-P1[1]*P1[1]-P1[2]*P1[2]-P1[3]*P1[3]) <<endl;
154 cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3] <<"};"<< " "<< sqrt(P2[0]*P2[0]-P2[1]*P2[1]-P2[2]*P2[2]-P2[3]*P2[3]) <<endl;
155 cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< " "<< sqrt(P3[0]*P3[0]-P3[1]*P3[1]-P3[2]*P3[2]-P3[3]*P3[3]) <<endl;
156 cout << "MAX====> " << maxprob << endl;
157
158
159 }
160
161 }
162 printf("MAXprob = %.10f\n",maxprob);
163*/
164
165
166 p->initializePhaseSpace(getNDaug(),getDaugs());
167 EvtVector4R D1 = p->getDaug(0)->getP4();
168 EvtVector4R D2 = p->getDaug(1)->getP4();
169 EvtVector4R D3 = p->getDaug(2)->getP4();
170
171 double P1[4], P2[4], P3[4];
172 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
173 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
174 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
175
176// P1[0] = 0.619505234038183; P1[1] = -0.091552038742807; P1[2] = -0.326369642554372; P1[3] = 0.145835944211028;
177// P2[0] = 0.537869938150145; P2[1] = 0.060139034763020; P2[2] = -0.082936309719630; P2[3] = -0.187328468377079;
178// P3[0] = 0.875173341522316; P3[1] = 0.176856301894967; P3[2] = 0.816309543635941; P3[3] = -0.223764071089403;
179// P1[0] = 0.958942320949657; P1[1] =-0.130886193017869 ; P1[2] = 0.458368295934040; P1[3] = -0.666871795529042;
180// P2[0] = 0.840895041827643; P2[1] = -0.279686661949455; P2[2] = -0.173310031669695; P2[3] = 0.595924907258905;
181// P3[0] = 0.205192165542905 ; P3[1] = 0.148385551896119; P3[2] = -0.008476135905980; P3[3] = 0.042367724031334;
182// P1[0] = ; P1[1] = ; P1[2] = ; P1[3] = ;
183// P2[0] = ; P2[1] = ; P2[2] = ; P2[3] = ;
184// P3[0] = ; P3[1] = ; P3[2] = ; P3[3] = ;
185 // P1[0] =0.61748841695829326248 ; P1[1] =-0.16431886625773378663 ; P1[2] = -0.015984626326897594106; P1[3] = -0.32623606810397931532;
186 // P2[0] =0.61165109530060246534 ; P2[1] = -0.16009556658150744801; P2[2] = -0.016518497076036600668; P2[3] = -0.32325478223947284873;
187 // P3[0] =0.73916048774110421071 ; P3[1] = 0.32441443283924120689; P3[2] = 0.032503123402934194774; P3[3] = 0.64949085034345221956;
188
189 double value;
190 int g0[6]={1,1,1,1,1,1};
191 int spin[6]={0,1,1,1,1,0};
192 int nstates=6;
193 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
194
195 setProb(value);
196
197 return ;
198
199}
200
201void EvtDsToKSKpi0::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 EvtDsToKSKpi0::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//------------base---------------------------------
213double EvtDsToKSKpi0::SCADot(double a1[4], double a2[4])
214{
215 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
216 return _cal;
217}
218
219double EvtDsToKSKpi0::barrier(int l, double sa, double sb, double sc, double r, double mass)
220{
221 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
222 if(q < 0) q = 1e-16;
223 double z;
224 z = q*r*r;
225 double sa0;
226 sa0 = mass*mass;
227 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
228 if(q0 < 0) q0 = 1e-16;
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/*
237double EvtDsToKSKpi0::Barrier(int l, double sa, double sb, double sc, double r2)
238{
239 double F;
240 double tmp = sa+sb-sc;
241 double q = 0.25*tmp*tmp/sa-sb;
242 if (q < 0) q = 1e-16;
243 double z = q*r2;
244 if (l==1) {
245 F = sqrt(2.0*z/(1.0+z));
246 }
247 else if (l==2) {
248 double z2 = z*z;
249 F = sqrt(13.0*z2/(9.0+3.0*z+z2));
250 } else {
251 F = 1.0;
252 }
253 return F;
254}
255*/
256void EvtDsToKSKpi0::calt1(double daug1[4], double daug2[4], double t1[4])
257{
258 double p, pq, tmp;
259 double pa[4], qa[4];
260 for(int i=0; i<4; i++) {
261 pa[i] = daug1[i] + daug2[i];
262 qa[i] = daug1[i] - daug2[i];
263 }
264 p = SCADot(pa,pa);
265 pq = SCADot(pa,qa);
266 tmp = pq/p;
267 for(int i=0; i<4; i++) {
268 t1[i] = qa[i] - tmp*pa[i];
269 }
270}
271void EvtDsToKSKpi0::calt2(double daug1[4], double daug2[4], double t2[4][4])
272{
273 double p, r;
274 double pa[4], t1[4];
275 calt1(daug1,daug2,t1);
276 r = SCADot(t1,t1)/3.0;
277 for(int i=0; i<4; i++) {
278 pa[i] = daug1[i] + daug2[i];
279 }
280 p = SCADot(pa,pa);
281 for(int i=0; i<4; i++) {
282 for(int j=0; j<4; j++) {
283 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
284 }
285 }
286}
287//-------------------prop--------------------------------------------
288
289double EvtDsToKSKpi0::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}
308double EvtDsToKSKpi0::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
309{
310 double widm = 0.;
311 double m = sqrt(sa);
312 double tmp = sb-sc;
313 double tmp1 = sa+tmp;
314 double q = 0.25*tmp1*tmp1/sa-sb;
315 if(q<0) q = 1e-16;
316 double tmp2 = mass2+tmp;
317 double q0 = 0.25*tmp2*tmp2/mass2-sb;
318 if(q0<0) q0 = 1e-16;
319 double z = q*r2;
320 double z0 = q0*r2;
321 double F = (1+z0)/(1+z);
322 double t = q/q0;
323 widm = t*sqrt(t)*mass/m*F;
324 return widm;
325}
326void EvtDsToKSKpi0::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
327{
328 double a[2], b[2];
329 a[0] = 1;
330 a[1] = 0;
331 b[0] = mass2-sa;
332 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
333 Com_Divide(a,b,prop);
334}
335
336void EvtDsToKSKpi0::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
337{
338 double a[2], b[2];
339 double tmp = sb-sc;
340 double tmp1 = sa+tmp;
341 double q2 = 0.25*tmp1*tmp1/sa-sb;
342 if(q2<0) q2 = 1e-16;
343
344 double tmp2 = mass2+tmp;
345 double q02 = 0.25*tmp2*tmp2/mass2-sb;
346 if(q02<0) q02 = 1e-16;
347
348 double q = sqrt(q2);
349 double q0 = sqrt(q02);
350 double m = sqrt(sa);
351 double q03 = q0*q02;
352 double tmp3 = log(mass+2*q0)+1.2926305904; // log(mpi0+mpip) = -1.2926305904;
353
354 double h = GS1*q/m*(log(m+2*q)+1.2926305904);
355 double h0 = GS1*q0/mass*tmp3;
356 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
357 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
358 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
359
360 a[0] = 1.0+d*width/mass;
361 a[1] = 0.0;
362 b[0] = mass2-sa+width*f;
363 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
364 Com_Divide(a,b,prop);
365}
366
367
368void EvtDsToKSKpi0::propagatorFlatte(double mass, double width, double sa, double sb, double sc, int r, double prop[2]){
369 double q, qKsK,qetapi;
370 double rhoab[2], rhoKsK[2];
371 q = 0.25*(sa+sb-sc)*(sa+sb-sc)/sa-sb;
372 qetapi=0.25*(sa+0.547862-0.13957039)*(sa+0.547862-0.13957039)/sa-0.547862*0.547862;
373 if(r == 0) qKsK = 0.25*sa - 0.49368*0.49368;
374 if(r == 1) qKsK = 0.25*(sa +sb-sc)*(sa +sb-sc)/sa-sb;
375 if(qetapi>0){
376 rhoab[0] = 2*sqrt(qetapi/sa);
377 rhoab[1] = 0;
378 }
379 if(qetapi<0){
380 rhoab[0] = 0;
381 rhoab[1] = 2*sqrt(-qetapi/sa);
382 }
383 if(qKsK>0){
384 rhoKsK[0] = 2*sqrt(qKsK/sa);
385 rhoKsK[1] = 0;
386 }
387 if(qKsK<0){
388 rhoKsK[0] = 0;
389 rhoKsK[1] = 2*sqrt(-qKsK/sa);
390 }
391 double a[2], b[2];
392 a[0] = 1;
393 a[1] = 0;
394 b[0] = mass*mass - sa + 0.341*rhoab[1] + 0.892*0.341*rhoKsK[1];
395 b[1] = - (0.341*rhoab[0] + 0.892*0.341*rhoKsK[0]);
396 Com_Divide(a,b,prop);
397}
398
399void EvtDsToKSKpi0::KPiSLASS(double sa, double sb, double sc, double prop[2]) {
400 const double m1430 = 1.441;
401 const double sa0 = 1.441*1.441; // m1430*m1430;
402 const double w1430 = 0.193;
403 const double Lass1 = 0.25/sa0;
404 double tmp = sb-sc;
405 double tmp1 = sa0+tmp;
406 double q0 = Lass1*tmp1*tmp1-sb;
407 if(q0<0) q0 = 1e-16;
408 double tmp2 = sa+tmp;
409 double qs = 0.25*tmp2*tmp2/sa-sb;
410 double q = sqrt(qs);
411 double width = w1430*q*m1430/sqrt(sa*q0);
412 double temp_R = atan(m1430*width/(sa0-sa));
413 if(temp_R<0) temp_R += math_pi;
414 double deltaR = -109.7 + temp_R;
415 double temp_F = atan(0.226*q/(2.0-3.819*qs)); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
416 if(temp_F<0) temp_F += math_pi;
417 double deltaF = 0.1 + temp_F;
418 double deltaS = deltaR + 2.0*deltaF;
419 double t1 = 0.96*sin(deltaF);
420 double t2 = sin(deltaR);
421 double CF[2], CS[2];
422 CF[0] = cos(deltaF);
423 CF[1] = sin(deltaF);
424 CS[0] = cos(deltaS);
425 CS[1] = sin(deltaS);
426 prop[0] = t1*CF[0] + t2*CS[0];
427 prop[1] = t1*CF[1] + t2*CS[1];
428}
429
430double EvtDsToKSKpi0::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
431 double pR[4], pD[4];
432 double temp_PDF, v_re;
433 temp_PDF = 0.0;
434 v_re = 0.0;
435 double B[2], s1, s2, s3, sR, sD;
436 for(int i=0; i<4; i++){
437 pR[i] = P1[i] + P2[i];
438 pD[i] = pR[i] + P3[i];
439 }
440 s1 = SCADot(P1,P1);
441 s2 = SCADot(P2,P2);
442 s3 = SCADot(P3,P3);
443 sR = SCADot(pR,pR);
444 sD = SCADot(pD,pD);
445 //int G[4][4];
446 //for(int i=0; i!=4; i++){
447 // for(int j=0; j!=4; j++){
448 // if(i==j){
449 // if(i==0) G[i][j] = 1;
450 // else G[i][j] = -1;
451 // }
452 // else G[i][j] = 0;
453 // }
454 //}
455 if(Ang == 0){
456 B[0] = 1;
457 B[1] = 1;
458 temp_PDF = 1;
459 }
460 if(Ang == 1){
461 B[0] = barrier(1,sR,s1,s2,3.0,mass);
462 B[1] = barrier(1,sD,sR,s3,5.0,1.9683);
463 //B[0] = Barrier(1,sR,s1,s2,9.0);
464 //B[1] = Barrier(1,sD,sR,s3,25.0);
465 double t1[4], T1[4];
466 calt1(P1,P2,t1);
467 calt1(pR,P3,T1);
468 temp_PDF = 0;
469 for(int i=0; i<4; i++){
470 temp_PDF += t1[i]*T1[i]*G[i][i];
471 }
472 }
473 if(Ang == 2){
474 B[0] = barrier(2,sR,s1,s2,3.0,mass);
475 B[1] = barrier(2,sD,sR,s3,5.0,1.9683);
476 // B[0] = Barrier(2,sR,s1,s2,9.0);
477 // B[1] = Barrier(2,sD,sR,s3,25.0);
478 double t2[4][4], T2[4][4];
479 calt2(P1,P2,t2);
480 calt2(pR,P3,T2);
481 temp_PDF = 0;
482 for(int i=0; i<4; i++){
483 for(int j=0; j<4; j++){
484 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
485 }
486 }
487 }
488 v_re = temp_PDF*B[0]*B[1];
489 return v_re;
490}
491
492
493void EvtDsToKSKpi0::calEva(double* Ks0, double* Kc, double* Pi0, double *mass1, double *width1, double *amp, double *phase,int* g0,int* spin, int* modetype, int nstates, double & Result)
494{
495 //double Ks0[4], Pic[4], Pi0[4];
496 //double rRes = 9.0;
497 double P12[4], P23[4], P13[4];
498 double cof[2], amp_PDF[2], PDF[2];
499 double snpi, sck, sks0;
500 double s12, s13, s23;
501 for(int i=0; i<4; i++){
502 P12[i] = Kc[i] + Ks0[i];
503 P13[i] = Pi0[i] + Ks0[i];
504 P23[i] = Kc[i] + Pi0[i];
505 }
506 sck = SCADot(Kc,Kc);
507 snpi = SCADot(Pi0,Pi0);
508 sks0 = SCADot(Ks0,Ks0);
509 s12 = SCADot(P12,P12);
510 s13 = SCADot(P13,P13);
511 s23 = SCADot(P23,P23);
512 double pro[2], temp_PDF, amp_tmp[2],temp_PDF1 ,temp_PDF2,pro1[2],pro2[2];
513 double mass1sq;
514 amp_PDF[0] = 0;
515 amp_PDF[1] = 0;
516 PDF[0] = 0;
517 PDF[1] = 0;
518 amp_tmp[0] = 0;
519 amp_tmp[1] = 0;
520 for(int i=0; i<nstates; i++) {
521 amp_tmp[0] = 0;
522 amp_tmp[1] = 0;
523 mass1sq = mass1[i]*mass1[i];
524 cof[0] = amp[i]*cos(phase[i]);
525 cof[1] = amp[i]*sin(phase[i]);
526 temp_PDF = 0;
527
528 if(modetype[i] == 23){
529 temp_PDF = DDalitz(Kc, Pi0, Ks0, spin[i], mass1[i]);
530 if(g0[i]==1) propagatorRBW(mass1sq,mass1[i],width1[i],s23,sck,snpi,rRes,spin[i],pro);
531 if(g0[i]==2) KPiSLASS(s23,sck,snpi,pro);
532 if(g0[i]==0){
533 pro[0] = 1;
534 pro[1] = 0;
535 }
536 amp_tmp[0] = temp_PDF*pro[0];
537 amp_tmp[1] = temp_PDF*pro[1];
538
539}
540 if(modetype[i] == 12){
541 temp_PDF = DDalitz(Ks0, Kc, Pi0, spin[i], mass1[i]);
542 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s12,sks0,sck,rRes,spin[i],pro);
543 if(g0[i]==2) KPiSLASS(s12,sks0,sck,pro);
544 if(g0[i]==3) propagatorFlatte(mass1[i],width1[i],s12,sks0,sck,1,pro);//Only for a0(980)
545 if(g0[i]==0){
546 pro[0] = 1;
547 pro[1] = 0;
548 }
549 amp_tmp[0] = temp_PDF*pro[0];
550 amp_tmp[1] = temp_PDF*pro[1];
551 }
552 if(modetype[i] == 13){
553 temp_PDF = DDalitz(Ks0, Pi0, Kc, spin[i], mass1[i]);
554 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s13,sks0,snpi,rRes,spin[i],pro);
555 if(g0[i]==2) KPiSLASS(s13,sks0,snpi,pro);
556 if(g0[i]==0){
557 pro[0] = 1;
558 pro[1] = 0;
559 }
560 amp_tmp[0] = temp_PDF*pro[0];
561 amp_tmp[1] = temp_PDF*pro[1];
562 }
563 Com_Multi(amp_tmp,cof,amp_PDF);
564 PDF[0] += amp_PDF[0];
565 PDF[1] += amp_PDF[1];
566}
567 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
568// cout<<"value = "<<value<<endl;
569 if(value <=0) value = 1e-20;
570 Result = value;
571// cout<<"value = "<<value<<" Result = "<<Result<<endl;
572}
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")
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
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
EvtDsToKSKpi0::decay
void decay(EvtParticle *p)
Definition EvtDsToKSKpi0.cc:128
EvtDsToKSKpi0::clone
EvtDecayBase * clone()
Definition EvtDsToKSKpi0.cc:38
EvtDsToKSKpi0::~EvtDsToKSKpi0
virtual ~EvtDsToKSKpi0()
Definition EvtDsToKSKpi0.cc:32
EvtDsToKSKpi0::getName
void getName(std::string &name)
Definition EvtDsToKSKpi0.cc: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