BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKppipipi0.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: EvtDToKppipipi0.cc
11// the necessary file: EvtDToKppipipi0.hh
12//
13// Description: D+ -> K+ pi+ pi- pi0, see DocDB-1221
14//
15// Modification history:
16//
17// Liaoyuan Dong Dec. 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/EvtDToKppipipi0.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
32EvtDToKppipipi0::~EvtDToKppipipi0() {}
33
34void EvtDToKppipipi0::getName(std::string& model_name){
35 model_name="DToKppipipi0";
36}
37
41
42void EvtDToKppipipi0::init(){
43 checkNArg(0);
44 checkNDaug(4);
45 checkSpinParent(EvtSpinType::SCALAR);
46 //cout << "Initializing EvtDToKppipipi0" << endl;
47
48 mK1270 = 1.272; mK1400 = 1.403;
49 GK1270 = 0.087; GK1400 = 0.174;
50 mKstr0 = 0.89555; mrhop = 0.77511;
51 GKstr0 = 0.0473; Grhop = 0.1491;
52 mKstrp = 0.89167; mrho0 = 0.77526;
53 GKstrp = 0.0514; Grho0 = 0.1474;
54 mOmega = 0.78266;
55 GOmega = 0.00868;
56 msigma = 0.526;
57 Gsigma = 0.535;
58 mK1650 = 1.65;
59 GK1650 = 0.158;
60 ma1 = 1.230;
61 Ga1 = 0.420;
62
63 rho[7] = 1.0;
64 phi[7] = 0.0; //K1400+
65 rho[0] = 0.3578145188;
66 phi[0] =-1.8624743029; //1S
67 rho[1] = 0.3468367648;
68 phi[1] = 0.1168108412; //1P
69 rho[2] = 0.4257117223;
70 phi[2] =-0.7433416190; //2S
71 rho[3] = 0.3042749891;
72 phi[3] = 4.3671838579; //2P
73 rho[4] = 1.3376376890;
74 phi[4] =-3.0855437414; //3S
75 rho[5] = 1.3376376890;
76 phi[5] =-3.0855437414; //3S
77 rho[6] = 0.5098861534;
78 phi[6] =-1.7222823562; //5S
79 rho[8] = 0.6422897385;
80 phi[8] = 2.5529952524; //Komega
81
82 modetype[0]=1;
83 modetype[1]=1;
84 modetype[2]=2;
85 modetype[3]=2;
86 modetype[4]=3;
87 modetype[5]=4;
88 modetype[6]=5;
89 modetype[7]=9;
90 modetype[8]=11;
91
92 //for (int i=0; i<9; i++) {
93 // cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
94 //}
95
96 mD = 1.86966;
97 rD = 5;
98 math_pi = 3.1415926;
99
100 pi = 3.1415926;
101
102 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
103 int EE[4][4][4][4] =
104 { { {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
105 {{0,0,0,0}, {0,0,0,0}, {0,0,0,1}, {0,0,-1,0}},
106 {{0,0,0,0}, {0,0,0,-1}, {0,0,0,0}, {0,1,0,0} },
107 {{0,0,0,0}, {0,0,1,0}, {0,-1,0,0}, {0,0,0,0} } },
108 { {{0,0,0,0}, {0,0,0,0}, {0,0,0,-1}, {0,0,1,0} },
109 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
110 {{0,0,0,1}, {0,0,0,0}, {0,0,0,0}, {-1,0,0,0}},
111 {{0,0,-1,0}, {0,0,0,0}, {1,0,0,0}, {0,0,0,0} } },
112 { {{0,0,0,0}, {0,0,0,1}, {0,0,0,0}, {0,-1,0,0}},
113 {{0,0,0,-1}, {0,0,0,0}, {0,0,0,0}, {1,0,0,0} },
114 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
115 {{0,1,0,0}, {-1,0,0,0}, {0,0,0,0}, {0,0,0,0} } },
116 { {{0,0,0,0}, {0,0,-1,0}, {0,1,0,0}, {0,0,0,0} },
117 {{0,0,1,0}, {0,0,0,0}, {-1,0,0,0}, {0,0,0,0} },
118 {{0,-1,0,0}, {1,0,0,0}, {0,0,0,0}, {0,0,0,0} },
119 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} } } };
120 for (int i=0; i<4; i++) {
121 for (int j=0; j<4; j++) {
122 G[i][j] = GG[i][j];
123 for (int k=0; k<4; k++) {
124 for (int l=0; l<4; l++) {
125 E[i][j][k][l] = EE[i][j][k][l];
126 }
127 }
128 }
129 }
130}
131
132void EvtDToKppipipi0::initProbMax() {
133 setProbMax(100000);//setProbMax(40000.0);
134}
135
136void EvtDToKppipipi0::decay( EvtParticle *p){
137 //-----------for max value------------------
138/* double maxprob = 0.0;
139 for(int ir=0;ir<=60000000;ir++){
140 p->initializePhaseSpace(getNDaug(),getDaugs());
141 EvtVector4R kp = p->getDaug(0)->getP4();
142 EvtVector4R pip = p->getDaug(1)->getP4();
143 EvtVector4R pim = p->getDaug(2)->getP4();
144 EvtVector4R pi0 = p->getDaug(3)->getP4();
145 double value;
146 double Kp[4],Pip[4],Pim[4],Pi0[4];
147 Kp[0] = kp.get(0); Pip[0] = pip.get(0); Pim[0] = pim.get(0); Pi0[0] = pi0.get(0);
148 Kp[1] = kp.get(1); Pip[1] = pip.get(1); Pim[1] = pim.get(1); Pi0[1] = pi0.get(1);
149 Kp[2] = kp.get(2); Pip[2] = pip.get(2); Pim[2] = pim.get(2); Pi0[2] = pi0.get(2);
150 Kp[3] = kp.get(3); Pip[3] = pip.get(3); Pim[3] = pim.get(3); Pi0[3] = pi0.get(3);
151 //double Prob = calPDF(Kp, Pip, Pim, Pi0);
152 calPDF(Kp, Pip, Pim, Pi0, value);
153 if(value>maxprob) {
154 maxprob=value;
155 std::cout << "Max PDF = " << ir << " prob= " << value << std::endl;
156 }
157 }
158 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
159 return;*/
160 //-----------------------------------------------
161 p->initializePhaseSpace(getNDaug(),getDaugs());
162 EvtVector4R kp = p->getDaug(0)->getP4();
163 EvtVector4R pip = p->getDaug(1)->getP4();
164 EvtVector4R pim = p->getDaug(2)->getP4();
165 EvtVector4R pi0 = p->getDaug(3)->getP4();
166
167 double Kp[4],Pip[4],Pim[4],Pi0[4];
168 Kp[0] = kp.get(0); Pip[0] = pip.get(0); Pim[0] = pim.get(0); Pi0[0] = pi0.get(0);
169 Kp[1] = kp.get(1); Pip[1] = pip.get(1); Pim[1] = pim.get(1); Pi0[1] = pi0.get(1);
170 Kp[2] = kp.get(2); Pip[2] = pip.get(2); Pim[2] = pim.get(2); Pi0[2] = pi0.get(2);
171 Kp[3] = kp.get(3); Pip[3] = pip.get(3); Pim[3] = pim.get(3); Pi0[3] = pi0.get(3);
172
173 double value;
174 calPDF(Kp, Pip, Pim, Pi0, value);
175 setProb(value);
176 return;
177
178}
179
180double EvtDToKppipipi0::calPDF(double Kp[], double Pip[], double Pim[], double Pi0[], double & Result) {
181 double cof[2], amp_tmp[2], amp_PDF[2], PDF[2];
182 double flag[3], mass_R[2], width_R[2];
183 double pro[2];
184 double sa[3], sb[3], sc[3], B[3];
185 double t1D[4], t1V[4], t1A[4];
186 double pS[4], pV[4], pD[4], pA[4];
187 double pro1[2], pro2[2],proKPi_S[2];
188
189 double rD2 = 25.0;
190 double rRes2 = 9.0;
191 double mass1[9] ={mKstr0, mKstr0, mKstrp, mKstrp, ma1, ma1, mK1270, mK1400, mOmega };
192 double mass2[9] ={ mrhop, mrhop, mrho0, mrho0, mrhop, mrhop, mrhop, mKstr0, mrho0 };
193 double width1[9]={GKstr0, GKstr0, GKstrp, GKstrp, Ga1, Ga1, GK1270, GK1400, GOmega };
194 double width2[9]={ Grhop, Grhop, Grho0, Grho0, Grhop, Grhop, Grhop, GKstr0, Grho0 };
195 double g0[9] ={ 1, 1, 1, 1, 1, 1, 1, 1, 1 };
196 double g1[9] ={ 1, 1, 1, 1, 1, 1, 1, 1, 1 };
197 double g2[9] ={ 0, 1, 0, 1, 0, 0, 0, 0, 1 };
198 double temp_PDF = 0;
199 PDF[0]=0;
200 PDF[1]=0;
201
202 for(int i=0; i<9; i++) {
203 flag[0] = g0[i]; flag[1] = g1[i];flag[2] = g2[i];
204 mass_R[0] = mass1[i]; mass_R[1] = mass2[i];
205 width_R[0] = width1[i]; width_R[1] = width2[i];
206
207 amp_tmp[0] = 0;
208 amp_tmp[1] = 0;
209
210 cof[0] = rho[i]*cos(phi[i]);
211 cof[1] = rho[i]*sin(phi[i]);
212
213 if(modetype[i] == 1){
214 double t1V1[4], t1V2[4], t2D[4][4];
215 temp_PDF = 0;
216 double pV1[4], pV2[4];
217 for(int ii=0; ii!=4; ii++){
218 pV1[ii] = Kp[ii] + Pim[ii];
219 pV2[ii] = Pip[ii] + Pi0[ii];
220 pD[ii] = pV1[ii] + pV2[ii];
221 }
222 sa[0] = SCADot(pV1,pV1);
223 sb[0] = SCADot(Kp,Kp);
224 sc[0] = SCADot(Pim,Pim);
225 sa[1] = SCADot(pV2,pV2);
226 sb[1] = SCADot(Pip,Pip);
227 sc[1] = SCADot(Pi0,Pi0);
228 sa[2] = SCADot(pD,pD);
229 sb[2] = sa[0];
230 sc[2] = sa[1];
231 if(flag[0] == 1){
232 propagatorRBW(mass_R[0]*mass_R[0],mass_R[0],width_R[0],sa[0],sb[0],sc[0],rRes2,1, pro1);
233 }
234 else if(flag[0] == 0){pro1[0] = 1; pro1[1] = 0;}
235 if(flag[1] == 1){
236 propagatorGS(mass_R[1]*mass_R[1],mass_R[1],width_R[1],sa[1],sb[1],sc[1],rRes2, pro2);
237 }
238 else if(flag[1] == 0){pro2[0] = 1; pro2[1] = 0;}
239 B[0] = Barrier(mass_R[0]*mass_R[0],1,sa[0],sb[0],sc[0],rRes2);
240 B[1] = Barrier(mass_R[1]*mass_R[1],1,sa[1],sb[1],sc[1],rRes2);
241 calt1(Kp,Pim,t1V1);
242 calt1(Pip,Pi0,t1V2);
243 if(flag[2] == 0){
244 for(int ii=0; ii!=4; ii++){
245 temp_PDF += (G[ii][ii])*t1V1[ii]*t1V2[ii];
246 }
247 B[2] = 1;
248 }
249 if(flag[2] == 1){
250 calt1(pV1,pV2,t1D);
251 for(int ii=0; ii!=4; ii++){
252 for(int j=0; j!=4; j++){
253 for(int k=0; k!=4; k++){
254 for(int l=0; l!=4; l++){
255 temp_PDF += E[ii][j][k][l]*pD[ii]*t1D[j]*t1V1[k]*t1V2[l]*
256 (G[ii][ii])*(G[j][j])*(G[l][l])*(G[k][k]);
257 }
258 }
259 }
260 }
261 B[2] = Barrier(mD*mD,1,sa[2],sb[2],sc[2],rD2);
262 }
263 if(flag[2] == 2){
264 calt2(pV1,pV2,t2D);
265 for(int ii=0; ii!=4; ii++){
266 for(int j=0; j!=4; j++){
267 temp_PDF += t2D[ii][j]*t1V1[ii]*t1V2[j]*(G[ii][ii])*(G[j][j]);
268 }
269 }
270 B[2] = Barrier(mD*mD,2,sa[2],sb[2],sc[2],rD2);
271 }
272 Com_Multi(pro1,pro2,pro);
273 amp_tmp[0] = temp_PDF*B[0]*B[1]*B[2]*pro[0];
274 amp_tmp[1] = temp_PDF*B[0]*B[1]*B[2]*pro[1];
275
276 }
277 else if(modetype[i] == 2){
278 double t1V1[4], t1V2[4], t2D[4][4];
279 temp_PDF = 0;
280 double pV1[4], pV2[4];
281 for(int ii=0; ii!=4; ii++){
282 pV1[ii] = Kp[ii] + Pi0[ii];
283 pV2[ii] = Pip[ii] + Pim[ii];
284 pD[ii] = pV1[ii] + pV2[ii];
285 }
286 sa[0] = SCADot(pV1,pV1);
287 sb[0] = SCADot(Kp,Kp);
288 sc[0] = SCADot(Pi0,Pi0);
289 sa[1] = SCADot(pV2,pV2);
290 sb[1] = SCADot(Pip,Pip);
291 sc[1] = SCADot(Pim,Pim);
292 sa[2] = SCADot(pD,pD);
293 sb[2] = sa[0];
294 sc[2] = sa[1];
295 if(flag[0] == 1){
296 propagatorRBW(mass_R[0]*mass_R[0],mass_R[0],width_R[0],sa[0],sb[0],sc[0],rRes2,1, pro1);
297 }
298 else if(flag[0] == 0){pro1[0] = 1; pro1[1] = 0;}
299 if(flag[1] == 1){
300 propagatorGS(mass_R[1]*mass_R[1],mass_R[1],width_R[1],sa[1],sb[1],sc[1],rRes2, pro2);
301 }
302 else if(flag[1] == 0){pro2[0] = 1; pro2[1] = 0;}
303 B[0] = Barrier(mass_R[0]*mass_R[0],1,sa[0],sb[0],sc[0],rRes2);
304 B[1] = Barrier(mass_R[1]*mass_R[1],1,sa[1],sb[1],sc[1],rRes2);
305 calt1(Kp,Pi0,t1V1);
306 calt1(Pip,Pim,t1V2);
307 if(flag[2] == 0){
308 for(int ii=0; ii!=4; ii++){
309 temp_PDF += (G[ii][ii])*t1V1[ii]*t1V2[ii];
310 }
311 B[2] = 1;
312 }
313 if(flag[2] == 1){
314 calt1(pV1,pV2,t1D);
315 for(int ii=0; ii!=4; ii++){
316 for(int j=0; j!=4; j++){
317 for(int k=0; k!=4; k++){
318 for(int l=0; l!=4; l++){
319 temp_PDF += E[ii][j][k][l]*pD[ii]*t1D[j]*t1V1[k]*t1V2[l]*
320 (G[ii][ii])*(G[j][j])*(G[l][l])*(G[k][k]);
321 }
322 }
323 }
324 }
325 B[2] = Barrier(mD*mD,1,sa[2],sb[2],sc[2],rD2);
326 }
327 if(flag[2] == 2){
328 calt2(pV1,pV2,t2D);
329 for(int ii=0; ii!=4; ii++){
330 for(int j=0; j!=4; j++){
331 temp_PDF += t2D[ii][j]*t1V1[ii]*t1V2[j]*(G[ii][ii])*(G[j][j]);
332 }
333 }
334 B[2] = Barrier(mD*mD,2,sa[2],sb[2],sc[2],rD2);
335 }
336 Com_Multi(pro1,pro2,pro);
337 amp_tmp[0] = temp_PDF*B[0]*B[1]*B[2]*pro[0];
338 amp_tmp[1] = temp_PDF*B[0]*B[1]*B[2]*pro[1];
339
340 }else if(modetype[i] == 3){
341 temp_PDF = 0;
342 double t2A[4][4];
343 for(int ii=0; ii!=4; ii++){
344 pV[ii] = Pip[ii] + Pi0[ii];
345 pA[ii] = pV[ii] + Pim[ii];
346 pD[ii] = pA[ii] + Kp[ii];
347 }
348 sa[0] = SCADot(pV,pV);
349 sb[0] = SCADot(Pip,Pip);
350 sc[0] = SCADot(Pi0,Pi0);
351 sa[1] = SCADot(pA,pA);
352 sb[1] = sa[0];
353 sc[1] = SCADot(Pim,Pim);
354 sa[2] = SCADot(pD,pD);
355 sb[2] = sa[1];
356 sc[2] = SCADot(Kp,Kp);
357
358 if(flag[0] == 1){
359 propagatorGS(mass_R[1]*mass_R[1],mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2, pro1);
360 }
361 else if(flag[0] == 0){pro1[0] = 1; pro1[1] = 0;}
362
363 if(flag[1] == 1){
364 propagatorRBW_a1(mass_R[0]*mass_R[0],mass_R[0], width_R[0], sa[1], sb[1], sc[1], rRes2, flag[2], pro2);
365 }
366 else if(flag[1] == 0){pro2[0] = 1; pro2[1] = 0;}
367
368 B[1] = Barrier(mass_R[1]*mass_R[1], 1,sa[0],sb[0],sc[0],rRes2);
369 B[2] = Barrier(mD*mD, 1,sa[2],sb[2],sc[2],rD2);
370 calt1(Pip,Pi0,t1V);
371 calt1(pA,Kp,t1D);
372 if(flag[2] == 0){
373 for(int ii=0; ii!=4; ii++){
374 for(int j=0; j!=4; j++){
375 temp_PDF += t1D[ii]*(G[ii][j] - pA[ii]*pA[j]/sa[1])*t1V[j]*G[ii][ii]*G[j][j];
376 }
377 }
378 B[0] = 1;
379 }
380 else if(flag[2] == 2){
381 calt2(pV,Pim,t2A);
382 for(int ii=0; ii!=4; ii++){
383 for(int j=0; j!=4; j++){
384 temp_PDF += t1D[ii]*t2A[ii][j]*t1V[j]*G[ii][ii]*G[j][j];
385 }
386 }
387 B[0] = Barrier(mass_R[0]*mass_R[0], 2,sa[1],sb[1],sc[1],rRes2);
388 }
389 Com_Multi(pro1,pro2,pro);
390 amp_tmp[0] = temp_PDF*B[0]*B[1]*B[2]*pro[0];
391 amp_tmp[1] = temp_PDF*B[0]*B[1]*B[2]*pro[1];
392
393 }else if(modetype[i] == 4){
394 temp_PDF = 0;
395 double t2A[4][4];
396 for(int ii=0; ii!=4; ii++){
397 pV[ii] = Pim[ii] + Pi0[ii];
398 pA[ii] = pV[ii] + Pip[ii];
399 pD[ii] = pA[ii] + Kp[ii];
400 }
401 sa[0] = SCADot(pV,pV);
402 sb[0] = SCADot(Pim,Pim);
403 sc[0] = SCADot(Pi0,Pi0);
404 sa[1] = SCADot(pA,pA);
405 sb[1] = sa[0];
406 sc[1] = SCADot(Pip,Pip);
407 sa[2] = SCADot(pD,pD);
408 sb[2] = sa[1];
409 sc[2] = SCADot(Kp,Kp);
410
411 if(flag[0] == 1){
412 propagatorGS(mass_R[1]*mass_R[1],mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2, pro1);
413 }
414 else if(flag[0] == 0){pro1[0] = 1; pro1[1] = 0;}
415
416 if(flag[1] == 1){
417 propagatorRBW_a1(mass_R[0]*mass_R[0],mass_R[0], width_R[0], sa[1], sb[1], sc[1], rRes2, flag[2], pro2);
418 }
419 else if(flag[1] == 0){pro2[0] = 1; pro2[1] = 0;}
420
421 B[1] = Barrier(mass_R[1]*mass_R[1], 1,sa[0],sb[0],sc[0],rRes2);
422 B[2] = Barrier(mD*mD, 1,sa[2],sb[2],sc[2],rD2);
423 calt1(Pim,Pi0,t1V);
424 calt1(pA,Kp,t1D);
425 if(flag[2] == 0){
426 for(int ii=0; ii!=4; ii++){
427 for(int j=0; j!=4; j++){
428 temp_PDF += t1D[ii]*(G[ii][j] - pA[ii]*pA[j]/sa[1])*t1V[j]*G[ii][ii]*G[j][j];
429 }
430 }
431 B[0] = 1;
432 }
433 else if(flag[2] == 2){
434 calt2(pV,Pip,t2A);
435 for(int ii=0; ii!=4; ii++){
436 for(int j=0; j!=4; j++){
437 temp_PDF += t1D[ii]*t2A[ii][j]*t1V[j]*G[ii][ii]*G[j][j];
438 }
439 }
440 B[0] = Barrier(mass_R[0]*mass_R[0], 2,sa[1],sb[1],sc[1],rRes2);
441 }
442 Com_Multi(pro1,pro2,pro);
443 amp_tmp[0] = temp_PDF*B[0]*B[1]*B[2]*pro[0];
444 amp_tmp[1] = temp_PDF*B[0]*B[1]*B[2]*pro[1];
445
446 }else if(modetype[i] == 5){
447 temp_PDF = 0;
448 double t2A[4][4];
449 for(int ii=0; ii!=4; ii++){
450 pV[ii] = Pim[ii] + Pi0[ii];
451 pA[ii] = pV[ii] + Kp[ii];
452 pD[ii] = pA[ii] + Pip[ii];
453 }
454 sa[0] = SCADot(pV,pV);
455 sb[0] = SCADot(Pim,Pim);
456 sc[0] = SCADot(Pi0,Pi0);
457 sa[1] = SCADot(pA,pA);
458 sb[1] = sa[0];
459 sc[1] = SCADot(Kp,Kp);
460 sa[2] = SCADot(pD,pD);
461 sb[2] = sa[1];
462 sc[2] = SCADot(Pip,Pip);
463 if(flag[0] == 1){
464 propagatorGS(mass_R[1]*mass_R[1],mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2, pro1);
465 }
466 else if(flag[0] == 0){pro1[0] = 1; pro1[1] = 0;}
467
468 if(flag[1] == 1){
469 propagatorRBW(mass_R[0]*mass_R[0],mass_R[0], width_R[0], sa[1], sb[1], sc[1], rRes2, flag[2], pro2);
470 }
471 else if(flag[1] == 0){pro2[0] = 1; pro2[1] = 0;}
472 B[1] = Barrier(mass_R[1]*mass_R[1],1,sa[0],sb[0],sc[0],rRes2);//rho
473 B[2] = Barrier(mD*mD,1,sa[2],sb[2],sc[2],rD2);
474 calt1(Pim,Pi0,t1V);
475 calt1(pA,Pip,t1D);
476 if(flag[2] == 0){
477 for(int ii=0; ii!=4; ii++){
478 for(int j=0; j!=4; j++){
479 temp_PDF += t1D[ii]*(G[ii][j] - pA[ii]*pA[j]/sa[1])*t1V[j]*G[ii][ii]*G[j][j];
480 }
481 }
482 B[0] = 1;
483 }
484 else if(flag[2] == 2){
485 calt2(pV,Kp,t2A);
486 for(int ii=0; ii!=4; ii++){
487 for(int j=0; j!=4; j++){
488 temp_PDF += t1D[ii]*t2A[ii][j]*t1V[j]*G[ii][ii]*G[j][j];
489 }
490 }
491 B[0] = Barrier(mass_R[0]*mass_R[0],2,sa[1],sb[1],sc[1],rRes2);
492 }
493 Com_Multi(pro1,pro2,pro);
494 amp_tmp[0] = temp_PDF*B[0]*B[1]*B[2]*pro[0];
495 amp_tmp[1] = temp_PDF*B[0]*B[1]*B[2]*pro[1];
496
497 }else if(modetype[i] == 6){
498 temp_PDF = 0;
499 double t2A[4][4];
500 for(int ii=0; ii!=4; ii++){
501 pV[ii] = Pip[ii] + Pim[ii];
502 pA[ii] = pV[ii] + Kp[ii];
503 pD[ii] = pA[ii] + Pi0[ii];
504 }
505 sa[0] = SCADot(pV,pV);
506 sb[0] = SCADot(Pip,Pip);
507 sc[0] = SCADot(Pim,Pim);
508 sa[1] = SCADot(pA,pA);
509 sb[1] = sa[0];
510 sc[1] = SCADot(Kp,Kp);
511 sa[2] = SCADot(pD,pD);
512 sb[2] = sa[1];
513 sc[2] = SCADot(Pi0,Pi0);
514 if(flag[0] == 1){
515 propagatorGS(mass_R[1]*mass_R[1],mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2, pro1);
516 }
517 else if(flag[0] == 0){pro1[0] = 1; pro1[1] = 0;}
518
519 if(flag[1] == 1){
520 propagatorRBW(mass_R[0]*mass_R[0],mass_R[0], width_R[0], sa[1], sb[1], sc[1], rRes2, 1, pro2);
521 }
522 else if(flag[1] == 0){pro2[0] = 1; pro2[1] = 0;}
523 B[1] = Barrier(mass_R[1]*mass_R[1],1,sa[0],sb[0],sc[0],rRes2);//rho
524 B[2] = Barrier(mD*mD,1,sa[2],sb[2],sc[2],rD2);
525 calt1(Pip,Pim,t1V);
526 calt1(pA,Pi0,t1D);
527 if(flag[2] == 0){
528 for(int ii=0; ii!=4; ii++){
529 for(int j=0; j!=4; j++){
530 temp_PDF += t1D[ii]*(G[ii][j] - pA[ii]*pA[j]/sa[1])*t1V[j]*G[ii][ii]*G[j][j];
531 }
532 }
533 B[0] = 1;
534 }
535 else if(flag[2] == 2){
536 calt2(pV,Kp,t2A);
537 for(int ii=0; ii!=4; ii++){
538 for(int j=0; j!=4; j++){
539 temp_PDF += t1D[ii]*t2A[ii][j]*t1V[j]*G[ii][ii]*G[j][j];
540 }
541 }
542 B[0] = Barrier(mass_R[0]*mass_R[0],2,sa[1],sb[1],sc[1],rRes2);
543 }
544 Com_Multi(pro1,pro2,pro);
545 amp_tmp[0] = temp_PDF*B[0]*B[1]*B[2]*pro[0];
546 amp_tmp[1] = temp_PDF*B[0]*B[1]*B[2]*pro[1];
547
548 } else if(modetype[i] == 7){
549 temp_PDF = 0;
550 double t2A[4][4];
551 for(int ii=0; ii!=4; ii++){
552 pV[ii] = Kp[ii] + Pim[ii];
553 pA[ii] = pV[ii] + Pi0[ii];
554 pD[ii] = pA[ii] + Pip[ii];
555 }
556 sa[0] = SCADot(pV,pV);
557 sb[0] = SCADot(Kp,Kp);
558 sc[0] = SCADot(Pim,Pim);
559 sa[1] = SCADot(pA,pA);
560 sb[1] = sa[0];
561 sc[1] = SCADot(Pi0,Pi0);
562 sa[2] = SCADot(pD,pD);
563 sb[2] = sa[1];
564 sc[2] = SCADot(Pip,Pip);
565 if(flag[0] == 1){
566 propagatorRBW(mass_R[1]*mass_R[1],mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2, 1, pro1);
567 }
568 else if(flag[0] == 0){pro1[0] = 1; pro1[1] = 0;}
569 if(flag[1] == 1){
570 propagatorRBW(mass_R[0]*mass_R[0],mass_R[0], width_R[0], sa[1], sb[1], sc[1], rRes2, flag[2], pro2);
571 }
572 else if(flag[1] == 0){pro2[0] = 1; pro2[1] = 0;}
573
574 B[1] = Barrier(mass_R[1]*mass_R[1],1,sa[0],sb[0],sc[0],rRes2);
575 B[2] = Barrier(mD*mD,1,sa[2],sb[2],sc[2],rD2);
576 calt1(Kp,Pim,t1V);
577 calt1(pA,Pip,t1D);
578 if(flag[2] == 0){
579 for(int ii=0; ii!=4; ii++){
580 for(int j=0; j!=4; j++){
581 temp_PDF += t1D[ii]*(G[ii][j] - pA[ii]*pA[j]/sa[1])*t1V[j]*G[ii][ii]*G[j][j];
582 }
583 }
584 B[0] = 1;
585 }
586 else if(flag[2] == 2){
587 calt2(pV,Pi0,t2A);
588 for(int ii=0; ii!=4; ii++){
589 for(int j=0; j!=4; j++){
590 temp_PDF += t1D[ii]*t2A[ii][j]*t1V[j]*G[ii][ii]*G[j][j];
591 }
592 }
593 B[0] = Barrier(mass_R[0]*mass_R[0],2,sa[1],sb[1],sc[1],rRes2);
594 }
595 Com_Multi(pro1,pro2,pro);
596 amp_tmp[0] = temp_PDF*B[0]*B[1]*B[2]*pro[0];
597 amp_tmp[1] = temp_PDF*B[0]*B[1]*B[2]*pro[1];
598
599 } else if(modetype[i] == 8){
600 temp_PDF = 0;
601 double t2A[4][4];
602 for(int ii=0; ii!=4; ii++){
603 pV[ii] = Kp[ii] + Pi0[ii];
604 pA[ii] = pV[ii] + Pim[ii];
605 pD[ii] = pA[ii] + Pip[ii];
606 }
607 sa[0] = SCADot(pV,pV);
608 sb[0] = SCADot(Kp,Kp);
609 sc[0] = SCADot(Pi0,Pi0);
610 sa[1] = SCADot(pA,pA);
611 sb[1] = sa[0];
612 sc[1] = SCADot(Pim,Pim);
613 sa[2] = SCADot(pD,pD);
614 sb[2] = sa[1];
615 sc[2] = SCADot(Pip,Pip);
616 if(flag[0] == 1){
617 propagatorRBW(mass_R[1]*mass_R[1],mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2, 1, pro1);
618 }
619 else if(flag[0] == 0){pro1[0] = 1; pro1[1] = 0;}
620 if(flag[1] == 1){
621 propagatorRBW(mass_R[0]*mass_R[0],mass_R[0], width_R[0], sa[1], sb[1], sc[1], rRes2, flag[2], pro2);
622 }
623 else if(flag[1] == 0){pro2[0] = 1; pro2[1] = 0;}
624
625 B[1] = Barrier(mass_R[1]*mass_R[1],1,sa[0],sb[0],sc[0],rRes2);
626 B[2] = Barrier(mD*mD,1,sa[2],sb[2],sc[2],rD2);
627 calt1(Kp,Pi0,t1V);
628 calt1(pA,Pip,t1D);
629 if(flag[2] == 0){
630 for(int ii=0; ii!=4; ii++){
631 for(int j=0; j!=4; j++){
632 temp_PDF += t1D[ii]*(G[ii][j] - pA[ii]*pA[j]/sa[1])*t1V[j]*G[ii][ii]*G[j][j];
633 }
634 }
635 B[0] = 1;
636 }
637 else if(flag[2] == 2){
638 calt2(pV,Pim,t2A);
639 for(int ii=0; ii!=4; ii++){
640 for(int j=0; j!=4; j++){
641 temp_PDF += t1D[ii]*t2A[ii][j]*t1V[j]*G[ii][ii]*G[j][j];
642 }
643 }
644 B[0] = Barrier(mass_R[0]*mass_R[0],2,sa[1],sb[1],sc[1],rRes2);
645 }
646 Com_Multi(pro1,pro2,pro);
647 amp_tmp[0] = temp_PDF*B[0]*B[1]*B[2]*pro[0];
648 amp_tmp[1] = temp_PDF*B[0]*B[1]*B[2]*pro[1];
649
650 } else if(modetype[i] == 9){
651 temp_PDF = 0;
652 double t2A[4][4];
653 for(int ii=0; ii!=4; ii++){
654 pV[ii] = Kp[ii] + Pim[ii];
655 pA[ii] = pV[ii] + Pip[ii];
656 pD[ii] = pA[ii] + Pi0[ii];
657 }
658 sa[0] = SCADot(pV,pV);
659 sb[0] = SCADot(Kp,Kp);
660 sc[0] = SCADot(Pim,Pim);
661 sa[1] = SCADot(pA,pA);
662 sb[1] = sa[0];
663 sc[1] = SCADot(Pip,Pip);
664 sa[2] = SCADot(pD,pD);
665 sb[2] = sa[1];
666 sc[2] = SCADot(Pi0,Pi0);
667 if(flag[0] == 1){
668 propagatorRBW(mass_R[1]*mass_R[1],mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2, 1, pro1);
669 }
670 else if(flag[0] == 0){pro1[0] = 1; pro1[1] = 0;}
671 if(flag[1] == 1){
672 propagatorRBW(mass_R[0]*mass_R[0],mass_R[0], width_R[0], sa[1], sb[1], sc[1], rRes2, flag[2], pro2);
673 }
674 else if(flag[1] == 0){pro2[0] = 1; pro2[1] = 0;}
675
676 B[1] = Barrier(mass_R[1]*mass_R[1],1,sa[0],sb[0],sc[0],rRes2);
677 B[2] = Barrier(mD*mD,1,sa[2],sb[2],sc[2],rD2);
678 calt1(Kp,Pim,t1V);
679 calt1(pA,Pi0,t1D);
680 if(flag[2] == 0){
681 for(int ii=0; ii!=4; ii++){
682 for(int j=0; j!=4; j++){
683 temp_PDF += t1D[ii]*(G[ii][j] - pA[ii]*pA[j]/sa[1])*t1V[j]*G[ii][ii]*G[j][j];
684 }
685 }
686 B[0] = 1;
687 }
688 else if(flag[2] == 2){
689 calt2(pV,Pip,t2A);
690 for(int ii=0; ii!=4; ii++){
691 for(int j=0; j!=4; j++){
692 temp_PDF += t1D[ii]*t2A[ii][j]*t1V[j]*G[ii][ii]*G[j][j];
693 }
694 }
695 B[0] = Barrier(mass_R[0]*mass_R[0],2,sa[1],sb[1],sc[1],rRes2);
696 }
697 Com_Multi(pro1,pro2,pro);
698 amp_tmp[0] = temp_PDF*B[0]*B[1]*B[2]*pro[0];
699 amp_tmp[1] = temp_PDF*B[0]*B[1]*B[2]*pro[1];
700
701 } else if(modetype[i] == 11){
702 //rho0 pi0
703 temp_PDF = 0;
704 double pV1[4], pV2[4], qV1[4], qV2[4];
705 for(int ii=0; ii!=4; ii++){
706 pV2[ii] = Pip[ii] + Pim[ii];
707 qV2[ii] = Pip[ii] - Pim[ii];
708 pV1[ii] = pV2[ii] + Pi0[ii];
709 qV1[ii] = pV2[ii] - Pi0[ii];
710 pD[ii] = pV1[ii] + Kp[ii];
711 }
712 for(int ii=0; ii!=4; ii++){
713 for(int j=0; j!=4; j++){
714 for(int k=0; k!=4; k++){
715 for(int l=0; l!=4; l++){
716 temp_PDF += E[ii][j][k][l]*pV1[ii]*qV1[j]*Kp[k]*qV2[l]*
717 (G[ii][ii])*(G[j][j])*(G[k][k])*(G[l][l]);
718 }
719 }
720 }
721 }
722 sa[0] = SCADot(pV2,pV2);
723 sb[0] = SCADot(Pip,Pip);
724 sc[0] = SCADot(Pim,Pim);
725 sa[1] = SCADot(pV1,pV1);
726 sb[1] = sa[0];
727 sc[1] = SCADot(Pi0,Pi0);
728 sa[2] = SCADot(pD,pD);
729 sb[2] = sa[1];
730 sc[2] = SCADot(Kp,Kp);
731 propagatorGS(mass_R[1]*mass_R[1],mass_R[1],width_R[1],sa[0],sb[0],sc[0],rRes2, pro1);
732 if(flag[0] == 1){
733 propagatorRBW(mass_R[0]*mass_R[0],mass_R[0], width_R[0], sa[1], sb[1], sc[1], rRes2, 1, pro2);
734 }
735 else if(flag[0] == 0){pro1[0] = 1; pro1[1] = 0;}
736 Com_Multi(pro1,pro2,pro);
737
738 B[1] = Barrier(mass_R[1]*mass_R[1],1,sa[0],sb[0],sc[0],rRes2);
739 B[0] = Barrier(mass_R[0]*mass_R[0],1,sa[1],sb[1],sc[1],rRes2);
740 B[2] = Barrier(mD*mD,1,sa[2],sb[2],sc[2],rD2);
741
742 amp_tmp[0] = temp_PDF*B[0]*B[1]*B[2]*pro[0];
743 amp_tmp[1] = temp_PDF*B[0]*B[1]*B[2]*pro[1];
744
745 //rho+ pi-
746 temp_PDF = 0;
747 //double pV1[4], pV2[4], qV1[4], qV2[4];
748 for(int ii=0; ii!=4; ii++){
749 pV2[ii] = Pip[ii] + Pi0[ii];
750 qV2[ii] = Pip[ii] - Pi0[ii];
751 pV1[ii] = pV2[ii] + Pim[ii];
752 qV1[ii] = pV2[ii] - Pim[ii];
753 pD[ii] = pV1[ii] + Kp[ii];
754 }
755 for(int ii=0; ii!=4; ii++){
756 for(int j=0; j!=4; j++){
757 for(int k=0; k!=4; k++){
758 for(int l=0; l!=4; l++){
759 temp_PDF += E[ii][j][k][l]*pV1[ii]*qV1[j]*Kp[k]*qV2[l]*
760 (G[ii][ii])*(G[j][j])*(G[k][k])*(G[l][l]);
761 }
762 }
763 }
764 }
765 sa[0] = SCADot(pV2,pV2);
766 sb[0] = SCADot(Pip,Pip);
767 sc[0] = SCADot(Pi0,Pi0);
768 sa[1] = SCADot(pV1,pV1);
769 sb[1] = sa[0];
770 sc[1] = SCADot(Pim,Pim);
771 sa[2] = SCADot(pD,pD);
772 sb[2] = sa[1];
773 sc[2] = SCADot(Kp,Kp);
774 propagatorGS(mass_R[1]*mass_R[1],mass_R[1],width_R[1],sa[0],sb[0],sc[0],rRes2, pro1);
775 if(flag[0] == 1){
776 propagatorRBW(mass_R[0]*mass_R[0],mass_R[0], width_R[0], sa[1], sb[1], sc[1], rRes2, 1, pro2);
777 }
778 else if(flag[0] == 0){pro1[0] = 1; pro1[1] = 0;}
779 Com_Multi(pro1,pro2,pro);
780
781 B[1] = Barrier(mass_R[1]*mass_R[1],1,sa[0],sb[0],sc[0],rRes2);
782 B[0] = Barrier(mass_R[0]*mass_R[0],1,sa[1],sb[1],sc[1],rRes2);
783 B[2] = Barrier(mD*mD,1,sa[2],sb[2],sc[2],rD2);
784
785 amp_tmp[0] += temp_PDF*B[0]*B[1]*B[2]*pro[0];
786 amp_tmp[1] += temp_PDF*B[0]*B[1]*B[2]*pro[1];
787
788 //rho- pi+
789 temp_PDF = 0;
790 //double pV1[4], pV2[4], qV1[4], qV2[4];
791 for(int ii=0; ii!=4; ii++){
792 pV2[ii] = Pim[ii] + Pi0[ii];
793 qV2[ii] = Pim[ii] - Pi0[ii];
794 pV1[ii] = pV2[ii] + Pip[ii];
795 qV1[ii] = pV2[ii] - Pip[ii];
796 pD[ii] = pV1[ii] + Kp[ii];
797 }
798 for(int ii=0; ii!=4; ii++){
799 for(int j=0; j!=4; j++){
800 for(int k=0; k!=4; k++){
801 for(int l=0; l!=4; l++){
802 temp_PDF += E[ii][j][k][l]*pV1[ii]*qV1[j]*Kp[k]*qV2[l]*
803 (G[ii][ii])*(G[j][j])*(G[k][k])*(G[l][l]);
804 }
805 }
806 }
807 }
808 sa[0] = SCADot(pV2,pV2);
809 sb[0] = SCADot(Pim,Pim);
810 sc[0] = SCADot(Pi0,Pi0);
811 sa[1] = SCADot(pV1,pV1);
812 sb[1] = sa[0];
813 sc[1] = SCADot(Pip,Pip);
814 sa[2] = SCADot(pD,pD);
815 sb[2] = sa[1];
816 sc[2] = SCADot(Kp,Kp);
817 propagatorGS(mass_R[1]*mass_R[1],mass_R[1],width_R[1],sa[0],sb[0],sc[0],rRes2, pro1);
818 if(flag[0] == 1){
819 propagatorRBW(mass_R[0]*mass_R[0],mass_R[0], width_R[0], sa[1], sb[1], sc[1], rRes2, 1, pro2);
820 }
821 else if(flag[0] == 0){pro1[0] = 1; pro1[1] = 0;}
822 Com_Multi(pro1,pro2,pro);
823
824 B[1] = Barrier(mass_R[1]*mass_R[1],1,sa[0],sb[0],sc[0],rRes2);
825 B[0] = Barrier(mass_R[0]*mass_R[0],1,sa[1],sb[1],sc[1],rRes2);
826 B[2] = Barrier(mD*mD,1,sa[2],sb[2],sc[2],rD2);
827
828 amp_tmp[0] += temp_PDF*B[0]*B[1]*B[2]*pro[0];
829 amp_tmp[1] += temp_PDF*B[0]*B[1]*B[2]*pro[1];
830
831 }
832
833 Com_Multi(amp_tmp,cof,amp_PDF);
834 PDF[0] += amp_PDF[0];
835 PDF[1] += amp_PDF[1];
836 }
837 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
838 Result = value;
839}
840
841void EvtDToKppipipi0::Com_Multi(double a1[2], double a2[2], double res[2])
842{
843 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
844 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
845}
846void EvtDToKppipipi0::Com_Divide(double a1[2], double a2[2], double res[2])
847{
848 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
849 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
850 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
851}
852double EvtDToKppipipi0::SCADot(double a1[4], double a2[4])
853{
854 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
855 return _cal;
856}
857double EvtDToKppipipi0::Barrier(double mass2, int l, double sa, double sb, double sc, double r2)
858{
859 double F;
860 double tmp = sa+sb-sc;
861 double q = fabs(0.25*tmp*tmp/sa-sb);
862 //if (q < 0) q = 1e-16;
863 double tmp2 = mass2+sb-sc;
864 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
865 //if (q0 < 0) q0 = 1e-16;
866 double z = q*r2;
867 double z0 = q0*r2;
868 if (l==1) {
869 F = sqrt((1.0+z0)/(1.0+z));
870 }
871 else if (l==2) {
872 double z2 = z*z; double z02 = z0*z0;
873 F = sqrt((9.0+3.0*z0+z02)/(9.0+3.0*z+z2));
874 } else {
875 F = 1.0;
876 }
877 return F;
878}
879void EvtDToKppipipi0::calt1(double daug1[4], double daug2[4], double t1[4])
880{
881 double p, pq, tmp;
882 double pa[4], qa[4];
883 for(int i=0; i<4; i++) {
884 pa[i] = daug1[i] + daug2[i];
885 qa[i] = daug1[i] - daug2[i];
886 }
887 p = SCADot(pa,pa);
888 pq = SCADot(pa,qa);
889 tmp = pq/p;
890 for(int i=0; i<4; i++) {
891 t1[i] = qa[i] - tmp*pa[i];
892 }
893}
894void EvtDToKppipipi0::calt2(double daug1[4], double daug2[4], double t2[4][4])
895{
896 double p, r;
897 double pa[4], t1[4];
898 calt1(daug1,daug2,t1);
899 r = SCADot(t1,t1)/3.0;
900 for(int i=0; i<4; i++) {
901 pa[i] = daug1[i] + daug2[i];
902 }
903 p = SCADot(pa,pa);
904 for(int i=0; i<4; i++) {
905 for(int j=0; j<4; j++) {
906 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
907 }
908 }
909}
910void EvtDToKppipipi0::propagator(double mass2, double mass, double width, double sx, double prop[2])
911{
912 double a[2], b[2];
913 a[0] = 1;
914 a[1] = 0;
915 b[0] = mass2-sx;
916 b[1] = -mass*width;
917 Com_Divide(a,b,prop);
918}
919double EvtDToKppipipi0::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
920{
921 double widm = 0.;
922 double m = sqrt(sa);
923 double tmp = sb-sc;
924 double tmp1 = sa+tmp;
925 double q = fabs(0.25*tmp1*tmp1/sa-sb);
926 //if(q<0) q = 1e-16;
927 double tmp2 = mass2+tmp;
928 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
929 //if(q0<0) q0 = 1e-16;
930 double z = q*r2;
931 double z0 = q0*r2;
932 double t = q/q0;
933 if(l == 0) {widm = sqrt(t)*mass/m;}
934 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
935 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
936 return widm;
937}
938double EvtDToKppipipi0::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
939{
940 double widm = 0.;
941 double m = sqrt(sa);
942 double tmp = sb-sc;
943 double tmp1 = sa+tmp;
944 double q = fabs(0.25*tmp1*tmp1/sa-sb);
945 //if(q<0) q = 1e-16;
946 double tmp2 = mass2+tmp;
947 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
948 //if(q0<0) q0 = 1e-16;
949 double z = q*r2;
950 double z0 = q0*r2;
951 double F = (1+z0)/(1+z);
952 double t = q/q0;
953 widm = t*sqrt(t)*mass/m*F;
954 return widm;
955}
956void EvtDToKppipipi0::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
957{
958 double a[2], b[2];
959 a[0] = 1;
960 a[1] = 0;
961 b[0] = mass2-sa;
962 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
963 Com_Divide(a,b,prop);
964}
965void EvtDToKppipipi0::propagatorRBWl1(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
966{
967 double a[2], b[2];
968 a[0] = 1;
969 a[1] = 0;
970 b[0] = mass2-sa;
971 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
972 Com_Divide(a,b,prop);
973}
974void EvtDToKppipipi0::propagatorRBW_a1(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
975{
976 double a[2], b[2];
977 int iii=int(sqrt(sa)*1000)-1;
978 double a1width[3000]={
979 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
980 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
981 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
982 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
983 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
984 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
985 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
986 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
987 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
988 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
989 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
990 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
991 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
992 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
993 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
994 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
995 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
996 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
997 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
998 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
999 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1000 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1001 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1002 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1003 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1004 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1005 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1006 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1007 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1008 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1009 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1010 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1011 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1012 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1013 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1014 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1015 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1016 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1017 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1018 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1019 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1020 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1021 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1022 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1023 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1024 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1025 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1026 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1027 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1028 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1029 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1030 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1031 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,
1032 0.000000,0.000001,0.000001,0.000001,0.000002,0.000002,0.000002,0.000003,
1033 0.000004,0.000004,0.000005,0.000006,0.000007,0.000008,0.000009,0.000010,
1034 0.000011,0.000012,0.000014,0.000015,0.000017,0.000019,0.000021,0.000023,
1035 0.000025,0.000027,0.000029,0.000032,0.000035,0.000038,0.000041,0.000044,
1036 0.000047,0.000050,0.000054,0.000058,0.000062,0.000066,0.000070,0.000075,
1037 0.000079,0.000084,0.000089,0.000094,0.000100,0.000105,0.000111,0.000117,
1038 0.000124,0.000130,0.000137,0.000143,0.000151,0.000158,0.000165,0.000173,
1039 0.000182,0.000190,0.000199,0.000207,0.000216,0.000225,0.000235,0.000245,
1040 0.000256,0.000266,0.000277,0.000288,0.000300,0.000311,0.000322,0.000335,
1041 0.000347,0.000360,0.000373,0.000385,0.000400,0.000415,0.000429,0.000442,
1042 0.000457,0.000473,0.000488,0.000504,0.000520,0.000539,0.000555,0.000572,
1043 0.000590,0.000608,0.000626,0.000646,0.000664,0.000684,0.000704,0.000725,
1044 0.000745,0.000766,0.000787,0.000809,0.000828,0.000854,0.000878,0.000901,
1045 0.000927,0.000952,0.000973,0.001001,0.001027,0.001048,0.001080,0.001104,
1046 0.001132,0.001159,0.001189,0.001219,0.001245,0.001277,0.001308,0.001338,
1047 0.001370,0.001404,0.001433,0.001468,0.001498,0.001533,0.001570,0.001600,
1048 0.001638,0.001678,0.001711,0.001745,0.001780,0.001825,0.001857,0.001898,
1049 0.001941,0.001972,0.002017,0.002065,0.002104,0.002146,0.002189,0.002234,
1050 0.002277,0.002319,0.002369,0.002410,0.002461,0.002511,0.002557,0.002605,
1051 0.002661,0.002704,0.002762,0.002807,0.002855,0.002910,0.002965,0.003020,
1052 0.003074,0.003127,0.003178,0.003228,0.003288,0.003351,0.003409,0.003471,
1053 0.003532,0.003598,0.003660,0.003720,0.003793,0.003854,0.003910,0.003972,
1054 0.004050,0.004108,0.004181,0.004254,0.004309,0.004380,0.004464,0.004533,
1055 0.004603,0.004679,0.004756,0.004811,0.004898,0.004974,0.005048,0.005142,
1056 0.005215,0.005279,0.005363,0.005449,0.005533,0.005604,0.005695,0.005783,
1057 0.005869,0.005971,0.006060,0.006142,0.006247,0.006332,0.006409,0.006502,
1058 0.006594,0.006713,0.006784,0.006889,0.006995,0.007079,0.007190,0.007303,
1059 0.007381,0.007487,0.007592,0.007710,0.007801,0.007910,0.008032,0.008149,
1060 0.008247,0.008378,0.008462,0.008559,0.008706,0.008843,0.008943,0.009091,
1061 0.009207,0.009308,0.009448,0.009555,0.009698,0.009810,0.009936,0.010020,
1062 0.010186,0.010320,0.010474,0.010611,0.010742,0.010840,0.011011,0.011167,
1063 0.011281,0.011395,0.011541,0.011714,0.011853,0.012046,0.012169,0.012277,
1064 0.012460,0.012617,0.012805,0.012922,0.013072,0.013234,0.013389,0.013561,
1065 0.013704,0.013917,0.014025,0.014239,0.014425,0.014600,0.014716,0.014958,
1066 0.015114,0.015325,0.015488,0.015630,0.015797,0.016035,0.016206,0.016404,
1067 0.016591,0.016842,0.016964,0.017199,0.017392,0.017557,0.017798,0.017987,
1068 0.018178,0.018337,0.018631,0.018829,0.019008,0.019221,0.019467,0.019698,
1069 0.019941,0.020166,0.020379,0.020585,0.020806,0.021040,0.021309,0.021482,
1070 0.021764,0.022046,0.022306,0.022478,0.022736,0.023049,0.023222,0.023536,
1071 0.023756,0.024090,0.024337,0.024557,0.024873,0.025099,0.025392,0.025682,
1072 0.025995,0.026291,0.026498,0.026927,0.027119,0.027377,0.027804,0.028135,
1073 0.028279,0.028682,0.028871,0.029355,0.029531,0.029956,0.030243,0.030592,
1074 0.030873,0.031246,0.031494,0.031771,0.032167,0.032515,0.032881,0.033211,
1075 0.033653,0.033988,0.034394,0.034639,0.035055,0.035569,0.035879,0.036211,
1076 0.036611,0.036932,0.037489,0.037779,0.038284,0.038723,0.039001,0.039574,
1077 0.039854,0.040274,0.040881,0.041117,0.041644,0.042055,0.042531,0.043030,
1078 0.043354,0.043832,0.044277,0.044956,0.045284,0.045828,0.046440,0.046800,
1079 0.047518,0.047727,0.048258,0.048850,0.049316,0.049992,0.050486,0.050987,
1080 0.051410,0.051928,0.052613,0.053110,0.053824,0.054351,0.055078,0.055654,
1081 0.056030,0.056763,0.057245,0.057832,0.058569,0.059292,0.060048,0.060569,
1082 0.061056,0.061869,0.062612,0.063186,0.063886,0.064655,0.065198,0.065815,
1083 0.066649,0.067577,0.068012,0.068967,0.069630,0.070181,0.070786,0.071989,
1084 0.072764,0.073466,0.074461,0.075093,0.075994,0.076834,0.077455,0.078709,
1085 0.079581,0.080408,0.080884,0.081965,0.082882,0.083658,0.084824,0.085513,
1086 0.086662,0.087602,0.088678,0.089492,0.090641,0.091369,0.092494,0.093484,
1087 0.094615,0.095385,0.096168,0.097668,0.098611,0.099630,0.100772,0.102020,
1088 0.103145,0.104110,0.105071,0.106604,0.107791,0.108451,0.109509,0.111356,
1089 0.112026,0.113921,0.114507,0.116071,0.117027,0.118213,0.120164,0.120701,
1090 0.122121,0.123894,0.124937,0.126134,0.127391,0.128882,0.130056,0.131649,
1091 0.133046,0.134275,0.135119,0.137072,0.138476,0.139612,0.140388,0.142734,
1092 0.143576,0.145445,0.147414,0.148856,0.149891,0.150963,0.152477,0.153717,
1093 0.155275,0.156859,0.158462,0.159257,0.161865,0.163182,0.164465,0.165538,
1094 0.167003,0.169257,0.171211,0.172093,0.173261,0.174639,0.176510,0.177684,
1095 0.179077,0.181041,0.182446,0.184769,0.184926,0.186741,0.188844,0.190884,
1096 0.191714,0.192254,0.193921,0.195917,0.196766,0.199052,0.200603,0.201808,
1097 0.202699,0.204636,0.205712,0.206849,0.208741,0.209424,0.211698,0.212753,
1098 0.215516,0.215857,0.217790,0.217774,0.220454,0.221821,0.223466,0.224494,
1099 0.225632,0.227231,0.229456,0.229581,0.231537,0.232263,0.233834,0.234725,
1100 0.237079,0.238015,0.239400,0.240193,0.241693,0.243787,0.244317,0.244971,
1101 0.246711,0.248615,0.249387,0.250905,0.252702,0.253535,0.254385,0.255375,
1102 0.256671,0.258405,0.259741,0.260875,0.262131,0.262920,0.264860,0.265893,
1103 0.266016,0.267727,0.270039,0.270689,0.271047,0.272313,0.272474,0.274724,
1104 0.275813,0.275937,0.278793,0.278783,0.281407,0.281351,0.282481,0.284226,
1105 0.284113,0.284999,0.285655,0.288361,0.287856,0.288893,0.290211,0.291708,
1106 0.291985,0.294298,0.294849,0.296796,0.296197,0.296851,0.298011,0.300368,
1107 0.299982,0.302378,0.304363,0.303711,0.304729,0.306789,0.306378,0.307372,
1108 0.308720,0.309509,0.309712,0.310782,0.311699,0.312668,0.312755,0.313675,
1109 0.315311,0.316640,0.317217,0.317403,0.318478,0.319916,0.321803,0.322678,
1110 0.323237,0.324343,0.324433,0.324493,0.324969,0.325894,0.328563,0.328721,
1111 0.328954,0.330640,0.328164,0.331267,0.331695,0.333772,0.333619,0.334351,
1112 0.334605,0.336434,0.337510,0.336535,0.337362,0.338799,0.340732,0.339896,
1113 0.342707,0.343471,0.342318,0.342431,0.344543,0.345611,0.345786,0.346590,
1114 0.346610,0.347761,0.348914,0.349558,0.350577,0.352128,0.350982,0.354134,
1115 0.352773,0.353213,0.352972,0.354927,0.355784,0.355778,0.355801,0.357040,
1116 0.358013,0.358432,0.360045,0.359743,0.360238,0.359850,0.362184,0.361580,
1117 0.363430,0.362333,0.364397,0.364472,0.364370,0.365303,0.366644,0.367777,
1118 0.368604,0.367631,0.368324,0.369782,0.371121,0.370653,0.370040,0.371649,
1119 0.370201,0.373362,0.373900,0.374159,0.374916,0.374503,0.376703,0.372802,
1120 0.376191,0.379596,0.377325,0.376363,0.379369,0.379791,0.378703,0.380177,
1121 0.381762,0.381335,0.381374,0.384668,0.381763,0.382746,0.384723,0.385089,
1122 0.386229,0.386702,0.387749,0.384423,0.384714,0.384181,0.388489,0.388618,
1123 0.388179,0.390092,0.389871,0.390496,0.391181,0.390679,0.392614,0.392269,
1124 0.393899,0.393466,0.391421,0.391090,0.395586,0.391776,0.396882,0.393254,
1125 0.394400,0.395749,0.398063,0.397138,0.397585,0.397288,0.397847,0.395375,
1126 0.400170,0.400007,0.401191,0.398513,0.401922,0.400477,0.404257,0.403271,
1127 0.400677,0.403913,0.403172,0.404727,0.403406,0.404404,0.405265,0.406389,
1128 0.405738,0.402173,0.407831,0.405895,0.409172,0.408934,0.405915,0.408486,
1129 0.407320,0.407437,0.405444,0.408400,0.410909,0.412427,0.409881,0.411021,
1130 0.413001,0.410369,0.414702,0.413372,0.413095,0.410972,0.416346,0.416095,
1131 0.414132,0.414344,0.416952,0.415197,0.417583,0.416582,0.416622,0.416895,
1132 0.416576,0.415551,0.417925,0.414838,0.417051,0.416831,0.420000,0.419132,
1133 0.418173,0.417645,0.419679,0.419866,0.419581,0.421531,0.420878,0.422737,
1134 0.421872,0.421304,0.425486,0.424434,0.420842,0.426753,0.422761,0.422178,
1135 0.422372,0.424173,0.425582,0.425080,0.425831,0.423551,0.422949,0.425784,
1136 0.427977,0.427948,0.426368,0.425138,0.425351,0.428643,0.428148,0.427488,
1137 0.431704,0.430167,0.429655,0.429584,0.425458,0.430728,0.429845,0.431145,
1138 0.429180,0.428874,0.430720,0.430024,0.432034,0.431359,0.431535,0.432995,
1139 0.432425,0.432454,0.433140,0.432574,0.433814,0.433348,0.432886,0.435472,
1140 0.436517,0.432681,0.436999,0.435182,0.434834,0.435478,0.438255,0.436650,
1141 0.434464,0.438530,0.434077,0.436471,0.434012,0.436822,0.437505,0.440135,
1142 0.438322,0.438032,0.439001,0.440270,0.438661,0.439233,0.439274,0.437945,
1143 0.443080,0.439191,0.438233,0.440415,0.441063,0.440926,0.440929,0.439731,
1144 0.443584,0.439729,0.441597,0.442615,0.444637,0.443180,0.440789,0.440261,
1145 0.442202,0.445081,0.445484,0.445415,0.445532,0.442806,0.444188,0.441073,
1146 0.444299,0.445897,0.445279,0.442830,0.445506,0.445272,0.447267,0.443522,
1147 0.445519,0.446459,0.446753,0.446377,0.446129,0.446383,0.448556,0.446593,
1148 0.445293,0.449199,0.447590,0.445968,0.447482,0.448474,0.449890,0.450004,
1149 0.447765,0.449274,0.450652,0.448210,0.449360,0.449577,0.448575,0.452112,
1150 0.448780,0.451393,0.450200,0.452018,0.451182,0.452050,0.451748,0.451377,
1151 0.451402,0.448810,0.452311,0.452909,0.452491,0.452418,0.454190,0.454420,
1152 0.452121,0.452307,0.456857,0.453506,0.454058,0.457203,0.454394,0.453596,
1153 0.452240,0.453692,0.456516,0.453753,0.455541,0.452702,0.456481,0.452226,
1154 0.454280,0.454855,0.456297,0.456482,0.454154,0.455387,0.454748,0.455764,
1155 0.457282,0.455487,0.454822,0.454257,0.457678,0.454225,0.458689,0.456123,
1156 0.457011,0.457386,0.458351,0.458638,0.456164,0.455884,0.458525,0.457575,
1157 0.458340,0.458912,0.457836,0.461734,0.457545,0.460755,0.460960,0.459226,
1158 0.458613,0.461078,0.460958,0.460337,0.460237,0.461190,0.460760,0.457911,
1159 0.461310,0.459657,0.461960,0.461040,0.459578,0.461650,0.461550,0.461251,
1160 0.461054,0.463082,0.461732,0.461324,0.462547,0.461261,0.461629,0.464067,
1161 0.462430,0.462525,0.464232,0.462921,0.463202,0.465558,0.462914,0.461698,
1162 0.463963,0.463040,0.464275,0.461940,0.462913,0.465261,0.461500,0.463679,
1163 0.463354,0.465205,0.464529,0.462220,0.464279,0.463427,0.465387,0.465288,
1164 0.464839,0.464926,0.466100,0.465531,0.466187,0.464647,0.466285,0.465461,
1165 0.464134,0.466783,0.466763,0.466183,0.467089,0.464497,0.466080,0.466109,
1166 0.468166,0.466984,0.465335,0.466721,0.466856,0.465113,0.468377,0.467904,
1167 0.464546,0.468787,0.465648,0.469841,0.469477,0.466311,0.468700,0.465183,
1168 0.466559,0.470433,0.468563,0.468109,0.466980,0.467567,0.467670,0.466991,
1169 0.467992,0.468784,0.469406,0.469652,0.468527,0.470460,0.467308,0.470693,
1170 0.469539,0.468000,0.469295,0.467038,0.471908,0.468829,0.470663,0.469266,
1171 0.468975,0.470222,0.468649,0.469507,0.472307,0.471611,0.470419,0.471181,
1172 0.471140,0.473187,0.471086,0.469801,0.472234,0.472131,0.468996,0.470229,
1173 0.471597,0.469625,0.472230,0.470164,0.468404,0.472264,0.471336,0.471597,
1174 0.472280,0.471256,0.473151,0.471863,0.474458,0.471956,0.473099,0.473956,
1175 0.471725,0.472809,0.473065,0.473180,0.470611,0.473614,0.474263,0.472792,
1176 0.473543,0.472656,0.469728,0.473431,0.474538,0.475322,0.474962,0.473598,
1177 0.474114,0.473486,0.472934,0.473252,0.477149,0.471719,0.476383,0.473076,
1178 0.473952,0.473104,0.472459,0.474433,0.474494,0.473588,0.473839,0.478113,
1179 0.472435,0.475571,0.475194,0.475626,0.474617,0.474520,0.474472,0.476437,
1180 0.474512,0.474497,0.474628,0.476203,0.475698,0.473907,0.477144,0.479000,
1181 0.475553,0.477481,0.473998,0.476672,0.477115,0.477114,0.476282,0.476152,
1182 0.477009,0.479854,0.474354,0.477645,0.477517,0.477111,0.474843,0.476173,
1183 0.477321,0.477384,0.477880,0.475726,0.476004,0.478204,0.475586,0.477973,
1184 0.477935,0.480640,0.478234,0.476349,0.477493,0.476994,0.479815,0.477771,
1185 0.476333,0.476325,0.478245,0.477284,0.479238,0.478339,0.478966,0.478012,
1186 0.479304,0.480148,0.476125,0.481267,0.479801,0.476720,0.478898,0.479284,
1187 0.479153,0.480157,0.478681,0.479712,0.478993,0.479943,0.478349,0.478930,
1188 0.478052,0.477173,0.479244,0.480454,0.479128,0.480530,0.477843,0.478369,
1189 0.478561,0.478639,0.479191,0.481763,0.481321,0.480979,0.479702,0.479777,
1190 0.479384,0.477571,0.481880,0.478615,0.481303,0.478783,0.479384,0.480517,
1191 0.481928,0.481199,0.479041,0.479188,0.481491,0.482840,0.478766,0.481941,
1192 0.481298,0.478105,0.482933,0.479744,0.483361,0.482332,0.482556,0.482057,
1193 0.483616,0.480599,0.482245,0.481091,0.480871,0.481938,0.480678,0.481851,
1194 0.482902,0.482158,0.480187,0.481772,0.484967,0.483094,0.482133,0.483929,
1195 0.483354,0.483382,0.483964,0.479941,0.481375,0.480255,0.482184,0.482541,
1196 0.482032,0.483484,0.479492,0.483305,0.481070,0.483573,0.485689,0.485767,
1197 0.484221,0.481365,0.482440,0.481507,0.483418,0.480978
1198 };
1199 double width_a1=a1width[iii];
1200 a[0] = 1;
1201 a[1] = 0;
1202 b[0] = mass2-sa;
1203 //b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
1204 b[1] = -mass*width_a1;
1205 Com_Divide(a,b,prop);
1206}
1207void EvtDToKppipipi0::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
1208{
1209
1210 double GS1 = 0.636619783;
1211 double GS2 = 0.01860182466;
1212 double GS3 = 0.1591549458;
1213 double GS4 = 0.00620060822;
1214 double a[2], b[2];
1215 double tmp = sb-sc;
1216 double tmp1 = sa+tmp;
1217 double q2 = fabs(0.25*tmp1*tmp1/sa-sb);
1218 //if(q2<0) q2 = 1e-16;
1219
1220 double tmp2 = mass2+tmp;
1221 double q02 = fabs(0.25*tmp2*tmp2/mass2-sb);
1222 //if(q02<0) q02 = 1e-16;
1223
1224 double q = sqrt(q2);
1225 double q0 = sqrt(q02);
1226 double m = sqrt(sa);
1227 double q03 = q0*q02;
1228 double tmp3 = log(mass+2*q0)+1.2760418309; // log(mass_2Pion) = 1.2760418309;
1229
1230 double h = GS1*q/m*(log(m+2*q)+1.2760418309);
1231 double h0 = GS1*q0/mass*tmp3;
1232 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
1233 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
1234 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
1235
1236 a[0] = 1.0+d*width/mass;
1237 a[1] = 0.0;
1238 b[0] = mass2-sa+width*f;
1239 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
1240 Com_Divide(a,b,prop);
1241}
1242void EvtDToKppipipi0::rhoab(double sa, double sb, double sc, double res[2]) {
1243 double tmp = sa+sb-sc;
1244 double q = 0.25*tmp*tmp/sa-sb;
1245 if(q>=0) {
1246 res[0]=2.0*sqrt(q/sa);
1247 res[1]=0.0;
1248 } else {
1249 res[0]=0.0;
1250 res[1]=2.0*sqrt(-q/sa);
1251 }
1252}
1253void EvtDToKppipipi0::rho4Pi(double sa, double res[2]) {
1254 double temp = 1.0-0.3116765584/sa; // 0.3116765584=0.13957*0.13957*16
1255 if(temp>=0) {
1256 res[0]=sqrt(temp)/(1.0+exp(9.8-3.5*sa));
1257 res[1]=0.0;
1258 } else {
1259 res[0]=0.0;
1260 res[1]=sqrt(-temp)/(1.0+exp(9.8-3.5*sa));
1261 }
1262}
1263void EvtDToKppipipi0::propagatorsigma500(double sa, double sb, double sc, double prop[2]) {
1264 double f = 0.5843+1.6663*sa;
1265 const double M = 0.9264;
1266 const double mass2 = 0.85821696; // M*M
1267 const double mpi2d2 = 0.00973989245;
1268 double g1 = f*(sa-mpi2d2)/(mass2-mpi2d2)*exp((mass2-sa)/1.082);
1269 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
1270 rhoab(sa,sb,sc,rho1s);
1271 rhoab(mass2,sb,sc,rho1M);
1272 rho4Pi(sa,rho2s);
1273 rho4Pi(mass2,rho2M);
1274 Com_Divide(rho1s,rho1M,rho1);
1275 Com_Divide(rho2s,rho2M,rho2);
1276 double a[2], b[2];
1277 a[0] = 1.0;
1278 a[1] = 0.0;
1279 b[0] = mass2-sa+M*(g1*rho1[1]+0.0024*rho2[1]);
1280 b[1] = -M*(g1*rho1[0]+0.0024*rho2[0]);
1281 Com_Divide(a,b,prop);
1282}
1283void EvtDToKppipipi0::Flatte_rhoab(double sa, double sb, double sc, double rho[2])
1284{
1285 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
1286 if(q>0) {
1287 rho[0]=2* sqrt(q/sa);
1288 rho[1]=0;
1289 }
1290 else if(q<0){
1291 rho[0]=0;
1292 rho[1]=2*sqrt(-q/sa);
1293 }
1294}
1295void EvtDToKppipipi0::propagator980(double mass, double sx, double *sb, double *sc, double prop[2])
1296{
1297 double unit[2]={1.0};
1298 double ci[2]={0,1};
1299 double rho1[2];
1300 Flatte_rhoab(sx,sb[0],sc[0],rho1);
1301 double rho2[2];
1302 Flatte_rhoab(sx,sb[1],sc[1],rho2);
1303 double gK_f980=0.69465, gPi_f980=0.165;
1304 double tmp1[2]={gK_f980,0};
1305 double tmp11[2];
1306 double tmp2[2]={gPi_f980,0};
1307 double tmp22[2];
1308 Com_Multi(tmp1,rho1,tmp11);
1309 Com_Multi(tmp2,rho2,tmp22);
1310 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
1311 double tmp31[2];
1312 Com_Multi(tmp3, ci,tmp31);
1313 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
1314 Com_Divide( unit,tmp4, prop);
1315}
1316void EvtDToKppipipi0::KPiSLASS(double sa, double sb, double sc, double prop[2]) {
1317 const double m1430 = 1.441;
1318 const double sa0 = 2.076481; // m1430*m1430;
1319 const double w1430 = 0.193;
1320 const double Lass1 = 0.25/sa0;
1321 double tmp = sb-sc;
1322 double tmp1 = sa0+tmp;
1323 double q0 = fabs(Lass1*tmp1*tmp1-sb);
1324 //if(q0<0) q0 = 1e-16;
1325 double tmp2 = sa+tmp;
1326 double qs = fabs(0.25*tmp2*tmp2/sa-sb);
1327 double q = sqrt(qs);
1328 double width = w1430*q*m1430/sqrt(sa*q0);
1329 double temp_R = atan(m1430*width/(sa0-sa));
1330 if(temp_R<0) temp_R += math_pi;
1331 double deltaR = -109.7*math_pi/180.0 + temp_R;
1332 double temp_F = atan(0.226*q/(2.0-3.8194*qs)); // 2.0*0.113 = 0.226; -33.8*0.113 = -3.8194
1333 if(temp_F<0) temp_F += math_pi;
1334 double deltaF = 0.1*math_pi/180.0 + temp_F;
1335 double deltaS = deltaR + 2.0*deltaF;
1336 double t1 = 0.96*sin(deltaF);
1337 double t2 = sin(deltaR);
1338 double CF[2], CS[2];
1339 CF[0] = cos(deltaF);
1340 CF[1] = sin(deltaF);
1341 CS[0] = cos(deltaS);
1342 CS[1] = sin(deltaS);
1343 prop[0] = t1*CF[0] + t2*CS[0];
1344 prop[1] = t1*CF[1] + t2*CS[1];
1345
1346}
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
a1width
double a1width[3000]
Definition EvtDsTopipipipi0pi0.cc:300
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
EvtDToKppipipi0::getName
void getName(std::string &name)
Definition EvtDToKppipipi0.cc:34
EvtDToKppipipi0::clone
EvtDecayBase * clone()
Definition EvtDToKppipipi0.cc:38
EvtDToKppipipi0::decay
void decay(EvtParticle *p)
Definition EvtDToKppipipi0.cc:136
EvtDToKppipipi0::~EvtDToKppipipi0
virtual ~EvtDToKppipipi0()
Definition EvtDToKppipipi0.cc:32
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
rb::Grhop
const double Grhop
Definition TConstant.h:19
b
const double b
Definition slope.cxx:9