BOSS 7.1.0
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
38EvtDecayBase* EvtDToKSpieta::clone(){
39 return new EvtDToKSpieta;
40}
41
42void EvtDToKSpieta::init(){
43 checkNArg(0);
44 checkNDaug(3);
45 checkSpinParent(EvtSpinType::SCALAR);
46 ma0 = 0.987; 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.26085;//K1430-kseta
58 phi[1] = 2.6200;
59
60 spin[0]=0;
61 spin[1]=0;
62
63 modetype[0]=23;
64 modetype[1]=13;
65
66 //cout << "Initializing EvtDToKSpieta" << endl;
67 //for (int i=0; i<2; i++) {
68 // cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
69 //}
70 mass_Ks=0.4977;
71 mass_Eta=0.547862;
72 mD = 1.86965;
73 rD = 5;
74 metap = 0.95778;
75 mkstr = 0.89594;
76 mk0 = 0.497611;
77 mass_Kaon = 0.493677;
78 mass_Pion = 0.13957;
79 mass_Pi0 = 0.1349768;
80 math_pi = 3.1415926;
81
82 pi = 3.1415926;
83 mpi = 0.13957;
84 g1 = 0.5468;
85 g2 = 0.23;
86
87 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
88 int EE[4][4][4][4] =
89 { { {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
90 {{0,0,0,0}, {0,0,0,0}, {0,0,0,1}, {0,0,-1,0}},
91 {{0,0,0,0}, {0,0,0,-1}, {0,0,0,0}, {0,1,0,0} },
92 {{0,0,0,0}, {0,0,1,0}, {0,-1,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,0}, {0,0,0,0}, {0,0,0,0} },
95 {{0,0,0,1}, {0,0,0,0}, {0,0,0,0}, {-1,0,0,0}},
96 {{0,0,-1,0}, {0,0,0,0}, {1,0,0,0}, {0,0,0,0} } },
97 { {{0,0,0,0}, {0,0,0,1}, {0,0,0,0}, {0,-1,0,0}},
98 {{0,0,0,-1}, {0,0,0,0}, {0,0,0,0}, {1,0,0,0} },
99 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
100 {{0,1,0,0}, {-1,0,0,0}, {0,0,0,0}, {0,0,0,0} } },
101 { {{0,0,0,0}, {0,0,-1,0}, {0,1,0,0}, {0,0,0,0} },
102 {{0,0,1,0}, {0,0,0,0}, {-1,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,0,0}, {0,0,0,0}, {0,0,0,0} } } };
105 for (int i=0; i<4; i++) {
106 for (int j=0; j<4; j++) {
107 G[i][j] = GG[i][j];
108 for (int k=0; k<4; k++) {
109 for (int l=0; l<4; l++) {
110 E[i][j][k][l] = EE[i][j][k][l];
111 }
112 }
113 }
114 }
115}
116
117void EvtDToKSpieta::initProbMax() {
118 setProbMax(13.00);//setProbMax(1680.0);
119}
120
121void EvtDToKSpieta::decay( EvtParticle *p){
122 //-----------for max value------------------
123/* double maxprob = 0.0;
124 for(int ir=0;ir<=60000000;ir++){
125 p->initializePhaseSpace(getNDaug(),getDaugs());
126 EvtVector4R ks0 = p->getDaug(0)->getP4();
127 EvtVector4R pic = p->getDaug(1)->getP4();
128 EvtVector4R pi0 = p->getDaug(2)->getP4();
129
130 double Ks0[4],Pic[4],Pi0[4];
131 Ks0[0] = ks0.get(0); Pic[0] = pic.get(0); Pi0[0] = pi0.get(0);
132 Ks0[1] = ks0.get(1); Pic[1] = pic.get(1); Pi0[1] = pi0.get(1);
133 Ks0[2] = ks0.get(2); Pic[2] = pic.get(2); Pi0[2] = pi0.get(2);
134 Ks0[3] = ks0.get(3); Pic[3] = pic.get(3); Pi0[3] = pi0.get(3);
135 double value;
136 calPDF(Ks0, Pic, Pi0, value);
137 if(value>maxprob) {
138 maxprob=value;
139 std::cout << "Max PDF = " << ir << " prob= " << value << std::endl;
140 }
141 }
142 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
143 return;*/
144 //-----------------------------------------------
145 p->initializePhaseSpace(getNDaug(),getDaugs());
146 EvtVector4R ks0 = p->getDaug(0)->getP4();
147 EvtVector4R pic = p->getDaug(1)->getP4();
148 EvtVector4R pi0 = p->getDaug(2)->getP4();
149
150 double Ks0[4],Pic[4],Pi0[4];
151 Ks0[0] = ks0.get(0); Pic[0] = pic.get(0); Pi0[0] = pi0.get(0);
152 Ks0[1] = ks0.get(1); Pic[1] = pic.get(1); Pi0[1] = pi0.get(1);
153 Ks0[2] = ks0.get(2); Pic[2] = pic.get(2); Pi0[2] = pi0.get(2);
154 Ks0[3] = ks0.get(3); Pic[3] = pic.get(3); Pi0[3] = pi0.get(3);
155
156 double value;
157 calPDF(Ks0, Pic, Pi0, value);
158 setProb(value);
159 return;
160}
161
162double EvtDToKSpieta::calPDF(double Ks0[], double Pic[], double Pi0[], double & Result) {
163 double cof[2], amp_tmp1[2], amp_tmp2[2], amp_tmp[2], amp_PDF[2], PDF[2];
164 double flag[3], mass_R[2], width_R[2];
165 double pro[2];
166 double sa[3], sb[3], sc[3], B[3];
167 double t1D[4], t1V[4], t1A[4];
168 double pS[4], pV[4], pD[4], pA[4];
169 double pro1[2], pro2[2],proKPi_S[2];
170
171 double rD2 = 25.0;
172 double rRes2 = 9.0;
173 double rRes = 9.0;
174 double mass1[2] ={ ma0, mKn_1430 };
175 double width1[2]={ Ga0, GKn_1430 };
176 double g0[2] ={ 3, 1 };
177 double temp_PDF = 0;
178// PDF[0]=0;
179// PDF[1]=0;
180 double P12[4], P23[4], P13[4];
181 double scpi, snpi, sks0;
182 double s12, s13, s23;
183 for(int i=0; i<4; i++){
184 P12[i] = Pic[i] + Ks0[i];
185 P13[i] = Pi0[i] + Ks0[i];
186 P23[i] = Pic[i] + Pi0[i];
187 }
188 scpi = SCADot(Pic,Pic);
189 snpi = SCADot(Pi0,Pi0);
190 sks0 = SCADot(Ks0,Ks0);
191 s12 = SCADot(P12,P12);
192 s13 = SCADot(P13,P13);
193 s23 = SCADot(P23,P23);
194 double mass1sq;
195 double Amp_KPiS[2];
196 amp_PDF[0] = 0;
197 amp_PDF[1] = 0;
198 PDF[0] = 0;
199 PDF[1] = 0;
200// amp_tmp[0] = 0;
201// amp_tmp[1] = 0;
202
203 for(int i=0; i<2; i++) {
204 amp_tmp[0] = 0;
205 amp_tmp[1] = 0;
206 mass1sq = mass1[i]*mass1[i];
207 cof[0] = rho[i]*cos(phi[i]);
208 cof[1] = rho[i]*sin(phi[i]);
209 temp_PDF = 0;
210 if(modetype[i] == 23){
211 temp_PDF = DDalitz( Pic, Pi0, Ks0, spin[i], mass1[i]);
212 if(g0[i]==1) propagatorRBW(mass1sq,mass1[i],width1[i],s23,scpi,snpi,rRes,spin[i],pro);
213 if(g0[i]==2) KPiSLASS(s23,scpi,snpi,pro);
214 if(g0[i]==3) propagatorFlatte(mass1[i],width1[i],s23,scpi,snpi,0,pro);
215 if(g0[i]==0){
216 pro[0] = 1;
217 pro[1] = 0;
218 }
219 amp_tmp[0] = temp_PDF*pro[0];
220 amp_tmp[1] = temp_PDF*pro[1];
221 }
222 if(modetype[i] == 12){
223 temp_PDF = DDalitz(Ks0, Pic, Pi0, spin[i], mass1[i]);
224 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s12,sks0,scpi,rRes,spin[i],pro);
225 if(g0[i]==2) {propagatorFlatte(mass1[i],width1[i],s12,sks0,scpi,1,pro);
226 pro[0]=pro[0]*(0.01+0.990*0.990)/(0.01+s12);
227 pro[1]=pro[1]*(0.01+0.990*0.990)/(0.01+s12);
228 }
229 if(g0[i]==3) propagatorFlatte(mass1[i],width1[i],s12,sks0,scpi,1,pro);//Only for a0(980)
230
231 if(g0[i]==0){
232 pro[0] = 1;
233 pro[1] = 0;
234 }
235 amp_tmp[0] = temp_PDF*pro[0];
236 amp_tmp[1] = temp_PDF*pro[1];
237 }
238 if(modetype[i] == 13){
239 temp_PDF = DDalitz(Ks0, Pi0, Pic, spin[i], mass1[i]);
240 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s13,sks0,snpi,rRes,spin[i],pro);
241 if(g0[i]==2) { // K*1430 Flatte
242 double skm2[2]={sks0, mass_Ks *mass_Ks};
243 double spi2[2]={snpi, mass_Eta *mass_Eta};
244 propagatorKstr1430(mass1[i],s13,skm2,spi2,pro);
245 };
246 if(g0[i]==3) propagatorFlatte(mass1[i],width1[i],s13,sks0,snpi,1,pro);//Only for a0(980)
247 if(g0[i]==0){
248 pro[0] = 1;
249 pro[1] = 0;
250 }
251 amp_tmp[0] = temp_PDF*pro[0];
252 amp_tmp[1] = temp_PDF*pro[1];
253 }
254 if(modetype[i] == 132){
255 KPiSLASS(s12,sks0,scpi,Amp_KPiS);
256 amp_tmp[0] = Amp_KPiS[0];
257 amp_tmp[1] = Amp_KPiS[1];
258 }
259//cout<<amp_tmp[0]<<" "<<amp_tmp[1]<<endl;
260 Com_Multi(amp_tmp,cof,amp_PDF);
261 PDF[0] += amp_PDF[0];
262 PDF[1] += amp_PDF[1];
263}
264
265
266// printf("PDF[0] : %.10f \n",PDF[0]);
267// printf("PDF[1] : %.10f \n",PDF[1]);
268 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
269 //printf("Value : %.10f \n",value);
270
271 Result = value;
272
273}
274
275void EvtDToKSpieta:: 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 EvtDToKSpieta:: 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 EvtDToKSpieta:: 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 EvtDToKSpieta:: 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 EvtDToKSpieta:: 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 EvtDToKSpieta:: 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 EvtDToKSpieta:: 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 EvtDToKSpieta:: 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 EvtDToKSpieta:: 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 EvtDToKSpieta:: 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 EvtDToKSpieta:: 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 EvtDToKSpieta:: 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}
453
454void EvtDToKSpieta:: PiPiSWASS(double sa, double sb, double sc, double prop[2]) {
455 double tmp = sb-sc;
456 double tmp2 = sa+tmp;
457 double qs = 0.25*tmp2*tmp2/sa-sb;
458 double q = sqrt(qs);
459 double a0 = -0.11/mass_Pion;
460 prop[0] = 1/(1+a0*a0*q*q);
461 prop[1] = a0*q/(1+a0*a0*q*q);
462}
463void EvtDToKSpieta:: KPiSLASS(double sa, double sb, double sc, double prop[2]) {
464 const double m1430 = 1.441;
465 const double sa0 = 1.441*1.441; // m1430*m1430;
466 const double w1430 = 0.193;
467 const double Lass1 = 0.25/sa0;
468 double tmp = sb-sc;
469 double tmp1 = sa0+tmp;
470 double q0 = Lass1*tmp1*tmp1-sb;
471 if(q0<0) q0 = 1e-16;
472 double tmp2 = sa+tmp;
473 double qs = 0.25*tmp2*tmp2/sa-sb;
474 double q = sqrt(qs);
475 double width = w1430*q*m1430/sqrt(sa*q0);
476 double temp_R = atan(m1430*width/(sa0-sa));
477 if(temp_R<0) temp_R += math_pi;
478 double deltaR = -1.915 + temp_R; //fiR=-109.7
479 double temp_F = atan(0.226*q/(2.0-3.819*qs)); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
480 if(temp_F<0) temp_F += math_pi;
481 double deltaF = 0.002 + temp_F; //fiF=0.1
482 double deltaS = deltaR + 2.0*deltaF;
483 double t1 = 0.96*sin(deltaF);
484 double t2 = sin(deltaR);
485 double CF[2], CS[2];
486 CF[0] = cos(deltaF);
487 CF[1] = sin(deltaF);
488 CS[0] = cos(deltaS);
489 CS[1] = sin(deltaS);
490 prop[0] = t1*CF[0] + t2*CS[0];
491 prop[1] = t1*CF[1] + t2*CS[1];
492}
493
494void EvtDToKSpieta:: Flatte_rhoab(double sa, double sb, double sc, double rho[2]){
495 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
496 if(q>0) {
497 rho[0]=2* sqrt(q/sa);
498 rho[1]=0;
499 }
500 else if(q<0){
501 rho[0]=0;
502 rho[1]=2*sqrt(-q/sa);
503 }
504}
505
506void EvtDToKSpieta:: propagatorKstr1430(double mass, double sx, double *sb, double *sc, double prop[2]) //K*1430 Flatte
507{
508 double unit[2]={1.0};
509 double ci[2]={0,1};
510 double rho1[2];
511 Flatte_rhoab(sx,sb[0],sc[0],rho1);
512 double rho2[2];
513 Flatte_rhoab(sx,sb[1],sc[1],rho2);
514 double gKPi_Kstr1430=0.2990, gEtaPK_Kstr1430=0.0529;
515 double tmp1[2]={gKPi_Kstr1430,0};
516 double tmp11[2];
517 double tmp2[2]={gEtaPK_Kstr1430,0};
518 double tmp22[2];
519 Com_Multi(tmp1,rho1,tmp11);
520 Com_Multi(tmp2,rho2,tmp22);
521 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
522 double tmp31[2];
523 Com_Multi(tmp3, ci,tmp31);
524 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
525 Com_Divide( unit,tmp4, prop);
526}
527double EvtDToKSpieta:: DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
528 double pR[4], pD[4];
529 double temp_PDF, v_re;
530 temp_PDF = 0.0;
531 v_re = 0.0;
532 double B[2], s1, s2, s3, sR, sD;
533 for(int i=0; i<4; i++){
534 pR[i] = P1[i] + P2[i];
535 pD[i] = pR[i] + P3[i];
536 }
537 s1 = SCADot(P1,P1);
538 s2 = SCADot(P2,P2);
539 s3 = SCADot(P3,P3);
540 sR = SCADot(pR,pR);
541 sD = SCADot(pD,pD);
542 int G[4][4];
543 for(int i=0; i!=4; i++){
544 for(int j=0; j!=4; j++){
545 if(i==j){
546 if(i==0) G[i][j] = 1;
547 else G[i][j] = -1;
548 }
549 else G[i][j] = 0;
550 }
551 }
552 if(Ang == 0){
553 B[0] = 1;
554 B[1] = 1;
555 temp_PDF = 1;
556 }
557 if(Ang == 1){
558 B[0] = barrier(1,sR,s1,s2,3.0,mass);
559 B[1] = barrier(1,sD,sR,s3,5.0,1.869);
560 // B[0] = Barrier(1,sR,s1,s2,9.0);
561 // B[1] = Barrier(1,sD,sR,s3,25.0);
562 double t1[4], T1[4];
563 calt1(P1,P2,t1);
564 calt1(pR,P3,T1);
565 temp_PDF = 0;
566 for(int i=0; i<4; i++){
567 temp_PDF += t1[i]*T1[i]*G[i][i];
568 }
569 }
570 if(Ang == 2){
571 B[0] = barrier(2,sR,s1,s2,3.0,mass);
572 B[1] = barrier(2,sD,sR,s3,5.0,1.869);
573 // B[0] = Barrier(2,sR,s1,s2,9.0);
574 // B[1] = Barrier(2,sD,sR,s3,25.0);
575 double t2[4][4], T2[4][4];
576 calt2(P1,P2,t2);
577 calt2(pR,P3,T2);
578 temp_PDF = 0;
579 for(int i=0; i<4; i++){
580 for(int j=0; j<4; j++){
581 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
582 }
583 }
584 }
585 v_re = temp_PDF*B[0]*B[1];
586 return v_re;
587}
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:121
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