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