BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToEta3pi.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: EvtDsToEta3pi.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/EvtDsToEta3pi.hh"
27#include <string>
28#include "EvtGenBase/EvtConst.hh"
29#include "EvtGenBase/EvtFlatte.hh"
30#include "EvtGenBase/EvtDecayTable.hh"
31using std::endl;
32
33EvtDsToEta3pi::~EvtDsToEta3pi() {}
34
35void EvtDsToEta3pi::getName(std::string& model_name){
36 model_name="DsToEta3pi";
37}
38
42
43void EvtDsToEta3pi::init(){
44
45 // check that there are 0 arguments
46 checkNArg(0);
47 checkNDaug(4);
48
49 checkSpinParent(EvtSpinType::SCALAR);
50 checkSpinDaughter(0,EvtSpinType::SCALAR);
51 checkSpinDaughter(1,EvtSpinType::SCALAR);
52 checkSpinDaughter(2,EvtSpinType::SCALAR);
53 checkSpinDaughter(3,EvtSpinType::SCALAR);
54
55 int mode = 0;
56 phi[0] = 2.4600e+00; modetype[0] = 1;
57 phi[1] = 0.0000e+00; modetype[1] = 2;
58 phi[2] = -1.5394e+00; modetype[2] = 3;
59 phi[3] = -6.0185e+00; modetype[3] = 4;
60 phi[4] = -6.0185e+00; modetype[4] = 5;
61 phi[5] = -6.1365e+00; modetype[5] = 6;
62 phi[6] = -6.1365e+00; modetype[6] = 7;
63 phi[7] = 1.3514e+00; modetype[7] = 8;
64 phi[8] = 2.4527e+00; modetype[8] = 9;
65 phi[9] = -2.0956e+00; modetype[9] = 10;
66 phi[10]= -2.0956e+00; modetype[10] = 11;
67
68 rho[0] = -1.5823e+00;
69 rho[1] = 10.000e+00;
70 rho[2] = 1.3750e+01;
71 rho[3] = -3.7257e-01;
72 rho[4] = -3.7257e-01;
73 rho[5] = -3.6293e+00;
74 rho[6] = -3.6293e+00;
75 rho[7] = -8.9184e+00;
76 rho[8] = -1.4292e+01;
77 rho[9] = 1.2518e+00;
78 rho[10]= 1.2518e+00;
79
80 mass1[0] = 9.9000e-01;
81 mass1[1] = 7.7526e-01;
82 mass1[2] = 5.2600e-01;
83 mass1[3] = 9.9000e-01;
84 mass1[4] = 9.9000e-01;
85 mass1[5] = 9.9000e-01;
86 mass1[6] = 9.9000e-01;
87 mass1[7] = 9.6500e-01;
88 mass1[8] = 5.2600e-01;
89 mass1[9] = 9.9000e-01;
90 mass1[10]= 9.9000e-01;
91
92 mass2[0] = 7.7526e-01;
93 mass2[1] = 1.3462e+00;
94 mass2[2] = 1.3462e+00;
95 mass2[3] = 1.4050e+00;
96 mass2[4] = 1.4050e+00;
97 mass2[5] = 1.8000e+00;
98 mass2[6] = 1.8000e+00;
99 mass2[7] = 1.8000e+00;
100 mass2[8] = 1.8000e+00;
101 mass2[9] = 1.4502e+00;
102 mass2[10]= 1.4502e+00;
103
104 width1[0] = 7.5000e-02;
105 width1[1] = 1.4910e-01;
106 width1[2] = 5.0000e-01;
107 width1[3] = 7.5000e-02;
108 width1[4] = 7.5000e-02;
109 width1[5] = 7.5000e-02;
110 width1[6] = 7.5000e-02;
111 width1[7] = 7.5000e-02;
112 width1[8] = 5.0000e-01;
113 width1[9] = 7.5000e-02;
114 width1[10]= 7.5000e-02;
115
116 width2[0] = 1.4910e-01;
117 width2[1] = 3.6900e-01;
118 width2[2] = 3.6900e-01;
119 width2[3] = 5.1000e-02;
120 width2[4] = 5.1000e-02;
121 width2[5] = 3.8600e-01;
122 width2[6] = 3.8600e-01;
123 width2[7] = 3.8600e-01;
124 width2[8] = 3.8600e-01;
125 width2[9] = 5.4000e-02;
126 width2[10]= 5.4000e-02;
127
128 //cout << "DsToEta3pi :"<< endl;
129 //for (int i=0; i<11; i++) {
130 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << " m1= " << mass1[i] << " m2= " << mass2[i] << " w1= " << width1[i] << " w2= " << width2[i] << endl;
131 //}
132
133 rRes = 3.0;
134 rRes2 = 9.0;
135 rD = 5.0;
136 rD2 = 25.0;
137 mass_Ks = 0.497611;
138 mass_Kaon = 0.49368;
139 mass_Pion = 0.13957;
140 mass_Pi0 = 0.1349766;
141 meta = 0.547862;
142 GS1 = 0.636619783;
143 GS2 = 0.01860182466;
144 GS3 = 0.1591549458; // 1/(2*math_2pi)
145 GS4 = 0.00620060822;
146 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
147 int EE[4][4][4][4] =
148 { { {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
149 {{0,0,0,0}, {0,0,0,0}, {0,0,0,1}, {0,0,-1,0}},
150 {{0,0,0,0}, {0,0,0,-1}, {0,0,0,0}, {0,1,0,0} },
151 {{0,0,0,0}, {0,0,1,0}, {0,-1,0,0}, {0,0,0,0} } },
152 { {{0,0,0,0}, {0,0,0,0}, {0,0,0,-1}, {0,0,1,0} },
153 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
154 {{0,0,0,1}, {0,0,0,0}, {0,0,0,0}, {-1,0,0,0}},
155 {{0,0,-1,0}, {0,0,0,0}, {1,0,0,0}, {0,0,0,0} } },
156 { {{0,0,0,0}, {0,0,0,1}, {0,0,0,0}, {0,-1,0,0}},
157 {{0,0,0,-1}, {0,0,0,0}, {0,0,0,0}, {1,0,0,0} },
158 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
159 {{0,1,0,0}, {-1,0,0,0}, {0,0,0,0}, {0,0,0,0} } },
160 { {{0,0,0,0}, {0,0,-1,0}, {0,1,0,0}, {0,0,0,0} },
161 {{0,0,1,0}, {0,0,0,0}, {-1,0,0,0}, {0,0,0,0} },
162 {{0,-1,0,0}, {1,0,0,0}, {0,0,0,0}, {0,0,0,0} },
163 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} } } };
164 for (int i=0; i<4; i++) {
165 for (int j=0; j<4; j++) {
166 G[i][j] = GG[i][j];
167 for (int k=0; k<4; k++) {
168 for (int l=0; l<4; l++) {
169 E[i][j][k][l] = EE[i][j][k][l];
170 }
171 }
172 }
173 }
174}
175
176void EvtDsToEta3pi::initProbMax() {
177 setProbMax(108000.0);
178}
179
180void EvtDsToEta3pi::decay( EvtParticle *p){
181
182// double maxprob = 0.0;
183// for(int ir=0;ir<=60000000;ir++){
184// p->initializePhaseSpace(getNDaug(),getDaugs());
185// EvtVector4R _pip1 = p->getDaug(0)->getP4();
186// EvtVector4R _pip2 = p->getDaug(1)->getP4();
187// EvtVector4R _pim = p->getDaug(2)->getP4();
188// EvtVector4R _eta = p->getDaug(3)->getP4();
189//
190// double _Pip1[4],_Pip2[4],_Pim[4],_Eta[4];
191// _Pip1[1]=_pip1.get(1);_Pip1[2]=_pip1.get(2);_Pip1[3]=_pip1.get(3);
192// _Pip2[1]=_pip2.get(1);_Pip2[2]=_pip2.get(2);_Pip2[3]=_pip2.get(3);
193// _Pim[1] = _pim.get(1);_Pim[2] = _pim.get(2);_Pim[3] = _pim.get(3);
194// _Eta[1] = _eta.get(1);_Eta[2] = _eta.get(2);_Eta[3] = _eta.get(3);
195//
196// double _prob;
197// int g0[11]={1,1,1,1,1,1,1,1,1,1,1};
198// int g1[11]={1,1,1,1,1,1,1,1,1,1,1};
199// int g2[11]={0,0,0,0,0,0,0,0,0,0,0};
200// int nstates=11;
201// calEvaMy(_Pip1, _Pip2, _Pim, _Eta, mass1, mass2, width1, width2, rho, phi, g0,g1,g2, modetype, nstates, _prob);
202// if(_prob>maxprob) {
203// maxprob=_prob;
204// printf("nrun = %d,maxprob = %.10f,prob = %.10f\n\n",ir,maxprob,_prob);
205// }
206// }
207
208
209 p->initializePhaseSpace(getNDaug(),getDaugs());
210 EvtVector4R pip1 = p->getDaug(0)->getP4();
211 EvtVector4R pip2 = p->getDaug(1)->getP4();
212 EvtVector4R pim = p->getDaug(2)->getP4();
213 EvtVector4R eta = p->getDaug(3)->getP4();
214
215 double Pip1[4],Pip2[4],Pim[4],Eta[4];
216 Pip1[1]=pip1.get(1);Pip1[2]=pip1.get(2);Pip1[3]=pip1.get(3);
217 Pip2[1]=pip2.get(1);Pip2[2]=pip2.get(2);Pip2[3]=pip2.get(3);
218 Pim[1] = pim.get(1);Pim[2] = pim.get(2);Pim[3] = pim.get(3);
219 Eta[1] = eta.get(1);Eta[2] = eta.get(2);Eta[3] = eta.get(3);
220
221 double prob;
222 int g0[11]={1,1,1,1,1,1,1,1,1,1,1};
223 int g1[11]={1,1,1,1,1,1,1,1,1,1,1};
224 int g2[11]={0,0,0,0,0,0,0,0,0,0,0};
225 int nstates=11;
226
227 calEvaMy(Pip1, Pip2, Pim, Eta, mass1, mass2, width1, width2, rho, phi, g0, g1, g2, modetype, nstates, prob);
228 setProb(prob);
229
230 return ;
231}
232
233void EvtDsToEta3pi::calEvaMy(double* Pip1, double* Pip2, double* Pim, double* Eta, double *mass1, double *mass2, double* width1, double* width2, double *amp, double *phase,int* g0, int* g1, int* g2, int* modetype, int nstates, double & Result){
234 double prho1[4],prho2[4],pa11[4],pa12[4],pa01[4],pa02[4],pD[4],pf11[4],pf12[4],pa0m[4],peta12951[4],peta12952[4],psigma1[4],psigma2[4];
235 Pip1[0] = sqrt(mass_Pion*mass_Pion + Pip1[1]*Pip1[1] + Pip1[2]*Pip1[2] + Pip1[3]*Pip1[3]);
236 Pip2[0] = sqrt(mass_Pion*mass_Pion + Pip2[1]*Pip2[1] + Pip2[2]*Pip2[2] + Pip2[3]*Pip2[3]);
237 Pim[0] = sqrt(mass_Pion*mass_Pion + Pim[1]*Pim[1] + Pim[2]*Pim[2] + Pim[3]*Pim[3]);
238 Eta[0] = sqrt(meta*meta + Eta[1]*Eta[1] + Eta[2]*Eta[2] + Eta[3]*Eta[3]);
239
240
241 for(int i=0;i!=4;i++) {
242 prho1[i]=Pim[i] +Pip1[i];
243 prho2[i]=Pim[i] +Pip2[i];
244 pa11[i] =prho1[i]+Pip2[i];
245 pa12[i] =prho2[i]+Pip1[i];
246 psigma1[i] =Pim[i]+Pip1[i];
247 psigma2[i] =Pim[i]+Pip2[i];
248 pa01[i] =Eta[i] +Pip1[i];
249 pa02[i] =Eta[i] +Pip2[i];
250 pa0m[i] =Eta[i] +Pim[i];
251 pf11[i] =Pim[i] +Pip1[i] +Eta[i];
252 pf12[i] =Pim[i] +Pip2[i] +Eta[i];
253 pD[i] =pa11[i] +Eta[i];
254 peta12951[i] = Pim[i] +Pip1[i] +Eta[i];
255 peta12952[i] = Pim[i] +Pip2[i] +Eta[i];
256 }
257 double cof[2], amp_tmp1[2], amp_tmp2[2], amp_tmp[2], amp_PDF[2], PDF[2];
258 double amp_Omg, phs_Omg;
259 amp_PDF[0] = 0;
260 amp_PDF[1] = 0;
261 PDF[0] = 0;
262 PDF[1] = 0;
263 amp_PDF[0] = 0;
264 amp_PDF[1] = 0;
265 PDF[0] = 0;
266 PDF[1] = 0;
267 double spion1,spion2,spim,seta,sa11,sa12,srho1,srho2,sa01,sa02,sD,sf11,sf12,sa0m,seta12951,seta12952,ssigma1,ssigma2;
268 double t1rho1[4],t2rho1[4][4],t1rho2[4],t2rho2[4][4],t1a01[4],t1a02[4],t2a11[4][4],t2a12[4][4],t1d01[4],t1d02[4],t1f11[4],t2f11[4][4],t1f12[4],t2f12[4][4],t1a11[4],t1a12[4],t1f11_2[4],t1f12_2[4];
269 double temp_PDF, tmp1, tmp2, temp[2],tmp3,tmp4,temp_PDF1;
270 double pro[2], pro0[2], pro1[2],pro2[2],pro3[2],pro4[2];
271 double t1D[4],t2D[4][4],B[3], t1Tm[4],A[2], Bc[3];
272
273 double mass1sq, mass2sq;
274
275 seta = SCADot(Eta,Eta);
276 spion1 = SCADot(Pip1,Pip1);
277 spion2 = SCADot(Pip2,Pip2);
278 spim = SCADot(Pim,Pim);
279
280 srho1 = SCADot(prho1,prho1);
281 srho2 = SCADot(prho2,prho2);
282 sa11 = SCADot(pa11,pa11);
283 sa12 = SCADot(pa12,pa12);
284
285 sa01 = SCADot(pa01,pa01);
286 sa02 = SCADot(pa02,pa02);
287 sa0m = SCADot(pa0m,pa0m);
288 sf11 = SCADot(pf11,pf11);
289 sf12 = SCADot(pf12,pf12);
290 sD = SCADot(pD,pD);
291 seta12951 = SCADot(peta12951,peta12951);
292 seta12952 = SCADot(peta12952,peta12952);
293 ssigma1 = SCADot(psigma1,psigma1);
294 ssigma2 = SCADot(psigma2,psigma2);
295 double seta2[2]={mass_Ks*mass_Ks,seta};
296 double spion12[2]={mass_Kaon*mass_Kaon,spion1};
297 double spion22[2]={mass_Kaon*mass_Kaon,spion2};
298 double spim2[2]={mass_Kaon*mass_Kaon,spim};
299
300 calt1(Pip1,Pim,t1rho1);
301 calt2(Pip1,Pim,t2rho1);
302 calt1(Pip2,Pim,t1rho2);
303 calt2(Pip2,Pim,t2rho2);
304
305 calt1(Pip1,Eta,t1a01);
306 calt1(Pip2,Eta,t1a02);
307
308 calt1(pa01,prho2,t1d01);
309 calt1(pa02,prho1,t1d02);
310
311 calt1(pa0m,Pip1,t1f11);
312 calt1(pa0m,Pip2,t1f12);
313
314 calt1(psigma1,Pip2,t1a11);
315 calt1(psigma2,Pip1,t1a12);
316
317 calt1(pa01,Pim,t1f11_2);
318 calt1(pa02,Pim,t1f12_2);
319
320 calt2(prho1,Pip2,t2a11);
321 calt2(prho2,Pip1,t2a12);
322
323 calt2(pa01,Pip1,t2f11);
324 calt2(pa02,Pip2,t2f12);
325
326 for(int i=0; i<nstates; i++){
327 int d=0;
328 cof[0] = amp[i]*cos(phase[i]);
329 cof[1] = amp[i]*sin(phase[i]);
330 mass1sq = mass1[i]*mass1[i];
331 mass2sq = mass2[i]*mass2[i];
332 temp_PDF = 0;
333 temp_PDF1 = 0;
334 amp_tmp[0] = 0;
335 amp_tmp[1] = 0;
336 amp_tmp1[0] = 0;
337 amp_tmp1[1] = 0;
338 amp_tmp2[0] = 0;
339 amp_tmp2[1] = 0;
340 tmp1 = 0;
341 tmp2 = 0;
342 tmp3 = 0;
343 tmp4 = 0;
344
345 if(modetype[i]==1)
346 {
347 if(g0[i]==1) propagatora0980(mass1[i],sa01,seta2,spion12,pro0);
348 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
349 if(g1[i]==1) propagatorGS(mass2sq,mass2[i],width2[i],srho2,spion2,spim,rRes2,pro1);
350 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; }
351 Com_Multi(pro0,pro1,pro);
352
353 for(int a=0; a<4; a++)
354 {
355 temp_PDF += G[a][a]*t1d01[a]*t1rho2[a];
356 }
357
358 B[1] = Barrier(1,srho2,spion2,spim,rRes2);
359 B[2] = Barrier(1,sD,sa01,srho2,rD2);
360 tmp1 = B[1]*B[2]*temp_PDF;
361 amp_tmp1[0] = tmp1*pro[0];
362 amp_tmp1[1] = tmp1*pro[1];
363
364 temp_PDF=0;
365 if(g0[i]==1) propagatora0980(mass1[i],sa02,seta2,spion22,pro0);
366 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
367 if(g1[i]==1) propagatorGS(mass2sq,mass2[i],width2[i],srho1,spion1,spim,rRes2,pro1);
368 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; }
369 Com_Multi(pro0,pro1,pro);
370
371 for(int a=0; a<4; a++)
372 {
373 temp_PDF += G[a][a]*t1d02[a]*t1rho1[a];
374 }
375
376 B[1] = Barrier(1,srho1,spion1,spim,rRes2);
377 B[2] = Barrier(1,sD,sa02,srho1,rD2);
378 tmp2 = B[1]*B[2]*temp_PDF;
379 amp_tmp2[0] = tmp2*pro[0];
380 amp_tmp2[1] = tmp2*pro[1];
381 }
382 if(modetype[i]==2)
383 {
384 if(g0[i]==1) propagatorGS(mass1sq,mass1[i],width1[i],srho1,spion1,spim,rRes2,pro0);
385 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
386 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],sa11,srho1,spion2,rRes2,g2[i],pro1);
387 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; }
388 Com_Multi(pro0,pro1,pro);
389 calt1(pa11,Eta,t1D);
390 B[0] = Barrier(1,srho1,spion1,spim,rRes2);
391 B[2] = Barrier(1,sD,sa11,seta,rD2);
392 if(g2[i]==0)
393 {
394 for(int a=0; a<4; a++)
395 {
396 for(int j=0; j<4; j++)
397 {
398 temp_PDF += t1D[a]*(G[a][j]-pa11[a]*pa11[j]/sa11)*t1rho1[j]*G[a][a]*G[j][j];
399 }
400 }
401 tmp1 = B[0]*B[2]*temp_PDF;
402 }
403// if(g2[i]==2)
404// {
405// for(int a=0; a<4; a++)
406// {
407// for(int j=0; j<4; j++)
408// {
409// temp_PDF += t1D[a]*t2a11[a][j]*t1rho1[j]*G[a][a]*G[j][j];
410// }
411// }
412// B[1] = Barrier(2,sa11,srho1,spion2,rRes2);
413// tmp1 = B[0]*B[1]*B[2]*temp_PDF;
414// }
415 amp_tmp1[0] = tmp1*pro[0];
416 amp_tmp1[1] = tmp1*pro[1];
417
418 temp_PDF=0;
419 if(g0[i]==1) propagatorGS(mass1sq,mass1[i],width1[i],srho2,spion2,spim,rRes2,pro0);
420 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
421 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],sa12,srho2,spion1,rRes2,g2[i],pro1);
422 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; }
423 Com_Multi(pro0,pro1,pro);
424 calt1(pa12,Eta,t1D);
425 B[0] = Barrier(1,srho2,spion2,spim,rRes2);
426 B[2] = Barrier(1,sD,sa12,seta,rD2);
427 if(g2[i]==0)
428 {
429 for(int a=0; a<4; a++)
430 {
431 for(int j=0; j<4; j++)
432 {
433 temp_PDF += t1D[a]*(G[a][j]-pa12[a]*pa12[j]/sa12)*t1rho2[j]*G[a][a]*G[j][j];
434 }
435 }
436 tmp2 = B[0]*B[2]*temp_PDF;
437 }
438// if(g2[i]==2)
439// {
440// for(int a=0; a<4; a++)
441// {
442// for(int j=0; j<4; j++)
443// {
444// temp_PDF += t1D[a]*t2a12[a][j]*t1rho2[j]*G[a][a]*G[j][j];
445// }
446// }
447// B[1] = Barrier(2,sa12,srho2,spion1,rRes2);
448// tmp2 = B[0]*B[1]*B[2]*temp_PDF;
449// }
450 amp_tmp2[0] = tmp2*pro[0];
451 amp_tmp2[1] = tmp2*pro[1];
452 }
453 if(modetype[i]==3)
454 {
455 if(g0[i]==1) propagatorsigma500(ssigma1,spion1,spim,pro0);
456 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
457 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],sa11,ssigma1,spion2,rRes2,1,pro1);
458 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; }
459// printf("pro0 = %.10f, %.10f\n", pro0[0], pro0[1]);
460// printf("pro1 = %.10f, %.10f\n", pro1[0], pro1[1]);
461 Com_Multi(pro0,pro1,pro);
462 calt1(pa11,Eta,t1D);
463 B[2] = Barrier(1,sD,sa11,seta,rD2);
464 B[1] = Barrier(1,sa11,ssigma1,spion2,rRes2);
465 for(int a=0; a<4; a++)
466 {
467 temp_PDF += t1D[a]*t1a11[a]*G[a][a];
468 }
469 tmp1 = B[1]*B[2]*temp_PDF;
470 amp_tmp1[0] = tmp1*pro[0];
471 amp_tmp1[1] = tmp1*pro[1];
472 temp_PDF=0;
473 if(g0[i]==1) propagatorsigma500(ssigma2,spion2,spim,pro0);
474 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
475 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],sa12,ssigma2,spion1,rRes2,1,pro1);
476 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; }
477// printf("pro0 = %f, %f\n", pro0[0], pro0[1]);
478// printf("pro1 = %f, %f\n", pro1[0], pro1[1]);
479 Com_Multi(pro0,pro1,pro);
480 calt1(pa12,Eta,t1D);
481 B[2] = Barrier(1,sD,sa12,seta,rD2);
482 B[1] = Barrier(1,sa12,ssigma2,spion1,rRes2);
483 for(int a=0; a<4; a++)
484 {
485 temp_PDF += t1D[a]*t1a12[a]*G[a][a];
486 }
487 tmp2 = B[1]*B[2]*temp_PDF;
488 amp_tmp2[0] = tmp2*pro[0];
489 amp_tmp2[1] = tmp2*pro[1];
490 }
491 if(modetype[i]==4)
492 {
493 if(g0[i]==1) propagatora0980(mass1[i],sa0m,spim2,seta2,pro0);
494 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; pro2[0]=1;pro2[1]=0;}
495 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],seta12951,sa0m,spion1,rRes2,0,pro1);
496 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; pro3[0]=1; pro3[1]=0;}
497 Com_Multi(pro0,pro1,pro);
498 calt1(pf11,Pip2,t1D);
499 tmp1 = 1;
500 amp_tmp1[0] = tmp1*pro[0];
501 amp_tmp1[1] = tmp1*pro[1];
502
503 temp_PDF=0;
504 if(g0[i]==1) propagatora0980(mass1[i],sa0m,spim2,seta2,pro0);
505 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; pro2[0]=1; pro2[1]=0;}
506 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],seta12952,sa0m,spion2,rRes2,0,pro1);
507 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; pro3[0]=1; pro3[1]=0;}
508 Com_Multi(pro0,pro1,pro);
509 calt1(pf12,Pip1,t1D);
510 tmp2 = 1;
511 amp_tmp2[0] = tmp2*pro[0];
512 amp_tmp2[1] = tmp2*pro[1];
513 }
514 if(modetype[i]==5)
515 {
516 if(g0[i]==1) propagatora0980(mass1[i],sa01,spion12,seta2,pro2);
517 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; pro2[0]=1;pro2[1]=0;}
518 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],seta12951,sa01,spim,rRes2,0,pro3);
519 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; pro3[0]=1; pro3[1]=0;}
520 Com_Multi(pro2,pro3,pro4);
521 calt1(pf11,Pip2,t1D);
522 tmp3 = 1;
523 amp_tmp1[0] = tmp3*pro4[0];
524 amp_tmp1[1] = tmp3*pro4[1];
525
526 temp_PDF1=0;
527 if(g0[i]==1) propagatora0980(mass1[i],sa02,spion22,seta2,pro2);
528 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; pro2[0]=1; pro2[1]=0;}
529 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],seta12952,sa02,spim,rRes2,0,pro3);
530 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; pro3[0]=1; pro3[1]=0;}
531 Com_Multi(pro2,pro3,pro4);
532 calt1(pf12,Pip1,t1D);
533 tmp4 = 1;
534 amp_tmp2[0] = tmp4*pro4[0];
535 amp_tmp2[1] = tmp4*pro4[1];
536 }
537 if(modetype[i]==6)
538 {
539 if(g0[i]==1) propagatora0980(mass1[i],sa0m,spim2,seta2,pro0);
540 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; pro2[0]=1;pro2[1]=0;}
541// if(g1[i]==0) propagatorRBW(mass2sq,mass2[i],width2[i],seta12951,sa0m,spion1,rRes2,0,pro1);
542 if(g1[i]==1) { pro1[0] = 1; pro1[1] = 0; pro3[0]=1; pro3[1]=0;}
543 Com_Multi(pro0,pro1,pro);
544 calt1(pf11,Pip2,t1D);
545 tmp1 = 1;
546 amp_tmp1[0] = tmp1*pro[0];
547 amp_tmp1[1] = tmp1*pro[1];
548 temp_PDF=0;
549 if(g0[i]==1) propagatora0980(mass1[i],sa0m,spim2,seta2,pro0);
550 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; pro2[0]=1; pro2[1]=0;}
551// if(g1[i]==0) propagatorRBW(mass2sq,mass2[i],width2[i],seta12952,sa0m,spion2,rRes2,0,pro1);
552 if(g1[i]==1) { pro1[0] = 1; pro1[1] = 0; pro3[0]=1; pro3[1]=0;}
553 Com_Multi(pro0,pro1,pro);
554 calt1(pf12,Pip1,t1D);
555 tmp2 = 1;
556 amp_tmp2[0] = tmp2*pro[0];
557 amp_tmp2[1] = tmp2*pro[1];
558 }
559 if(modetype[i]==7)
560 {
561
562 if(g0[i]==1) propagatora0980(mass1[i],sa01,spion12,seta2,pro2);
563 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; pro2[0]=1;pro2[1]=0;}
564 //if(g1[i]==0) propagatorRBW(mass2sq,mass2[i],width2[i],seta12951,sa01,spim,rRes2,0,pro3);
565 if(g1[i]==1) { pro1[0] = 1; pro1[1] = 0; pro3[0]=1; pro3[1]=0;}
566 Com_Multi(pro2,pro3,pro4);
567 calt1(pf11,Pip2,t1D);
568 tmp3 = 1;
569 amp_tmp1[0] = tmp3*pro4[0];
570 amp_tmp1[1] = tmp3*pro4[1];
571 temp_PDF1=0;
572 if(g0[i]==1) propagatora0980(mass1[i],sa02,spion22,seta2,pro2);
573 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; pro2[0]=1; pro2[1]=0;}
574 //if(g1[i]==0) propagatorRBW(mass2sq,mass2[i],width2[i],seta12952,sa02,spim,rRes2,0,pro3);
575 if(g1[i]==1) { pro1[0] = 1; pro1[1] = 0; pro3[0]=1; pro3[1]=0;}
576 Com_Multi(pro2,pro3,pro4);
577 calt1(pf12,Pip1,t1D);
578 tmp4 = 1;
579 amp_tmp2[0] = tmp4*pro4[0];
580 amp_tmp2[1] = tmp4*pro4[1];
581 }
582 if(modetype[i]==8)
583 {
584 if(g0[i]==1) propagator980(mass1[i],srho1,spion12,spim2,pro0);
585 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
586 //if(g1[i]==0) propagatorRBW(mass2sq,mass2[i],width2[i],seta12951,srho1,seta,rRes2,0,pro1);
587 if(g1[i]==1) { pro1[0] = 1; pro1[1] = 0; }
588 Com_Multi(pro0,pro1,pro);
589// printf("pro0 = %f, %f\n", pro0[0], pro0[1]);
590// printf("pro1 = %f, %f\n", pro1[0], pro1[1]);
591// printf("pro = %f, %f\n", pro[0], pro[1]);
592 tmp1 = 1;
593 amp_tmp1[0] = tmp1*pro[0];
594 amp_tmp1[1] = tmp1*pro[1];
595 temp_PDF=0;
596 if(g0[i]==1) propagator980(mass1[i],srho2,spion22,spim2,pro0);
597 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
598 //if(g1[i]==0) propagatorRBW(mass2sq,mass2[i],width2[i],seta12952,srho2,seta,rRes2,0,pro1);
599 if(g1[i]==1) { pro1[0] = 1; pro1[1] = 0; }
600 Com_Multi(pro0,pro1,pro);
601// printf("pro0 = %f, %f\n", pro0[0], pro0[1]);
602// printf("pro1 = %f, %f\n", pro1[0], pro1[1]);
603// printf("pro = %f, %f\n", pro[0], pro[1]);
604 calt1(pf12,Pip1,t1D);
605 tmp2 = 1;
606 amp_tmp2[0] = tmp2*pro[0];
607 amp_tmp2[1] = tmp2*pro[1];
608 }
609 if(modetype[i]==9)
610 {
611 if(g0[i]==1) propagatorsigma500(ssigma1,spion1,spim,pro0);
612 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
613 //if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],seta12951,ssigma1,seta,rRes2,0,pro1);
614 if(g1[i]==1) { pro1[0] = 1; pro1[1] = 0; }
615 Com_Multi(pro0,pro1,pro);
616 calt1(pf11,Pip2,t1D);
617 tmp1 = 1;
618 amp_tmp1[0] = tmp1*pro[0];
619 amp_tmp1[1] = tmp1*pro[1];
620 temp_PDF=0;
621 if(g0[i]==1) propagatorsigma500(ssigma2,spion2,spim,pro0);
622 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
623 //if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],seta12952,ssigma2,seta,rRes2,0,pro1);
624 if(g1[i]==1) { pro1[0] = 1; pro1[1] = 0; }
625 Com_Multi(pro0,pro1,pro);
626 calt1(pf12,Pip1,t1D);
627 tmp2 = 1;
628 amp_tmp2[0] = tmp2*pro[0];
629 amp_tmp2[1] = tmp2*pro[1];
630 }
631
632 if(modetype[i]==10)
633 {
634 if(g0[i]==1) propagatora0980(mass1[i],sa0m,seta2,spim2,pro0);
635 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0;pro2[0]=1;pro2[1]=0;}
636 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],sf11,sa0m,spion1,rRes2,1,pro1);
637 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0;pro3[0]=1;pro3[1]=0;}
638 Com_Multi(pro0,pro1,pro);
639 calt1(pf11,Pip2,t1D);
640 B[1] = Barrier(1,sf11,sa0m,spion1,rRes2);
641 B[2] = Barrier(1,sD,sf11,spion2,rD2);
642 for(int a=0; a<4; a++)
643 {
644 temp_PDF += t1D[a]*t1f11[a]*G[a][a];
645 }
646 tmp1 = B[1]*B[2]*temp_PDF;
647 amp_tmp1[0] = tmp1*pro[0];
648 amp_tmp1[1] = tmp1*pro[1];
649
650 temp_PDF=0;
651 if(g0[i]==1) propagatora0980(mass1[i],sa0m,seta2,spim2,pro0);
652 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; pro2[0]=1; pro2[1]=0;}
653 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],sf12,sa0m,spion2,rRes2,1,pro1);
654 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0;pro3[0]=1;pro3[1]=0; }
655 Com_Multi(pro0,pro1,pro);
656 calt1(pf12,Pip1,t1D);
657 B[1] = Barrier(1,sf12,sa0m,spion2,rRes2);
658 B[2] = Barrier(1,sD,sf12,spion1,rD2);
659 for(int a=0; a<4; a++)
660 {
661 temp_PDF += t1D[a]*t1f12[a]*G[a][a];
662 }
663 tmp2 = B[1]*B[2]*temp_PDF;
664 amp_tmp2[0] = tmp2*pro[0];
665 amp_tmp2[1] = tmp2*pro[1];
666
667 }
668
669 if(modetype[i]==11)
670 {
671 if(g0[i]==1) propagatora0980(mass1[i],sa01,seta2,spion12,pro2);
672 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0;pro2[0]=1;pro2[1]=0;}
673 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],sf11,sa01,spim,rRes2,1,pro3);
674 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0;pro3[0]=1;pro3[1]=0;}
675 Com_Multi(pro2,pro3,pro4);
676 calt1(pf11,Pip2,t1D);
677 Bc[1] = Barrier(1,sf11,sa01,spim,rRes2);
678 Bc[2] = Barrier(1,sD,sf11,spion2,rD2);
679 for(int a=0; a<4; a++)
680 {
681 temp_PDF1 += t1D[a]*t1f11_2[a]*G[a][a];
682 }
683 tmp3 = Bc[1]*Bc[2]*temp_PDF1;
684 amp_tmp1[0] = tmp3*pro4[0];
685 amp_tmp1[1] = tmp3*pro4[1];
686 //printf("amp_tmp1 = %f , %f\n", amp_tmp1[0], amp_tmp1[1]);
687
688 temp_PDF1=0;
689 if(g0[i]==1) propagatora0980(mass1[i],sa02,seta2,spion22,pro2);
690 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; pro2[0]=1; pro2[1]=0;}
691 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],sf12,sa02,spim,rRes2,1,pro3);
692 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0;pro3[0]=1;pro3[1]=0; }
693 Com_Multi(pro2,pro3,pro4);
694 calt1(pf12,Pip1,t1D);
695 Bc[1] = Barrier(1,sf12,sa02,spim,rRes2);
696 Bc[2] = Barrier(1,sD,sf12,spion1,rD2);
697 for(int a=0; a<4; a++)
698 {
699 temp_PDF1 += t1D[a]*t1f12_2[a]*G[a][a];
700 }
701 tmp4 = Bc[1]*Bc[2]*temp_PDF1;
702 amp_tmp2[0] = tmp4*pro4[0];
703 amp_tmp2[1] = tmp4*pro4[1];
704 //printf("amp_tmp2 = %f , %f\n", amp_tmp2[0], amp_tmp2[1]);
705 }
706 amp_tmp[0] = amp_tmp1[0]+amp_tmp2[0];
707 amp_tmp[1] = amp_tmp1[1]+amp_tmp2[1];
708 Com_Multi(amp_tmp,cof,amp_PDF);
709 PDF[0] += amp_PDF[0];
710 PDF[1] += amp_PDF[1];
711 }
712
713 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
714 if(value <=0) {
715 value = 1e-20;
716 }
717 Result = value;
718}
719
720void EvtDsToEta3pi::Com_Multi(double a1[2], double a2[2], double res[2])
721{
722 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
723 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
724}
725void EvtDsToEta3pi::Com_Divide(double a1[2], double a2[2], double res[2])
726{
727 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
728 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
729 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
730}
731double EvtDsToEta3pi::SCADot(double a1[4], double a2[4])
732{
733 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
734 return _cal;
735}
736double EvtDsToEta3pi::Barrier(int l, double sa, double sb, double sc, double r2)
737{
738 double F;
739 double tmp = sa+sb-sc;
740 double q = 0.25*tmp*tmp/sa-sb;
741 if (q < 0) q = 1e-16;
742 double z = q*r2;
743 if (l==1) {
744 F = sqrt(2.0*z/(1.0+z));
745 }
746 else if (l==2) {
747 double z2 = z*z;
748 F = sqrt(13.0*z2/(9.0+3.0*z+z2));
749 } else {
750 F = 1.0;
751 }
752 return F;
753}
754void EvtDsToEta3pi::calt1(double daug1[4], double daug2[4], double t1[4])
755{
756 double p, pq, tmp;
757 double pa[4], qa[4];
758 for(int i=0; i<4; i++) {
759 pa[i] = daug1[i] + daug2[i];
760 qa[i] = daug1[i] - daug2[i];
761 }
762 p = SCADot(pa,pa);
763 pq = SCADot(pa,qa);
764 tmp = pq/p;
765 for(int i=0; i<4; i++) {
766 t1[i] = qa[i] - tmp*pa[i];
767 }
768}
769
770void EvtDsToEta3pi::calt2(double daug1[4], double daug2[4], double t2[4][4])
771{
772 double p, r;
773 double pa[4], t1[4];
774 calt1(daug1,daug2,t1);
775 r = SCADot(t1,t1);
776 for(int i=0; i<4; i++)
777 {
778 pa[i] = daug1[i] + daug2[i];
779 }
780 p = SCADot(pa,pa);
781 for(int i=0; i<4; i++)
782 {
783 for(int j=0; j<4; j++)
784 {
785 t2[i][j] = t1[i]*t1[j] - 1.0/3*r*(G[i][j]-pa[i]*pa[j]/p);
786 }
787 }
788}
789
790void EvtDsToEta3pi::propagator(double mass2, double mass, double width, double sx, double prop[2])
791{
792 double a[2], b[2];
793 a[0] = 1;
794 a[1] = 0;
795 b[0] = mass2-sx;
796 b[1] = -mass*width;
797 Com_Divide(a,b,prop);
798}
799double EvtDsToEta3pi::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
800{
801 double widm = 0.;
802 double m = sqrt(sa);
803 double tmp = sb-sc;
804 double tmp1 = sa+tmp;
805 double q = 0.25*tmp1*tmp1/sa-sb;
806 if(q<0) q = 1e-16;
807 double tmp2 = mass2+tmp;
808 double q0 = 0.25*tmp2*tmp2/mass2-sb;
809 if(q0<0) q0 = 1e-16;
810 double z = q*r2;
811 double z0 = q0*r2;
812 double t = q/q0;
813 if(l == 0) {widm = sqrt(t)*mass/m;}
814 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
815 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
816 return widm;
817}
818double EvtDsToEta3pi::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
819{
820 double widm = 0.;
821 double m = sqrt(sa);
822 double tmp = sb-sc;
823 double tmp1 = sa+tmp;
824 double q = 0.25*tmp1*tmp1/sa-sb;
825 if(q<0) q = 1e-16;
826 double tmp2 = mass2+tmp;
827 double q0 = 0.25*tmp2*tmp2/mass2-sb;
828 if(q0<0) q0 = 1e-16;
829 double z = q*r2;
830 double z0 = q0*r2;
831 double F = (1+z0)/(1+z);
832 double t = q/q0;
833 widm = t*sqrt(t)*mass/m*F;
834 return widm;
835}
836void EvtDsToEta3pi::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
837{
838 double a[2], b[2];
839 a[0] = 1;
840 a[1] = 0;
841 b[0] = mass2-sa;
842 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
843 Com_Divide(a,b,prop);
844}
845
846
847void EvtDsToEta3pi::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
848{
849 double a[2], b[2];
850 double tmp = sb-sc;
851 double tmp1 = sa+tmp;
852 double q2 = 0.25*tmp1*tmp1/sa-sb;
853 if(q2<0) q2 = 1e-16;
854
855 double tmp2 = mass2+tmp;
856 double q02 = 0.25*tmp2*tmp2/mass2-sb;
857 if(q02<0) q02 = 1e-16;
858
859 double q = sqrt(q2);
860 double q0 = sqrt(q02);
861 double m = sqrt(sa);
862 double q03 = q0*q02;
863 double tmp3 = log(mass+2*q0)+1.2760418309; // log(mass_2Pion) = 1.2760418309;
864
865 double h = GS1*q/m*(log(m+2*q)+1.2760418309);
866 double h0 = GS1*q0/mass*tmp3;
867 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
868 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
869 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
870
871 a[0] = 1.0+d*width/mass;
872 a[1] = 0.0;
873 b[0] = mass2-sa+width*f;
874 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
875 Com_Divide(a,b,prop);
876}
877void EvtDsToEta3pi::rhoab(double sa, double sb, double sc, double res[2]) {
878 double tmp = sa+sb-sc;
879 double q = 0.25*tmp*tmp/sa-sb;
880 if(q>=0) {
881 res[0]=2.0*sqrt(q/sa);
882 res[1]=0.0;
883 } else {
884 res[0]=0.0;
885 res[1]=2.0*sqrt(-q/sa);
886 }
887}
888void EvtDsToEta3pi::rho4Pi(double sa, double res[2]) {
889 double temp = 1.0-0.3116765584/sa; // 0.3116765584=0.13957*0.13957*16
890 if(temp>=0) {
891 res[0]=sqrt(temp)/(1.0+exp(9.8-3.5*sa));
892 res[1]=0.0;
893 } else {
894 res[0]=0.0;
895 res[1]=sqrt(-temp)/(1.0+exp(9.8-3.5*sa));
896 }
897}
898void EvtDsToEta3pi::propagatorsigma500(double sa, double sb, double sc, double prop[2]) {
899 double f = 0.5843+1.6663*sa;
900 const double M = 0.9264;
901 const double mass2 = 0.85821696; // M*M
902 const double mpi2d2 = 0.00973989245;
903 double g1 = f*(sa-mpi2d2)/(mass2-mpi2d2)*exp((mass2-sa)/1.082);
904 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
905 rhoab(sa,sb,sc,rho1s);
906 rhoab(mass2,sb,sc,rho1M);
907 rho4Pi(sa,rho2s);
908 rho4Pi(mass2,rho2M);
909 Com_Divide(rho1s,rho1M,rho1);
910 Com_Divide(rho2s,rho2M,rho2);
911 double a[2], b[2];
912 a[0] = 1.0;
913 a[1] = 0.0;
914 b[0] = mass2-sa+M*(g1*rho1[1]+0.0024*rho2[1]);
915 b[1] = -M*(g1*rho1[0]+0.0024*rho2[0]);
916 Com_Divide(a,b,prop);
917}
918void EvtDsToEta3pi::propagatora0_980(double mass2, double mass, double width, double sa, double sb, double sc, int charge, double prop[2]){
919 double q, qKsK;
920 double rho[2], rhoKsK[2];
921 rhoab(sa,sb,sc,rho);
922 if(charge == 0) {
923 qKsK = 0.25*sa-0.2437199424;
924 }
925 else if(charge == 1) {
926 double tmp1 = sa+3.899750596e-03;
927 qKsK = 0.25*tmp1*tmp1/sa-0.247619692996;
928 }
929 if(qKsK>0){
930 rhoKsK[0] = 2*sqrt(qKsK/sa);
931 rhoKsK[1] = 0;
932 }
933 else if(qKsK<0){
934 rhoKsK[0] = 0;
935 rhoKsK[1] = 2*sqrt(-qKsK/sa);
936 }
937 double a[2], b[2];
938 a[0] = 1;
939 a[1] = 0;
940 b[0] = mass2-sa+0.341*rho[1]+0.304172*rhoKsK[1]; // 0.304172=0.892*0.341
941 b[1] = -0.341*rho[0]-0.304172*rhoKsK[0];
942 Com_Divide(a,b,prop);
943}
944void EvtDsToEta3pi::Flatte_rhoab(double sa, double sb, double sc, double rho[2])
945{
946 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
947 double b[2];
948 if(q>0) {
949 rho[0]=2* sqrt(q/sa);
950 rho[1]=0;
951 }
952 else if(q<0){
953 rho[0]=0;
954 rho[1]=2*sqrt(-q/sa);
955 }
956}
957void EvtDsToEta3pi::propagator980(double mass, double sx, double *sb, double *sc, double prop[2])
958{
959 double unit[2]={1.0};
960 double ci[2]={0,1};
961 double rho1[2];
962 Flatte_rhoab(sx,sb[0],sc[0],rho1);
963 double rho2[2];
964 Flatte_rhoab(sx,sb[1],sc[1],rho2);
965
966 double gK_f980=0.69466, gPi_f980=0.165;
967 double tmp1[2]={gK_f980,0};
968 double tmp11[2];
969 double tmp2[2]={gPi_f980,0};
970 double tmp22[2];
971 Com_Multi(tmp1,rho1,tmp11);
972 Com_Multi(tmp2,rho2,tmp22);
973 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
974 double tmp31[2];
975 Com_Multi(tmp3, ci,tmp31);
976 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
977 Com_Divide( unit,tmp4, prop);
978}
979
980void EvtDsToEta3pi::propagatora0980(double mass, double sx, double *sb, double *sc, double prop[2])
981{
982 double unit[2]={1.0};
983 double ci[2]={0,1};
984 double rho1[2];
985 Flatte_rhoab(sx,sb[0],sc[0],rho1);
986 double rho2[2];
987 Flatte_rhoab(sx,sb[1],sc[1],rho2);
988
989 double gKK_a980=0.892*0.341, gPiEta_a980=0.341;
990 double tmp1[2]={gKK_a980,0};
991 double tmp11[2];
992 double tmp2[2]={gPiEta_a980,0};
993 double tmp22[2];
994 Com_Multi(tmp1,rho1,tmp11);
995 Com_Multi(tmp2,rho2,tmp22);
996 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
997 double tmp31[2];
998 Com_Multi(tmp3, ci,tmp31);
999 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
1000 Com_Divide( unit,tmp4, prop);
1001}
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")
g1
TF1 * g1
Definition DataBase/tau_mode.c:66
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:252
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
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
EvtDsToEta3pi::~EvtDsToEta3pi
virtual ~EvtDsToEta3pi()
Definition EvtDsToEta3pi.cc:33
EvtDsToEta3pi::clone
EvtDecayBase * clone()
Definition EvtDsToEta3pi.cc:39
EvtDsToEta3pi::getName
void getName(std::string &name)
Definition EvtDsToEta3pi.cc:35
EvtDsToEta3pi::decay
void decay(EvtParticle *p)
Definition EvtDsToEta3pi.cc:180
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
charge
float charge
Definition prepare_data.cxx:28
b
const double b
Definition slope.cxx:9