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