BOSS 7.1.0
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 Oct.15, 2022 Module created
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
39EvtDecayBase* EvtD0ToKpipi0::clone(){
40 return new EvtD0ToKpipi0;
41}
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.32950;//rho1450-pipi0
53 rho[2] = 0.34447;//K*892n-Kpi
54 rho[3] = 0.37795;//K*982p-Kpi0
55 rho[4] = 0.15531;//K0*1430n-Kpi
56 rho[5] = 0.27334;//K0*1430p-Kpi0
57 rho[6] = 0.34413;//K2*1430n-Kpi
58 rho[7] = 0.34723;//K2*1430p-Kpi0
59 rho[8] = 0.38826;//K*1680n-Kpi
60 rho[9] = 0.86076;//KpiSW
61
62 phi[0] = 0;//rho770
63 phi[1] = -3.5053; //rho1450-pipi0
64 phi[2] = 0.12764; //K*892n-Kpi
65 phi[3] = 9.1910; //K*982p-Kpi0
66 phi[4] = -5.4878; //K0*1430n-Kpi
67 phi[5] = 4.4317; //K0*1430p-Kpi0
68 phi[6] = -0.69265;//K2*1430n-Kpi
69 phi[7] = -3.5656; //K2*1430p-Kpi0
70 phi[8] = -4.3046; //K*1680n-Kpi
71 phi[9] = 8.8127; //KpiSW
72
73 spin[0] = 1;
74 spin[1] = 1;
75 spin[2] = 1;
76 spin[3] = 1;
77 spin[4] = 0;
78 spin[5] = 0;
79 spin[6] = 2;
80 spin[7] = 2;
81 spin[8] = 1;
82 spin[9] = 0;
83
84 modetype[0] = 23;
85 modetype[1] = 23;
86 modetype[2] = 12;
87 modetype[3] = 13;
88 modetype[4] = 12;
89 modetype[5] = 13;
90 modetype[6] = 12;
91 modetype[7] = 13;
92 modetype[8] = 12;
93 modetype[9] = 132;
94
95 //std::cout << "Initializing EvtD0ToKpipi0" << std::endl;
96 //for (int i=0; i<10; i++) {
97 // cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
98 //}
99
100 mass_Ks=0.4977;
101 mass_Eta=0.547862;
102 mD = 1.86965;
103 rD = 5;
104 metap = 0.95778;
105 mkstr = 0.89594;
106 mk0 = 0.497611;
107 mass_Kaon = 0.493677;
108 mass_Pion = 0.13957;
109 mass_Pi0 = 0.1349768;
110 math_pi = 3.1415926;
111
112 pi = 3.1415926;
113 mpi = 0.13957;
114 g1 = 0.5468;
115 g2 = 0.23;
116
117 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
118 for (int i=0; i<4; i++) {
119 for (int j=0; j<4; j++) {
120 G[i][j] = GG[i][j];
121 }
122 }
123}
124
125void EvtD0ToKpipi0::initProbMax() {
126 setProbMax(144.5);//setProbMax;
127}
128
129void EvtD0ToKpipi0::decay( EvtParticle *p){
130 //-----------for max value------------------
131 /* double maxprob = 0.0;
132 for(int ir=0;ir<=60000000;ir++){
133 p->initializePhaseSpace(getNDaug(),getDaugs());
134 EvtVector4R ks0 = p->getDaug(0)->getP4();
135 EvtVector4R pic = p->getDaug(1)->getP4();
136 EvtVector4R pi0 = p->getDaug(2)->getP4();
137
138 double Ks0[4],Pic[4],Pi0[4];
139 Ks0[0] = ks0.get(0); Pic[0] = pic.get(0); Pi0[0] = pi0.get(0);
140 Ks0[1] = ks0.get(1); Pic[1] = pic.get(1); Pi0[1] = pi0.get(1);
141 Ks0[2] = ks0.get(2); Pic[2] = pic.get(2); Pi0[2] = pi0.get(2);
142 Ks0[3] = ks0.get(3); Pic[3] = pic.get(3); Pi0[3] = pi0.get(3);
143 double value;
144 calPDF(Ks0, Pic, Pi0, value);
145 if(value>maxprob) {
146 maxprob=value;
147 std::cout << "Max PDF = " << ir << " prob= " << value << std::endl;
148 }
149 }
150 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
151 return;*/
152 //-----------------------------------------------
153 p->initializePhaseSpace(getNDaug(),getDaugs());
154 EvtVector4R ks0 = p->getDaug(0)->getP4();
155 EvtVector4R pic = p->getDaug(1)->getP4();
156 EvtVector4R pi0 = p->getDaug(2)->getP4();
157
158 double Ks0[4],Pic[4],Pi0[4];
159 Ks0[0] = ks0.get(0); Pic[0] = pic.get(0); Pi0[0] = pi0.get(0);
160 Ks0[1] = ks0.get(1); Pic[1] = pic.get(1); Pi0[1] = pi0.get(1);
161 Ks0[2] = ks0.get(2); Pic[2] = pic.get(2); Pi0[2] = pi0.get(2);
162 Ks0[3] = ks0.get(3); Pic[3] = pic.get(3); Pi0[3] = pi0.get(3);
163
164 double value;
165 calPDF(Ks0, Pic, Pi0, value);
166 setProb(value);
167 return;
168}
169
170double EvtD0ToKpipi0::calPDF(double Ks0[], double Pic[], double Pi0[], double & Result) {
171 double cof[2], amp_tmp1[2], amp_tmp2[2], amp_tmp[2], amp_PDF[2], PDF[2];
172 double flag[3], mass_R[2], width_R[2];
173 double pro[2];
174 double sa[3], sb[3], sc[3], B[3];
175 double t1D[4], t1V[4], t1A[4];
176 double pS[4], pV[4], pD[4], pA[4];
177 double pro1[2], pro2[2],proKPi_S[2];
178
179 double rD2 = 25.0;
180 double rRes2 = 9.0;
181 double rRes = 9.0;
182 double mass1[10] ={ mrho_770, mrho_1450, mKn892, mKp892, mK0_1430, mK0_1430, mK2_1430, mK2_1430, mKn_1680, 1.0};
183 double width1[10]={ Grho_770, Grho_1450, GKn892, GKp892, GK0_1430, GK0_1430, GK2_1430, GK2_1430, GKn_1680, 1.0};
184 double g0[10] ={ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
185 double temp_PDF = 0;
186 double P12[4], P23[4], P13[4];
187 double scpi, snpi, sks0;
188 double s12, s13, s23;
189 for(int i=0; i<4; i++){
190 P12[i] = Pic[i] + Ks0[i];
191 P13[i] = Pi0[i] + Ks0[i];
192 P23[i] = Pic[i] + Pi0[i];
193 }
194 scpi = SCADot(Pic,Pic);
195 snpi = SCADot(Pi0,Pi0);
196 sks0 = SCADot(Ks0,Ks0);
197 s12 = SCADot(P12,P12);
198 s13 = SCADot(P13,P13);
199 s23 = SCADot(P23,P23);
200 double mass1sq;
201 double Amp_KPiS[2];
202 amp_PDF[0] = 0;
203 amp_PDF[1] = 0;
204 PDF[0] = 0;
205 PDF[1] = 0;
206
207 for(int i=0; i<10; i++) {
208 amp_tmp[0] = 0;
209 amp_tmp[1] = 0;
210 mass1sq = mass1[i]*mass1[i];
211 cof[0] = rho[i]*cos(phi[i]);
212 cof[1] = rho[i]*sin(phi[i]);
213 temp_PDF = 0;
214 if(modetype[i] == 23){
215 temp_PDF = DDalitz( Pic, Pi0, Ks0, spin[i], mass1[i]);
216 if(g0[i]==1) propagatorRBW(mass1sq,mass1[i],width1[i],s23,scpi,snpi,rRes,spin[i],pro);
217 if(g0[i]==2) KPiSLASS(s23,scpi,snpi,pro);
218 if(g0[i]==3) propagatorFlatte(mass1[i],width1[i],s23,scpi,snpi,0,pro);
219 if(g0[i]==0){
220 pro[0] = 1;
221 pro[1] = 0;
222 }
223 amp_tmp[0] = temp_PDF*pro[0];
224 amp_tmp[1] = temp_PDF*pro[1];
225 }
226 if(modetype[i] == 12){
227 temp_PDF = DDalitz(Ks0, Pic, Pi0, spin[i], mass1[i]);
228 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s12,sks0,scpi,rRes,spin[i],pro);
229 if(g0[i]==2) {propagatorFlatte(mass1[i],width1[i],s12,sks0,scpi,1,pro);
230 pro[0]=pro[0]*(0.01+0.990*0.990)/(0.01+s12);
231 pro[1]=pro[1]*(0.01+0.990*0.990)/(0.01+s12);
232 }
233 if(g0[i]==3) propagatorFlatte(mass1[i],width1[i],s12,sks0,scpi,1,pro);//Only for a0(980)
234
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] == 13){
243 temp_PDF = DDalitz(Ks0, Pi0, Pic, spin[i], mass1[i]);
244 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s13,sks0,snpi,rRes,spin[i],pro);
245 if(g0[i]==2) { // K*1430 Flatte
246 double skm2[2]={sks0, mass_Ks *mass_Ks};
247 double spi2[2]={snpi, mass_Eta *mass_Eta};
248 propagatorKstr1430(mass1[i],s13,skm2,spi2,pro);
249 };
250 if(g0[i]==3) propagatorFlatte(mass1[i],width1[i],s13,sks0,snpi,1,pro);//Only for a0(980)
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] == 132){
259 KPiSLASS(s12,sks0,scpi,Amp_KPiS);
260 amp_tmp[0] = Amp_KPiS[0];
261 amp_tmp[1] = Amp_KPiS[1];
262 }
263 Com_Multi(amp_tmp,cof,amp_PDF);
264 PDF[0] += amp_PDF[0];
265 PDF[1] += amp_PDF[1];
266}
267
268
269 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
270
271 Result = value;
272
273}
274
275void EvtD0ToKpipi0:: Com_Multi(double a1[2], double a2[2], double res[2])
276{
277 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
278 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
279}
280void EvtD0ToKpipi0:: Com_Divide(double a1[2], double a2[2], double res[2])
281{
282 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
283 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
284 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
285}
286//------------base---------------------------------
287double EvtD0ToKpipi0:: SCADot(double a1[4], double a2[4])
288{
289 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
290 return _cal;
291}
292
293double EvtD0ToKpipi0:: barrier(int l, double sa, double sb, double sc, double r, double mass)
294{
295 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
296 if(q < 0) q = 1e-16;
297 double z;
298 z = q*r*r;
299 double sa0;
300 sa0 = mass*mass;
301 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
302 if(q0 < 0) q0 = 1e-16;
303 double z0 = q0*r*r;
304 double F = 0.0;
305 if(l == 0) F = 1;
306 if(l == 1) F = sqrt((1+z0)/(1+z));
307 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
308 return F;
309}
310
311double EvtD0ToKpipi0:: Barrier(int l, double sa, double sb, double sc, double r2)
312{
313 double F;
314 double tmp = sa+sb-sc;
315 double q = 0.25*tmp*tmp/sa-sb;
316 if (q < 0) q = 1e-16;
317 double z = q*r2;
318 if (l==1) {
319 F = sqrt(2.0*z/(1.0+z));
320 }
321 else if (l==2) {
322 double z2 = z*z;
323 F = sqrt(13.0*z2/(9.0+3.0*z+z2));
324 } else {
325 F = 1.0;
326 }
327 return F;
328}
329
330//------------------spin-------------------------------------------
331void EvtD0ToKpipi0:: calt1(double daug1[4], double daug2[4], double t1[4])
332{
333 double p, pq, tmp;
334 double pa[4], qa[4];
335 for(int i=0; i<4; i++) {
336 pa[i] = daug1[i] + daug2[i];
337 qa[i] = daug1[i] - daug2[i];
338 }
339 p = SCADot(pa,pa);
340 pq = SCADot(pa,qa);
341 tmp = pq/p;
342 for(int i=0; i<4; i++) {
343 t1[i] = qa[i] - tmp*pa[i];
344 }
345}
346void EvtD0ToKpipi0:: calt2(double daug1[4], double daug2[4], double t2[4][4])
347{
348 double p, r;
349 double pa[4], t1[4];
350 calt1(daug1,daug2,t1);
351 r = SCADot(t1,t1)/3.0;
352 for(int i=0; i<4; i++) {
353 pa[i] = daug1[i] + daug2[i];
354 }
355 p = SCADot(pa,pa);
356 for(int i=0; i<4; i++) {
357 for(int j=0; j<4; j++) {
358 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
359 }
360 }
361}
362//-------------------prop--------------------------------------------
363void EvtD0ToKpipi0:: propagator(double mass2, double mass, double width, double sx, double prop[2])
364{
365 double a[2], b[2];
366 a[0] = 1;
367 a[1] = 0;
368 b[0] = mass2-sx;
369 b[1] = -mass*width;
370 Com_Divide(a,b,prop);
371}
372double EvtD0ToKpipi0:: wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
373{
374 double widm = 0.;
375 double m = sqrt(sa);
376 double tmp = sb-sc;
377 double tmp1 = sa+tmp;
378 double q = 0.25*tmp1*tmp1/sa-sb;
379 if(q<0) q = 1e-16;
380 double tmp2 = mass2+tmp;
381 double q0 = 0.25*tmp2*tmp2/mass2-sb;
382 if(q0<0) q0 = 1e-16;
383 double z = q*r2;
384 double z0 = q0*r2;
385 double t = q/q0;
386 if(l == 0) {widm = sqrt(t)*mass/m;}
387 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
388 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
389 return widm;
390}
391double EvtD0ToKpipi0:: widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
392{
393 double widm = 0.;
394 double m = sqrt(sa);
395 double tmp = sb-sc;
396 double tmp1 = sa+tmp;
397 double q = 0.25*tmp1*tmp1/sa-sb;
398 if(q<0) q = 1e-16;
399 double tmp2 = mass2+tmp;
400 double q0 = 0.25*tmp2*tmp2/mass2-sb;
401 if(q0<0) q0 = 1e-16;
402 double z = q*r2;
403 double z0 = q0*r2;
404 double F = (1+z0)/(1+z);
405 double t = q/q0;
406 widm = t*sqrt(t)*mass/m*F;
407 return widm;
408}
409void EvtD0ToKpipi0:: propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
410{
411 double a[2], b[2];
412 a[0] = 1;
413 a[1] = 0;
414 double q=0.25*(sa+sb-sc)*(sa+sb-sc)/sa-sb;
415
416 b[0] = mass2-sa;
417 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
418 Com_Divide(a,b,prop);
419}
420
421
422void EvtD0ToKpipi0:: propagatorFlatte(double mass, double width, double sa, double sb, double sc, int r, double prop[2]){
423 double q, qKsK,qetapi;
424// double qKsK,qetapi;
425 double rhoab[2], rhoKsK[2];
426 q = 0.25*(sa+sb-sc)*(sa+sb-sc)/sa-sb;
427 qetapi=0.25*(sa+0.547862-0.13957039)*(sa+0.547862-0.13957039)/sa-0.547862*0.547862;
428 if(r == 0) qKsK = 0.25*sa - 0.49368*0.49368;
429 if(r == 1) qKsK = 0.25*(sa+sb-sc)*(sa+sb-sc)/sa - sb;
430 if(qetapi>0){
431 rhoab[0] = 2*sqrt(qetapi/sa);
432 rhoab[1] = 0;
433 }
434 if(qetapi<0){
435 rhoab[0] = 0;
436 rhoab[1] = 2*sqrt(-qetapi/sa);
437 }
438 if(qKsK>0){
439 rhoKsK[0] = 2*sqrt(qKsK/sa);
440 rhoKsK[1] = 0;
441 }
442 if(qKsK<0){
443 rhoKsK[0] = 0;
444 rhoKsK[1] = 2*sqrt(-qKsK/sa);
445 }
446 double a[2], b[2];
447 a[0] = 1;
448 a[1] = 0;
449 b[0] = mass*mass - sa + 0.341*rhoab[1] + 0.892*0.341*rhoKsK[1];
450 b[1] = - (0.341*rhoab[0] + 0.892*0.341*rhoKsK[0]);
451 Com_Divide(a,b,prop);
452}
453void EvtD0ToKpipi0:: PiPiSWASS(double sa, double sb, double sc, double prop[2]) {
454 double tmp = sb-sc;
455 double tmp2 = sa+tmp;
456 double qs = 0.25*tmp2*tmp2/sa-sb;
457 double q = sqrt(qs);
458 double a0 = -0.11/mass_Pion;
459 prop[0] = 1/(1+a0*a0*q*q);
460 prop[1] = a0*q/(1+a0*a0*q*q);
461}
462void EvtD0ToKpipi0:: KPiSLASS(double sa, double sb, double sc, double prop[2]) {
463 const double m1430 = 1.441;
464 const double sa0 = 1.441*1.441; // m1430*m1430;
465 const double w1430 = 0.193;
466 const double Lass1 = 0.25/sa0;
467 double tmp = sb-sc;
468 double tmp1 = sa0+tmp;
469 double q0 = Lass1*tmp1*tmp1-sb;
470 if(q0<0) q0 = 1e-16;
471 double tmp2 = sa+tmp;
472 double qs = 0.25*tmp2*tmp2/sa-sb;
473 double q = sqrt(qs);
474 double width = w1430*q*m1430/sqrt(sa*q0);
475 double temp_R = atan(m1430*width/(sa0-sa));
476 if(temp_R<0) temp_R += math_pi;
477 double deltaR = -1.915 + temp_R; //fiR=-109.7
478 double temp_F = atan(0.226*q/(2.0-3.819*qs)); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
479 if(temp_F<0) temp_F += math_pi;
480 double deltaF = 0.002 + temp_F; //fiF=0.1
481 double deltaS = deltaR + 2.0*deltaF;
482 double t1 = 0.96*sin(deltaF);
483 double t2 = sin(deltaR);
484 double CF[2], CS[2];
485 CF[0] = cos(deltaF);
486 CF[1] = sin(deltaF);
487 CS[0] = cos(deltaS);
488 CS[1] = sin(deltaS);
489 prop[0] = t1*CF[0] + t2*CS[0];
490 prop[1] = t1*CF[1] + t2*CS[1];
491}
492
493void EvtD0ToKpipi0:: Flatte_rhoab(double sa, double sb, double sc, double rho[2]){
494 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
495 if(q>0) {
496 rho[0]=2* sqrt(q/sa);
497 rho[1]=0;
498 }
499 else if(q<0){
500 rho[0]=0;
501 rho[1]=2*sqrt(-q/sa);
502 }
503}
504
505void EvtD0ToKpipi0:: propagatorKstr1430(double mass, double sx, double *sb, double *sc, double prop[2]) //K*1430 Flatte
506{
507 double unit[2]={1.0};
508 double ci[2]={0,1};
509 double rho1[2];
510 Flatte_rhoab(sx,sb[0],sc[0],rho1);
511 double rho2[2];
512 Flatte_rhoab(sx,sb[1],sc[1],rho2);
513 double gKPi_Kstr1430=0.2990, gEtaPK_Kstr1430=0.0529;
514 double tmp1[2]={gKPi_Kstr1430,0};
515 double tmp11[2];
516 double tmp2[2]={gEtaPK_Kstr1430,0};
517 double tmp22[2];
518 Com_Multi(tmp1,rho1,tmp11);
519 Com_Multi(tmp2,rho2,tmp22);
520 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
521 double tmp31[2];
522 Com_Multi(tmp3, ci,tmp31);
523 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
524 Com_Divide( unit,tmp4, prop);
525}
526double EvtD0ToKpipi0:: DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
527 double pR[4], pD[4];
528 double temp_PDF, v_re;
529 temp_PDF = 0.0;
530 v_re = 0.0;
531 double B[2], s1, s2, s3, sR, sD;
532 for(int i=0; i<4; i++){
533 pR[i] = P1[i] + P2[i];
534 pD[i] = pR[i] + P3[i];
535 }
536 s1 = SCADot(P1,P1);
537 s2 = SCADot(P2,P2);
538 s3 = SCADot(P3,P3);
539 sR = SCADot(pR,pR);
540 sD = SCADot(pD,pD);
541 int G[4][4];
542 for(int i=0; i!=4; i++){
543 for(int j=0; j!=4; j++){
544 if(i==j){
545 if(i==0) G[i][j] = 1;
546 else G[i][j] = -1;
547 }
548 else G[i][j] = 0;
549 }
550 }
551 if(Ang == 0){
552 B[0] = 1;
553 B[1] = 1;
554 temp_PDF = 1;
555 }
556 if(Ang == 1){
557 B[0] = barrier(1,sR,s1,s2,3.0,mass);
558 B[1] = barrier(1,sD,sR,s3,5.0,1.864);
559 double t1[4], T1[4];
560 calt1(P1,P2,t1);
561 calt1(pR,P3,T1);
562 temp_PDF = 0;
563 for(int i=0; i<4; i++){
564 temp_PDF += t1[i]*T1[i]*G[i][i];
565 }
566 }
567 if(Ang == 2){
568 B[0] = barrier(2,sR,s1,s2,3.0,mass);
569 B[1] = barrier(2,sD,sR,s3,5.0,1.864);
570 double t2[4][4], T2[4][4];
571 calt2(P1,P2,t2);
572 calt2(pR,P3,T2);
573 temp_PDF = 0;
574 for(int i=0; i<4; i++){
575 for(int j=0; j<4; j++){
576 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
577 }
578 }
579 }
580 v_re = temp_PDF*B[0]*B[1];
581 return v_re;
582}
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:129
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