BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKpipi0.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: EvtD0ToKpipi0.cc
11// the necessary file: EvtD0ToKpipi0.hh
12//
13// Description: D0 -> K- pi+ pi0, see https://indico.ihep.ac.cn/event/17067/contributions/116851/attachments/63347/73471/charm_9.6.pdf
14//
15// Modification history:
16//
17// Liaoyuan Dong Sat Apr 29 03:10:55 CST 2023
18//
19//------------------------------------------------------------------------
20#include "EvtGenBase/EvtPatches.hh"
21#include "EvtGenBase/EvtParticle.hh"
22#include "EvtGenBase/EvtGenKine.hh"
23#include "EvtGenBase/EvtPDL.hh"
24#include "EvtGenBase/EvtReport.hh"
25#include "EvtGenModels/EvtD0ToKpipi0.hh"
26#include "EvtGenBase/EvtComplex.hh"
27#include "EvtGenBase/EvtDecayTable.hh"
28#include <stdlib.h>
29#include <iostream>
30#include <cmath>
31using namespace std;
32
33EvtD0ToKpipi0::~EvtD0ToKpipi0() {}
34
35void EvtD0ToKpipi0::getName(std::string& model_name){
36 model_name="D0ToKpipi0";
37}
38
42
43void EvtD0ToKpipi0::init(){
44 checkNArg(0);
45 checkNDaug(3);
46 checkSpinParent(EvtSpinType::SCALAR);
47
48 mrho_770 = 0.76519, mrho_1450 = 1.465, mKn892 = 0.89555, mKp892 = 0.89167, mK0_1430 = 1.466, mK2_1430 = 1.432, mKn_1680 = 1.717;
49 Grho_770 = 0.150, Grho_1450 = 0.310, GKn892 = 0.0473, GKp892 = 0.0514, GK0_1430 = 0.174, GK2_1430 = 0.109, GKn_1680 = 0.322;
50
51 rho[0] = 1; //rho770-pipi0
52 rho[1] = 1.5285; //rho1450-pipi0
53 rho[2] = 0.37736; //K*892p-Kpi
54 rho[3] = 0.34628; //K*982n-Kpi0
55 rho[4] = 0.34568; //K0*1430n-Kpi
56 rho[5] = 0.29595; //K0*1430p-Kpi0
57 rho[6] = 0.33434; //K2*1430n-Kpi
58 rho[7] = 0.34835; //K2*1430p-Kpi0
59 rho[8] = 0.89301; //K*1680n-Kpi
60 rho[9] = 0.33549; //K*1680p-Kpi
61 rho[10] = 0.90975; //res-S
62 rho[11] = 0.22014; //res-Kpip
63 rho[12] = 0.82374; //res-Kpi0p
64
65 phi[0] = 0; //rho770
66 phi[1] = -3.6546; //rho1450-pipi0
67 phi[2] = 2.9582; //K*892p-Kpi
68 phi[3] = 0.17931; //K*982n-Kpi0
69 phi[4] = -5.8211; //K0*1430n-Kpi
70 phi[5] = 4.3925; //K0*1430p-Kpi0
71 phi[6] = 5.6749; //K2*1430n-Kpi
72 phi[7] = -3.6710; //K2*1430p-Kpi0
73 phi[8] = -4.9602; //K*1680n-Kpi
74 phi[9] = -5.2493; //K*1680p-Kpi
75 phi[10] = -2.093; //res-S
76 phi[11] = -6.6026; //res-Kpip
77 phi[12] = 10.376; //res-Kpi0p
78
79 spin[0] = 1;
80 spin[1] = 1;
81 spin[2] = 1;
82 spin[3] = 1;
83 spin[4] = 0;
84 spin[5] = 0;
85 spin[6] = 2;
86 spin[7] = 2;
87 spin[8] = 1;
88 spin[9] = 1;
89 spin[10] = 0;
90 spin[11] = 1;
91 spin[12] = 1;
92
93 modetype[0] = 23;
94 modetype[1] = 23;
95 modetype[2] = 13;
96 modetype[3] = 12;
97 modetype[4] = 12;
98 modetype[5] = 13;
99 modetype[6] = 12;
100 modetype[7] = 13;
101 modetype[8] = 12;
102 modetype[9] = 13;
103 modetype[10] = 12;
104 modetype[11] = 12;
105 modetype[12] = 13;
106
107/*
108 std::cout << "EvtD0ToKpipi0 (May 01, 2023) ==> Initialization" << std::endl;
109 for (int i=0; i<13; i++) {
110 cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
111 }
112*/
113 mass_Ks=0.4977;
114 mass_Eta=0.547862;
115 mD = 1.86965;
116 rD = 5;
117 metap = 0.95778;
118 mkstr = 0.89594;
119 mk0 = 0.497611;
120 mass_Kaon = 0.493677;
121 mass_Pion = 0.13957;
122 mass_Pi0 = 0.1349768;
123 math_pi = 3.1415926;
124 pi = 3.1415926;
125 mpi = 0.13957;
126
127 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
128 for (int i=0; i<4; i++) {
129 for (int j=0; j<4; j++) {
130 G[i][j] = GG[i][j];
131 }
132 }
133}
134
135void EvtD0ToKpipi0::initProbMax() {
136 setProbMax(161.0);//setProbMax;
137}
138
139void EvtD0ToKpipi0::decay( EvtParticle *p){
140 //-----------for max value------------------
141 // double maxprob = 0.0;
142 // for(int ir=0;ir<=60000000;ir++){
143 // p->initializePhaseSpace(getNDaug(),getDaugs());
144 // EvtVector4R ks0 = p->getDaug(0)->getP4();
145 // EvtVector4R pic = p->getDaug(1)->getP4();
146 // EvtVector4R pi0 = p->getDaug(2)->getP4();
147
148 // double Ks0[4],Pic[4],Pi0[4];
149 // Ks0[0] = ks0.get(0); Pic[0] = pic.get(0); Pi0[0] = pi0.get(0);
150 // Ks0[1] = ks0.get(1); Pic[1] = pic.get(1); Pi0[1] = pi0.get(1);
151 // Ks0[2] = ks0.get(2); Pic[2] = pic.get(2); Pi0[2] = pi0.get(2);
152 // Ks0[3] = ks0.get(3); Pic[3] = pic.get(3); Pi0[3] = pi0.get(3);
153 // double value;
154 // calPDF(Ks0, Pic, Pi0, value);
155 // if(value>maxprob) {
156 // maxprob=value;
157 // std::cout << "Max PDF = " << ir << " prob= " << value << std::endl;
158 // }
159 // }
160 // std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
161 // return;
162 //-----------------------------------------------
163 p->initializePhaseSpace(getNDaug(),getDaugs());
164 EvtVector4R ks0 = p->getDaug(0)->getP4();
165 EvtVector4R pic = p->getDaug(1)->getP4();
166 EvtVector4R pi0 = p->getDaug(2)->getP4();
167
168 double Ks0[4],Pic[4],Pi0[4];
169 Ks0[0] = ks0.get(0); Pic[0] = pic.get(0); Pi0[0] = pi0.get(0);
170 Ks0[1] = ks0.get(1); Pic[1] = pic.get(1); Pi0[1] = pi0.get(1);
171 Ks0[2] = ks0.get(2); Pic[2] = pic.get(2); Pi0[2] = pi0.get(2);
172 Ks0[3] = ks0.get(3); Pic[3] = pic.get(3); Pi0[3] = pi0.get(3);
173
174 //Ks0[0] =0.665612; Pic[0] =0.655979; Pi0[0] = 0.566459;
175 //Ks0[1] =0.121309; Pic[1] =0.185548; Pi0[1] = -0.140599;
176 //Ks0[2] =0.289703; Pic[2] =-0.499909;Pi0[2] = 0.017034;
177 //Ks0[3] =-0.31393; Pic[3] =0.355657; Pi0[3] = 0.022672;
178
179 double value;
180 calPDF(Ks0, Pic, Pi0, value);
181
182 setProb(value);
183 return;
184}
185
186double EvtD0ToKpipi0::calPDF(double Ks0[], double Pic[], double Pi0[], double & Result) {
187 double cof[2], amp_tmp1[2], amp_tmp2[2], amp_tmp[2], amp_PDF[2], PDF[2];
188 double flag[3], mass_R[2], width_R[2];
189 double pro[2];
190 double sa[3], sb[3], sc[3], B[3];
191 double t1D[4], t1V[4], t1A[4];
192 double pS[4], pV[4], pD[4], pA[4];
193 double pro1[2], pro2[2],proKPi_S[2];
194
195 double rD2 = 25.0;
196 double rRes2 = 9.0;
197 double rRes = 9.0;
198 double mass1[13] ={ mrho_770, mrho_1450, mKp892, mKn892, mK0_1430, mK0_1430, mK2_1430, mK2_1430, mKn_1680, mKn_1680, 1.414, 1.414, 1.414};
199 double width1[13]={ Grho_770, Grho_1450, GKp892, GKn892, GK0_1430, GK0_1430, GK2_1430, GK2_1430, GKn_1680, GKn_1680, 0.232, 0.232, 0.232};
200 double g0[13] ={ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0};
201 double temp_PDF = 0;
202 double P12[4], P23[4], P13[4];
203 double scpi, snpi, sks0;
204 double s12, s13, s23;
205 for(int i=0; i<4; i++){
206 P12[i] = Pic[i] + Ks0[i];
207 P13[i] = Pi0[i] + Ks0[i];
208 P23[i] = Pic[i] + Pi0[i];
209 }
210 scpi = SCADot(Pic,Pic);
211 snpi = SCADot(Pi0,Pi0);
212 sks0 = SCADot(Ks0,Ks0);
213 s12 = SCADot(P12,P12);
214 s13 = SCADot(P13,P13);
215 s23 = SCADot(P23,P23);
216 double mass1sq;
217 double Amp_KPiS[2];
218 amp_PDF[0] = 0;
219 amp_PDF[1] = 0;
220 PDF[0] = 0;
221 PDF[1] = 0;
222
223 for(int i=0; i<13; i++) {
224 amp_tmp[0] = 0;
225 amp_tmp[1] = 0;
226 mass1sq = mass1[i]*mass1[i];
227 cof[0] = rho[i]*cos(phi[i]);
228 cof[1] = rho[i]*sin(phi[i]);
229 temp_PDF = 0;
230 if(modetype[i] == 23){
231 temp_PDF = DDalitz( Pic, Pi0, Ks0, spin[i], mass1[i]);
232 if(g0[i]==1) propagatorRBW(mass1sq,mass1[i],width1[i],s23,scpi,snpi,rRes,spin[i],pro);
233 if(g0[i]==2) KPiSLASS(s23,scpi,snpi,pro);
234 if(g0[i]==3) propagatorFlatte(mass1[i],width1[i],s23,scpi,snpi,0,pro);
235 if(g0[i]==0){
236 pro[0] = 1;
237 pro[1] = 0;
238 }
239 amp_tmp[0] = temp_PDF*pro[0];
240 amp_tmp[1] = temp_PDF*pro[1];
241 }
242 if(modetype[i] == 12){
243 temp_PDF = DDalitz(Ks0, Pic, Pi0, spin[i], mass1[i]);
244 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s12,sks0,scpi,rRes,spin[i],pro);
245 if(g0[i]==2) {propagatorFlatte(mass1[i],width1[i],s12,sks0,scpi,1,pro);
246 pro[0]=pro[0]*(0.01+0.990*0.990)/(0.01+s12);
247 pro[1]=pro[1]*(0.01+0.990*0.990)/(0.01+s12);
248 }
249 if(g0[i]==3) propagatorFlatte(mass1[i],width1[i],s12,sks0,scpi,1,pro);//Only for a0(980)
250
251 if(g0[i]==0){
252 pro[0] = 1;
253 pro[1] = 0;
254 }
255 amp_tmp[0] = temp_PDF*pro[0];
256 amp_tmp[1] = temp_PDF*pro[1];
257 }
258 if(modetype[i] == 13){
259 temp_PDF = DDalitz(Ks0, Pi0, Pic, spin[i], mass1[i]);
260 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s13,sks0,snpi,rRes,spin[i],pro);
261 if(g0[i]==2) { // K*1430 Flatte
262 double skm2[2]={sks0, mass_Ks *mass_Ks};
263 double spi2[2]={snpi, mass_Eta *mass_Eta};
264 propagatorKstr1430(mass1[i],s13,skm2,spi2,pro);
265 };
266 if(g0[i]==3) propagatorFlatte(mass1[i],width1[i],s13,sks0,snpi,1,pro);//Only for a0(980)
267 if(g0[i]==0){
268 pro[0] = 1;
269 pro[1] = 0;
270 }
271 amp_tmp[0] = temp_PDF*pro[0];
272 amp_tmp[1] = temp_PDF*pro[1];
273 }
274 if(modetype[i] == 132){
275 KPiSLASS(s12,sks0,scpi,Amp_KPiS);
276 amp_tmp[0] = Amp_KPiS[0];
277 amp_tmp[1] = Amp_KPiS[1];
278 }
279 Com_Multi(amp_tmp,cof,amp_PDF);
280 PDF[0] += amp_PDF[0];
281 PDF[1] += amp_PDF[1];
282 }
283
284
285 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
286
287 Result = value;
288
289}
290
291void EvtD0ToKpipi0:: Com_Multi(double a1[2], double a2[2], double res[2])
292{
293 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
294 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
295}
296void EvtD0ToKpipi0:: Com_Divide(double a1[2], double a2[2], double res[2])
297{
298 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
299 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
300 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
301}
302//------------base---------------------------------
303double EvtD0ToKpipi0:: SCADot(double a1[4], double a2[4])
304{
305 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
306 return _cal;
307}
308
309double EvtD0ToKpipi0:: barrier(int l, double sa, double sb, double sc, double r, double mass)
310{
311 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
312 if(q < 0) q = 1e-16;
313 double z;
314 z = q*r*r;
315 double sa0;
316 sa0 = mass*mass;
317 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
318 if(q0 < 0) q0 = 1e-16;
319 double z0 = q0*r*r;
320 double F = 0.0;
321 if(l == 0) F = 1;
322 if(l == 1) F = sqrt((1+z0)/(1+z));
323 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
324 return F;
325}
326
327double EvtD0ToKpipi0:: Barrier(int l, double sa, double sb, double sc, double r2)
328{
329 double F;
330 double tmp = sa+sb-sc;
331 double q = 0.25*tmp*tmp/sa-sb;
332 if (q < 0) q = 1e-16;
333 double z = q*r2;
334 if (l==1) {
335 F = sqrt(2.0*z/(1.0+z));
336 }
337 else if (l==2) {
338 double z2 = z*z;
339 F = sqrt(13.0*z2/(9.0+3.0*z+z2));
340 } else {
341 F = 1.0;
342 }
343 return F;
344}
345
346//------------------spin-------------------------------------------
347void EvtD0ToKpipi0:: calt1(double daug1[4], double daug2[4], double t1[4])
348{
349 double p, pq, tmp;
350 double pa[4], qa[4];
351 for(int i=0; i<4; i++) {
352 pa[i] = daug1[i] + daug2[i];
353 qa[i] = daug1[i] - daug2[i];
354 }
355 p = SCADot(pa,pa);
356 pq = SCADot(pa,qa);
357 tmp = pq/p;
358 for(int i=0; i<4; i++) {
359 t1[i] = qa[i] - tmp*pa[i];
360 }
361}
362void EvtD0ToKpipi0:: calt2(double daug1[4], double daug2[4], double t2[4][4])
363{
364 double p, r;
365 double pa[4], t1[4];
366 calt1(daug1,daug2,t1);
367 r = SCADot(t1,t1)/3.0;
368 for(int i=0; i<4; i++) {
369 pa[i] = daug1[i] + daug2[i];
370 }
371 p = SCADot(pa,pa);
372 for(int i=0; i<4; i++) {
373 for(int j=0; j<4; j++) {
374 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
375 }
376 }
377}
378//-------------------prop--------------------------------------------
379void EvtD0ToKpipi0:: propagator(double mass2, double mass, double width, double sx, double prop[2])
380{
381 double a[2], b[2];
382 a[0] = 1;
383 a[1] = 0;
384 b[0] = mass2-sx;
385 b[1] = -mass*width;
386 Com_Divide(a,b,prop);
387}
388double EvtD0ToKpipi0:: wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
389{
390 double widm = 0.;
391 double m = sqrt(sa);
392 double tmp = sb-sc;
393 double tmp1 = sa+tmp;
394 double q = 0.25*tmp1*tmp1/sa-sb;
395 if(q<0) q = 1e-16;
396 double tmp2 = mass2+tmp;
397 double q0 = 0.25*tmp2*tmp2/mass2-sb;
398 if(q0<0) q0 = 1e-16;
399 double z = q*r2;
400 double z0 = q0*r2;
401 double t = q/q0;
402 if(l == 0) {widm = sqrt(t)*mass/m;}
403 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
404 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
405 return widm;
406}
407double EvtD0ToKpipi0:: widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
408{
409 double widm = 0.;
410 double m = sqrt(sa);
411 double tmp = sb-sc;
412 double tmp1 = sa+tmp;
413 double q = 0.25*tmp1*tmp1/sa-sb;
414 if(q<0) q = 1e-16;
415 double tmp2 = mass2+tmp;
416 double q0 = 0.25*tmp2*tmp2/mass2-sb;
417 if(q0<0) q0 = 1e-16;
418 double z = q*r2;
419 double z0 = q0*r2;
420 double F = (1+z0)/(1+z);
421 double t = q/q0;
422 widm = t*sqrt(t)*mass/m*F;
423 return widm;
424}
425void EvtD0ToKpipi0:: propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
426{
427 double a[2], b[2];
428 a[0] = 1;
429 a[1] = 0;
430 double q=0.25*(sa+sb-sc)*(sa+sb-sc)/sa-sb;
431
432 b[0] = mass2-sa;
433 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
434 Com_Divide(a,b,prop);
435}
436
437
438void EvtD0ToKpipi0:: propagatorFlatte(double mass, double width, double sa, double sb, double sc, int r, double prop[2]){
439 double q, qKsK,qetapi;
440 // double qKsK,qetapi;
441 double rhoab[2], rhoKsK[2];
442 q = 0.25*(sa+sb-sc)*(sa+sb-sc)/sa-sb;
443 qetapi=0.25*(sa+0.547862*0.547862-0.13957039*0.13957039)*(sa+0.547862*0.547862-0.13957039*0.13957039)/sa-0.547862*0.547862;
444 if(r == 0) qKsK = 0.25*sa - 0.49368*0.49368;
445 if(r == 1) qKsK = 0.25*(sa+sb-sc)*(sa+sb-sc)/sa - sb;
446 if(qetapi>0){
447 rhoab[0] = 2*sqrt(qetapi/sa);
448 rhoab[1] = 0;
449 }
450 if(qetapi<0){
451 rhoab[0] = 0;
452 rhoab[1] = 2*sqrt(-qetapi/sa);
453 }
454 if(qKsK>0){
455 rhoKsK[0] = 2*sqrt(qKsK/sa);
456 rhoKsK[1] = 0;
457 }
458 if(qKsK<0){
459 rhoKsK[0] = 0;
460 rhoKsK[1] = 2*sqrt(-qKsK/sa);
461 }
462 double a[2], b[2];
463 a[0] = 1;
464 a[1] = 0;
465 b[0] = mass*mass - sa + 0.341*rhoab[1] + 0.892*0.341*rhoKsK[1];
466 b[1] = - (0.341*rhoab[0] + 0.892*0.341*rhoKsK[0]);
467 Com_Divide(a,b,prop);
468}
469void EvtD0ToKpipi0:: PiPiSWASS(double sa, double sb, double sc, double prop[2]) {
470 double tmp = sb-sc;
471 double tmp2 = sa+tmp;
472 double qs = 0.25*tmp2*tmp2/sa-sb;
473 double q = sqrt(qs);
474 double a0 = -0.11/mass_Pion;
475 prop[0] = 1/(1+a0*a0*q*q);
476 prop[1] = a0*q/(1+a0*a0*q*q);
477}
478void EvtD0ToKpipi0:: KPiSLASS(double sa, double sb, double sc, double prop[2]) {
479 const double m1430 = 1.441;
480 const double sa0 = 1.441*1.441; // m1430*m1430;
481 const double w1430 = 0.193;
482 const double Lass1 = 0.25/sa0;
483 double tmp = sb-sc;
484 double tmp1 = sa0+tmp;
485 double q0 = Lass1*tmp1*tmp1-sb;
486 if(q0<0) q0 = 1e-16;
487 double tmp2 = sa+tmp;
488 double qs = 0.25*tmp2*tmp2/sa-sb;
489 double q = sqrt(qs);
490 double width = w1430*q*m1430/sqrt(sa*q0);
491 double temp_R = atan(m1430*width/(sa0-sa));
492 if(temp_R<0) temp_R += math_pi;
493 double deltaR = -1.915 + temp_R; //fiR=-109.7
494 double temp_F = atan(0.226*q/(2.0-3.819*qs)); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
495 if(temp_F<0) temp_F += math_pi;
496 double deltaF = 0.002 + temp_F; //fiF=0.1
497 double deltaS = deltaR + 2.0*deltaF;
498 double t1 = 0.96*sin(deltaF);
499 double t2 = sin(deltaR);
500 double CF[2], CS[2];
501 CF[0] = cos(deltaF);
502 CF[1] = sin(deltaF);
503 CS[0] = cos(deltaS);
504 CS[1] = sin(deltaS);
505 prop[0] = t1*CF[0] + t2*CS[0];
506 prop[1] = t1*CF[1] + t2*CS[1];
507}
508
509void EvtD0ToKpipi0:: Flatte_rhoab(double sa, double sb, double sc, double rho[2]){
510 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
511 if(q>0) {
512 rho[0]=2* sqrt(q/sa);
513 rho[1]=0;
514 }
515 else if(q<0){
516 rho[0]=0;
517 rho[1]=2*sqrt(-q/sa);
518 }
519}
520
521void EvtD0ToKpipi0:: propagatorKstr1430(double mass, double sx, double *sb, double *sc, double prop[2]) //K*1430 Flatte
522{
523 double unit[2]={1.0};
524 double ci[2]={0,1};
525 double rho1[2];
526 Flatte_rhoab(sx,sb[0],sc[0],rho1);
527 double rho2[2];
528 Flatte_rhoab(sx,sb[1],sc[1],rho2);
529 double gKPi_Kstr1430=0.2990, gEtaPK_Kstr1430=0.0529;
530 double tmp1[2]={gKPi_Kstr1430,0};
531 double tmp11[2];
532 double tmp2[2]={gEtaPK_Kstr1430,0};
533 double tmp22[2];
534 Com_Multi(tmp1,rho1,tmp11);
535 Com_Multi(tmp2,rho2,tmp22);
536 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
537 double tmp31[2];
538 Com_Multi(tmp3, ci,tmp31);
539 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
540 Com_Divide( unit,tmp4, prop);
541}
542double EvtD0ToKpipi0:: DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
543 double pR[4], pD[4];
544 double temp_PDF, v_re;
545 temp_PDF = 0.0;
546 v_re = 0.0;
547 double B[2], s1, s2, s3, sR, sD;
548 for(int i=0; i<4; i++){
549 pR[i] = P1[i] + P2[i];
550 pD[i] = pR[i] + P3[i];
551 }
552 s1 = SCADot(P1,P1);
553 s2 = SCADot(P2,P2);
554 s3 = SCADot(P3,P3);
555 sR = SCADot(pR,pR);
556 sD = SCADot(pD,pD);
557 int G[4][4];
558 for(int i=0; i!=4; i++){
559 for(int j=0; j!=4; j++){
560 if(i==j){
561 if(i==0) G[i][j] = 1;
562 else G[i][j] = -1;
563 }
564 else G[i][j] = 0;
565 }
566 }
567 if(Ang == 0){
568 B[0] = 1;
569 B[1] = 1;
570 temp_PDF = 1;
571 }
572 if(Ang == 1){
573 B[0] = barrier(1,sR,s1,s2,3.0,mass);
574 B[1] = barrier(1,sD,sR,s3,5.0,1.864);
575 double t1[4], T1[4];
576 calt1(P1,P2,t1);
577 calt1(pR,P3,T1);
578 temp_PDF = 0;
579 for(int i=0; i<4; i++){
580 temp_PDF += t1[i]*T1[i]*G[i][i];
581 }
582 }
583 if(Ang == 2){
584 B[0] = barrier(2,sR,s1,s2,3.0,mass);
585 B[1] = barrier(2,sD,sR,s3,5.0,1.864);
586 double t2[4][4], T2[4][4];
587 calt2(P1,P2,t2);
588 calt2(pR,P3,T2);
589 temp_PDF = 0;
590 for(int i=0; i<4; i++){
591 for(int j=0; j<4; j++){
592 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
593 }
594 }
595 }
596 v_re = temp_PDF*B[0]*B[1];
597 return v_re;
598}
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
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
EvtD0ToKpipi0::~EvtD0ToKpipi0
virtual ~EvtD0ToKpipi0()
Definition EvtD0ToKpipi0.cc:33
EvtD0ToKpipi0::clone
EvtDecayBase * clone()
Definition EvtD0ToKpipi0.cc:39
EvtD0ToKpipi0::decay
void decay(EvtParticle *p)
Definition EvtD0ToKpipi0.cc:139
EvtD0ToKpipi0::getName
void getName(std::string &name)
Definition EvtD0ToKpipi0.cc:35
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